Coexistence in an Inhomogeneous Environment

We examine the two-dimensional extension of the model of Kessler and Sander of competition between two species identical except for dispersion rates. In this class of models, the spatial inhomogeneity of reproduction rates gives rise to an implicit cost of dispersal, due to the tendency to leave favorable locations. Then, as in the Hamilton-May model with its explicit dispersal cost, the tradeoff between dispersal case and the beneficial role of dispersal in limiting fluctuations, leads to an advantage of one dispersal rate over another, and the eventual extinction of the disadvantaged species. In two dimensions we find that while the competition leads to the elimination of one species at high and low population density, at intermediate densities the two species can coexist essentially indefinitely. This is a new phenomenon not present in either the one-dimensional form of the Kessler-Sander model nor in the totally connected Hamilton-May model, and points to the importance of geometry in the question of dispersal.


Introduction
A central question in ecology is how species determine their migration rate. Migration entails a number of possible benefits to the population, among them reducing overcrowding areas and the discovery of new resources to exploit. There are also significant costs associated with migration. Especially in patchy environments, there may be a significant risk of death (due to predation or other exogenous causes) in trying to cross over unsuitable areas to find a colonizable site. Presumably, migration rates are under evolutionary control, and so can in principle be adjusted to some optimal level. There have been a large number of studies in this field, both from an analytical and a numerical perspective.
The current work is an extension of the work of Kessler and Sander [1] reconciling two seemingly contradictory approaches to this question of optimal dispersal. On the one hand, there has been extensive study of the rate equations approach to modeling. This approach had its beginnings in a result of Hastings [2], who showed both in a two patch system and in a continuous onedimensional system that the species with the slower dispersal rate is stable with respect to the introduction of a faster dispersing species, always driving the faster species to extinction, as long as the two patches were not identical. This is due to the fact that the overall effect of dispersal was to take individuals from favorable areas to less suitable places, leaving the slow species in control of the favorable site. The higher production rate of progeny at the favorable site eventually leads to a takeover of the poorer site as well. Thus, dispersal in a inhomogeneous environment carries with it an implicit cost. This finding was reinforced by Dockery and coworkers [3]. Examining the continuous one-dimensional version of the Hastings model, they studied the reaction-diffusion system: Here, f and s represent the local population density of the fast (f) and slow (s) migrating species, with D f wD s . The inhomogeneous distribution of resources is modeled by a spatially varying growth rate, a(x), common to the two species. The competition is of logistic type, with each species competing equally against both fellow members of the same species and members of the other variety. The Hastings result is that the pure slow state with f~0 is linearly stable, and the pure fast state, s~0 is unstable. Dockery, et al. were able to prove the stronger result that, independent of the exact form of a(x) and of the initial conditions, the pure slow state is the only global attractor and the faster species always went to extinction. The implication of both these works is that any mutation resulting in a lower level of migration should fix in the population, driving the migration rate eventually to zero.
This strong result is in sharp contrast to what occurs in the individual based models originated by Hamilton and May [4]. Looking at a system with an infinite site of patches fully connected by dispersal, and with an explicit cost of dispersal, they showed that the evolutionary stable strategy was a large level of dispersal. This optimal level of dispersal approached unity as the dispersal cost was lowered to zero, and surprisingly, remained finite even when the dispersal cost approached 100% of the dispersers.
Kessler and Sander [1] reconciled these results by showing that the demographic fluctuations inherent in the individual based model, and absent in the rate equation approach, were responsible for the very different results of the two models. By adding a noise term to the differential equation, Eq. (1), they showed that fluctuations have a negative effect and that dispersal, by reducing such fluctuations by uniformizing the population in space, has a beneficial effect. Thus, the tradeoff between the dispersal cost and the fluctuation reduction benefit gives rises to a nonzero optimal value of dispersal rate. This result underlies the decrease of the optimal dispersal rate with increasing population density seen in the Hamilton-May-Comins [5] extension of their model to arbitrary population density, since increasing population density decreases fluctuations and so reduces the beneficial aspects of dispersal.
To verify their findings, Kessler and Sander studied an individual based version of the one-dimensional Dockery model. The system has two species of individuals, fast and slow, residing on sites of a one-dimensional lattice x i~0 ,1, . . . L with no-flux boundary conditions. The system is designed to reproduce a spatially discretized version of the Dockery model in the deterministic limit. All individuals have a probability of a(x i )dt of reproducing in some small time interval dt, depending on their location x i . Furthermore, individuals at x i have a concentrationdependent probability of death, (f i zs i )bdt, where f i and s i are the numbers of fast and slow individuals sharing the particular site. They also move to a nearest-neighbor site with probability D f ,s dt, depending on their identity. Kessler and Sander chose to examine a distribution of resources corresponding to two ''oases'' of food separated by a desert: They set b~1 N0 , so that in the absence of diffusion, the total population at each site would saturate at the local carrying capacity N 0 a(x i ). The results of this model showed that, for a given level of inhomogeneity g, when N 0 was large, the rate equations correctly predict that the slower species wins. On the other hand, for N 0 small, the ability of the faster species to use diffusion to minimize fluctuations led to its victory, against the predictions of the rate equations. One surprise was that the boundary between large and small N 0 was surprisingly large, in the range of hundreds per lattice site, against the naive expectation that the rate equations would be reliable for N 0 's of order 10 or so. The interplay between inhomogeneity and N 0 was consistent with the observation that in the Hamilton-May-Comins [5] extension of the model, wherein an arbitrary carrying capacity is allowed, the   evolutionary stable dispersal rate decreases with N 0 , and with increasing dispersal cost.
The model we examine herein is the natural two-dimensional generalization of the Kessler-Sander model. The system consists of a square two-dimensional lattice, x i ,y i~0 ,1, . . . L. The dynamics is exactly the same as in the one dimensional case, with the dispersal move now occurring to one of the four nearest-neighbor sites. We choose for the inhomogeneity a ''checkerboard'' pattern with oases in the corners and the center and deserts in-between, see Fig. 1. In all our simulations, we have set L~50.

