Complex eco-evolutionary dynamics induced by the coevolution of predator–prey movement strategies

The coevolution of predators and prey has been the subject of much empirical and theoretical research that produced intriguing insights into the interplay of ecology and evolution. To allow for mathematical analysis, models of predator–prey coevolution are often coarse-grained, focussing on population-level processes and largely neglecting individual-level behaviour. As selection is acting on individual-level properties, we here present a more mechanistic approach: an individual-based simulation model for the coevolution of predators and prey on a fine-grained resource landscape, where features relevant for ecology (like changes in local densities) and evolution (like differences in survival and reproduction) emerge naturally from interactions between individuals. Our focus is on predator–prey movement behaviour, and we present a new method for implementing evolving movement strategies in an efficient and intuitively appealing manner. Throughout their lifetime, predators and prey make repeated movement decisions on the basis of their movement strategies. Over the generations, the movement strategies evolve, as individuals that successfully survive and reproduce leave their strategy to more descendants. We show that the movement strategies in our model evolve rapidly, thereby inducing characteristic spatial patterns like spiral waves and static spots. Transitions between these patterns occur frequently, induced by antagonistic coevolution rather than by external events. Regularly, evolution leads to the emergence and stable coexistence of qualitatively different movement strategies within the same population. Although the strategy space of our model is continuous, we often observe the evolution of discrete movement types. We argue that rapid evolution, coexistent movement types, and phase shifts between different ecological regimes are not a peculiarity of our model but a result of more realistic assumptions on eco-evolutionary feedbacks and the number of evolutionary degrees of freedom.


