A damage spreading transition in a stochastic host–pathogen system

One of the leading proposals for solving the biodiversity problem is the Janzen–Connell hypothesis, suggesting that the abundance of a species is limited by a host-specific exploiter. Motivated by this model, here we analyze a spatially explicit host–pathogen system, looking for coexistence conditions under stochastic dynamics. Above the standard extinction transition associated with the failure of the pathogen to invade, we report another, damage spreading transition, marking the point where macroscopic clusters of host individuals disappear. Beyond its practical significance, this transition is apparently a generic landmark along the axis of decreasing stochasticity, if the deterministic dynamics support cycles or quasicycles.


Introduction
Seeking to understand the forces that shape biodiversity is one of the most important issues of modern science [1], both for practical reasons (such as attempts to conserve various species of fauna and flora) and as a fundamental theoretical issue. Within this framework, one particularly disturbing problem is the co-occurrence of many competing species in a relatively small area, despite the competitive exclusion principle [2][3][4][5].
When many species are competing for a single limiting resource, the theory predicts that only one of them, the fittest, will survive. In general, models show that the number of coexisting species is bounded from above by the number of limiting resources (niches). For example, niche partitioning explains quite nicely the coexistence of a few bird species, each having a different beak size and different diet. However, in some cases this explanation appears to be unsatisfactory [6]. In particular, it is quite difficult to explain the coexistence of several hundred tree species in a km 2 of tropical forest, or the large number of plankton species that coexist within small regions of a freshwater lake, as reflecting the superiority of each species in its unique niche [6,7].
One popular solution to this apparent paradox was suggested years ago by Janzen [8] and Connell [9] (JC). These authors hypothesized the existence of host-specific predators or pathogens for every tree species, say, in a tropical forest. Such a pathogen makes the areas directly surrounding the parent tree inhospitable for seedlings/trees of the same species. In the absence of this pathogen, the fittest species will win the Darwinian competition game, selecting out all other taxa. The pathogen dictates a sparse matrix of the fittest trees in the forest, leaving room for the second-best competitor with its own parasite, and so on.
Empirical support for the JC hypothesis comes from studies demonstrating that recruitment patters are skewed away from reproductive adult trees [10][11][12], and from those showing the density dependence of the parasite/herbivore species [13]. However, we are not familiar with studies of the persistence of such a spatially explicit host-pathogen system, identifying the conditions under which both the victim and the exploiter survive in the long run.
Here we consider a toy model for the JC mechanism or, in general, for a birth-death-infection process in a spatial domain. It turns out that, beyond its relevance to the biodiversity problem, our numerical study reveals an interesting feature of the stochastic dynamics: the appearance of a damage spreading transition (DST) [14,15]. In terms of stochastic processes and infection models we interpret this DST as a milestone along the path connecting the strongly stochastic regime (where history-to-history variations are large and the use of partial differential equations is insufficient) to the quasideterministic limit.

The SIB model
We have simulated the dynamics of a d-dimensional square lattice (d = 1, 2), where each cell is in one out of three states. A S cell contains a susceptible tree, an I cell is occupied with an infected tree and B denotes a 'bare soil' cell. An infected tree may either die (the cell becomes bare soil) or infect one of its nearest-neighbor susceptible cells. A susceptible cell may recruit one of its nearest-neighbor B's. The process rates are These rates may take any positive values. One of them may be used in order to set the overall timescale, so the system is characterized by two parameters. In the well-mixed (zerodimensional) limit, the corresponding deterministic differential equations are In the limit k 3 = 0 (no recruitment) this process is simply the Kermack-McKendrick susceptible-infected-recovered (SIR) epidemic model [16]. The other extreme, k 3 → ∞ (i.e. every bare soil site instantly becomes a tree) corresponds to the susceptible-infected-susceptible (SIS) dynamics [17]. The stochastic dynamics of each models in a spatial domain supports an extinction transition (ET) at some critical value of infection rate k 1c . This transition is in the directed percolation (DP) universality class for SIS [18] and in the dynamic percolation universality class for SIR [19]. Away from the transition point, the description of the system by deterministic rate equations becomes more and more accurate as k 1 increases, in a manner analyzed recently in [20,21].
Since the sum of S, I and B sites is equal to the total volume of the system N , B = N − S − I and the mean-field equations can be reduced to a two-dimensional (2D) system. The well mixed system admits a coexistence fixed point at S * = k 2 /k 1 and I * = (N − S * )(k 3 /[k 1 + k 3 ]). This is the only stable solution as long as k 1 > k 2 (or, equivalently, as long as the primary reproductive number R 0 = k 1 /k 2 is greater than one).
If the community dynamics of forest trees is subject to the JC mechanism, there is one superior species that will take over the system if all pathogens are removed. A customary method (see e.g. [3]) to model this predominant species dynamics is to allow its seedlings to replace all other residents of the forest. Accordingly, in our model (equation (1)) S and I corresponds, respectively, to healthy and infected individuals of the superior species, while a 'bare soil' site (B) may contain an individual of any other species, as all of them may be invaded by the offspring of the superior.
Here we utilize this approach to study the host-pathogen dynamics of the superior species only, neglecting the complex dynamics of the other species in the community. As we shall see, the conditions under which this top-level host-pathogen system is persistent are not trivial. In future works we intend to iterate this method and to introduce more species (and their pathogens) in descending order of competitive strength, to model the whole community dynamics.

