A model of trophodynamics

Nonlinear Analysis 35 (1999) 37 – 57

A model of trophodynamics Brad Lackey∗; 1 Department of Mathematical Sciences, University of Alberta, Edmonton, Alberta, Canada T6G 2G1 Received 22 April 1997; accepted 5 May 1997

Keywords: Trophodynamics; Second-order dierential equations; Method of equivalence

1. Introduction Our goal is to construct and analyse a mathematical model of population dynamics treated as the dynamics of trophic species. Here, we will even consider a single population as a producer or consumer of some basic object. Following the tradition in trophodynamics, we will call these objects basal species, even though they may not be species in the biological sense. In the following section, we will consider some dierent models of a single species population, but in each case the resulting equations will be of the form dx = Q(x; t) − N; dt dN = G(x; N; t); dt where x is some measure of the basal species (usually log density) and N is the density of the population in question. Note that this model is a second-order ordinary dierential equation in x, and the population is related to dx=dt. The relationship between trophodynamics and secondorder dierential equations was rst realized in [2] in the context of Volterra production. There, the basal species is produced at a constant rate by the population, and hence the authors call it a surrogate of biomass. In our context, the basal species dynamics and relationship with the population may be quite general, so this term is not appropriate. ∗ 1

E-mail: [email protected]. Research supported by the Izaak W. Killam Trust.

0362-546X/98/$19.00 ? 1998 Elsevier Science Ltd. All rights reserved. PII: S 0 3 6 2 - 5 4 6 X ( 9 8 ) 0 0 0 9 7 - 2

38

B. Lackey / Nonlinear Analysis 35 (1999) 37 – 57

Yet, the use of second-order dierential equations is the crucial feature. Traditionally, a single population would be modelled by a rst-order dierential equation, such as the logistic equation or some relative thereof. Unfortunately, dynamical systems of this type only allow monotonic approach to equilibrium. Yet, laboratory experiments supply examples where the population density exhibits damped oscillation about the carrying capacity. Perhaps the clearest account of this is F.M. Williams’ chemostat experiments with Selenastrum and Chlorella in [24]. We will take this example up again later. Following our single population examples, we will note the multispecies analog of our trophodynamic equations. The analysis of these equations can be quite complicated, and will be reported elsewhere. We then turn to computing invariants of a second-order dierential equation. We assert that this is the most important aspect of trophodynamic modelling. In de ning the trophic species, often complicated identi cations are made; thus, what the variables of the dynamical equations describe is not always clear. One should be concerned about changing the de nition of the variables, thus changing the structure of the dierential equation, and even dubious about any conclusions that draw explicitly upon a particular choice of variable. Naturally, if one changes variables in a reasonable way (say via a dieomorphism) then one can track the modi cations that occur. Yet, if given two second-order equations in dierent variables, it is not at all clear when they describe the same trophodynamic system. Thus, we are lead to computing the invariants. That is, computing the objects which do not change (or change in a predictable fashion) under a change of variable. Any qualitative conclusions drawn from a trophodynamic system must be phrased in terms of the invariants. Our nal section is dedicated to this: we will nd that stability, social interactions, and logistic growth are all describable as invariant conditions, and thus do not depend on how we formulate the trophic species. We also oer two appendices reviewing the mathematical construct necessary for our analyses. 2. A single trophic population We rst consider models in which a single population consumes (or produces) one basal species. Let us denote the density of the population by N and the density of the basal species as b. We start with presuming dynamics of the form db = Q(b; t) + P(b; N; t); dt dN = G(b; N; t): dt The function Q models the dynamics of the basal species on its own, hence we have taken this function independent of N . We use P to represent the production or reduction of the basal species by the population. The function G governs the population dynamics, which we allow to depend on both b and N . As it stands, this is a general planar dynamical system. To continue in any meaningful way, we need to make some

B. Lackey / Nonlinear Analysis 35 (1999) 37 – 57

39

simplifying assumptions; we will concentrate on the structure of P. Let us note, though, that a traditional choice for G is logistic: G = g(b; t)N − r(b; t)N 2 . By allowing the parameters g; r to depend on b, we may model populations whose growth rate (and carrying capacity) depend on, say, the amount of food available or chemical density present. Yet, the assumption of logistic population growth will not simplify our later analysis, so we will not presume it. This case of logistic growth will appear again, so we merely point it out as a historically successful model. Turning to the structure of P, we will consider two cases: density-independent consumption=production and density-dependent consumption=production. The former is contained in the latter (by choosing the dependence to be constant). We consider the two cases separately to illustrate the two stages in constructing the second-order differential equation which models our system. 2.1. Density-independent consumption=production We assume that each member of the population produces (or consumes) the basal species at a xed rate, independent of all other factors. Thus, the production function must take the form P(b; N; t) = kN . Here, k is a constant whose sign distinguishes consumption from production. We rescale the basal species by setting b = −k x. Thus, the equation governing the dynamics of the rescaled basal species is dx = Q(x; t) − N: dt We would point out that the function Q(x; t) is not the same as Q(b; t), but since we have not stipulated any special form for it, we will not introduce new notation to describe this change. Also, it is important to recognize that if k¿0 (the production case), then x takes negative values. Again, not introducing any new notation for the change in the function G, the general model for density-independent production=consumption is dx = Q(x; t) − N; dt dN = G(x; N; t): dt To illustrate the properties of this system, let us now consider an extended example of a population which produces a chemical or organic waste product which is self-inhibitory. A well-studied example of such behaviour is Chlorella; this algal population produces an autotoxin, named chlorellin, which inhibits its own growth. An excellent exposition of this is in Ch. 6 of [20]. As well, extensive research on the eects of chlorellin on Chlorella can be found in [16–18]. Example 2.1. Consider a population (with density N ) which produces a chemical or organic waste toxic to itself. Let us denote the density of this autotoxin as b. We will assume, as above, the only source for the toxin is from the population, and this

