Spatial Pattern Switching Enables Cyclic Evolution in Spatial Epidemics

Infectious diseases often spread as spatial epidemic outbreak waves. A number of model studies have shown that such spatial pattern formation can have important consequences for the evolution of pathogens. Here, we show that such spatial patterns can cause cyclic evolutionary dynamics in selection for the length of the infectious period. The necessary reversal in the direction of selection is enabled by a qualitative change in the spatial pattern from epidemic waves to irregular local outbreaks. The spatial patterns are an emergent property of the epidemic system, and they are robust against changes in specific model assumptions. Our results indicate that emergent spatial patterns can act as a rich source for complexity in pathogen evolution.


Introduction
Recent studies show that in spatial models, evolutionary dynamics of infectious diseases change with respect to predictions from non-spatial model, which implicitly assume complete mixing of individuals [1][2][3][4][5]. Spatial epidemic waves have been observed for several diseases [6][7][8][9]. Starting with Rand et al. [10], a number of authors have analyzed a spatial pathogen-host model in which a host population with local reproduction is infected by a lethal pathogen that is transmitted through direct local contact [4,[11][12][13]. In this class of models, local extinction of host populations caused by the pathogen is balanced by re-colonization of empty space by reproduction of uninfected hosts. Evolution towards increased transmissibility is limited, because pathogens that are too infectious exhaust their local host population before enough new hosts are born for the pathogen to persist.
Spatial pattern formation can also act to limit the duration of the infectious period. We recently examined a spatial epidemic model in which infection of hosts leads to waning immunity instead of host death [14], and we found that, after the system selforganizes into epidemic waves, natural selection is directed toward increasing outbreak frequency. Outbreak frequency is optimal for infections of relatively short duration so that pathogens evolve towards short lasting infections. The mechanistic explanation for competition for frequency between waves is that, when two infection waves collide, this typically is followed by local extinction of the pathogen. Subsequently, the pathogen with the higher outbreak frequency will be the first to reinvade the susceptible host population. Frequency selection also occurs in spatial competition in chemical reactions [15], autocatalytic hypercycles [16] and parasitoid-host systems [1,17].
Our previous analysis [14] adopted a standard evolutionary approach in that ecological dynamics were operating on a faster timescale than mutation and evolution. However, for many pathogens, in particular viruses and bacteria, this assumption need not be true, as even within a single host the pathogen can accumulate many mutants as a by-product of the large mutation rate, large copy numbers, and short generation time [18,19]. In this paper we investigate the effects of increasing mutation rate in our spatial disease model (see Model section for a description of the model). It will turn out that co-occurrence of pathogens with a substantial difference in infection period can trigger a novel type of spatial patterns, which can temporarily reverse the selection pressure in the system. As a result, the model exhibits so-called evolutionary cycling in the length of the infection period.

Model
In our spatial epidemic model, hosts are situated in a regular square contact network. Hosts can be in three states: susceptible (S), infected (I), or resistant to infection (R). There are three possible state transitions, namely infection, acquiring resistance, and loss of resistance. Infection is a local contact process, where each infected host can infect its eight direct neighbors, and this occurs at rate b. Acquiring resistance occurs after a fixed time t i since infection (i.e. t i can be interpreted as the duration of the infectious period). Loss of resistance occurs after a fixed duration t r . Note that only infection is a spatial process. Because the total population size of hosts is constant, infection is both density and frequency dependent. All time units are relative to the resistant period (which is set at t r = 1). In this paper, we investigate evolution of the length of the infectious period. Upon each new infection, the infectious period can, with equal probability, increase or decrease with a fixed mutation size step 6Dt i with mutation probability m. As a consequence, strains with different infectious period will co-occur and compete for susceptible hosts.
We assume full cross-resistance, and no co-infection or superinfection between strains, so that after infection by a particular strain a host cannot be infected by any other strain.
In our earlier paper [14], we used a probabilistic cellular automaton model with a small fixed time step. We tested decreasing the time step until the behavior approached a continuous time process. However, such a method requires quite small time steps, which implies a large waste of computational power, as most time is used for updating hosts that do nothing within the time step. Instead, here we adopt a much faster continuous time method by using an event based updating procedure. For each infected individual, we schedule a next infection event, using an exponential distribution with rate parameter b. A scheduled infection event will only cause a new infection if, at the time of the event, the individual is still infectious, and if a randomly chosen neighbor individual is susceptible to infection. Acquiring resistance and loss of resistance of individuals are also scheduled as events. All events are efficiently managed in a binary heap data structure [20]. The next scheduled event is always located at the top of the binary tree, and updating the tree is an efficient process (scaling with the logarithm of the number of events). The method is vey flexible, as e.g. other probability distributions for the infectious period, such as lognormal or gamma distributions, can be implemented without loss of computational speed. Furthermore, using larger neighborhoods for the infection process does not substantially increase the computation time, as the neighborhood is only evaluated at the time of a potential infection. A non-spatial analogue of the model can be investigated by implementing infection at a global scale, that is, an infected individual can infect any other individual in the system with equal probability, thus effectively removing the spatial component from the model. A deterministic continuous approximation of the mean field dynamics [21] would consist of a set of three (SIR type) delay differential equations. Typical single run simulations of the spatial model, as reported in this paper, take around one hour to run on a Pentium computer. For large simulations several nodes of the Dutch national computer cluster LISA were used overnight. A C-code for the updating procedure can be obtained from the corresponding author upon request.

Rock-Scissor-Paper Dynamics
In Fig. 1A through 1C (for movies see supplemental Videos S1, S2, S3), we investigate pair-wise competition between three pathogen genotypes that differ in infection period, with infection periods of t I = 0.4, t I = 0.6 and t I = 0.8, respectively. We observe that these three pathogen types have a cyclical dominance structure, where each infection period can outcompete exactly one other infection period, resulting in a so-called ''rock-scissorpaper'' dynamics. In each pair-wise competition two distinct genotypes are initialized each on one side of a 2016100 field. Infection and resistance is randomly initialized in large blocks, in order to speed up the pattern formation. The middle vertical row is kept susceptible for the first 10 time units, so that the patterns can fully develop before the competition is started. Fig. 1A, t I = 0.4 against t I = 0.6, and Fig. 1B, t I = 0.6 against t I = 0.8, both show that the genotype with the shorter infection period wins the competition. This outcome is surprising, as a shorter infection period is a disadvantage at an individual level, as it generates fewer secondary infections. This selection for shorter infection period was also reported in our previous paper [14], and it can be explained by selection for higher outbreak frequency. At the spatial location where the two genotypes co-occur, the disease locally goes extinct after each wave of infection. The genotype that subsequently can cause the next disease outbreak faster will reinvade the area more quickly and it will gradually increase its domain of dominance. This selection for shorter infection period depends on the spatial pattern formation, and particularly on the epidemic waves that cause local extinction of the disease. In the non-spatial analogue of our model a shorter infection period is always a competitive disadvantage, as it locally generates less secondary infections per infected individual (i.e. it has a reduced reproduction number R 0 ; see ref. [22]). The spatial patterns cancel this selection at the local individual level, because after each epidemic wave both disease genotypes locally will go extinct.
However, in Fig. 1C, for a two-fold difference in infection period of t I = 0.4 against t I = 0.8, the genotype with the longer infection period wins. This reversal of selection is accompanied by a local change in the spatial pattern; that is, at the interface between the two genotypes the dynamics now switches to a finescaled irregular outbreak pattern instead of regular epidemic waves. Here, the epidemic waves from the fast genotype run into remaining resistance fragments of the slow genotype. As a consequence, the waves break up and become irregular, and, most importantly, the disease no longer locally goes extinct (Fig. 1D). In the absence of local extinction, selection will favor the disease genotype with the longer infection period, because it locally generates more secondary infections (i.e. it has a larger R 0 , see ref. [20]).
Summing up, we have two opposite directions of selection in our system, depending on the spatial pattern at the competition interface. When the difference between infection periods is relatively small, the spatial pattern consists of epidemic waves, and both genotypes will go extinct in between successive epidemic waves. In this spatial regime, selection is for increased outbreak frequency, because the higher frequency genotype will be able to reinvade the area more quickly. In contrast, when the difference between infection periods is large, the spatial pattern at the competition interface will be irregular (i.e. turbulent), and local extinction of both genotypes between successive epidemic waves does no longer occur. In this spatial regime, selection is for increased infection period, because the genotype with longer infection period will be able to locally cause more secondary infections. In a way, the situation is reminiscent of a life history trade-off, or in particular ''r'' versus ''K'' selection theory [23,24]. Here, the epidemic waves create an unstable ''r'' strategy environment, selecting for fast regeneration and growth, whereas the local irregular outbreak pattern creates a (relatively) stable ''K'' strategy environment, selecting for increased local competition strength. In the next paragraph we will quantify the exact

Author Summary
Parasites are commonly believed to evolve to make as many infections as possible. In large scale simulations of disease spread, however, natural selection can instead act to maximize outbreak frequency. Here, pathogens that cause short infections can be rewarded for their prudence by a rapid subsequent outbreak. Very surprisingly, in a simple spatial model for disease spread so-called evolutionary cycling can be observed, that is, the infection period indefinitely keeps evolving up and down. The reversal in the direction of selection is triggered by a switch in the spatial patterns from epidemic waves to irregular local outbreaks. This finding offers an alternative mechanism to explain why infectious diseases show so much variation in the lengths of their infectious periods.
boundaries for these two modes of selection, depending on the infection periods of both competitors.

Evolutionary Cycling
In Fig. 2A, a pair-wise competition plot is shown for infection periods up to t I = 1.0. For small differences in infection period, there is selection for increased outbreak frequency, leading to an evolutionary stable attractor (ESS) around t i = 0.2 (i.e. in the ESS, R 0 = t I b = 6.4). However, if the difference in infection period is large, the fine-scaled local irregular outbreak pattern develops, and selection favors the genotype with the longer infection period. As a consequence, for a small mutation rate of m = 0.01, an evolving  population will converge to the ESS, because at any time the difference between competing genotypes will be small (Fig. 2B, dotted line). However, for a larger mutation rate of m = 0.05, the evolutionary dynamics are qualitatively different, and they converge to large amplitude cycling of the infection period (Fig. 2B, solid line). In the spatial pattern of the cyclic evolutionary attractor (Fig. 3A) it appears that distinct areas differ substantially in infection period. In Fig. 3B, after 20 time units, the areas of irregular outbreaks have invaded the nearby epidemic wave areas, whereas the other areas have decreased in infection period (see supplemental Video S4). In Fig. 3C the local direction of selection for this 20 time unit period is plotted, showing a strong dependence of the direction of selection on the local spatial pattern. In areas of irregular outbreaks the infection period increases whereas in areas of epidemic waves it decreases over time.

Phase Transition and Hysteresis
When the mutation rate is slowly increased, the evolutionary dynamics show a sudden shift from the ESS to the cyclic evolutionary attractor (Fig. 4A). For increasing mutation rate, the distribution of mutants around the ESS gradually widens (due to quasispecies dynamics, see ref. [19]), until the first local irregular outbreaks can develop when neighboring genotypes differ enough in infection period. Hereafter, the dynamics quickly shifts to the cyclic attractor. If we now gradually reduce the mutation rate (Fig. 4B), it turns out that the cyclic evolutionary attractor can be maintained for quite small mutation rates, before the system falls back to the ESS. There is thus a considerable region, for 0.027#m#0.036, where the system displays bi-stability between the ESS and the cyclic evolutionary attractor. This bi-stability can be explained by the fact that the cyclic evolutionary attractor reinforces itself, as it is associated with large differences in infection period across the field which induces the irregular outbreak pattern. It should be noted that, due to the stochastic nature of the model, the hysteresis is only a transient phenomenon; that is, allowing long enough simulation time, within the bi-stable region the dynamics will spend time in both alternative attractors. Increasing the system size will, within the bi-stable region, promote the cyclic evolutionary attractor, as the first irregular outbreaks can develop anywhere in the spatial domain, and increasing the size of the domain will thus increase the chance of the irregular outbreak pattern to develop. Once the irregular outbreak pattern has developed somewhere in the field, it will expand to the rest of the system. Also, on a small spatial domain, the cyclic evolutionary attractor is more easily lost, because persistence of the attractor depends on different spatial regions being in different phases of the attractor. On a small field global synchronization occurs more easily, and this will result in the system falling back to the ESS (i.e. t i = 0.2) attractor.

How Robust is the Evolutionary Cycling?
In the previous section, the first local irregular outbreak pattern is generated by increasing the mutation rate. One could argue that the mutation rate that is necessary to generate this pattern is quite large, and hence the evolutionary cycling might be considered unlikely to play a role in real epidemics. However, it turns out that the irregular outbreak pattern can also originate from various other sources. For instance, cyclic evolutionary dynamics can be observed for a small mutation rate of m = 0.001 (or even smaller) if only a small region of the field has some variance in the length of the resistant period. In Fig. 5 such a small inhomogeneous area is introduced in the middle of the field. Within this inhomogeneous region, the disease will develop a small scale spatial pattern in which there is no local extinction. As a consequence, within this region, the disease will evolve to maximum length of the infection period (Fig. 5, dotted line). Subsequently, the long infection periods in middle of the field will seed the irregular outbreak pattern and the resulting evolutionary cyclic dynamics in the rest of the field (Fig. 5, solid line). The mutation rate in this case will set the timescale of the evolutionary cycle, but even at very small mutation rates the evolutionary cycling dynamics persist. The evolutionary cycling here is much more regular than in Fig. 2B, because now the maximal infection period remains continuously present in the system, whereas before it occasionally was lost and had to reemerge from de novo mutation.
The cyclic evolutionary attractor can also be induced for small mutation rates if some local movement of individuals is included in the model, or if occasional large effect mutations are considered (results not shown). Furthermore, the occurrence of the cyclic evolutionary attractor does not depend on the specific assumptions of the model. We imposed a maximum to the length of the infection period, but such an upper limit could also be obtained implicitly by e.g. implementing a trade-off between the length of the infection period and the infectivity of a genotype. We extensively tested robustness of our results against various model assumptions, for instance by changing the contact network topology (using 4 or 6 infection neighbors), changing the boundary conditions (reflecting boundaries and empty boundary conditions), and testing various statistical distributions for the duration of the infectious period and the resistant period (e.g. log-normal and gamma distributions). The key property that is necessary for the selection for shorter infection periods to occur is recurrent local extinction of the disease. This coincides with a local dynamics that is unstable (i.e. the mean field dynamics converges to large amplitude oscillations or even extinction), inducing the spatiotemporal pattern of epidemic waves. If instead, the local (and mean field) dynamics of the disease converges to a stable endemic equilibrium, selection will favor genotypes with maximal infection period.

Discussion
In the spatial epidemic model we present here, the local oscillating dynamics give rise to epidemic waves and recurrent local extinction of the disease. These spatial patterns enable selection for a reduced infection period, as this will increase the outbreak frequency. However, the system can switch to a local fine-scaled irregular outbreak pattern, which emerges at the interface between pathogens with large enough difference in infection period. Here, local extinction of the disease no longer occurs, and selection is in the direction of increased infection period. These opposite directions of selection can lead to large amplitude cyclic evolutionary dynamics.
Cyclic evolution was reported before for non-spatial systems [25][26][27]. It can result from co-evolution between species or traits. In contrast, in the model we present here, we observe cyclic dynamics in a single evolving trait, namely the length of the infection period. In our system, a change in the spatial pattern acts as a switch for the direction of selection. Interestingly, both alternative spatial patterns cause selection in a direction that promotes the switch to the other pattern, resulting in a continual cycling between the two patterns and directions of selection. Alternative spatial patterns can also act to induce bi-stability, in the case that each of the two patterns causes a direction of selection that enforces the current pattern [4,17]. Note that the observed cyclic evolutionary dynamics cannot be explained with currently popular spatial analytic tools, such as spatial adaptive dynamics [28], spatial moment approximation techniques [29] or spatial inclusive fitness measurements [30]. The fine-scaled irregular outbreak pattern is only a transient pattern that occurs at the interface between genotypes, and the local selection always favors the disease genotype with longer infection period. However, as we demonstrated, this local selection can be overruled by local extinction and recolonization, which can cause selection to act in an opposite direction. In this paper, we use numerical computational methods, such as constructing the pair-wise competition plot of Fig. 2A. We feel that the current emphasis on spatial approximation techniques acts to underestimate the potential (non-linear) effects of spatial patterns on disease dynamics. In particular, spatial patterns can undergo sudden changes, such as the transition from epidemic waves to local irregular outbreaks that is observed in this paper, and such transitions are often hard (or even impossible) to predict with the current approximation techniques. Maybe these analytic tools can be improved to include, or at least predict, such transitions; for instance by developing 'early warning signals' for these shifts in spatial pattern [31], but this is not easy. In the meantime, we want to advocate the complementary use of numerical computational methods that incorporate the full non-linearity of the system, as such methods can increase both our qualitative and quantitative understanding of the impact of spatial pattern formation on disease dynamics. Spatial epidemic waves have been reported for many infectious diseases [6][7][8][9], and recurrent local extinction is occurring for many epidemic diseases [32]. We have shown that these spatio-temporal patterns can have profound effects on selection for disease properties. Most notably, epidemic waves and recurrent local extinction can induce selection for a short infection period, which is a property of many diseases, and which is otherwise hard to explain without invoking strong trade-off assumptions [33][34][35]. We want to emphasize that the spatial patterns can induce selection for properties that emerge at a scale beyond that of the individual or the direct neighborhood of individuals [36]. Selection for increased outbreak frequency in epidemic waves is an intriguing example where a property that appears at a scale of outbreak centers overrides the local individual selection for maximizing secondary infections.

Supporting Information
Video S1 This video shows competition between two pathogen types which differ in infection period. Infection period 0.4 (green) is winning against period 0.6 (yellow) because it causes a higher outbreak frequency. The video corresponds to Fig. 1A  Video S2 This video shows competition between infection period 0.6 (yellow) and period 0.8 (orange). Again, the pathogen with the shorter infection period and higher outbreak frequency wins. The video corresponds to Fig. 1B  Video S3 This video shows competition between infection period 0.8 (orange) and period 0.4 (green). Here, turbulence develops at the interface between the two types, and the type with the longer infection period slowly wins. Note that green is increasing initially, before the turbulence has developed. The video corresponds to Fig. 1C  Video S4 This video shows cyclic evolution in the length of the infection period. The direction of selection depends on the local spatial pattern. In some areas the infection period is increasing and in other areas it is decreasing, depending on the local spatial pattern. The video corresponds to the period between Figs. 3A and 3B. Found at: doi:10.1371/journal.pcbi.1001030.s004 (6.00 MB AVI)