Introduction
Predator-prey coevolution has fascinated biologists for decades (Cott 1940;Pimentel 1961;Levin & Udovic 1977;Dawkins & Krebs 1979). In the ecological arena, predator-prey interactions can lead to complex non-equilibrium dynamics (Turchin 2003). On top of these ecological predator-prey interactions, an evolutionary arms race may occur, where adaptive changes in the prey population impose new selective pressures on the predator population, and vice versa. Experimental findings suggest that the ecological and the evolutionary dynamics can be intertwined in an intricate manner (Yoshida et al. 2003(Yoshida et al. , 2007Becks et al. 2010). In natural systems, it is a major challenge to unravel this complexity (Hendry 2019). It is therefore no surprise that theoretical models have played a crucial role for the understanding of the ecology and evolution of predator-prey interactions (Fussmann et al. 2007, Govaert et al. 2019. Predator-prey coevolution models have traditionally been based on the frameworks of population genetics (Nuismer et al. 2005;Kopp & Gavrilets 2006;Cortez & Weitz 2014;Yamamichi & Ellner 2016), quantitative genetics (Gavrilets 1997;Mougi & Iwasa 2010;Cortez 2018), and adaptive dynamics Marrow et al. 1996;Flaxman & Lou 2009), where each approach encompasses a wide range of models. The approaches differ in their assumptions on the nature of genetic variation (discrete vs continuous), the occurrence and distribution of mutations, and the interaction of alleles within and across loci, but they have in common that they strive for analytical tractability. To achieve this, highly simplifying assumptions need to be made on the traits that are the target of selection. For example, predators are often characterized by a one-dimensional attack strategy, prey by a one-dimensional avoidance strategy, and prey capture rates are assumed to be maximal when the predator attack strategy matches the prey avoidance strategy ( Van der Laan & Hogeweg 1995, Marrow et al. 1996, Gavrilets 1997, Nuismer et al. 2005, Kopp & Gavrilets 2006, Yamamichi & Ellner 2016, see Abrams 2000 for a general overview). Such simplification allows for an elegant and seemingly general characterisation of coevolutionary outcomes, but the question arises whether it captures the essence of predator-prey interactions, which in natural systems are mediated by complex behavioural action and reaction patterns. More recent studies have explored a multidimensional trait space (Gilman et al. 2012;Débarre et al. 2014), demonstrating that an increasing number of traits may have a destabilizing effect, as there are more possibilities for instability. However, these models also assume that the predation rate is determined by the match between unspecified traits of the two species.
Behavioural traits play an important role in predator-prey interactions. Experiments have demonstrated that predators and prey show strong behavioural responses to each other (Gilliam & Fraser 1987;Savino & Stein 1989;Ehlinger 1990;Hammond et al. 2007;Simon et al. 2019). These responses differ across species and spatial scales, and they are likely the product of strategies that incorporate various information sources from the environment. In the literature, behavioural interactions between predators and prey are often discussed in game-theoretical terms as 'behavioural response races' (Sih 1984(Sih , 2005 or 'predator-prey shell games' (Mitchell & Lima 2002). These games usually play out in space, where resources, prey and predators are heterogeneously distributed among different patches. Prey have to balance foraging for a resource with predator avoidance, while predators are faced with the task of predicting prey behaviour. A full dynamical analysis of such interactions is a forbidding task (Flaxman & Lou 2009). Therefore, analytical approaches have to use short-cuts, such as the assumption that at all times individual predators or prey behave in such a way that they maximize their fitness under their given local circumstances (Iwasa 1982;Gilliam & Fraser 1987;Abrams 2007). It is often doubtful whether such short-cuts are realistic. For example, it is unlikely that evolution will fine-tune behaviour to such an extent that it is locally optimal under all circumstances and when conditions rapidly change (McNamara & Houston 2009;McNamara & Weissing 2010). It is therefore important to complement analytical theory with simulation-based approaches that can make more realistic assumptions on the ecological setting and the (evolving) traits governing the behaviour of predator and prey.
Several such simulation models have been developed (Huse et al. 1999;Kimbrell & Holt 2004;Flaxman et al. 2011;Patin et al. 2020). All these models are individual-based, meaning that both the ecology and the evolution of predator-prey interactions reflect the fate of individual agents. Such an approach is natural because selection acts on individual characteristics, while population-level phenomena are aggregates and/or emerging features of these characteristics. Additionally, individual-based models enable the implementation of mechanisms, both at the level of individual behaviour and with regard to the environmental setting. The existing individual-based models have shown that they can validate findings of analytical models (Huse et al. 1999, cf. Iwasa 1982 and theoretical expectations such as the ideal free distribution (Flaxman et al. 2011). As illustrated by the recent study of Patin and colleagues (2020), such models can unravel the importance of random movement, memory use, and other factors that are difficult to study within the scope of analytical models.
Virtually all individual-based models consider individual interactions at a coarsegrained spatial scale (for an exception see Kimbrell & Holt 2004). The environment is assumed to be structured in discrete patches and individuals have the task of choosing a patch that provides an optimal balance between resource abundance and safety. Within patches, predator-prey interactions are governed by patch-level population dynamics and not based on single individuals moving and behaving in space. Accordingly, these models do not capture interactions that take place at the individual level.
Here, we consider a model where individuals move and interact in more fine-grained space, where only few individuals co-occur at the same location. Predators and prey evolve situation-dependent movement strategies that determine the likelihood of finding food resources and avoiding predation. The implementation of movement strategies is a crucial ingredient of our model. Instead of assuming that movement is directly guided by fitness expectations (as, for example, in Iwasa 1982;Abrams 2007), the strategies in our model are different realisations of an inherited proximate mechanism that uses environmental inputs to evaluate the 'suitability' of environmental situations. We assume that predator and prey individuals continually scan their environment and judge the suitability of each movement direction on basis of the local densities of resources, prey and predators in these directions; subsequently, they move in the direction with the highest suitability score. The movements made by individuals determine their survival and foraging success in case of prey, and their prey capture rate in case of predators. These in turn affect the number of offspring produced. On an evolutionary timescale, successful individuals transmit their evaluation strategy to many offspring, subject to some mutation. Over the generations, successful strategies will spread, thus selecting evaluation mechanisms, that use proximate environmental cues to make adaptive decisions under the local circumstances. The evolution of the resulting movement strategies in turn has the potential to change the spatial pattern of resource densities and abundances and, hence, the nature of the trophic interactions.
We would like to stress that our model is not intended to mimic any particular biological system. Instead, our goal is to obtain conceptual insights into the coevolution of predators 1 3 and prey, when movement decisions are the target of selection and commonly made simplifications on the environment and trait interactions are relaxed. For this purpose, we have kept the other aspects of our model, including our assumptions on the evolving movement strategies, as simple as possible. With this conceptual model, we want to gain insight into the following questions: What kinds of movement strategies do evolve and which information sources are used by predators and prey? How do predator-prey interactions reflect and shape the resource landscape? Are the evolved populations monomorphic, or do different movement strategies coexist in the same population? How fast is evolutionary change in relation to ecological change; does the interplay of ecology and evolution lead to novel eco-evolutionary patterns and dynamics?

Ecological interactions
To aid intuition, consider the prey to be herbivores feeding on a resource, henceforth called grass, and predators feeding on herbivores. All individuals live in an environment consisting of a grid of 512*512 cells with wrap-around boundaries, such that individuals leaving the grid on one side reappear on the diametrically opposed side. A cell can host one or several herbivores and predators, but larger concentrations are unlikely as the number of grid cells is an order of magnitude larger than the typical number of individuals. Ecological interactions occur in discrete time steps. A time step (we think of a day) contains a grass growth phase, a movement phase and a foraging phase with predator-prey interactions. Grass grows at a constant rate of 0.01 per time step up to a maximum density of 1. Next, herbivores and predators move between cells based on their inherited movement strategy as described below. Herbivores visiting a cell deplete the grass and gain the corresponding amount of energy. Multiple herbivores occupying the same cell share the amount of grass on that cell. If a predator encounters a herbivore on the same cell, the predator succeeds to capture the herbivore with probability 0.5, in which case the herbivore is killed and the predator gains one unit of energy. If several predators co-occur with several herbivores (which only rarely ever happens), the successful predators kill and consume all the herbivores present in the cell. The killed herbivores are equally distributed between successful predators.

Movement strategies
Movement strategies are based on the evaluation of nearby cells (see Fig. 1). For each evaluated grid cell, a 'suitability score' S is calculated. S is the weighted sum S = w g G + w h H + w p P, where G, H and P are the grass density, herbivore density, and predator density in the cell, respectively. The weighing factors w g , w h , and w p are individual properties; they are genetically encoded and transmitted from parent to offspring.
Grass density represents the total amount of grass in a given cell, while for herbivore and predator densities we convoluted the otherwise discrete presence-absence values of agents via a Gaussian filter in a neighbourhood distance of one (Lindeberg 1994). This yields continuous values of herbivore and predator densities that are diffused around the actual positions of the agents, much like we would expect from olfactory cues or similar. An individual can thus sense the presence of other individuals even when they are outside, but close to, the individual's movement range. The density furthermore indicates cells where individuals could move to in the next timestep. Individuals evaluate all cells in their movement range, and move to the cell with the highest suitability. We apply a certain amount of noise to the calculated suitability scores. This adds the possibility for stochastic movements, particularly if weights are small. When the weights are large, this noise becomes trivial to the comparison between cells.
In the simulations shown, the movement range has a radius of one for the herbivores (9 cells, Fig. 1b), and a radius of two for the predators (25 cells). If predators have the same movement radius as their prey, herbivores can reliably escape predation by moving away from high predator densities. Only when predators can move further than herbivores do both parties need to predict the behaviour of the other party, and interactions become more intricate.

Evolution of the evaluation mechanism
We consider haploid parthenogenetic populations with discrete, non-overlapping generations. Population sizes vary throughout the simulation, as herbivore numbers are diminished by predation and herbivores and predators produce offspring in relation to consumed resources. Herbivores produce on average 0.1 offspring per unit of resource consumed. The a Individuals evaluate all cells in their movement range as to their 'suitability' and move to the cell of highest suitability. The suitability S of each cell is the weighted sum of the local grass density G, the local herbivore density H, and the predator density P, where the weighing factors w g , w h , and w p are genetically determined and hence evolvable. b Herbivores and predators move on the same rectangular grid, but their movement range per time step can be different. The plot illustrates the movement range of herbivores (blue, radius = 1) and predators (red, radius = 2) for our standard configuration 1 3 predators are supposed to be generalists; their expected number of offspring is a baseline value of 0.6, with 0.6 offspring added per prey caught. The realized number of offspring of an individual is drawn from a Poisson distribution based on the calculated expected value. At the start of a new generation, offspring are placed at random in the 'dispersal range' around the position of their parent. In the simulations shown below we used a dispersal radius of 1 for both species (9 cells, Fig. 1b). The above values for the food-to-offspring conversion rates and the dispersal range are used in our default scenario, but other parameter settings will be discussed as well.
Each individual has three gene loci with alleles w g , w h , and w p that correspond to the weights encoding the evaluation of environmental suitability (Fig. 1a) and, hence, determine the individual movement strategy. The movement strategy determines the types of habitat most likely visited and therefore the individual's intake rate and, in the case of the herbivore, the individual's probability of escaping predation. At the end of a generation (after 100 time steps), surviving individuals produce offspring depending on their total food intake. Each offspring inherits the genetic parameters of its parent, subject to rare mutations. A mutation occurs with probability µ = 0.001 per locus, in which case the original value is changed by an amount drawn from a Cauchy distribution with location 0 and scale parameter 0.001. At the beginning of the simulation, the weights of predators and prey were initialised with a draw from this distribution, implying that most weights started close to zero.

Results
We conducted many hundreds of long-term simulations that shared the characteristic features that we will demonstrate for one exemplary simulation run. Over a wide range of parameter combinations (see below), extinction of the predator or prey population is virtually assured. We therefore chose parameters that allow for extended coexistence between predators and prey, still permitting for extinction to occur as a result of the ecological and evolutionary dynamics. A video of such a simulation can be accessed under https:// youtu. be/ cLUCE x6Mlnk, and the program can be easily run locally by download from our digital resources or the github repository. Figure 2 shows a typical simulation run of our model. Over a period of 25,000 generations, the interacting movement strategies of predators and prey induce characteristic population dynamics (Fig. 2a) that appear regular for a while but then switch to a new dynamical state. The snapshots in Fig. 2b illustrate that spatial patterns underlie these states, which in turn reflect the dominant movement strategies in the predator and prey populations (Fig. 2c). It is apparent that the phase shifts in the population dynamics and spatial patterns reflect evolutionary changes in the movement strategies. Sometimes, these shifts occur within a few generations. Evolution is stochastic, highly dynamic and continually produces new strategies and ecological patterns.

Overview: transitions, pattern formation, and polymorphism
Between generation 35,000 and 38,000, the herbivore and predator populations show fairly regular oscillations (Fig. 3a). On the landscape level ( Fig. 2b.1), this manifests as the repeated expansion and depletion of herbivores and predators across the landscape. Figure 2c.1 illustrates the underlying movement strategies of predators and prey in generation 37,000. For ease of comparison, the weights of the three inputs were normalized by dividing each weight by the sum of the absolute values of all three weights. The purple and the cyan dots indicate that two movement strategies coexist in the prey population: both give a very low weight to grass density and a (different) negative weight to prey density. Hence, herbivores care surprisingly little about resource availability and move primarily to avoid conspecifics. The two morphs differ in the way predator density is considered. As indicated by the colours of the dots, prey morph 1 weighs predator density negatively (with a Fig. 2 Eco-evolutionary pattern formation. a Population size of predators (red) and herbivores (blue). b Landscape snapshots at three time points of the simulation (indicated by arrows in panel A). Grass density is shown in green, herbivore density in blue and predator density in red. Other colors emerge from additive color mixing, yellow for example signifies areas of high grass and predator density, purple that herbivores and predators occupy the same area. Black areas correspond to empty cells with low resources. c Movement strategies of herbivores (purple-cyan-blue) and predators (red-yellow-black) at each of the three snapshots, measured by the relative magnitude of each the three weights (absolute values sum to one). Shown is a subsample of 100 individuals per population, with jitter noise added around the true position (0.08 in x and y).The relative magnitude of the weight for predator density is color-coded for herbivores from −0.7 (= purple) to 0.0 (pink/cyan) and 0.7 (= blue), and for predators from −1.0 (= red) to 0.0 (= yellow/grey) and 0.2 (= black). In the panel of generation 37,000, the herbivore population consists of a predator-averse (-0.25, morph 1), and a predator-prone morph (+ 0.03, morph 2). Also the predator population consists a predatoraverse (−0.99, morph a), and a predator-prone morph (+ 0.15, morph b) relative weight of ~ 0.25), while prey morph 2 weighs predator density slightly positively (~ 0.03). In other words, morph 1 avoids areas with high predator density, while morph 2 is, surprisingly, attracted to such areas. We will elaborate on this counterintuitive finding and the dynamics between the morphs in the next section. Also the predator population is dimorphic: morph a (red dots) mainly considers (and strongly avoids) conspecifics (with a relative weight of ~ 0.99), while morph b (black dots) primarily moves according to large prey densities, while being split on its preference for grass, otherwise having a slight positive preference for conspecifics (~ 0.15).
After generation 38,000, the polymorphism in the herbivore population is lost, and a new predator morph emerges that has a strong preference for high grass densities (Fig. 2c,  generation 45,000). These predators aggregate on high grass patches and exclude herbivores from such areas. As a consequence, stable spatial patterns form in the distributions of predators, herbivores and resources, and the population dynamics are stabilized (Fig. 2b.2).

Fig. 3
Polymorphism and trait cycles. a Population size of predators (red) and herbivores (blue) between generations 35,000 and 35,500. The shown dynamic continues until generation ~ 38,000. b Relative frequency of the predator-avoiding morph 1 (pink) and the predator-prone morph 2 (cyan) in the herbivore population. c Relative frequency of the predator-avoiding (red) and the predator-prone (black) morph in the predator population The other predator morph strongly avoids conspecifics and roams also within the herbivore aggregations.
After generation 47,000, the spatial organisation is lost and population cycles emerge again. The population cycles are repeatedly interrupted by brief periods of stability (generation 52,300; 52,650; 53,500), and cycles differ in amplitude as a function of the involved movement types. In generation 52,600, we observe spiral wave patterns. Both populations are monomorphic at this point. Predators track low grass densities, allowing them to follow the herbivores, as these leave tracks of depleted grass in their way. The herbivores primarily avoid conspecifics, but also have a preference for predator densities and a small positive preference for grass. The reaction norms of these morphs are depicted in figure S3.
In summary, population dynamics and spatial patterns both reflect the evolved movement strategies of herbivores and predators. For example, dynamic spatial patterns such as spirals (Fig. 2b.3) are produced by predators chasing after herbivores, as these expand across the landscape. In contrast, spatial aggregations as in Fig. 2b.2 reflect a 'sit-and-wait' strategy of the predator, where individual predators prefer high grass densities and tend to remain there; the predators can thus monopolize high grass patches and consume those herbivores that are eventually lured in by the high grass density. The population dynamics of the entire simulation as well as the evolution of weights is shown in figures S1 and S2 of the supplementary material. For a more dynamical depiction of the spatial dynamics, we refer the reader to the animation videos and executable of our model in the digital supplement.

Polymorphism and trait cycles
We now take a closer look at the oscillations that dominate the simulation between generations 35,000 and 38,000, and the polymorphisms that occur during this phase. The densities of both populations are subject to oscillations with a period length of ca. 8 generations and a quarter phase lag between herbivores and predators (Fig. 3a). As we have seen in Fig. 2c.1, both populations are polymorphic at this point: some herbivores are 'predatoraverse' in that they avoid high predator densities (morph 1), while other herbivores are 'predator-prone' in that they have a slight positive preference for high predator densities (morph 2). Figure 3b shows how the relative frequencies of the two morphs fluctuate over a period of 500 generations (pink: predator-averse morph 1; cyan: predator-prone morph 2). Likewise, Fig. 3c shows the fluctuations in the two predator morphs: the predator-averse morph a that strongly avoids locations with high predator densities and the predator-prone morph b that has a weak positive preference for high predator densities (and, in addition, is strongly attracted by high prey densities).
The four morphs exhibit stochastic but regular oscillations in relative frequency, with a period of ca. 50 generations. Each of the two morphs in the two populations is adapted to a morph in the other population, and increases when this morph is common: When the predator-averse herbivore morph 1 is common, the predator-averse predator morph b increases, followed by the increase of the predator-prone herbivore morph 2, which when common induces the increase of the predator-prone predator morph a. The predator-averse herbivore morph 1, which avoids high predator densities, is adapted to escape predation from predators with a positive preference for predator densities (morph b), but suffers from predation by the predator-averse predator morph a. In turn, morph a does not efficiently capture predator-prone herbivores, because it moves away from the patches with high predator density that are attractive for the herbivore morph 2. Thus, even though the behaviour of the predator-prone herbivore morph 2 seems counterintuitive, it is an effective adaptation against predators within the framework of our model, where individuals need to predict each other's location in the next timestep.
In addition to the large-scale oscillations, the morph frequencies show epicycles influenced by the oscillating population dynamics. These do not show a consistent pattern and emerge from the repeated expansions and contractions of both populations across the landscape, as well as stochastic effects in the spatial distributions of individuals.

Evolutionary transitions and their effect on ecological patterns
We now take a closer look at the evolutionary transitions from one kind of pattern to another one. For this purpose, we will focus on a 400-generation segment between generation 52,400 and 52,800 (which includes snapshot 3 of Fig. 2). In this period, the simulation undergoes three shifts. At first, predator and prey population sizes are stable, and the landscape configuration is stable as well, with herbivores forming loose aggregations ( Fig. 4a.1). Around generation 52,450, the situation changes, the spatial structure is lost Fig. 4 Effect of mutations on evolutionary transitions. a Population size of predators (red) and herbivores (blue) between generations 52,400 and 52,800. (B, C, D) Evolution of the three weighing factors determining the movement strategy in the predator population: b weighing factor w G (= the weight given to grass density); c weighing factor w H (= the weight given to herbivore density); d weighing factor w P (= the weight given to predator density). In most of the generations shown, the predator population is polymorphic for one or more weighing factors. The relative frequencies of the coexisting trait values within each generation are encoded by a color gradient from 0.0 (= white) to 0.3 (= red) and 1.0 (= blue). Weight values are shown on a tanh-transformed scale. (1) and (2) show snapshots of the landscape snapshots at generation 52,400 and 52,800 and predator and prey populations begin to oscillate. After generation 52,600, a rapid shift occurs that stabilizes the population dynamics again, only to fall back into oscillations after generation 52,700, with a now higher amplitude than before. Indeed, both the predator and prey populations come close to extinction during this phase.
The phase transitions can be traced back to mutations of movement strategies in the predator and prey population. The rapid shift in generation 52,610 for example is caused by a mutation in the herbivore population, changing predator preference from positive to strongly negative. The other two shifts are induced by mutational changes in the predator population. Figure 4 shows the distribution of all three weights in the predator population during this time. The first shift (in generation 52,445) occurs when two new predator mutations appear that both have a negative grass weight (i.e., an aversion of high grass densities). During the following oscillations, the previous two grass-prone morphs disappear. In addition, the other two weights of the movement strategy of the predator shift to values close to zero, implying that in this time period the predators are only guided by their aversion for high grass densities. The predator population then remains largely static for ~ 100 generations, after which the above-mentioned adaptation in the prey population occurs, which leads to counteradaptations in the predator population: First, the predator density weight mutates back to be negative, then the grass weight mutates to be near-neutral, and finally the herbivore density weight mutates to a more positive value, thereby producing again population oscillations and spatial dynamics (Fig. 4.2). Thus, from generation 52,610 to generation 52,700, three mutations occur in the predator population that ultimately lead to the phase shift.
Several things should be noted here: First, evolutionary changes occur on a similar time scale as ecological changes, and the occurrence and spread of new movement strategies shapes the ecological dynamics, which in turn determine the success of the different strategies. Ecological and evolutionary processes are thus strongly intertwined. Second, the predator population displays a high level of polymorphism during this period. The existence of these polymorphisms allows for subsequent adaptive change. Third, phase shifts can require several antecedent mutations, that can then lead to an abrupt change of the ecological dynamics. And finally, evolution can both be stabilizing and destabilizing, and adaptations in the predator population not necessarily lead to an improvement of population-level fitness.

Model sensitivity and parameter settings
We have deliberately focussed on a single simulation run in order to explore the dynamics in great detail. In the presented simulation, dispersal occurs locally in a range of one around the parental individual. Herbivores produce offspring at a conversion rate of 0.1. Predators produce offspring at a conversion rate of 0.6, and have a baseline food intake of 0.6. For the same parameter settings, replicate simulations can strongly differ from each other ( fig. S7 in the supplementary material). However, the basal elements discussed above (periods of stasis followed by periods of oscillations, and vice versa; abrupt transitions; polymorphism) were observed in all simulation runs based on our default parameter values. The duration of the various phases can be very different, with some replicates maintaining oscillations or remaining static for long periods of time. Extinction occurs frequently in our simulations, with four simulations out of ten running beyond generation 100,000, and one going extinct around generation 10,000 (figure S7).

3
Different parameter settings have a predictable effect on the eco-evolutionary patterns generated by our simulations. For example, offspring dispersal has a strong influence on the simulations. By default, we considered local offspring dispersal (dispersal radius = 1). When offspring are dispersed more widely (dispersal radius = 10), phase shifts still occur and spatial patterns emerge, but extinction occurs more regularly than before ( fig S8 in the  supplementary material). When offspring are distributed randomly across the landscape, extinction occurs within few generations. The spatial structure emerging from local reproduction thus promotes stability.
The conversion rates of herbivores and predators influence the stability and overall level of the population dynamics. If the conversion rate of herbivores is increased from 0.1 to 0.2, their average population size increases, but extinction becomes much more frequent, with few replicates reaching generation 10,000 ( fig S9 in the supplementary material). If the conversion rate or the baseline food intake of the predator is reduced, the overall abundance of predators decreases and the ratio of prey to predators becomes larger, oscillations dampen, and spatial patterns like spiral waves or rapid expansions vanish. As a consequence, the feedback between landscape structure and movement strategy evolution is diminished and the simulations tend to produce stable ecological and evolutionary dynamics, although polymorphism remains common (figure S10 in the supplement). In the digital supplement, we provide an executable of our model that allows the interested reader to explore other parameter settings such as lower conversion rates, random movement or more complex movement strategies controlled by recursive networks.

Discussion
We introduced a new method to model evolvable movement strategies and applied this method to the antagonistic coevolution of movement decisions in predators and prey. The movement strategies are characterized by heritable parameters that determine how an individual evaluates the environmental cues that determine its movement decisions. Movement is thus not guided directly by fitness expectations (as in Iwasa 1982;Abrams 2007), but by an inherited mechanism that evaluates the 'suitability' of the available options on the basis of proximate environmental information and then bases the next move on the comparison of suitabilities. Suitability judgements are not necessarily aligned with fitness expectations (which may not be well-defined in a highly dynamic setting), but they are 'adaptive' in the sense that they are shaped by natural selection. As movement affects foraging success and predator-prey encounter rates, the evaluation mechanism is an important determinant of lifetime reproductive success and hence subject to natural selection.
In our model of antagonistic coevolution, selection pressures vary strongly in space and time, leading to rapid evolution and rich spatial dynamics. Although our model is still very simple, we observe a range of phenomena that do not occur in models with coarser spatial scales and with fewer evolutionary degrees of freedom. The populations in our model do not evolve stable movement patterns but instead exhibit intricate ecological and evolutionary dynamics, including regular frequency cycles between movement strategies and the spontaneous advent of novel strategies and counter-strategies. Regularly, qualitatively different movement strategies coexist as polymorphisms for extended periods of time. For example, we observed the evolutionary emergence of sit-and-wait predators, while other predators were chasing their prey (Fig. 2c.2). The movement strategies in herbivores and predators determine the spatial pattern of resource depletion and predator-prey encounters.
Over the generations, these patterns are fluent, as they change with the evolution of the underlying movement strategies. Rapid transitions between patterns (e.g. between static spots and spiral waves) can occur; these reflect coevolutionary changes rather than changes in external conditions. We will now discuss some of these findings in more detail and in the context of other models of predator-prey coevolution.
Previous work on the coevolution of predator-prey movement has produced interesting spatial dynamics, such as evolutionarily optimal strategies leading to 'predator-prey chases' across habitats (Abrams 2007) and the coupling of ecological dynamics across habitats via evolving conditional strategies (Flaxman et al. 2011). However, in these studies 'movement' corresponds to a choice between a small number of densely populated habitats. Within habitat patches, the interactions of predators and their prey is not modelled at the individual level but by patch-level dynamic equations (but see Patin et al. 2020). In contrast to these habitat choice models, we consider movement in a fine-grained spatial environment. From ecological models, it is known that in such fine-grained environments the interplay of diffusion and predator-prey interactions can induce a diversity of spatial patterns, including rotating spirals and static stripes or spots, depending on the parameters of the ecological interaction (Hassell et al. 1991;Comins & Hassell 1996;Alonso et al. 2002;Banerjee 2015). Hence, patterns may shift abruptly when ecological parameters change, either externally as in the publications above, or intrinsically through evolutionary processes as in our model. Repeated switches between patterns can thus occur even though the extrinsic parameters have not changed. Further, evolution of movement strategies can produce qualitatively novel patterns unknown from ecological models that assume random dispersal. These eco-evolutionary patterns are a product of the feedback between the movement-generated distribution of individuals in space and the evolution of movement strategies based on the local conditions that individuals experience.
Behavioural polymorphisms commonly occur in our simulations and can either be persistent (Fig. 3) or fleeting (Fig. 4); either fluctuate or remain stable over time. Polymorphisms are also predicted by analytical models of coevolution (Senthilnathan & Gavrilets 2021), but they occur much more frequently in mechanistic individual-based models (Botero et al. 2010;Long & Weissing 2020). The behavioural dimorphisms emerging in the herbivore and predator population (Fig. 3) show that systematic behavioural variation can emerge in two species due to coevolution. The existence of behavioural polymorphism can have important ecological and evolutionary implications (Sih et al. 2012, Wolf & Weissing 2010. In predator-prey systems, intraspecific variation can, for example, be crucial for the sustained persistence of both species (Senthilnathan & Gavrilets 2021). Evolution in our model is either dominated by the advent of novel mutations, or by frequency-dependent oscillations between different morphs present in the population. The latter produce trait cycles reminiscent of population genetics models (Nuismer et al. 2005;Kopp & Gavrilets 2006;Cortez & Weitz 2014), where oscillations occur between the frequencies of different alleles present in the population. This is exactly the pattern we see in Fig. 3, where two morphs of predators and prey oscillate in frequency, with each morph being adapted to one specific morph of the other species. We thus recover the matching alleles assumption of population genetics models, but without assuming this interaction a priori. Instead it naturally emerges from the ecological interactions of our model. Mutation-limited evolution is much more erratic, as it depends on the stochastic occurrence of one or several sequential mutations (Fig. 4). The newly arising morphs frequently induce phase shifts that transform the ecological dynamics. We did not observe trait cycles as described in quantitative genetics (Gavrilets 1997;Mougi & Iwasa 2010;Cortez 2018) or adaptive dynamics approaches (Dieckmann et al. 1995;Marrow et al. 1996), where oscillations occur due to systematic shifts in a continuous trait that either corresponds to the mean value of a normal distribution (quantitative genetics) or to the trait value of a monomorphic resident population (adaptive dynamics). In line with these modelling frameworks, we also assume that mutational step sizes are typically small. However, the Cauchy distribution of mutational step sizes allows for rare mutations of large effect, which may explain why continuous variation plays a less prominent role in the evolutionary processes observed here (see Wolf et al. 2008). We kept our model as simple as possible, in order to demonstrate that not much structure is required for obtaining the ecological and evolutionary patterns described above. It would be interesting to study the implications of features such as sexual reproduction and different modes of inheritance, or a spatially heterogeneous resource distribution, but this is beyond the scope of our study. Here we only discuss the implementation of movement strategies by three weighing factors (w g , w h , and w p ) in our model. This way, there are three 'evolutionary degrees of freedom' in our model, giving larger scope to nonequilibrium dynamics than traditional single-trait models (Leimar, 2009;McNamara & Weissing 2010;Débarre et al. 2014). While previous work on coevolution in multidimensional phenotype space is still framed in terms of simple phenotype matching rules (Gilliam et al. 1987;Débarre et al. 2014), the interactions in our model are mediated by movement strategies. How these strategies interact is, however, an emergent property of the model. The nature of this interaction not only depends on the evolved strategies in conspecifics and antagonists, but also on the local environmental conditions, in which individuals encounter each other. Trade-offs for herbivores between resource acquisition and predation avoidance come about naturally in this case and do not need to be assumed ad hoc. The emergent interactions between individuals allows for much more rapid and unpredictable evolution due to the emergence of 'surprising' (and sometimes counterintuitive) strategies and counter-strategies, the advent of novel forms of behaviour, and by allowing for behaviour that is (at least partly) stochastic and unpredictable. We anticipate that, quite generally, models with more evolutionary degrees have much richer eco-evolutionary dynamics than most conventional models.
Having said this, we are fully aware that our behavioural model is still unrealistically simple and well-behaved. Our model could be extended quite naturally by exploring different modes of inheritance, and by basing the calculation of environmental suitability on a more complex algorithm, such as an evolving artificial neural network (ANN) (Huse et al. 1999;Enquist & Ghirlanda 2005;Morales et al. 2005). Evolving regulatory networks (including ANNs) have a number of important features (Wagner 1979, Van den Berg & Weissing 2015, Van Gestel & Weissing 2016, such as the emergence of cryptic variation (since the same phenotypic strategy can be encoded by very different networks), which allows for much faster evolution in the face of environmental change. Perhaps most importantly, network models tend to have an intricate genotype-phenotype map, implying that small mutations can have large and unexpected implications at the phenotypic level. It is therefore not surprising that ANN-based pilot simulations on predator-prey coevolution exhibit even richer eco-evolutionary dynamics than occurring in the present model.