40

B. Lackey / Nonlinear Analysis 35 (1999) 37 – 57

Fig. 1. Behaviour of an Equilibrium. A typical sink (top left) yields an S-shaped curve (top right) as the population approaches equilibrium. (Here: M = 13 ,  = 4,  = 3.) A spiral (bottom left) yields damped oscillatory behaviour (bottom right) as the population approaches equilibrium. (Here: M = 2,  = 1,  = 3.)

production occurs at constant rate k. For simplicity, we will assume that the toxin dissipates in proportion to its density. Thus, the basal species dynamics is modelled by db = −b + kN: dt As for the population, we take a constant growth rate and a death rate in proportion to the density of the toxin. Against traditional logic, we do not take any density-dependent regulatory terms. Thus, we have dN = N − bN: dt As suggested, we rescale the toxin density by b = −k x. Consequently, −x represents a total production of the toxin, having the units of population density times a unit time interval. Now, the system simplify to dx = −x − N; dt dN = (M + x)N: dt Recall that x¡0, so that although the second equation appears to indicate that N will grow without bound, this is not the case. Actually, the stationary points of this system are simple to classify; there are precisely two such points: (x; N ) = (0; 0) and (−M; M ). The point (0; 0) is easily veri ed to be an unstable saddle point. The other point (−M; M ) is stable; it is a sink if 2 ≥ 4M and a spiral if 2 ¡4M . In the laboratory, it would be quite dicult to measure the toxin density. Thus, if we ignore the x-component of the solution curve, and just consider the N -component, the two

B. Lackey / Nonlinear Analysis 35 (1999) 37 – 57

41

cases (sink and spiral) yield dierent qualitative behaviour of the population. The sink exhibits the standard S-shaped curve as the population approaches its carrying capacity N = M . In the spiral case, the population grows past M , levels o and begins to die back, dropping below the equilibrium, having a damped oscillatory behaviour converging to the carrying capacity. 2.2. Density-dependent consumption=production Now, we consider the case where the population consumes (or produces) the basal species in a density-dependent way. As in the previous case, will suppose that each member produces independently of the population density. Yet now, we will allow the production to depend upon the density of the basal species (but not on time explicitly). That is, we consider a basal species dynamics of the form db = Q(b; t) − k(b)N dt We have taken the negative sign here by convention only. Thus, consumption is modelled with k(b)¿0 and production by k(b)¡0. In order to proceed, we assume that k(b) 6= 0 in our domain. As k is nonvanishing, let us take any function f(b) so that df=db = 1=k. Now, we compute Q(b; t) dx df db = · =− − N: dt db dt f(b) We are to solve b = f−1 (x) and substitute this in for b in the above equation. Yet f is invertible if and only if it is monotone, and this is guaranteed by the assumption that k(b) 6= 0. Again, we do not introduce any new notation, and just write Q(x; t) for −Q( f−1 (x); t)=x. As an aside, we note that normalized Gompertzian growth ts nicely into this model. Let us take b to represent the population of the basal species as a proportion of its maximum, so 0 ≤ b ≤ 1. Taking x = log b (hence x ≤ 0), and Q = −x, the basal species will have its independent growth modelled by dx = −x dt −t

which has solution x = x(0)e−t . Thus, b = e k x(0)e . As the above equation has x(0)¡0, this is precisely the Gompertzian growth formula for the normalized b. We conclude, the most general model of density-dependent consumption has the form dx = Q(x; t) − N; dt dN = G(x; N; t): dt

42

B. Lackey / Nonlinear Analysis 35 (1999) 37 – 57

Example 2.2. We now consider the model as applied to F.M. Williams’ chemostat experiment. This actually has two parts: rst, we model the batch experiment, where there is no in ux of nutrients into the system, and then treat the chemostat as a perturbation which introduces nutrients and washes out a proportion of the population. We assume that the cell density of Selenastrum and Chlorella increases in direct proportion to the quantity of nitrate nutrients they consume. See [24] for a justi cation of this. In a batch experiment, there is no addition of nutrients, so assuming consumption in proportion to the level of nutrient available, the model for nutrient=cell density is db = −kbN; dt dN = bN : dt It is simpler to work directly with these equations. We have (d=dt)(b + kN ) = 0, so b = c −kN for some constant c. Therefore, we can reduce the system to just the single dierential equation dN = (c − kN )N; dt and, thus, the population density exhibits logistic growth. We model the activation of the chemostat as a perturbation of the above system: db =  − kbN; dt dN = bN − N: dt Here,  is the dierence of the input nutrients and washed our nutrients per unit time, and  is the proportion of the population washed out by the ow. Properly,  should be a function of N , but this goes against our modelling criterion that Q is not such. We will assume that our culture density is suciently large when the chemostat is turned on that there is no signi cant change in  as N approaches the steady state, hence may be assumed constant. In this case , which would in general depend on N as well, may be assumed constant for the same reason. Oddly enough, Williams’ experiments satis ed this criterion as smaller populations tended to washout of the chemostat before cell division could occur. Again, it is easier to work directly with this system. The only stationary point is (b; N ) = (=; =k) which is stable. In the case  ≥ (=2)2 , the equilibrium is a sink, and hence the population exhibits an S-shaped approach to its carrying capacity. If ¡(=2)2 then the point is a spiral, and the population exhibits damped oscillation. Example 2.3. We consider modelling a plants chemical defense against herbivores using Rhoades’ mechanism [1, 14, 19]. In this example, N will represent the density of our herbivore. We will have two competing basal species: Q representing the leaf

