ALTERNATIVE MODELS FOR CYCLIC LEMMING DYNAMICS

Many natural population growths and interactions are affected by seasonal changes, suggesting that these natural population dynamics should be modeled by nonautonomous differential equations instead of autonomous differential equations. Through a series of carefully derived models of the well documented high-amplitude, large-period fluctuations of lemming populations, we argue that when appropriately formulated, autonomous differential equations may capture much of the desirable rich dynamics, such as the existence of a periodic solution with period and amplitude close to that of approximately periodic solutions produced by the more natural but mathematically daunt ing nonautonomous models. We start this series of models from the Barrow model, a well formulated model for the dynamics of food-lemming interaction at Point Barrow (Alaska, USA) with sufficient experimental data. Our work suggests that an autonomous system can indeed be a good approximation to the moss-lemming dynamics at Point Barrow. This, together with our bifurcation analysis, indicates that neither seasonal factors (expressed by time dependent moss growth rate and lemming death rate in the Barrow model) nor the moss growth rate and lemming death rate are the main culprits of the observed multi-year lemming cycles. We suspect that the main culprits may include high lemming predation rate, high lemming birth rate, and low lemming self-limitation rate.

1. Introduction.Pioneer works on resource-consumer dynamics include the well known works of Lotka (1925) and Volterra (1926) which introduced the classical Lotka-Volterra predator-prey model that arguably forms the foundation of mathematical ecology.One of the most frequently used resource-consumer models is the Rosenzweig-MacArthur (1963) model, which produces two generic asymptotic behaviors -equilibria and limit cycles.Bazykin (1974) added a self-limitation term to the Rosenzweig-MacArthur model to account for the rather ubiquitous density dependent mortality rate (see also Bazykin et al. 1998).All these models produce oscillatory solutions that seem to mimic the fluctuating populations observed in nature.
Indeed, large-scale high-amplitude oscillations in populations of small rodents such as voles and lemmings have been a constant inspiration to numerous influential and thought provoking articles since the pioneering work of Elton (Elton 1924, Hanski et al. 2001).Lemming is a mouse like arctic rodent characterized by a small, short body that is about 13 cm (about 5 in) long, with a very short tail.Lemmings live in extensive burrows near the water, feed on vegetation, and build nests out of hair, grass, moss, and lichen.The female produces several broods a year, each of which contains about five young.Many researchers believe that such 2 HAO WANG AND YANG KUANG oscillations are more or less the exhibitions of the usual characteristics of resourceconsumer interactions.More specifically, ecologists tend to believe that the cause of such oscillations is either an interaction between lemmings and their food supply (Turchin and Batzli 2001), or an interaction between lemmings and their many predators (Hanski et al. 2001).
All key reproductive events for lemmings take place during winter.Empirical research is difficult on organisms that live under snow due to the extreme low temperature.However, ecologists have managed to obtain collections of data, even though they are very limited in many aspects.A family of two and three dimensional autonomous and nonautonomous differential equation models incorporating various biological mechanisms are described and their numerical dynamics are compared by Turchin and Batzli (2001).They arrived at the conclusion that their three dimensional nonautonomous differential equations Barrow model (named after the brown lemmings at Point Barrow, Alaska, USA) is the most appropriate one in their family of models.The objective of this paper is to show that when appropriately formulated, a two dimensional autonomous differential equation model (the Bazykin predator-prey model) can capture much of the desirable dynamics exhibited by their Barrow model.Specifically, we show that the period and amplitude of the stable periodic solution of our autonomous lemming-moss model are close to that of the approximately periodic solutions of the Barrow model.
The Barrow model is formulated to describe the apparent seasonal interactions among lemmings, moss and vascular plants.It takes the following form (Turchin and Batzli 2001). (1.1) Here V =vascular plant (v.plant) density, M =moss density and H=lemming density.τ = τ (t) = t − t ∈ [0, 1) is the season indicator.Summer (two months from melt-off in mid-June to first heavy frosts in mid-August) is represented by the interval of 5/6 ≤ τ < 1.In summer, Winter is represented by the interval of 0 ≤ τ < 5/6.In winter, Description of parameters and their value ranges are given in Table 1.At τ = 0 (transition between summer and winter), V(t) suffers a sudden 90% reduction due to the first heavy frost.Observe that G s is the summer lemming death rate in summer divided by the conversion rate R. We will call it the modified summer lemming death rate, or simply, the lemming summer death rate.Similar statements are true for G w .Together or alone, they are ambiguously referred as lemming death rate.Different forms (linear and logistic) are used in vascular plants and moss regrowth terms.Detailed justification and background information can be found in Turchin and Batzli (2001).Seasonal effect is taken into account, especially in resource (vascular plants and moss) regrowth term (Oksanen 1990).From the shape indicated by field data and the solution profiles generated by some plausible lemming-moss models, Turchin et al. (2000) suggest that the predator nature of the lemming is the main driving force of the observed dramatic lemming population cycles.In other words, the lemming-moss interaction generates lemming population cycles.    1 Or equivalently, the moss dynamics controls and induces lemming population oscillations.This hypothesis is hence referred as bottom-up regulation.Figure 1 represents a typical Barrow model simulation result.This periodic Barrow model, while embodies many realistic features, it is somewhat awkward to describe, simulate and almost impossible to analyze mathematically.A natural and interesting question is whether some simpler and plausible autonomous differential equation models can reasonably approximate its rich and realistic dynamics.Our answer to this question is positive.
In the next section, we formulate an intermediate two dimensional nonautonomous moss-lemming model resulted from canceling vascular plant variable and adding self-limitation term in the Barrow model.We then replace the season-dependent moss growth rate and the modified lemming death rate by suitable continuous time dependent functions.We carry out some simple mathematical analysis for this intermediate model.In section three, we formulate a two-dimensional autonomous model by simply choosing moss growth rate and lemming death rate to be their naturally weighted mean values, respectively.This model is mathematically tractable and a global qualitative analysis is attempted.In addition, we present bifurcation diagrams for several key parameters to gain additional insights on how these parameters affect the amplitudes and periods of those periodic solutions.We then proceed to the conclusion section of this paper: comparing the simpler autonomous model with those non-autonomous models in four key aspects expressed by their oscillatory solution profiles.We provide a brief discussion describing the implications of our findings in the contexts of some specific existing biological observations and theoretical statements.The paper ends with an ad hoc procedures for finding out if an autonomous differential equation model can be an acceptable approximation to a nonautonomous one.

