Advent of extreme events in predator populations

We study the dynamics of a ring of patches with vegetation–prey–predator populations, coupled through interactions of the Lotka–Volterra type. We find that the system yields aperiodic, recurrent and rare explosive bursts of predator density in a few isolated spatial patches from time to time. Further, the global predator biomass also exhibits sudden uncorrelated occurrences of large deviations from the mean as the coupled system evolves. The maximum value of the predator population in a patch, as well as the maximum value of the predator biomass, increases with coupling strength. These trends are further corroborated by fits to Generalized Extreme Value distributions, where the location and scale factor of the distribution increases markedly with coupling strength, indicating the crucial role of coupling interactions in the generation of extreme events. These results indicate how occurrences of extremely large predator populations can emerge in coupled population dynamics, and in a more general context they suggest a generic class of deterministic nonlinear systems that can naturally exhibit extreme events

www.nature.com/scientificreports/ trophic system composed of interacting vegetation (denoted by u), prey (denoted by v) and predator (denoted by w) populations are given by the following functions f (u, v, w), g (u, v, w) and h (u, v, w): Here the consumer-resources and predator-prey interactions are modelled using the Holling type II term and the Lotka-Volterra term, described by functions f 1 (u, v) and f 2 (v, w) respectively. The Holling type II functional form f 1 (u, v) = uv/(1 + ku) represents the coupling of vegetation (u) and prey (v) with coupling strength α 1 . The Lotka-Volterra functionial form f 2 (v, w) = vw , represents the coupling of prey (v) and predator (w) with coupling strength α 2 . The parameters a, b and c represent intrinsic growth rates of trophic species u, v and w. Further the model allows the predator w to maintain a low equilibrium level w = w ⋆ even when the population of its usual prey is low. This situation arises when alternative food sources are available for the predator (such as red squirrels for the Canadian lynx in this case) and is mathematically modelled by linearizing the predator growth rate about the equilibrium level as (w − w ⋆ ) . In this work we chose the parameters to be a = 1.0 , b = 1.0 , c = 10.0 , w ⋆ = 0.006 , α 1 = 0.5 , α 2 = 1.0 and k = 0.05 , consistent with empirical data as reported in Ref. 15 . Now we expand the scope of the model above to mimic a collection of such population patches, consisting of vegetation, preys and predators. The populations in a local patch interact with other nearby patches in such a way that the predator of one patch can attack prey in neighbouring patches, namely the patches are coupled via cross-predation. Figure 1 shows a schematic diagram of the interaction of a patch with nearby patches. Specifically, we consider a ring of such patches, with each patch indexed spatially by i, i = 1, . . . N , where N is the total number of patches in the ring. The populations of vegetation, prey and predator in patch i is represented by u i , v i and w i respectively.
The form of the predator-prey interaction between nearest neighbouring patches is of the Lotka-Volterra type and is given by the following set of dynamical equations: The coupling constant C reflects the strength of interaction among patches, and in this work we focus on this crucial parameter. Note that in formal terms this constitues a conjugately-coupled system 16 .
So the scenario this model captures is the following: predators are usually more mobile and can reach other patches to find prey. So the predators of neighbouring patches influence the population of prey, as the prey (for instance, the snowshoe hare in this particular example) can be devoured by predators (the Candian lynx in this example) of neighbouring patches. As a consequence the population of predators depends on the population of prey of the neighbouring patches, as they can make forays into neighbouring areas to consume prey. This scenario is modelled by the inter-patch coupling in the standard Lotka-Volterra form.

Spatial distribution and temporal evolution of the population densities in the network
The first quantity of relevance in this system is the local population density in the patches, and their temporal fluctuations. We will focus on the deviation of the local population densities from the mean value, i.e., we will look for the emergence of extreme events in the local patches, as evident in explosive population densities at specific spatial locations. Figure 2 shows the time evolution of the vegetation, prey and predator populations of one representative patch from the network at high coupling strength C = 1.0 . We observe that, while the population densities are mostly confined to low values, the evolution is punctuated by sudden boosts to very high values. For instance, local predator population densities can shoot up more than 10 standard deviations away from the mean value. This is evident in the lower right panel of Fig. 2, where one can see instances where w exceeds the 10σ threshold. The instants at which these large fluctuations occur are relatively rare and completely uncorrelated in time and space. These large deviations are also clearly evident through the space-time plot of the evolution of predator populations in different patches in the network, shown in Fig. 3, where the extreme values of predator populations are visible as sparse bright randomly located dots in the figure.
Note that if consider a lower threshold (e.g. µ + 4σ ), the vegetation, prey and predator all exhibit events that exceed the threshold, with prey and predator populations yielding above-threshold events even when uncoupled. However significantly, such occurrences in uncoupled patches are strictly periodic, and stem from the intrinsic pulse-like solutions of the constituent patches.
In Fig. 4 we display the spatial distributions of the predator populations w i for the patches i = 1, . . . N at a representative instant of time, for the case of uncoupled patches and coupled patches. It is clearly seen that there are a few patches in the coupled case that grow explosively and have a predator population much larger than the mean value. These results qualitatively suggest that the coupling of population patches give rise to extreme predator populations at a small number of spatial locations, analogous to an extreme event in space. Such catastrophic events are rare in space, and occurs only at a couple of sites, but signal significant damage as they may entail serious control costs.