B. Lackey / Nonlinear Analysis 35 (1999) 37 – 57

43

biomass, and T representing the tannin content. We will suppose that Q has normalized Gompertzian growth. Rhoades’ plant response mechanism is encoded in the relation T = d · Qg where 0 ≤ g ≤ 1. Thus for g = 1, we have the tannin-leaf ratio is constant, T=Q = d, indicating no response to defoliation. Yet if g = 0, then T=Q = d=Q, a reduction on Q yields a dramatic increase in tannin per leaf. As Q and T are related, we may choose either to represent our consumed basal species, here we will choose T . As Q satis es normalized Gompertzian growth (in the absence of the herbivore), we have d log Q = − log Q − (consumption term): dt If we assume that each herbivore consumes the leaves in direct proportion to the number of leaves available, then consumption is modelled by dQ=dt = −kQN , or equivalently (d log Q)=dt = − kN . So the dynamics of the leaf biomass is d log Q = − log Q − kN: dt We can easily convert this to tannin dynamics by recognizing that log T = log d + g log Q. If we write x = (1=gk) log T for log tannin density, we have a basal species equation dx  log d = −x − − N: dt gk The dynamics of our herbivore may be quite complicated. Yet an important ratio in determining the energy the herbivores derive from the leaves is T=Q. One can easily verify that the log of this ratio is directly related to x, log(T=Q) = (g − 1)kx + (log d)=g. 2.3. Converting to a second-order equation As we have shown above, our assumptions lead dynamics to the form dx = Q(x; t) − N; dt dN = G(x; N; t): dt The remainder of our model is based on the concept of the ‘excess producer’. We de ne this as y = Q(x; t) − N , so that dx=dt = y. Thus, y can be viewed as producing the basal species, in excess of what is being consumed by the population. In both of the examples above, the stable equilibrium occurs as y → 0. The dynamics of the excess producer is easily computed: dy @Q dx @Q dN = + − dt @x dt @t dt @Q @Q = y+ − G(x; Q(x; t) − y; t): @x @t

44

B. Lackey / Nonlinear Analysis 35 (1999) 37 – 57

Fig. 2. The Excess Producer, y, is de ned to be the ‘population’ whose eect would be to produce the basal species. At each time step, the basal species would grow independently to the height of the new bar, the population consumes N . What remains is the eect of the excess producer. Note that the density of the excess producer may be positive or negative.

We conclude that our dynamics are expressible as dx = y; dt dy + F(x; y; t) = 0: dt That is, we have constructed a second-order dierential equation in x, whose solutions model the growth of both the basal species and the population. In fact, if we consider the structure of F(x; y; t) = G(x; Q(x; t) − y; t) − (@Q=@x)y − (@Q=@t), we see that if G is linear in N (corresponding to exponential growth) then so is F in y. In general, if G is a polynomial of degree n in N , then F is also in y. In particular, if G is logistic growth, then so is F with possibly a source term due to the independent dynamics modelled with Q. The concept of the excess producer has a second bene t. Regardless of the complicated transformations that may go into a change in x; N , the excess producer always changes as the derivative of x. Thus, we have replaced the task of understanding how two rst-order equations change with both variables to one with only one second-order equation and one changing variable (after all, the changes in y are completely determined by those of x). One may be concerned, and rightfully so, that we have lost some information. But this is not the case: N is completely recoverable from y, and the transformational character of F takes both the transformations of x and N into account. 3. The multispecies model We now consider the analog of our simple model above in the case of several interacting populations. As above, we will separate the basal species from the populations; let us denote the densities of these by b1 ; : : : ; bm and N 1 ; : : : ; N r , respectively. As in the case above, we will assume that each member of population a consumes (or produces)

B. Lackey / Nonlinear Analysis 35 (1999) 37 – 57

45

basal species j independent of the population densities. Thus, the preliminary form of our model is db j = Q j (b; t) − (K1j (b)N 1 + · · · + Krj (b)N r ): dt We have suppressed the indices of functional variables for readability; of course, Q j (b; t) means Q j (b1 ; : : : ; bm ; t). Naturally, we have some conditions on the structure of the matrix Kaj (b). First we require that over our domain, the rank of K is constant. This is more that just a technical assumption, but let us postpone this discussion until the example below. As K has constant rank, we can nd matrices A(b) and B(b) such that   I 0 AKB−1 = n 0 0 where In is the n × n identity matrix and n is the rank of K. The matrices A and B are not unique, and we will discuss elsewhere the freedom in their choice, and what role this freedom plays in the theory. Our second condition is merely a integrability criterion: we assume that @Akj (b) @Alj (b) = : @bl @b k This guarantees that we can nd a collection of functions fj (b) so that @fj (b) = Akj (b): @b k These functions serve that same role as in the single population case. We de ne xj = fj (b). Then dx @f db = = A(Q − KN ) dt @b dt = AQ − AKB

−1



BN = AQ −

In 0

 0 BN: 0

We have suppressed the indices in the above equations for ease of reading. Considerj

j

ing the right side of the equation, we de ne x = ( wx  ) and AQ =( QR ) where j = 1; : : : ; n j

and  = n + 1; : : : ; m. Similarly, we de ne BN =( Nz ) where j has the same range and  = n + 1; : : : ; r. Then our basal species equations split into two component equations:        d x Q(x; w; t) I 0 N = − n R(x; w; t) 0 0 z dt w  =

Q(x; w; t) R(x; w; t)



 −

N 0

 :

46

B. Lackey / Nonlinear Analysis 35 (1999) 37 – 57