2.
A nonautonomous lemming-moss model.Since moss is the main food supply for lemmings, it is thus natural to reduce the Barrow model to a two dimension one by canceling vascular plant variable.Following Bazykin (1974), we include self-limitation term in the reduced system to account for various possible mechanisms that may introduce additional lemming density dependent mortality.Indeed, Chitty (1996) argues for such self-limitation effect in lemming dynamics.This leads to the following nonautonomous lemming-moss model with self-limitation. (2.1) The self-limitation rate E depends on many factors, possibly includes lemming behavioral changes such as from normal to more aggressive behavior as lemming population density increases.Chitty (1996) suggests that behavioral changes can be key factors inducing lemming cycles.As we will see, mathematically, the selflimitation term actually keeps lemming's population away from very low level, which makes the model more realistic Simulation results (see Figure 2) strongly suggest that the above nonautonomous lemming-moss model and Barrow model produce qualitatively similar dynamics.
To facilitate a systematic qualitative study of the above model, we would like to convert it into a model with fewer and more familiar parameters and variables.To this end, we let   1 and E = 0.01.
In reality, u(τ ), d(τ ) change continuously as season progresses.We thus try to use two continuous functions u(t), d(t) to replace the discontinuous functions employed by Turchin and Batzli (2001) in their Barrow model.This simple and natural modification enables a good amount of mathematical analysis.The new system takes the form of (2.3) In the following, we let to approximate the discontinuous function u(τ ) and d(τ ) below u(τ ) = 0 (0 ≤ τ ≤ 5/6), or u s (5/6 ≤ τ < 1), Using the median values of u s , d s , d w , we can visually and statistically compare u(t), d(t) with u(τ (t)) and d(τ (t)) (see Figure 3).
We now in a position to start our mathematical analysis for the nonautonomous model (2.   1 and l = 0.01.Since x = 0 and y = 0 are solutions of model (2.3), we see that the positive cone As usual, it is easy to show that (0, 0) is a saddle, stable along y axis and unstable along x axis.
The linearized system of (2.3) at the origin takes the form of The conclusion of the theorem follows immediately.
Our next analytical result says that model (2.3) is dissipative, which is equivalent to say that all nonnegative solutions are eventually uniformly bounded.
3. An autonomous lemming-moss model.We are now in a position to formulate an autonomous version of the lemming-moss interaction model.The method is simply replacing the season dependent parameters by their weighted averages.This yields where ū = u s /6 (mean 2), d = (5d w +d s )/6 (mean 8.975) are the weighted averages.This model is called the Bazykin predator-prey model (see Bazykin et al., 1998 andKuznetsov, 2004, p 325).Short of chaotic dynamics, the Bazykin model is capable of generating rich and complicated.For example, it may have three positive steady states and several limit cycles, and it may undergo several codimension 2 bifurcations and generate a large homoclinic loop (Kuznetsov, 2004, p. 325).Indeed, many mathematical questions on the qualitative properties of the model remain open (Hwang and Kuang, 2006).If l = 0 in (3.1), then it is the well studied Rosenzweig-MacArthur Holling type II predator-prey model (Kuang, 1990).Figure 5 is a solution of model (3.1) with the same initial condition and parameter values as the solution of (2.3) depicted in Figure 4.By comparing the simulation result with that of model (2.3) (Figure 6), we see that the period is slightly longer.However, their amplitudes are comparable in log scale.
A systematical bifurcation and global phase plane analysis has been carried out by Bazykin (Bazykin et al. 1998) and the graphical description of the 11 cases can also be found in (Kuznetsov, 2004, p. 331).The lemming-moss model (3.1) dynamics (with the parameters given in Table 1) seems to exhibit dynamics as depicted by case 1 (the interior steady state is unique and a sink, no positive limit cycle) or case 3 (the interior steady state is unique and a spiral source, there exists at least one and most likely a unique limit cycle).
From the bifurcation diagrams, we see that the following statements hold.

Select u s (which is 6ū) as the bifurcation variable and choose median values
for other parameters.If 6 < u s < 24 (special region for lemming), the system always behaves like case 3 and limit cycle exists all the time.
In the following, we analyze model (3.1) with b ∈ [0.5, 1], varying values of l, and other parameters taking the median values given in Table 1.For such parameters, (3.1) has steady states (0, 0), (K, 0) and appear to have an unique positive steady state (x * , y * ).As in the case of model (2.3), the nonnegative cone is invariant and model (3.1) is dissipative.It is also clear that (0, 0) is a saddle with y-axis as its stable manifold and x-axis as its unstable manifold.(K, 0) is a saddle with x-axis as its stable manifold.For convenience, bars above u and d will be dropped from now on.
The following simple lemma asserts that if the moss density at the positive steady state is high enough (higher than half of the difference of the carrying capacity and the Michaelis-Menten constant c), then the positive steady state is locally asymptotically stable.This is well known when l = 0 (Rosenzweig and MacArthur 1963).The variational matrix at (x * , y * ) is , in which case we see that the trace of J(x * , y * ) is negative and the determinant of J(x * , y * ) is positive.Hence, (x * , y * ) is locally asymptotically stable.
The next theorem provides an explicit condition for the unique positive steady state to be locally asymptotically stable.As to be expected, the median parameter values given in Table 1 fail to satisfy the condition.
Then the positive steady state (x * , y * ) is locally asymptotically stable.

Proof.
The positive steady state (x * , y * ) satisfies which is equivalent to We have H(K) < 0 and We see that ).This shows that x * > (K − c)/2 and by Lemma 3.1 we conclude that (x * , y * ) is locally asymptotically stable.
Our main mathematical theorem below suggests that low lemming consumption rate of moss and very high intraspecific competition among lemmings can stabilize the system dynamics.Again, the median parameter values given in Table 1 fail to satisfy the condition.Theorem 3.2.Assume model (3.1) has a unique positive steady state and l > a/(4c).Then the positive steady state (x * , y * ) is globally asymptotically stable.