The replica method: application to the susceptible-infected-susceptible process
In what follows we will study the synchronization process of a system and its replica when both are subject to an identical microscopic noise. In this section we introduce this method and apply it to the contact (SIS) process in one dimension (1D). The contrast between the replica synchronization for SIS (considered here) and for SIB (in the next section) emphasizes the surprising emergence of a new critical point, associated with the DST.
The SIS process, simulated here, is the limit k 3 → ∞ of equation (1). Every site is either susceptible or infected. In any elementary step an infected site is chosen at random and recovered (becomes S) with probability 1/(1 + λ), where λ = k 1 /(k 1 + k 2 ). With probability λ/(1 + λ) it attempts to infect one of its (randomly chosen) nearest neighbors, so if this neighbor was susceptible, it becomes infected too, and the process iterates. This process is also known as the contact process [18].
When such a stochastic process is simulated in a spatial domain, the results differ from the mean-field predictions. The phase transition point is above its mean-field value k 1 = k 2 (λ = 1/2), as the clustering of infected sites lowers the effective infection rate. The equilibrium density of infected sites is zero below λ c 3.297 [18]. This ET belongs to the DP universality class [18], so the spatial correlation length ξ ⊥ (measured by monitoring the number of infection events per site, when the infection is ignited with a single seed [20]) and the mean survival time of the infection, T (the average lifetime of a single seed ignited infection) both diverge at the transition. Since the critical exponents of DP are known with a very high degree of accuracy Above the transition, the average survival time and spatial extant of an infection are both infinite. To measure the correlation length and the lifetime of a fluctuation, we adopt a damage spreading approach. First, the simulation runs for a long time, until a 'typical' configuration is obtained. Then, the configuration is replicated, and a minimal excitation is introduced (one susceptible site becomes infected). From now on, the two replica dynamics are simulated under the same microscopic noise: a site is picked at random in replica A, and the same site is picked in replica B. If the selected site is, say, infected, a new random number is drawn in order to decide if this site will recover or will try to infect the site to its left or to its right, and the same decision is made for both replicas. Replica dynamics differ from each other only if a process on replica A cannot take place on replica B or vice versa. For example, when a selected infected site tries to infect the site on its right, it may do so in replica A (e.g.) when the target site is susceptible, but not in replica B where the same site is already infected. In such cases, the 'distance' between replicas (the number of sites in different states) decreases by one unit after this step. Once the distance reaches zero (at which the two replicas are synchronized) it will stay zero forever.
Using this technique, we can monitor the lifetime T (average time until synchronization) and the spatial extent ξ ⊥ (calculated using the number of infection events per site, only counting infections that happen in replica B and not in A) of an excitation. In fact, the standard seed method-tracking a process that starts with a single seed (below the transition) until it dies-is a special case of the damage spreading technique, since below the transition the typical state is simply the empty state. As seen in figures 1(a) and (b), the ξ ⊥ and T diverge above the transition with the same critical exponents, and the curves of ξ a reasonable numerical accuracy, to the same critical value of λ. The noise is slightly larger above the transition, as the preparation of a 'typical' state becomes more difficult close to the transition.
Far above the transition, for large values of λ, the lifetime and the spatial extent of the damage decay monotonically. In our model, λ dictates not only the infection rate but also the diffusion constant of the epidemic. Accordingly, for large values of λ, the system approaches its deterministic limit, as it is made of 'effective patches', each containing a large number of sites. In the deterministic limit, the SIS process satisfies the Fisher equation [20], and admits a uniform steady state solution where the amplitude of a damage is healed like t −1 .