Global maximum of predator populations
Now we quantitatively estimate the maximum vegetation, prey and predator population densities in a patch, denoted by u max , v max and w max respectively, attained in the course of the network dynamics 17,18 . We estimate this by finding the global maximum of u i , v i and w i , for i = 1, . . . N (where N is sufficiently large), sampled over a time interval T (where T is much longer than the intrinsic oscillation period). This will help us gauge the www.nature.com/scientificreports/ magnitude of the extreme event and its relation to coupling strengths. Here we present results with T = 50 , with no loss of generality. Figure 5 shows the maximum u max , v max and w max , for a wide range of coupling strengths of the population patches. In the figure we depict the scaled values of the maxima, where u max , v max and w max are divided by their values in the uncoupled case. These scaled quantities help us assess the increase in the maxima in coupled networks, compared to that obtained in uncoupled patches, allowing us to specifically gauge the coupling-induced effects on the emergent maximum population densities. It is evident from the simulation results that the magnitude of the maximum vegetation and prey populations (i.e. u max and v max denoted by the red and green curves) do not increase much as coupling strength increases. However, the maximum predator population increases very significantly with coupling strength, with w max at C = 1 exceeding over six-fold the value obtained for uncoupled patches.
Further, we also explored different α 2 values in Eq. (1), reflecting different intra-patch predation rates, where the dynamics of the uncoupled system exhibited very low peaks in the population densities, specifically an order of magintude smaller than that occurring in the uncoupled system displayed in Fig. 2. Interestingly, there too we observed clear evidence of extreme events in predator populations arising under coupling. This suggests that the presence of significant peaks in the coupled system is not necessary for the advent of extreme events.
Additionally note that increasing the strength of the predator-prey interaction in a single patch does not give rise to extreme events in the population dynamics. Rather, as this coupling strength increases, the maximum of the predator population decreases and the dynamics becomes more regular. This lends further support to the importance of inter-patch coupling in the generation of large deviations in predator populations.

temporal evolution of biomass
The next quantity of interest is the total biomass of the vegetation, prey and predators, denoted by b u , b v and b w respectively. This represents a collective dynamical quantity, and is given at an instant of time t as follows: www.nature.com/scientificreports/ where N is the system size. Again, we will examine the presence of large excursions from mean-values, as such explosive growth indicate extreme events in time of a global quantity. Figure 6 shows the biomass of the vegetation ( b u ), prey ( b v ) and predators ( b w ). The uncoupled case is shown alongside as a reference. It is clear that when the patches are uncoupled, b u , b v and b w do not experience any large fluctuations. Further, the biomass of vegetation and prey for the coupled case also stays bounded within the 4σ threshold. However, interestingly, the predator biomass in coupled patches occasionally builds up to extreme values, crossing the 4σ threshold. This is also corroborated through a comparison of the maximum values of the biomass of the vegetation, prey and predator in the coupled network vis-a-vis the uncoupled patches. The maximum prey biomass is almost unchanged on coupling, and the maximum vegetation biomass in a coupled network exhibits less than 20% change from uncoupled values. On the other hand, the maximum predator biomass in a coupled network exhibits a three-fold increase compared to uncoupled patches (cf. Fig. 7). So coupling has a very significant effect on the predator biomass, and one finds clear evidence of the emergence of extreme events in time for this global quantity.
In order to ascertain that the extreme values are uncorrelated and aperiodic we examine the time intervals between successive extreme events in the predator biomass evolution. Figure 8 shows the return map of the intervals between extreme events and it is clearly shows no regularity. The probability distribution of the intervals is also Poisson distributed and so the extreme predator population buildups are uncorrelated aperiodic events.
Also note that spatial heterogeneity plays a role in the emergence of extreme events. If one considers a homogeneous system, the biomass does not exhibit any large deviations from the mean. Only a system comprised of a spatially inhomogeneous distribution of vegetation, predator and prey, yields extreme values of the global biomass. This suggests the importance of spatial heterogeneity in the sudden emergence of large biomass from time to time.