Proof.
Notice that both the origin and the boundary steady state (K, 0) are saddles and (K, 0) has its unstable manifold pointed inward to the positive cone {(x, y) : x > 0, y > 0}.The conclusion of the theorem follows easily from the Poncare-Bendixson theorem and the dissisipativity of the system once we show   We shall show that a (x + c) 2 < l x for 0 < x < K.This is equivalent to say that l > ax/(x + c) 2 ≡ L(x) for 0 < x < K. Notice that L (x) = a(c − x)/(x + c) 3 , we see that L(x) achieves its maximum at c with the value of a/(4c).This shows that 4cl > a implies that l > ax/(x + c) 2 ≡ L(x) for 0 < x < K.In other words, we have shown that ∆(x) < 0 for 0 < x < K, which implies that model (3.1) has no positive limit cycle when 4cl > a.
The following theorem follows naturally from the dissipativity and analytical nature of model (3.1) (Kuang 1989).Theorem 3.3.Whenever (x * , y * ) is a spiral source, there exists at least one stable positive limit cycle in the strip {(x, y) : 0 < x < K, y > 0}.

Conclusion: comparison of autonomous and nonautonomous models.
Comparing solutions of models (2.3) and (3.1) in Figure 6, we can see that the main dynamical features (oscillatory nature, approximate periods and amplitudes) To gain a close look at how amplitude and period of the autonomous model change as some of the parameters vary within the range given in Table 1, we turn to some amplitude bifurcation diagrams and period diagrams.Amplitude bifurcation diagrams with respect to four key parameters are shown in Figure 7. Amplitude bifurcation diagrams for moss growth rate u and the modified lemming death rate d have similar behaviors, neither has any bifurcation point.In fact, both parameters have only very small influence on amplitudes of the cycles in their ecologically meaningful ranges.This seems to explain why amplitudes of autonomous model (3.1) are very close to those of nonautonomous model (2.3).In contrast, both bifurcation diagrams for conversion rate b and self-limitation rate l have a (Hopf) bifurcation point, respectively.
Period diagrams for the same four parameters are plotted in Figure 8. Period as a function of u and d is a smooth and slow changing one.Hence this explains why periods of autonomous model (3.1) are close to those of nonautonomous model (2.3).On the other hand, the period curves are discontinuous at the bifurcation point for both b and l, respectively.This coincides well with the observations from the amplitude bifurcation diagrams.The period changes fast when b or l varies near its bifurcation point.In this neighborhood of parameters, the periods of autonomous model (3.1) can be significantly different from that of nonautonomous one.
Table 3 contains the comparisons of the lemming variable of all four models appeared in this paper with real data from Gilg et al. (2003).The features compared are 1) the approximate periods, 2) the order of magnitudes from lemming low density level to its high density level, 3) the number of years from low lemming density level to high lemming density level vs the number of years from high lemming density level to low lemming density level, 4) and the peak shape.According to these key dynamical features, the autonomous moss-lemming model (3.1) actually provides a slightly better approximation than the Barrow model, if we take into account each and every period.
Our work clearly suggests that autonomous system can indeed be a good approximation to the moss-lemming dynamics at Point Barrow.This together with our bifurcation analysis indicate that neither seasonal factor (expressed by time dependent moss growth rate and lemming death rate in Barrow model), nor the moss growth rate and lemming death rate are the main culprits of the observed multi-year lemming cycles.Main culprits may include high lemming predation rate a (since conversion rate is more or less constant, this is equivalent to high lemming birth rate) and low self-limitation rate l as they enable cyclic population dynamics.
Indeed, Hoffmann (1958) suggested that reproductive changes play relatively minor roles in inducing cyclic population dynamics.However, Hoffmann (1958) thought that important factors causing population oscillations must include mortality changes.Our result suggests that this thought is not well grounded theoretically.Our conclusion that seasonal effect does not play a key role in inducing cyclic lemming dynamics is in agreement with a similar statement of Krebs (1964) made on lemming cycles at a different location.In addition, our findings confirm the field observation of Erlinge et al. (2000) who stated that higher litter sizes (equivalent to higher birth rates) may trigger lemming cycles.However, our conclusion that low self-limitation rate promote cyclic dynamics contradicts Chitty's (1996) hypothesis that claims high level of aggressiveness (equivalent to higher self-limitation rate) induces lemming cycles.
We end our paper by recaping our ad hoc but effective procedure for determining whether a time dependent parameter can be replaced by a constant in an oscillating system: plot amplitude and period diagrams in a biologically reasonable range for that parameter.If there is no bifurcation point and the period curve is continuous, then we can replace that time dependent parameter by its weighted average.Biologically, this suggests that parameter is not a key one for generating the observed oscillatory solutions.
Lemming vs vascular+moss in Barrow model