Damage spreading transition: the SIB process
What happens when the same damage spreading technique is applied to the SIB system? Strictly speaking, unlike the SIS model, SIB supports two absorbing states, when all sites are S and when all sites are B. Thus it does not satisfy the conditions of the Grassberger-Jansen conjecture [18]. However, starting with a susceptible system and adding one infected tree, the probability of  DST ET Figure 2. The ET and the DST are indicated by dotted lines for two-dimensional, discrete time SIB dynamics. The ET is shifted to the lower value of λ = 2.1 for q = 0.9. As seen in the main panel, the lifetime of an excitation diverges at the ET, drops above it, and then rises again and diverges (for an infinite system) at the DST (λ 3). For a finite system, above the DST the lifetime of an excitation diverges like exp(−κ(λ)N ). A plot of κ versus λ above the transition appears in inset (a), and points to the DST point where κ → 0. Inset (b) shows Z = T −1/1.455 versus λ below the DST transition, where the lifetime of an excitation is finite, suggesting that the DST belongs to the DP class, in agreement with [23].
reaching the pure B state decreases exponentially with the size of the system, and may be neglected in the thermodynamic limit. Accordingly, one would expect a DP ET at some λ c [22], as indeed was observed in simulations, see figure 2. Note that the results presented in figure 2 are for a 2D system. In 1D the infected fronts isolate the bare soil from the susceptible trees and (at least with nearest neighbor dynamics) the infection fails to become endemic for any finite λ.
The SIB process simulation takes place as follows. In an elementary step, a site is chosen at random; if it is a susceptible site (a healthy tree) it tries to recruit one of its nearest neighbors with probability q (the results depicted in figure 2 were obtained with q = 0.9), and the target site becomes S if it was in B state. When the site chosen is in I state, the procedure is identical to the one described in the last section for SIS.
Surprisingly, in the SIB model the damage lifetime is not a monotonically decreasing function of λ above the transition. Instead, as demonstrated in figure 2, the healing time diverges at the ET transition point, starts to drop with λ, but then rises and diverges again at λ DST , the point marking what we label a DST. Above this point there is a finite chance that the two replicas never synchronize, i.e. that the damage never heals. An analysis of small samples (figure 2, inset) shows that below the transition the lifetime of an excitation (time to synchronization) approaches a finite value, independent of the sample size, while above the transition this lifetime grows exponentially with system size.
This DST transition appears to mark a milestone along the road to the deterministic limit. As explained, an increase in λ is equivalent to an increase in the diffusion constant and hence to a decrease in the relative importance of the fluctuations associated with demographic noise. Indeed, for large values of λ our simulations (not presented here) show well-defined, circular propagating fronts of infection, as observed in the deterministic model. Like fire fronts, these propagating waves of infection separate the inner bare soil region from the outer healthy forest. Thus, an elementary excitation in a replica (a healthy tree in the bare soil zone) may be the originator of a process that never synchronizes.
Another aspect of the DST may be understood by looking at the steady state properties of a finite size system. If a system supports an attractive periodic manifold (e.g. a limit cycle) in the deterministic limit, one can easily understand the appearance of the DST as reflecting a zero mode associated with the phase along the cycle. Demographic stochasticity and the discreteness of the reactants not only allow for healing but also provide a mechanism for damage spreading. Accordingly, as the system approaches its deterministic limit it has to be above the damage spreading transition. Conversely, right above the ET the effect of stochasticity is strong and the 'phase' does not manifest itself in the dynamics, leaving a region between the ET and the DST where the damage heals.
The SIB system is a little bit more subtle. Although its deterministic equations support an attractive fixed point, this point is a focus, and under demographic noise it becomes unstable [24], displaying McKane-Newman quasicycles [25,26]. The amplitude of these quasicycles is on the order of √ N , where N is the number of sites in a fully connected system, so their relative width (e.g. S/N , the fraction of sites in the susceptible state) shrinks to zero in the deterministic limit N → ∞. In the coarse grained sense the deterministic system supports an attractive fixed point, with no phase. However, the fluctuations δS are proportional to the square root of N , so as long as the individuals are O(1), cycles, phase and the DST do exist in the deterministic limit.
The damage spreading transition was suggested, originally, as an indication for the existence of different dynamical phases, such as those characterizing glasses. Later studies revealed that the DST point is non-universal, and that two physically equivalent microscopic dynamics, such as those found in a 'heat bath' versus Metropolis method simulation, may yield different results [27]. It turns out, however, that in our system the damage spreading transition is followed by a significant change in the spatial organization of the system, as depicted in figure 3. In the region between the ET and the DST the system supports a macroscopic cluster of susceptible (healthy) trees, with cluster size that scales with the volume of the system N , while all other clusters yield statistics with exponential cutoffs. At the transition point the macroscopic cluster disappears, and the cluster statistics fit a power-law. Above the DST there is no longer a macroscopic cluster, but the cluster statistics again admit an exponential cutoff.
This phenomena-the disappearance of the infinite cluster and the power-law statistics at the transition-is well known from the theory of percolation transitions. Here, however, the measured critical exponents differ significantly from those of 2D percolation. For example, the exponent that characterizes the cluster statistics power-law at the transition is 1.77, while percolation theory predicts a much higher value, τ = 187/91 = 2.05. The mass of the infinite cluster grows like (λ DST − λ) β , whereas for 2D percolation β = 5/36. Instead, we have measured β 0.22.

Discussion
While (perhaps) non-universal, for any specific microscopic dynamics of a real host-pathogen process the DST point denotes a shift in the spatial deployment of the forest. Apparently, this shift is required if the JC mechanism is to work. Below the DST, when there is an infinite cluster of healthy superior trees, the clusters of the second-best species will be small and isolated, and since this second species has its own pathogen, its population will be extinction-prone. It appears that the JC mechanism may support a fair amount of species richness only if the single-species host-pathogen systems are above the DST. This dictates some constraints on the cluster statistics and dynamics in the forest, constraints that may be tested against empirical datasets like those analyzed in [28].
To conclude, we report here on a damage spreading transition in a quite generic victim-host model that supports quasicycles in its deterministic limit. Our results suggest that this transition indicates that the phase, and not only the amplitude, starts to play a significant role in the dynamics, a necessary step on the way toward the deterministic limit. The JC mechanism is apparently applicable only above the transition, when the infinite cluster disappears.