Results
Figs. 2, 3, 4 display some snapshots from the two dimensional simulation. The color represents the relative abundance of fast and slow on a site, with all fast being pure red and all slow being pure blue. Fig. 2 displays a system with the parameters: g~0:4,N 0~4 00. Fig. 2(a) displays the initial state of the system for all the runs, namely only slows in the left half of the system and only fasts in the right half. In Fig. 2(b)-Fig. 2(e) it is possible to see Coexistence in an Inhomogeneous Environment PLOS ONE | www.plosone.org that initially the fasts predominate in the deserts and the slows in the oases. Eventually, for this set of parameters, with its relatively low population density, the fast dispersing species encroach on the oases and drive the slows toward extinction, as seen in n Fig. 2(f). Fig. 3 displays a system with the same degree of inhomogeneity, g~0:4 but a larger population density N 0~7 00. The system starts as in Fig. 2(a). Again, the fasts fare better in the deserts, but the eventual advantage accrues to the slows, as expected due to the higher population density. By Fig. 2(d), after 4500 steps, the victory of the slows is essentially complete, with the slows in the ''deserts'' making a valiant but futile last stand, in accord with the statement that the slows win at sufficiently large N 0 .
The above results are consistent with the picture that emerged from Kessler and Sander's simulations of the one-dimensional model. However, we show in Fig. 4 snapshots of a system with the same inhomogeneity g~0:4 and an intermediate value of the population density, N 0~6 00. Again, the system starts as in Fig. 2(a). In this figure it is clearly seen that despite the much longer run times compared with the above two cases, there is essentially no change in the ratio between the fasts and slows after 15000 time steps. This behavior leads to the speculation that the state displayed in this figure will be preserved also at infinite time, with the system having settled into a state of stable coexistence. Fig. 5 display the percentage of fast individuals in the entire system as a function of time for the three cases presented in the snapshots. For g~0:4, N 0~4 00 the percentage of fasts rises quickly to 100%. While g~0:4, N 0~7 00 the percentage of fasts decreases to zero. It should be remarked that in this case the victory of the slows in clear from the spatial snapshot long before the total number of fasts has dropped significantly. For g~0:4, N 0~6 00 the percentage of fasts fluctuates around 67.2%, with no discernible secular trend.
To investigate further, we scan the g, N 0 parameter space, the results of which are summarized in Fig. 6. The blue color represents parameters in which the slows win, the red color represents parameters in which the fasts win. The green color represents parameters with apparent coexistence. The overall trend of the dependence on g and N 0 are as in the one dimensional model with higher inhomogeneity and higher N 0 both favoring the slows. The new feature is of course the intermediate coexistence phase, which is seen to shrink in extent as g increases, at least for g §0: 4.
In order to clarify if the apparent coexistence phase is real or only the result of the finite running time of the simulation, we measured the time until victory as the coexistence region is approached from either above (high N 0 ) or below. Figs. 7 and 8 display the time (averaged over five independent runs) until one species wins as a function of N 0 . Fig. 7 displays the average time until the fasts win in the small N 0 regime and Fig. 8 displays the average time until the slows win in the higher N 0 regime. The time till extinction of the loser grows markedly as N 0 approaches the coexistence regime, and is consistent with a power-law divergence: Likewise, the time to extinction of the losing fast species grows as N 0 approaches the coexistence area from above, again consistent with a power-law divergence: The parameters of Eqs. (4,5) were found by least-squares fitting. The graphs lead to the conclusion that even after an infinite time (strictly speaking, only for an infinitely large system; see below) there will be no victory of either population because the winning times diverge as the coexistence region boundary is crossed.
It is interesting to note that the exponents characterizing the divergence are different on the two sides of the transition. In addition, the exponents decrease with increasing g. Graphs of the exponents A f , A s as a function of g are presented in Fig. 9. The parameters and their standard deviation are displayed in Table 1, 2. To point out the contrast between the coexistence dynamics we see in the two dimensional system and the behavior of the one dimensional simulation, we have measured the average extinction time as a function of N 0 for a fixed g~0:6 in the one dimensional system. The data is presented in Fig. 10. Here the extinction time   Table 2. Parameters and their error ranges for the time to extinction of the losing slow species. rises as the border between the fast dominated and the slow dominated phases is approached, but it shows no sign of diverging.
Up till now, we have focused only on one specific pattern of inhomogeneity in the two dimensional system, namely the ''checkerboard''. We now investigate a ''striped'' pattern of inhomogeneity, which is a uniform extension in the second dimension of the one dimensional a(x) (Fig. 11): We present snapshots from the simulation in Fig. 12 for g~0:4, N 0~6 00, which showed coexistence in the checkerboard geometry. As always, initially the fasts dominate the ''deserts'' and the slows the ''oases'', but as in the checkerboard geometry, the system settles into a coexistence state. Thus the coexistence is tied to the presence of spatial inhomogeneity, but the form of the inhomogeneity in evidently less crucial.
As in the checkerboard patterned system, there is a correlation between the local growth rate and the fraction of fasts in this coexistence state. In the striped system, it is easier to quantify this behavior, by averaging over the direction transverse to the stripes. A we see in Fig. 13, the correlation is surprisingly weak, with only a *15% variation between the fast/slow balance between the deserts and the oases.
The presence of coexistence in the striped 50|50 system is a bit nonintuitive, given its absence in the one-dimension system, which is just a extremely ''thin'' (50|1 version of the same system. The explanation is that the critical parameter underlying coexistence must be the width. We show in Fig. 14 the average time to victory for systems of width 1, 10 and 50. We see that as the width increases, the peak moves toward lower N 0 and increases in height. The maximum time to victory for the 50|10 system is 29000 steps, which is much larger than that of the 50|1 system, but nevertheless is finite. We conjecture that this should be the case for any finite width system, with the maximum time to victory diverging with the system size. For the 50|50 system, as was the case for the same sized checkerboard system, such a maximum time is evidently much larger than we can afford to simulate.

Summary and Discussion
In conclusion, we have shown that spatial inhomogeneity in a two-dimensional system leads to the coexistence of two different dispersal strategies, for some intermediate range of densities. This feature is absent in both one-dimensional systems and in fully-   connected systems. The transition to coexistence is marked by the power-law divergence of the extinction time of the underperforming strategy, as the carrying capacity is varied. More study is required to clarify the nature of this coexistence phase, and the critical behavior associated with the transition.
The phenomenon of coexistence in an inhomogeneous environment is not surprising from the perspective of niche theory. The ''oases'', with their high reproduction rate, are evidently favorable for the slows and the ''deserts'' for the fast dispersers. However, given the density-dependent nature of the effective fitness, this is a simplistic view. In one dimension, one species may exhibit a temporary advantage, however in the fullness of time, one species always drives its competitor to extinction. Evidently, however, in the two dimensional environment, this is no longer true near the region where the advantage changes between species in a one dimensional patterned environment. The question of coexistence evidently requires a more refined understanding.
Even with a standard density-independent fitness landscape, the interplay of dispersal and niches raises interesting questions worthy of further study. Dispersal will at a minimum tend to wash out the niche boundaries, with some degree of mixing at the niche boundaries. The question that arises is whether there is some degree of dispersal that wipes out completely the underlying niche structure of the environment and leads to domination by a single ''compromise'' strategy. In our system, the fairly small difference between the population balance in the ''oases'' and the ''deserts'' means that dispersal effects are strong enough to wash out most, but not all, of the differences between the niches. Given the absence of a clear boundary between the two niches, it is surprising that in the two dimensional system, coexistence is still possible. Comparing the present model to one with density-independent fitness differences should be very instructive.