j We again de ne the excess producers as y j = Q j (x; w; t) − N . Using the governing equation for the consumers, we can nd dynamical equations for y and z. Generally, these equations are complicated, but the nal form of our system is

dw = R (x; w; t); dt dxj = y j; dt dy j + F j (w; x; y; z; t) = 0; dt dz  + H  (w; x; y; z; t) = 0: dt In this large system there are essentially three dimensions of interest: n, m, and r. The system of second-order equations is n-dimensional, but there are parameters in each equations. In the basal equations, the w species are not produced, but enter the dierential equations by feedback; these species form an (m−n)-dimensional parameter space. As well, the z populations are non-producers, but do modify the population equations of the excess producers. Although these variables now have complicated structure, we might view the excess producers as intermediate species, whereas the z populations are top species. These non-producers form an (r−n)-dimensional parameter space. Perhaps the most interesting case is when n = m = r. Then we have a Volterra– Hamilton system which has been well studied [3,4]. The detailed analysis of this general case will be taken up elsewhere; here, we will focus on the case of the examples above: n = m = r = 1. Example 3.1. We will take a brief aside to discuss the construction of the A and B matrices in the case that K is a matrix of constants. First, let us note that the rank of K is number of linearly independent columns, say, n. We can take a permutation matrix P so that the matrix K1 = KP has these as its rst n columns. Let A be the  matrix which implements the row reduction of K1 . Thus AK1 =

In ∗ 0 0

. Now, let C

be the matrix which implements the column reduction of the result. Thus,   In 0 = AK1 C = AK(PC): 0 0 We take B = (CP)−1 . This procedure of row and column reduction indicates a form of trophic equivalence in our model. For instance, suppose that one row of K is a multiple of another. This is stating that the two basal species in question are being consumed by each population with precisely the same ratio. Trophically, we would identify them as a single species as their densities always appear in the same proportion. This fact is re ected in the model; row reduction will eliminate one of the rows in favour of the other, replacing the duplicate with a row of zeroes. After the reduction takes place, one of the species

B. Lackey / Nonlinear Analysis 35 (1999) 37 – 57

47

(actually a weighted average of the two) is treated as a consumed species (represented by an x) whereas the other becomes a nonconsumed species (represented by a w). In general, if a collection of rows is linearly dependent, then a number of the rows are selected as representing the entire system trophically, and the remainder are relegated to nonconsumed species (which enter the dynamics through feedback only). The same argument works with the populations. If one column is a multiple of another, then the two represented populations consume the basal species in exactly the same proportion, so they are identi ed trophically. Mathematically, row reduction will not change the fact that the columns have a common ratio, hence when the column reduction is applied, one column is selected as a consumer (an N ), whereas the other becomes a column of zeros (corresponding to a z).

4. Equivalence and invariants Consider the problem of analysing a single second-order dierential equation, x00 + F(t; x; x0 ) = 0. Under a transformation x = x (x), then the dierential equation will change to x 00 + F = 0. The objects which do not change (or change in a predictable manner) are called invariants of the system. Naturally, the invariants will not only depend on the structure of the dierential equation, but also on the class of transformations that we consider. The study of invariants of dierential systems is called the KCC theory after its pioneers Kosambi, Cartan, and Chern. There are several dierent techniques used in KCC theory. Here we will use Cartan’s Method of Equivalence. Generally speaking, two systems are said to be equivalent (with respect to the allowable transformations) if there is a transformation which takes one system into the other. Naturally, the invariants of the systems must coincide. The method of equivalence is a procedure to determine the invariants of a system by focusing on conditions for systems to be equivalent. For systems of second-order dierential equations, the appropriate references are [6, 8, 12, 13]. The implementation of the method of equivalence is the exterior dierential calculus. See Appendix A for brief review and citations for more detailed treatises. Speci cally, we want to consider the space of variables to be (t; x; y) where y = x0 . The key to treating y as an independent variable, and maintaining its meaning, is to require that on any solution curve of the dierential equation, the dierential form dx − y dt vanishes. This idea naturally extends to describe the dierential equation itself as the curve for which dy + F(t; x; y) dt vanishes as well as the previous form. This description is necessary and sucient to characterize the solution, for if (x(t); y(t)) is a curve for which both forms vanish, then dx = y; dt dy + F(t; x; y) = 0; dt which is just another way to write the solution of the second-order equation.

48

B. Lackey / Nonlinear Analysis 35 (1999) 37 – 57

The class of transformations we will consider must preserve these forms in the following sense. A solution curve in the coordinates (t; x ; y)  must also have the form d x − y dt vanishing. So any transformation from (t; x; y) to (t; x ; y)  must preserve the vanishing of these forms. That is the transformation must act as follows: d x − y dt = p · (dx − y dt) where p is some non-vanishing function. Identically, if two equations are equivalent then the solution curves must coincide. Thus, dy + F dt = q · (dx − y dt) + r · (dy + F dt) : Here, r is non-vanishing and q is arbitrary. Thus, the class of transformations we consider are all those which act on the two forms above in the described way. We now apply the method of equivalence to the coframe {dt; dx − y dt; dy + F dt}. See Appendix B for more details of the method. First, we note that since we are only taking x = x (x), the coframe {dt; dx; dy} transforms as      dt 1 0 0 dt       d x  =  0 a 0   dx  : dy 0 b a dy We can translate between the two coframes by      dt 1 0 0 dt      dx − y dt −y 1 0 dx =     : dy + F dt F 0 1 dy Thus by the overdetermined algorithm (see Appendix B, Overdetermined Algorithm), we have all the components of the matrix       1 0 0 1 0 0 1 0 0 1 0 0       p 1 0 0  =  −py  0 p 0   −y 1 0   0 a a r 0 q r F 0 1 qy − rF qa − br 0 − ab2 1a a a2 are invariants under allowable transformations. We may therefore reduce the number of unknown parameters in our structure by setting these invariants to global constants whenever possible: say r = a, p = a, and q = b. Then our matrix of invariants is   1 0 0   1 0:  −ay by − aF 0 1 Note that y will vanish at equilibria and F at in ection points of our solution. Thus, since we must have a 6= 0, we cannot normalize our parameters further. Nonetheless,

