Stoichiometric plant-herbivore models and their interpretation

The purpose of this note is to mechanistically formulate a mathematically tractable model that specifically deals with the dynamics of plant- herbivore interaction in a closed phosphorus (P)-limiting environment. The key to our approach is the employment of the plant cell P quota and the Droop equation for its growth. Our model takes the simple form of a system of two autonomous ordinary differential equations. It can be shown that our model includes the LKE model (Loladze, Kuang and Elser (2000)) as a special case. Our study reveals that the details of ecological stoichiometry models really matter for quantitative predictions of plant-herbivore dynamics, especially at intermediate ranges of the carrying capacity.

1. Introduction. All organisms are composed of chemical elements such as carbon, nitrogen, and phosphorus. Although the relative abundance of these chemical constituents is known to vary considerably among species and across trophic levels, most ecological studies have until very recently ignored the sources and consequences of this chemical heterogeneity. From theoretical perspectives, Lotka (1925) and other early workers highlighted potential complications raised by having multiple currencies in ecological dynamics, but most subsequent work has focused instead on the dynamic implications of single currency (e.g., energy or carbon) models. However, rapidly accumulating evidence suggests that the dynamic implications of chemical heterogeneity among species deserve much more study than the subject has yet received. This body of research, which is to date chiefly empirical in nature, places major emphasis on the consequences of chemical heterogeneity among species for consumer-resource dynamics and nutrient recycling in ecosystems. Such multiple currency considerations enable simultaneous assessment of both food quantity and food quality. We refer to this approach as "ecological stoichiometry" (Elser et al. 2000, Sterner andElser 2002).
To complement and take advantage of the fast-growing empirical study of ecological stoichiometry, a variety of stoichiometry-based population models have been proposed and studied in recent years. These models vary from simple phenomenological two-dimensional resource-consumer models to more mechanistically formulated systems consisting of dozens of ordinary differential equations. Simple phenomenological two-dimensional resource-consumer models such as that of Loladze et al. (2000Loladze et al. ( , 2004 are mathematically tractable but lack convincing mechanistic basis on the resource dynamics. In contrast, the more mechanistic but complex model of Kooijman (2000) is mathematically intractable. Other models, such as those of Andersen (1997) and Grover (2002), are of intermediate complexity and deal with specific settings such as an open environment or assume that the stoichiometries of all species are constant. In all these models, plant-herbivore interactions may shift from a (+, −) type to an unusual (−, −) class. This leads to dynamics with multiple equilibria, where bistability and deterministic extinction of the herbivore are possible. The most noteworthy dynamics is the birth of bistability as a result of large values of the carrying capacity (K), which divides the plant-herbivore phase plane into two regions: one region with low-density but good-quality plants that sustain high-densities of herbivores, the other region with high density but lowquality plants that can sustain only low densities of herbivores (Loladze et al. 2000; see also Van de Koppel et al. 1996). In general, expressing plant-herbivore interactions in stoichiometrically realistic terms reveals qualitatively new dynamical behavior.
The purpose of this note is to formulate a simple, mathematically tractable model that provides a more mechanistic interpretation of the dynamics of plant-herbivore interactions in a phosphorus (P )-limited environment. The key to our approach is the employment of variability in the P content of the plants, using the Droop equation for the plant's growth. Our model takes the simple form of a system of two autonomous ordinary differential equations. It can be shown that the model of Loladze, Kuang, and Elser (2000), which we shall henceforth call the LKE model, is simply a special case of our model. To aid our model formulation and its comparison with the LKE model, it is convenient to recall here the main LKE model assumptions. They are A1. The total mass of phosphorus P t in the entire system is fixed; i.e., the system is closed for phosphorus.
A2. Phosphorus to carbon ratio (P :C) in the plant varies, but it never falls below a minimum q (mg P/mg C); the herbivore maintains a constant P:C ratio, denoted by θ (mg P/mg C).
A3. All phosphorus in the system is divided into two pools: phosphorus in the herbivore and phosphorus in the plant.
Assumption (A3) essentially assumes that free phosphorus is immediately taken by the plant. The LKE model takes the relatively simple form x is the density of plant (in milligrams of carbon per liter, mg C/l); y is the density of herbivore (mg C/l); b is the intrinsic growth rate of plant (day −1 ); d is the specific loss rate of herbivore that includes metabolic losses (respiration) and death (day −1 ); e is a constant production efficiency (yield constant); K is the plant's constant carrying capacity that depends on some external factors such as light intensity; f (x) is the herbivore's ingestion rate, which may be a Holling type II functional response.
Note that The left-hand side of (1.1) is a logistic equation, where (P − θy)/q is the carrying capacity of the plant determined by phosphorus availability. The right-hand side shows that it can be viewed as Droop's equation (Droop 1973), where q is the minimal phosphorus content of the plant and (P − θy)/x is its actual phosphorus content.
2. Model formulation. Because our main purpose is to derive a mathematically tractable stoichiometry-based plant-herbivore model in a closed environment, we assume that the total amount of phosphorus in the system, P t remains constant. This is equivalent to assumption (A1) of the LKE model. We consider first that only phosphorus (P) is limiting the growth of both plants and herbivores. Later, we add the possibility that a second factor (such as carbon) limits their growth. Since the stoichiometry of herbivores is relatively stable compared to the stoichiometry of plants, we assume that the P:C ratio of the herbivore is a constant θ, as in assumption (A2) of the LKE model. However, in the following, we will not assume (A3), which is rather restrictive and not necessary for our model formulation.
If we let P p , P z and P f be the phosphorus in plant, phosphorus in herbivore, and the free phosphorus respectively, then P t = P p + P z + P f . Let x = x(t) be the plant density, y = y(t) be the herbivore density, and Q = Q(t) be the plant's cell quota for P ; then P p = Qx and P z = θy. Hence In the following, we let q be the plant's minimal cell quota for P , µ m be the plant's true maximal growth rate, D be its death rate, and f (x) be the herbivore's ingestion rate (functional response). We use a variable-internal-stores model based on the Droop equation that relates growth rate to the internal cell quota (Droop, 1973(Droop, , 1974. This approach has been systematically studied by Grover (1991Grover ( , 1997 in the context of plant competition and used by many others (Andersen (1997), Ducobu et al. (1998)). We then have the following equation for the plant growth: Following the notation of Loladze et al. (2000), we let e be the herbivore's yield constant, which measures the conversion rate of ingested plant into its own biomass when the plants are P rich (when Q ≥ θ), and d be the specific loss rate of herbivore that includes metabolic losses and death. If the plants are P poor (when Q < θ), then the conversion rate suffers a reduction, and it becomes eQ/θ. This approach follows the Liebig's (1840) minimum principle and is used in Loladze et al.'s (2000) model formulation. We have the following growth equation for herbivore: Finally, we need an equation governing the dynamics of Q, the plant's cell quota for P . We assume that Q's recruitment comes proportionally from the free phosphorus (αP f ) and its depletion because of cell growth is µ m (Q − q). This results in the following simple equation Since Q(0) ≥ q, mathematically, this ensures that Q(t) ≥ q for all t > 0. Since the cell metabolic process operates in a much faster pace than the growth of total biomass of either species, the quasi-steady-state argument allows us to approximate Q(t) by the solution of which takes the form of This together with (2.2) yields

Substituting (2.8) into (2.7) yields
Substituting the above into (2.3) and applying some straightforward simplification yields (2.10) The above equation can be rewritten as To compare this equation to that of the LKE model, we rewrite the above equation So far, we considered phosphorus only. At this point, we now introduce the possibility that carbon is also a potentially limiting factor. For simplicity, we assume that if only carbon acquisition would limit the growth of the plant, then its population dynamics can be described by the classic logistic equation. Applying Liebigs minimum principle to phosphorus limitation versus carbon limitation of the plant and the herbivore, we arrive at the following simple plant-herbivore model: We note that, compared to the LKE model, this simple two-dimensional model has a stronger mechanistic basis. The model parameters, such as α and µ m , can be directly obtained from physiological measurements.
3. Some special cases. The LKE model can be interpreted as a special case of our mechanistic formulation. The assumption (A3) of the LKE model stipulates that all phosphorus is in the herbivores and plants; that is, the concentration of freely available phosphorus (P f ) is zero. This is tantamount to saying that the phosphorus uptake rate of the plants is extremely efficient. That is, α = ∞ or, equivalently α −1 = 0. If, in addition, we assume that the plant death rate D is negligibly small compared to its maximal growth rate, we may approximate the value of (µ m − D)/µ m as 1. Indeed, if we assume α → ∞ and (µ m − D)/µ m ≈ 1, our mechanistic model simplifies to the following form: From equation (2.8), we notice that Hence, as α tends to ∞, we see that Q tends to (P t − θy)/x. Notice further that This demonstrates that the LKE model can be interpreted as a special case of our mechanistic model, under the assumption that all phosphorus is taken up by the plants and herbivores. Another important implication of our model (2.13) is that the plant dynamics reduces to a form that resembles the classical logistic equation when the herbivore is absent and (A3) holds It should be pointed out here that we did not assume the population suffers from a crowding effect explicitly. However, this crowding effect is implicitly provided by the fact that the total nutrient in the system (here P ) is fixed, and individuals have to compete for this resource. Observe that instead of the often-assumed carrying capacity of the form P t /q, here the carrying capacity has the expression of This says that although theoretically the environment may accommodate P t /q plants, the actual upper limit the plant biomass can attain is [(µ m − D)/µ m ]P t /q, somewhat less than that. The reason the maximal carrying capacity P t /q cannot be reached in practice is that the death toll in a population keeps the population below its potential maximum. In other words, it says that a population with a relatively low death rate will likely amass more biomass than a population with a relatively high death rate.
4. Dynamics along a productivity gradient. To end this short note, we provide some snapshots of the dynamics of model (2.13), where we choose the Monod function as the functional response of the herbivores, In our snapshot series, we varied the carrying capacity based on carbon (the parameter K). This might reflect, for instance, an increase in light supply, yielding a higher photosynthetic carbon assimilation rate of the plants (see de Koppel et al. (1996) for an empirical and theoretical investigation on herbivory patterns along a productivity gradient). The parameter values are listed in Table 1 (except for D, which is µ m − b = 0), where we used the ones given in Loladze, Kuang, and Elser (2000), which are estimated from Andersen (1997) and Urabe and Sterner (1996). We used MATLAB release 12 for the computer simulations. In Figure 1, solutions of model (2.13) are compared to that of the LKE model for various values of carrying capacity K. Both models demonstrate similar qualitative changes in the dynamics. As carrying capacity increases, the stable positive steady state loses its stability through a Hopf bifurcation, and gives rise to a limit cycle. Further increasing the carrying capacity will collapse the cyclic behavior through a heteroclinic bifurcation and returns the dynamics to a simple steady-state behavior where the herbivore has gone extinct because of the low food quality of the abundant plants. Observe that in the cases when K is low or very high, our model and the LKE model produce almost identical solutions ( Fig. 1A and 1D).    Table 1. really matter. The bifurcation diagrams for both models (see Figure 2) further confirm these observations.