Generalized extreme value distribution
Lastly we analyse the distribution of the maximum size of the predator populations, as well as the maximum size of the global predator biomass. In probability theory and statistics, the generalized extreme value (GEV) distribution is a family of continuous probability distributions developed within extreme value theory, and the GEV distribution is often used as an approximation to model the maxima of long finite sequences of random variables 19 . So in order to obtain a quantitative measure of the extreme events generated for different coupling strengths, we fit these probability distributions to the GEV distribution, given by:  www.nature.com/scientificreports/ Here ζ is the shape parameter, µ is the location parameter and σ is the scale parameter. Note that the scale parameter is the most relevant parameter in this case as it determines the spread of the distribution (see Fig. 9). So larger scale parameters indicate larger deviations from the mean. Figure 10 shows the location and scale parameters obtained by best-fit to the Generalized Extreme Value distribution of the maximum predator population in a patch w max , and the maximum predator biomass b w . Clearly the location and scale parameters increase monotonically with coupling strength, for coupling strengths higher than ∼ 0.4 , with the rise being approximately linear at high C. Increasing location parameters indicate that the average predator population in a patch and the average predator biomass increase almost linearly with coupling strength. Increasing scale parameters suggest that the distribution becomes increasingly spread out, and the tail  www.nature.com/scientificreports/ of the distribution extends to larger values. So more extreme predator populations can be expected to occur in the patches, from time to time, when the coupling between the patches is stronger.

Effect of increasing range of coupling
Finally, we have also explored the effect of increasing spatial range of interactions, with predator-prey interactions occuring over 2k neighbouring patches, given by the following set of dynamical equations: Here k reflects the spatial range of the coupling, and determines the number of neighbouring patches over which the prey-predator interactions occur. Specifically the coupling range 2k = k ′ , where k ′ is the degree of the underlying ring connection network. The case of nearest neighbour interactions considered earlier is the limiting case of this model with k = 1 . Also note that for a homogeneous system, the inter-patch coupling term is the same for all k, as the system reduces to a set of identical population patches given by dynamical equations where the second and third equations in Eqns. 1 are simply augmented by the product of variables v and w weighted by the coupling constant C.
(4)  www.nature.com/scientificreports/ Figure 11 shows the time evolution of the vegetation, prey and predator populations of one representative patch from the network at high coupling strength C = 1.0 for different spatial ranges of coupling interactions. It is evident that while the population densities are mostly confined to low values, the evolution of the predator population is punctuated by sudden boosts to more than 10 standard deviations away from the mean value, for the case of coupling to 4 nearest neighbour patches (cf. lower left panel of Fig. 11). Interesting however, these large population buildups are less intense when coupling occurs over a larger spatial range (cf. right panel of Fig. 11). So, while increasing coupling strength monotonically increases the occurence of extreme events (cf. Fig. 5), large deviations are more enhanced when the spatial range of coupling is smaller, i.e. strong interactions over shorter ranges give rise to more extreme deviations from mean values. Figure 12 shows the maximum u max , v max and w max , with respect to coupling range, varying from nearest neighbour connections to near global coupling. As before, in the figure we depict the scaled values of the maxima, where u max , v max and w max are divided by their values in the uncoupled case, as these scaled quantities help us assess the increase in the maxima in coupled networks, compared to that obtained in uncoupled patches, allowing us to specifically gauge the coupling-induced effects on the emergent maximum population densities. It is evident from the simulation results that the magnitude of the maximum vegetation and prey populations (i.e. u max and v max denoted by the red and green curves) do not increase much as coupling range increases. However, the maximum predator population increases very significantly with coupling range, with w max at k ′ = 2 exceeding over 6-fold the value obtained for uncoupled patches. However again, as noticed earlier, unlike the effect of increasing coupling strength (cf. Fig. 5), increasing coupling range does not enhance the maximum value monotonically. Rather, the extreme values are most pronounced for nearest neighbour coupling (i.e. k ′ = 2 ). As coupling range further increases, the maxima saturate to a value approximately twice that of the uncoupled case. Similar trends are also observed for the biomass of vegetation, prey and predator for varing coupling ranges.
So we find that increasing the spatial range of coupling does not lead to monotonic enhancement of extreme events. Rather, there is an optimal size of coupling neighbourhood that maximizes the effect, and increasing coupling neighbourhoods further leads to reduction of the magnitude of extreme events. This indicates that the advent of large deviations arise from a subtle and non-trivial interplay of coupling range and coupling strength, and largest deviations from mean values emerge when there are strong interactions over small number of neighbouring patches.

conclusions
In summary, we have studied the population dynamics of a ring of patches with vegetation, preys and predators, where the patches are coupled through interactions of the Lotka-Volterra type. We find that this system yields extreme events in the predator population in the patches, with bursts of explosive predator population growth   www.nature.com/scientificreports/ Our results then are important, both from the view-point of general models of coupled nonlinear systems, as well as for the more specific implications it may potentially hold for population dynamics. So first we have demonstrated how a deterministic system, with generic Lotka-Volterra type of interactions, can give rise to extreme events in space and time. Such examples are uncommon, and so they are significant in the context of general complex systems. Secondly, in the specific context of population dynamics, our model system suggests how predator population densities can grow explosively in certain patches. Though relatively rare, the magnitude of these extreme bursts of predator population density is so huge, that the damage or ensuing cost to control the event, is considerable. Further the biomass of predators can also grow extremely large at certain points in time, and this is of significance due to the catastrophic effects large predator populations can have on the ecosystem as a whole. Lastly, interestingly, the Lotka-Volterra class of interactions has also been used extensively in economic theory 20 , and so our results may have some bearing on extreme events in the financial context.