B. Lackey / Nonlinear Analysis 35 (1999) 37 – 57

we have the invariants ay    dt 1     d x − y dt  =  0 0 dy + F dt

49

and by − aF. Our coframe now transforms by   0 0 dt   a 0   dx − y dt  : b a dy − F dt

Now, we turn to the structure equations of this coframe. As indicated in Appendix B our next step is to de ne the lifted one forms on (t; x; y; a; b) space      1 1 0 0 dt       2  =  0 a 0   dx − y dt  : 3 0 b a dy − F dt We take the exterior derivative of these forms, and relate the results to the original frame and the Maurer–Cartan forms of the transformation group, which is any matrix of forms of the form   0 0 0   0 ! 0 : 0  ! The rst of these is simple: d1 = 0. Now, d2 = da ∧ (dx − y dt) − a dy ∧ dt   da b 1 ∧ 2 − a 3 − 2 ∧ 1 = a  a a da − b 1 ∧ 2 + 1 ∧ 3 : = a So, we see our modi ed Maurer–Cartan form is ! = da=a − b 1 (mod 2 ). Finally, our last structural equation is   da b 1 ∧ 3 − 2 + db ∧ 2 + a dF ∧ dt − b dy ∧ dt d3 = a a a   b @F 1 ∧ 3 (mod 2 ): = ! ∧ 3 + 2 − a @y We may make this torsion vanish by selecting b = (a=2)@F=@y. This indicates that we need to modify our lifted coframe as     dt 1 0 0 1  dx − y dt     :  2  =  0 a 0     1 @F 3 0 0 a (dx − y dt) dy − F dt + 2 @y 

50

B. Lackey / Nonlinear Analysis 35 (1999) 37 – 57

Continuing our exterior dierentiation, we nd that our new lifted form satis es   1 @2 F da 1 @F d3 = − 1 −  2 ∧ 3 a 2 @y 2a @y 2 !  2 1 @2 F @F 1 @2 F 1 @2 F 1 @F − F − 1 ∧ 2 : + + y + 4 @y 2 @x@y 2 @y 2 2 @t@y @x Our remaining torsion does not depend on the group parameter a, therefore is an invariant function of (t; x; y) which we cannot normalize. As a note, the reduction of b has given us the two (group dependent) invariants v = ay and  = a(F − 12 (@F=@y)y). However, the modi ed Maurer–Cartan form has been completely determined, !=

da 1 @F 1 @2 F − 1 − 2 : a 2 @y 2a @y 2

As we cannot reduce the group parameters any further, and we have uniquely determined our Maurer–Cartan form, we prolongate (see Appendix B, Prolongation). Thus, our new space consists of variables (t; x; y; a) which has the coframe (1 ; 2 ; 3 ; !). As each of these forms is uniquely determined, there is no structure group, and we have an invariant coframe. The structure equations for this coframe are d1 d2 d3 d!

= 0; = ! ∧ 2 + 1 ∧ 3 ; = ! ∧ 3 + P 1 ∧ 3 ; = S 1 ∧ 2 + T 2 ∧ 3 :

The following are invariants of the system. v = ay;   1 @F y ;  =a F− 2 @y  2 1 @2 F @F 1 @2 F 1 @2 F 1 @F − F 2 + − ; + y P= 4 @y 2 @x@y 2 @y 2 @t@y @x   1 @2 F @3 F @3 F @3 F S = + F − y − ; 2a @y 2 @y3 @x@y 2 @t@y 2 1 @3 F T = 2 3: 2a @y 5. Qualitative structure of the invariants We begin our discussion with the last of the invariants T . Proposition 5.1. T = 0 if and only if G(x; N; t) is a quadratic in N .

B. Lackey / Nonlinear Analysis 35 (1999) 37 – 57

51

Proof. We have T=

1 @3 F −1 @3 G = 2 ; 2a2 @y3 2a @N 3

and the result follows. Although this result is simple, the fact that T is an invariant tells us more. If we start with an equation with a quadratic growth equation for the population, and rescale the basal species in any way whatsoever (even nonlinear scaling), then the new population growth equation will also be quadratic. This statement is rather important in modelling as well. Linear and quadratic terms are used to represent reproduction and growth limitations of a species. Higher-order terms are generally used to model social interactions within the population [3, 4]. Thus, the invariant T is a precise measure of social interaction, or a lack thereof. Moreover, as T is an invariant, two systems can be equivalent only if their social interactions have a constant positive ratio (which is then 2a2 ). Another interesting use of invariants is the following. Proposition 5.2. The zeroth-order invariants  and v are collinear (for each x; t) if and only if F = y + Ky 2 (where ; K may be functions of x; t). Proof. Fix a x and t and use Taylor’s theorem to write F = c0 + c1 y + k(y). The function k(y) satis es limy→0 (k(y)=y) = 0. Then    1 1 dk  = a c0 + c1 y + k − y : 2 2 dy Now, we have that  = 12 v. Thus,   1 1 1 dk y = y: c0 + c1 y + k − 2 2 dy 2 Evaluating this equation at y = 0 show that c0 = 0. Then dividing the equation by y and letting y → 0 yields c1 =  and k = 12 (dk=dy)y. Yet, as k(0) = 0, it is easily recognized that the solution to this dierential equation is k(y) = Ky 2 . Thus, F = y + Ky 2 . Conversely, if F is of this form, it is easy to verify that  = 12 v. An interesting exercise is to extend this argument for Volterra–Hamilton systems. One can show that  j = 12 v j for several species implies not only quadratic dynamics without a source, but also that the linear growth terms all have the same rate. That is, we have that F j = y j + Kklj y k y l . This is called the presymbiotic condition, and is used heavily in modelling growth and development of colonial animals [3, 4]. Finally, we would like to cite a stability result due to Kosambi. If one introduces a covariant derivative along solution curves using the connection form ! above, then one discovers that the P invariant plays the same role in second-order dierential equations as the ag curvature in Riemannian and Finsler geometry. Thus sign of this invariant