Figure 1 .
Figure 1.Barrow model simulation result with parameters taking the median values shown in Table 1 Lemming vs moss in lemming-moss model

Figure 2 .
Figure 2. Nonautonomous lemming-moss model simulation result with parameters taking the median values shown in Table1and E = 0.01.
3) where u(t) and d(t) are continuous periodic functions in t, with period 1 and 0 ≤ u(t) ≤ u s , d s ≤ d(t) ≤ d w .

Figure 4 .
Figure 4.A typical solution of the nonautonomous lemmingmoss model (2.3) with parameters taking the median values shown in Table1and l = 0.01.

Figure 5 .
Figure 5.A typical solution of the autonomous lemming-moss model (3.1) with parameters taking the median values shown in Table1and l = 0.01.

Figure 7 .
Figure 7. Bifurcation Diagrams and limit cycles for the autonomous system

Theorem 3 . 1 .
Assume model (3.1) has a unique positive steady state and Period diagram for u s (= 6ū)

Figure 8 .
Figure 8. Period Diagrams for the autonomous system

Table 2 .
Parameters in lemming-moss models are represented by the red continuous curves.Using mean values of u s , d s , d w , we can compare u(t), d(t) with u(τ ), d(τ ) statistically.The standard/average error of u(t) with respect to u(τ (t)) is 1.538 and the relative error is 1.538/u s = 0.128.The standard/average error of d(t) with respect to d(τ (t)) is 0.401 and the relative error is 0.401/(d w − d s ) = 0.141.

Table 3 .
Comparison of all four lemming models We notice that both the amplitude and period in autonomous system are a little larger in average than those in model (2.3).Moreover, these features are all closely resemble those of the Barrow model.