52

B. Lackey / Nonlinear Analysis 35 (1999) 37 – 57

(or the eigenvalues in the Hamilton–Volterra case) determine if solution curves are converging or diverging. Thus, we have the following result. Theorem 5.3. (Kosambi [13]). If P¿0; then the solutions of the second-order dierential equation are Lyaponov stable.

6. Conclusions We have seen in several examples above, we may model a variety of ecological situations with essentially the same mathematical construction. Traditionally, single species and multispecies models are dynamical systems ( rst-order dierential equations). Yet, there is good evidence that almost every ecological system relies on complex biochemical interactions [11,21]. As well, even very controlled single species experiments show that many population control their growth through allelochemic production [20]. The concept of trophic systems, extends in a natural way. If we wish to model these systems, we include the chemicals produced by the populations as species as well. Some basic assumptions about the means of production allow us to rewrite the systems of rst-order equations as systems of second-order equations. Mathematically, the treatment of second-order dierential equations is much more complicated that of rst-order equations and, as is common, the theory is much richer. One can show [15] that away from singularities, all rst-order dierential equations are equivalent, and thus have no invariants. Yet for second-order equations, we do have invariants, as well as a unique connection [8] and thus an entire geometry [4]. Mathematical biology (or any mathematical science) is a constant exchange of information from scientist to mathematician. Given any phenomenon (and enough time), one can cook up some mathematical construction that models it well enough. Yet, the real goal of modelling is to present a mathematical object that models a range of natural behaviour, and is tractable enough to draw meaningful conclusions. We present the use of second-order dierential equations here for this reason. Biologists recognize that production and consumption of other species, of chemical, of minerals, or even of available space, is ubiquitous in all theories.

Appendix A. Exterior dierential calculus Perhaps the rst systematic use of dierential one-forms was by Pfa in the early 1800s to study partial dierential equations. Frobenius introduced the exterior derivative (under the name bilinear covariant) in 1877. By the early 1920s, Cartan virtually rewrote dierential geometry using the exterior calculus; and since then nearly every text in dierential geometry or Lie groups dedicates a sizable proportion to the exterior dierential calculus. To the non-geometer, perhaps the most accessible treatise on differential forms is Flanders’ classic text [9]. Another readable text, more geared toward applied mathematics and physics, is [5]. As well [22] give an excellent (and readable)

B. Lackey / Nonlinear Analysis 35 (1999) 37 – 57

53

account of the use of dierential forms in Lie groups, especially the Maurer–Cartan forms used heavily here. A.1. Vector elds versus dierential forms The rst step to understanding dierential geometry is through the simple statement directions ↔ directional derivative. If we have coordinates (x1 ; : : : ; xn ) and in these coordinates a vector eld v = (v1 ; : : : ; vn ), then the familiar directional derivative formula v[f] = v j

@f @ = v j j (f) @x j @x

identi es the eld v with a partial dierential operator v j (@=@x j ). By comparing components, we see that a coordinate system, (x j ) induces a basis for vectors {@=@x j }. Suppose we have a dierentiable function y =  (x j ), then the chain rule induces a linear transformation by ∗ (@=@x j ) = (@y =@x j )(@=@y ). This map is called the derivative (or Jacobian or push-forward) along . Dierential forms (often called one-forms) are the natural duals to vector elds. To each @=@x j we associate something denoted dx j . This object is determined by a pairing dx j (@=@xk ) = kj which we extend linearly to all vector elds. We then take {dx1 ; : : : ; dxn } to be a basis of the one-forms in the same way that the @=@x j form a basis for vector elds. Thus dierential forms are linear functions on vector elds, yielding dierentiable functions. We can de ne a transformation on forms akin to the derivative map above by ∗ !(v) = !(∗ v). In coordinates this appears as ∗ (dy ) = (@y =@x j )dx j , so ∗ maps from form on the y-space back onto form on the x-space, justifying its name: the pull-back. One construction we will use frequently is pull-backs of coframes. Given an open subset U ⊆ Rn with coordinates (x j ), at each point p ∈ U we can consider the span of {dx j }. This is an n-dimensional vector space which we call the cotangent space at p. A collection of one-forms {!1 ; : : : ; !n } is a coframing of U if these form span the cotangent space at each point of U . If  : U → U is a dieomorphism, then its Jacobian is invertible everywhere. Thus every coframe of U , say {! j }, gets pulled-back to a coframe of U by ! j = ∗ ! j . A.2. Exterior algebra On any vector space, one may construct an exterior algebra. We will show one way to do this for dierential forms. To any two one-forms, ! and , we de ne an object (called a two-form) by the formula ! ∧ (v; w) = !(v)(w) − !(w)(v). Thus, a two-form is a bilinear function on pairs of vector elds (again yielding a function) with the special property that (v; w) = −(w; v). This property is equivalent to saying that (v; v) = 0 for all vector elds v. Higher-order forms are de ned analogously. A k-form is a multilinear function on k vector elds (yielding a function) which vanishes whenever any two of the vector elds coincide. Forms of any degree share the pull-back with one-forms. If is a k-form, then we de ne ∗ (v1 ; : : : ; vk ) = (∗ v1 ; : : : ; ∗ vk ).

54

B. Lackey / Nonlinear Analysis 35 (1999) 37 – 57

A.3. Exterior calculus A special property of forms, which is not shared with vector elds, is a calculus. Given a function f, we can de ne the exterior derivative of f as the one-form df = (@f=@x j )dx j . Similarly, if ! = ak dxk then we de ne the exterior derivative of this one-form to be the two form d! = (@ak =@x j )dx j ∧ dxk . We can de ne the exterior derivative of higher-order forms in a similar way. The exterior derivative has many useful properties; we list three that we use consistently here. 1. d(f!) = df ∧ ! + f d! for a function f and form !; 2. d(! ∧ ) = d! ∧  + (−1)p ! ∧ d for a p-form ! and any form ; 3. d(∗ !) = ∗ (d!) for any form ! and smooth map . Appendix B. Cartan’s method of equivalence  Cartan’s equivalence problem is the following: given two coframings {! j } and E. {! j } of open subsets U and U of Rn , determine when there is a dieomorphism  : U → U so that ∗ ! j = ! j . This problem was solved completely by Cartan in [7] (for an English translation of the appropriate sections of this paper, see the Appendix of [10]). Cartan called this the restricted equivalence problem; today, we call the equivalence problem of an e-structure. In application, our coframe is only determined up to some transformation group. For instance, in our second-order dierential equation, we may use any coframe which describes the same solution curve. We have seen above that any two of these coframes dier by a transformation of the form   1 0 0   0 a 0 0 b a where a is a nonvanishing function of (t; x; y) and b is an arbitrary such function. If we call G the group of these matrices, then we say the equivalence problem of a second-order dierential equation has a G-structure. The appropriate question to ask for a G-structure is: does there exist a dieomorphism so that ∗ ! j = gkj ! k for some suitable selected groups element g(x). The Method of Equivalence is Cartan’s method to eliminate parameters of the transformation group, reducing it to a single element, {e}. That is, the method reduces the G-structure to an e-structure, to which we may apply Cartan’s solution. Generally speaking, the method is to compute invariants (which will be functions of the group parameters). If we then demand that an invariant be some suitable value, this allows us to express some group parameters in terms of others, eectively reducing the structure group to a subgroup. As we are really interested in the invariants, we will concentrate on the invariants that arise as we reduce to an e-structure and not pursue Cartan’s answer to the equivalence problem.

B. Lackey / Nonlinear Analysis 35 (1999) 37 – 57

55

For detailed modern treatments of Cartan’s Method of Equivalence, see [10, 15]. B.1. Lifted coframes The method works almost exclusively on a lifted coframe. That is, if G is our transformation group and {! j } our coframe, then we de ne  j = gkj ! k for each g ∈ G. These one-forms live on the space U × G, and although they do not span the cotangent space, we call then the lifted coframe. The utility of this idea follows from the following lemma. Lemma B.1. (Gardner [10]). There exists a dieomorphism  : U → U such that ˆ : U ×G → ∗ ! j = gkj ! k for some g(x) if and only if there is a dieomorphism  ∗ j j ˆ  =  . U × G such that  It might appear that we have solved our problem: we now have an equivalence problem in terms of the ’s which have no transformation group involved. Yet, recall that { j } does not form a coframe over the now larger domain U × G. The transformation group is still present; it appears in the indeterminacy of completing these forms to a coframe. Nonetheless, we have translated our problem into a more symmetric form. Equivalence now has the form g! = g! (suppressing the ∗ ). As this form is symmetric in the barred and unbarred objects, it suces to analyse one such lifted coframe; the other will have analogous equations. B.2. Overdetermined algorithm Suppose we are given two competing coframes, {! j } and { j }, of U . Let us denote the structure groups of these two coframes as G and H ; and since these are coframes, we may write ! = A for some matrix A (all these objects have corresponding barred counterparts). If we have an equivalence, then g! = g! and h  = h, and thus  g! = g! = gA = gAh−1 h = gAh−1 h : −1  we must Yet, g! = gAh h .  As g! and h  form a coframe for any selection of g;  h, −1 −1 −1 have that gAh = gAh . That is, the components of gAh are all invariants! Clearly, these terms will depend on our group parameters. We may set them to some appropriate value, which determines relationships between the group parameters. This reduces the number of free variables in the group, thus reducing it to a subgroup. As well as reducing the size of the structure group, we may as well determine some invariants under allowable transformation. These ‘zeroth-order’ invariants are often useful in describing some qualitative structure of the system in question. For instance, in our second-order dierential equation above, the collinearity of the zerothorder invariants describes logistic growth. In the multispecies case, one can check that the symbiotic condition also follows from this criterion, [3].

56

B. Lackey / Nonlinear Analysis 35 (1999) 37 – 57

B.3. Structure equations Now, we consider the true power of the exterior dierential calculus. If we have an equivalence, ∗  = , then we may exterior dierentiate both sides of this equation.  = ∗ d .  Thus and equivalence of the lifted coframe serves Thus we have d = d(∗ ) as an equivalence of it exterior derivative as well. We are lead to compute these equations, which are called the structure equations of the coframing. We compute, d j = dgkj ∧ !k + gkj d!k : Now, as {! j } is a coframe of U , we may express d! j in terms of ! j ∧ !k . Thus we can write the above equation as d j = kj ∧  k + Tklj  k ∧ l ; where kj = dglj (g−1 )lk are the Maurer–Cartan forms of the structure group G. These forms for a coframing, for the transformation group G. Yet, there are many other choices of frames which can serve the same role. There are two important (and competing) ideas here. We wish to determine the matrix of forms as much as possible, so we ‘absorb’ as many of the terms in the second term in the above equation into the rst. This has the eect of replacing the Maurer–Cartan forms with another collection of forms; if these new forms are to be a coframing of G, though, they must satisfy the structure equation of the transformation group. Thus, we may only absorb terms in a way compatible with the group. For instance, above we found Maurer–Cartan forms of the form   0 0 0   0 ! 0 : 0  ! Thus, when we absorb terms, we are required to absorb them into the 2-2 and 3-3 entries equally, into the 3-2 entry arbitrarily, and into no other places.  Thus both d and  As opposed to this, we have seen that d = d∗  = ∗ (d ). are preserved under equivalence. Now, is exactly the same in both barred and unbarred coframes as the structure groups are identical, so the ‘torsion’ terms Tklj must be invariants! Usually, these terms depend on the group parameters, so we may use them to reduce the structure group as always. It would appear that leaving as much information in the T terms as possible allows us to achieve our goal (of eliminating the structure group) more quickly. Yet, recall that by constructing the lifted coframe, we have reformulated the structure group as an indeterminacy in completing the lifted coframe into an honest coframe of U × G. Yet this indeterminacy is apparent in the structure equations: it is the choice of . So if we can completely determine , then we use it to complete the coframe on U × G for which we have an e-structure. Notice that in our second-order dierential equation above, we have done both. We rst used the torsion to reduce the structure group by nding that b = a 12 (@F=@y). Then

B. Lackey / Nonlinear Analysis 35 (1999) 37 – 57

57

with the smaller structure group, we absorbed terms until ! was uniquely determined, thus we achieving an e-structure on the enlarged space (t; x; y; a). References [1] P.L. Antonelli, Mathematical Essays on Growth and the Emergence of Form, Section 2D, Univ. Alberta Press, 1985. [2] P.L. Antonelli, P. Auger, R. Bradbury, Corals and star sh waves on the great barrier reef: analytical trophodynamics and 2-patch aggregation methods, Nonlinear Times and Digest 2 (1995) 261–289, or Math. Comput. Modelling, Special issue on Aggregation Methods in Population Dynamics, to appear. [3] P.L. Antonelli, R.H. Bradbury, Volterra–Hamilton Models in the Ecology and Evolution of Colonial Organisms, World Scienti c, Singapore, 1996. [4] P.L. Antonelli, R.S. Ingarden, M. Matsumoto, The Theory of Sprays and Finsler Spaces with Applications in Physics and Biology, D. Reidel and Kluwer Academic Press, 1993. [5] W.L. Burke, Applied Dierential Geometry, Cambridge, 1985.  Cartan, Observations sur le memoir precedent, Math. Zeitschrift 37 (1933) 619– 622. [6] E.  Cartan, Les sous-groupes des groupes continus de transformations, Ann. de l’Ecole  [7] E. Normale (1908) 57–194. [8] S.S. Chern, Sur la geometrie d’un systeme d’equations deierentialles du second ordre, Bull. Sci. Math. 63 (1939) 206– 212. [9] H. Flanders, Dierential Forms with Applications to the Physical Sciences, Dover, New York, 1989. [10] R.B. Gardner, The Method of Equivalence and Its Applications, CBMS 58, SIAM, Philadelphia, PA, 1989. [11] J.B. Harborne, Introduction to Ecological Biochemistry, Academic Press, New York, 1988. [12] D.D. Kosambi, Parallelism and Path-Spaces, Math. Zietschrift 37 (1933) 608– 618. [13] D.D. Kosambi, Systems of dierential equations of the second order, Quart. J. Math. (Oxford) 6 (1935) 1–12. [14] D.S. Mulyk, Allelochemically-based interspeci c interactions, Math. Comput. Modelling 13 (6) (1990), In: P.L. Antonelli (Ed.), Special issue on Proc. Internat. Workshop on Population Dynamics of Outbreaks, 1990, pp. 31–34. [15] P.J. Olver, Equivalence, Invariants, and Symmetry, Cambridge, 1995. [16] R. Pratt, Studies on Chlorella vulgaris, V, Some properties of the growth inhibitor formed by Chlorella cells, Amer J. Bot. 27 (1940) 52–56. [17] R. Pratt, Studies on Chlorella vulgaris, IX, In uence on the growth of Chlorella of continuous removal of chlorellin form the solution, Amer J. Bot. 31 (1944) 418– 421. [18] R. Pratt, Studies on Chlorella vulgaris, XI, Relation between surface tension and accumulation of chlorellin, Amer J. Bot. 35 (1948) 634 – 637. [19] D.F. Rhoades, Evolution of plant chemical defense against herbivores, In: G.A. Rosenthal, D.H. Janzen (Ed.), Herbivores, Their Interaction with Secondary Plant Metabolites, Academic Press, New York, 1979, pp. 3–54. [20] E.L. Rice, Allelopathy, Second Edition, Academic Press, 1984. [21] G.A. Rosenthal, D.H. Janzen, Herbivores, Their Interaction with Secondary Plant Metabolites, Academic Press, New York, 1979. [22] D.H. Sattinger, O.L. Weaver, Lie Groups and Algebras with Applications to Physics, Geometry, and Mechanics, Springer, Berlin, 1986. [23] M. Spivak, Calculus on Manifolds, W.A. Benjamin, New York, 1965. [24] F.M. Williams, Dynamics of microbial populations, Systems Anal. Simulation Ecology I (1971) 197– 267.