- Email: [email protected]

Division of Muthematics

and St&sties,

Canberra, A.C.T.,

Australia

Received ? January 1980; revised 17 July 1981

ABSTRACT A set of N coupled nonlinear differential equations for a plantation of N competing plants is analysed. They describe the observed increasing variability with age resulting from competition between plants with given initial variability. The observed variation in skewness from zero to negative, then to positive, with age is explained in terms of the linearized equations with random initial conditions and the nonlinear equations with periodic initial conditions. Negative correlation between neighboring equations are simply the backward diffusion equation

1.

plants is predicted. The linearized in 2-dimensional discrete space.

INTRODUCTION

In a recent paper [I] some new equations describing the growth of a plantation were introduced. They consist of N coupled nonlinear differential equations for a stand of N plants. Numerical simulation of the equations produced DBH (diameter at breast height) distributions whose skewness reversed its sign with time in close agreement with data on Tagetus patufa [2] and Pinus radiata [6]. The distributions also become bimodal when competition was intense, in agreement with the data of [2]. In the present paper we begin a mathematical analysis of the new equations, providing some insight into the first of these striking competition phenomena. The analysis produces a spatial pattern like that observed by several authors (see Sections 3 and 4). The equations in [I] can be written

(1.1) where A,(t) is the area at time t of a disc centered at x, S, is the area of a subset of A, obtained by partitioning the union of the discs, as described below, and u and K are constants for time scales small enough or large enough to permit the neglect of seasonal variations in climate. There is one MATHEMATICAL

BIOSCIENCES

59117-32 (1982)

QElsevier Science Publishing Co., Inc., 1982 52 Vanderbilt Ave., New York. NY 10017

17 0025-5564/82/030017+

16$02.75

18

D. J. GATES

equation for every point x in some set A. Each x specifies the location of a plant, and A is taken here as a square lattice of points. A disc, called a zone of influence (ZOI), represents a region over which the isolated plant would take up resources (light, water, nutrients), and the reduced area S, is the corresponding area for a plant in competition with its neighbors, and depends on the set {A}X of ZOI areas of those neighbors. The equations in (1.1) were derived from the simple energy conservation equations

dW, _

--k(E,-YW,), dt where W, is the weight of the plant at x, E, is the energy it absorbs per unit time, -yW, is the energy it requires for maintenance per unit time, y is a constant measuring the maintenance requirement per unit weight of plant, and k is a constant giving the conversion efficiency of energy into weight. The equation would not apply if E, were excessively large, but usually this is not the case. Then (1.1) is derived in [I] from empirical relationships expressing W, and E, in terms of A, and S,. The constants K and Q in (1.1) are proportional to k and ky respectively, the proportionality constants depending on plant shape and supply rate of energy. By definition the basal area of B, of the stem (trunk) of a plant is given in the model by

(1.2)

4 =A,/P,

where p is a constant whose value we shall not need. Given ~1,this determines the Ax’s completely. Other directly measurable quantities for plants are also taken as functions of A, and S,. The area S, is defined as the area of the set of points r satisfying

and (r-xla-

R", G (r-yl”-

R,

forall

yfx,

(1.3)

where 1~ (YG cc is a parameter depending on plant shape and structure. The equality in (1.3) is just the common boundary curve between S, and SY. If a = 1, say, this curve is a hyperbola and the area of intersection of the discs at x and y is shared fairly evenly. If a =2, the curve is the common chord of these two discs, and the larger disc gets the larger share of the intersection. If (r = 00, the larger disc gets all of the intersection. In general, (Ymeasures the degree of domination of larger plants over smaller ones. The relations (1.3)

GROWTH AND COMPETITION IN PLANTATIONS

19

imply that the S, are a partition of the area of the union of the discs, and this corresponds to a partitioning or sharing of the available resources among the competing plants. The relations were derived from a fundamental viewpoint in 131. Some elementary properties of (1.1) can be obtained immediately. Since S, G Ax, there is no growth unless K10,

(1.4)

i.e., all plants have sufficient resources for maintenance, as discussed in [I]. If total cowr, i.e. canopy closure and root system closure, has been achieved, then (assuming periodic boundary conditions) x

S,=

Na2,

(1.9

XEA

where Q is the spacing in the square lattice. Then the mean zone area

satisfies -

$+oA=na2,

(1.7)

which has solution a(r)=$+[a(o)-+rJ’.

(1.8)

Thus A(f)-$

as

t-too,

(1.9)

which implies that, ultimately, a balance is achieved between mean supply ~a* of resources and mean maintenance requirements ox (see [ 11). There is a net rate of growth (dA/dt > 0) if tca2 > ox,

and if this occurs beyond canopy closure, where A> sa2/2 K >

a condition

which we henceforth

7W/2.

assume.

(1.10) at least, then (1.11)

20 2.

D.

EXPLICIT

EQUATION

GATES

FOR a = 2

The sumulations mentioned above are time consuming and give limited insight into the time development. The immediate hindrance to a deeper analysis of this process is the complexity of the function S,({A},) which gives the reduced areas in terms of the disc areas {A}, of plants near the plant at x. For general a this function involves integrals which cannot be evaluated in closed form. However, in the case where (Y= 2 and all area is covered (e.g. a closed canopy), a dramatic simplification occurs, especially if the discs are small enough so that only the eight nearest neighbors of a given plant can compete with it. In this case the reduced area S, is polygon with at most eight sides, or more particularly, a rectangle with some of its corners cut off at 45”, as illustrated in Figure 1. The basis vectors in the 4 lattice directions are a,, where la, (= a and a is the lattice spacing. The area S’ of the rectangle is given by

where D, is the distance from x of the common chord between the discs at x and x+a,. One finds that

D, =

R2, - R& 2a

+ a2

(2.2)

’

FIG. I. The reduced area .S, of the plant at x is a rectangle if a = 2 and total cover is achieved.

with some comers

missing

21

GROWTH AND COMPETITION IN PLANTATIONS so

that

If the disc centered at x+a, + a2 is sufficiently large, a comer is lost from the upper right hand of the rectangle (Figure 1). It has side Rz+a’-R:,,

Ri+a,+az+a2-Rifa,

‘+

2a

2a

-a

=~(A,+A,+.,+,*-A,+,,-A,..I)

(2.4) if and only if this is positive. Thus the loss in area from this comer is

Combining (2.3) with (2.5) gives an expression substitution in the differential equations (1.1) gives ldA,

for &((A),)

which on

a

---+;Ax=a2--&AA,

where a, is the difference operator

a,j(x)=

f(x+da,)-f(x-fa,)9

(2.7)

(2.8) Gxj= max(2f,O)

(2.9)

and

h=a:+a;

(2.10)

is the discrete Laplacian. Note that only second differences a,? and a,a2, so defined, have meaning for the AX’s, The equations (2.6) are the basis of the following analysis.

D. J. GATES

22

Since total cover is assumed in these equations, (1.5) applies [as may be verified by direct substitution from (2.6)], and (1.8) consequently holds. It should be noted however that further terms may need to be added to the right side of (2.6) for large t if more distant neighbors interact with each other. Terms such as (AA,.)* and A*A, then appear. The terms linear in A, in (2.6) are essentially discrete space versions of the terms of the diffusion equation with time reversed (see Section 3). The quadratic terms are essentially discrete versions of those that arise in the equations of small deformations of elastic plates [4], as are the further terms mentioned in the preceeding paragraph. Further mechanical analogies are exploited in Section 5. The explicit form of S, in (2.6) enables it to be calculated more accurately and efficiently than it was in [I], where simulations of the equations were performed. This approach is recommended if more thorough simulations are attempted, although it is confined to the (Y=2 case. 3.

LINEAR

APPROXIMATION

Suppose the canopy is closed, so that (2.6) applies, but competition not yet produced large differences in plant sizes; say

I4 -Ayl+~

,%+A,

for allx,y,

has

(3.1)

where Ec 1. Then the terms in the bracket [ ] in (2.6) are of order .? and may be neglected in comparison with the remaining terms of order E. The resulting equations may be written

4 _ -$u,, -_

(3.2)

dt

where (3.3) Equation (3.2) is a backward-evolving with diffusion coefficient D=

diffusion

K/V.

equation

in discrete space (3.4)

Consequently, initial spatial variation in the U,‘S tends to be enhanced as time proceeds, rather than diffusing away. This was very evident in the simulation in [I]. The model in [I] [Equation (2.9)] implies that K is proportional to the rate of supply of plant resources and to the conversion efficiency of resource into plant matter. These factors therefore contribute through D to the rate of increase of variability in the plantation.

GROWTH AND COMPETITION IN PLANTATIONS

Equation

23

(3.2) has the general solution

where s = (s,, s,), the S, are integers (positive and negative), and I,, is the modified Bessel function of order n [5]. The summation is over all integer pairs, and the solution is applicable for initial conditions u, (0) which ensure convergence of the summation-for example th.e periodic initial conditions used in [I]. Note that our zero of time does not correspond to the time of planting, but to an arbitrary time that exceeds the time of canopy closure, when the equations (2.6) become applicable. Competition effects can be illustrated by following the growth of a plant of one size in a plantation of otherwise identical plants of a different size:

u,(O) =

U,

for

x=10,

i rl;

for

x+0,

(3.6)

(i=1,2).

(3.7)

where IU, -U,lKu, Then (3.5) becomes u,(t)=Vz+(U,-Lr,)(-1)‘+ke2~D’I,(Dt)l~(D~),

(3.8)

where x=(j,k)a, and we have used the properties Z_,(x)=Z,(x),Z,(-X) = (- l)“Z,,( x), and (4.2). All the functions I,( ZA) are positive and increasing for t >O, so the deviations from U, in (3.8) grow at least as fast as e2Dr. These deviations decrease as Ix/ increases, tending to 0 as 1x1--t co. They have the same sign as U, - U,if j + k is even, and the opposite sign if y + k is odd, so that the signs have a chessboard pattern. This pattern of increasing deviation is broadly what one would expect on physical grounds. If U,> U2, then the plant at 0 has an advantage which it continues to exploit. The growth of its four nearest neighbors is suppressed by competition with it, while the four next nearest neighbors (the diagonal neighbors) gain some growth advantage because they suffer less competition than normal (i.e. U,)from the nearest neighbors; and so on for more distant neighbors. The assumption (3.7) is important here, for if U,Z+U,,then (3.2) and hence (3.8) are not valid, and we would expect more than the four nearest neighbors of 0 to be suppressed. Negative correlations between the weights of neighboring plants in row crops were observed by Yoda et al. [8], between DBHs of neighboring trees in Ponderosa Pine forests by Cooper [9], and between yields of tea plants (Camellia sinensis and the yields of their four nearest neighbors on a square

24

D. J. GATES

lattice by Cannel1 et al. [lo]; further references given in Table 6 of Ford and Diggle [ 111. 4.

RANDOM

INITIAL

to related phenomena

are

CONDITIONS

Now suppose, as in [l], that the initial conditions are random, the ~~(0)‘s being identically distributed with mean p, Then the U, (t)‘s have mean

= II,

(4.1)

using the identity

2

qx>=ex.

/= -00

(4.2)

Thus (3.3) yields

EA,(t)=$+P (P/g9

(4.3)

where cc,,,= EA.(O). This result can be deduced directly from the differential equations (l.l), subject only to the canopy closure condition ES,(t)= u2, and not requiring (Y=2 or any linearization of the equations. Equation (4.3) shows that EA,(t)+

as

t-+00,

(4.4)

so that an overall balance is ultimately achieved between average maintenance requirement oEA, and average supply KU’ of resource (see [l]). This is essentially the same result as (1.9). Returning to (3.2), we have

4 =

DE

4u:-

x u,u,+,, i=l

=2DE(

u; + I(;+,,-2u,u,+,,)

=2DE{

(ux - ~,+a~)~}

20,

(4.5)

GROWTH AND COMPETITION IN PLANTATIONS

25

where we have used the translation invariance of the distribution The variance is therefore nondecreasing as t increases. We also have

of the u,‘s.

(4.6) If nearest neighbors are negatively [see also (4.20) below], then

correlated

initially,

as suggested by (3.8)

Var uX( t)T e4D’Var u,(O)

(4.7)

for small t. Then (3.3) implies that -zO’var U,( t)

VarA,(t)=e

> e(4D-22a)‘VarAX(0)

(4.8)

for small t. Recalling from (1.11) that D > a/2, we deduce that VarA,( t) also increases with r. This provides a theoretical explanation of the increasing variability induced by competition reported in [l] and widely elsewhere. Such an increase is perhaps fairly clear on qualitative grounds also. For the general central moments we similarly obtain, putting v, = u, - p,

(4.9) If the u, are initially

correlated

thus:

V xtfl,

(4.10)

=au,+e,,

where the E,‘S have mean zero and are independent

of the u,‘s, then

E{4r)l=

eZnD’(1- a)E{ o,(O)“}

(4.11)

M,(r)=e

n(2D--o)t(1- cY)M”(O)

(4.12)

and

for small r, where M, is the nth central moment of the A,‘s. Since we expect initially (i.e. at the beginning of canopy closure) that correlations will be small or negative, we deduce that M,( 1)/M,,(O) is increasing for small r. This further demonstrates the increasing variability induced by competition.

26

D. J. GATES

These results are valid only to first order in t. One can obtain results to higher order if one assumes that the A,(O)‘s are independent. This assumption is not really justified, because one would expect competition prior to canopy closure to result in some degree of correlation. The results are at least indicative. From (3.5) we deduce that the n th cumulant of the U,(C) is

IC~(t) = Kn(0)e2nD' i /=

x

I,(-Dt)”

-00

(4.13) i-

To calculate the variance, we note [5, p. 376, 9.6.341 that (4.14)

e

and [5, p. 376, 9.6.161 257 dsez'coss=

/0

(4.15)

27rZ,(2t),

so that ; js

Consequently

(4.16)

Z,(t)‘=2z,(2t)-Z,(t)2. -00

from (4.13) and (4.16)

VarA,( t)= {VarA,(O)}

e(4O- *o)r{21,(2Dt)-

I,( Dt)2}2.

(4.17)

Since 2Z,(2 Dt) - I,( DT)2 is increasing, we again have an exponential growth in variance, now valid for somewhat larger 1, but subject to the restriction (3.1). The skewness SK(t) of u,( 1) and of AX(t) obtained from (4.13) is given by

SK(t)= SK(o)

which is a decreasing

function

c

jz-oo

I,(+)’ (4.18)

{210(2D~)-II,,(Dt)2}3’

of t. Its small t behavior is I-3D2r2+O(r3).

(4.19)

27

GROWTH AND COMPETITION IN PLANTATIONS

I

[.\

\

,1’

t-0 I \

\

SK(O)-

t

l.._J” \ canopy

closure

FIG. 2. Sketch of prediction (-) by (4.19) and observations and simulations (-----) of variation of skewness SK(r) with time t.

This behavior is consistent with the observations and simulations from [I] as indicated in Figure 2. Correlations between different plants, again assuming initial independence, are readily obtained from (3.5):

zJ2Dr)zk(2Dr)

Corr(A,,A,+,)=(-l)J+k

{21,(2Dt)-

z(ot)}’

(4.20)

’

where x=( j, k)a. As in (39, the correlation again has a chessboard pattern of signs, where nearest neighbors are negatively correlated, next neighbors positively correlated, and so on in an alternating pattern, The small t behavior of (4.20) is (-Dt)‘+“/(j!k!), which shows the decrease in magnitude increases. 5.

PERIODIC

PATTERNS

of the correlations

(4.21) as the separation

OF PLANTS

The main phenomenon unexplained by the analysis here and in [l] is the increase in skewness through to positive values beyond canopy closure (Figure 2). The initial drop in skewness to negative values was explained in [7]. The reversal of this trend for small t beyond canopy closure is explained by (4.18), but under the questionable assumption of initial independence. We know from the simulations in [l] that the full differential equations (2.6) correctly model the complete behavior of the skewness, while the linearized equations (3.2) fail to explain the later behavior. We now show, for the rather artificial case of a periodic pattern of plant sizes, that this large t behavior of the skewness can be derived analytically.

28

D. J. GATES

l1

l2

l1

i

‘2

‘3

‘2

‘3

l1

‘2

l1

‘2

l2

-3

l2

‘3

FIG. 3. Periodic pattern of three plant sizes which yields analytically a variation in skewness (Figure 5) similar to the observations and simulations based on random initial conditions.

We assume that initially there are plants of 3 different sizes, with ZOIs of areas A,, A *, and A3 arranged in the pattern of Figure 3. The three different sizes are sufficient to produce a nonzero skewness in the AX’s or R,‘s. The mean of the AI’s is p=:(A,+2A,+A3), and the central moments

(5.1)

(n = 2,3) (5.2)

are conveniently

expressed in terms of the variables W= A, -2A,

+ A,

(5.3)

and Z=A,-A,,

(5.4)

giving w= +2z2 CL==

,(j

9

w(wz+12z2) CL3 =

64

’

(5.5)

Thus the sign of the skewness is just the sign of W, as is clear from (5.3). In fact W and Z are alternative measures of the skewness and standard deviation, respectively. The equations (2.6) take different forms depending on whether W>O or WCO. These two cases, together with the W=O case, are illustrated in

.......... ..... 583 29

GROWTH AND COMPETITION IN PLANTATIONS

w

w=o

w>o

FIG. 4. Three different patterns (br two different and one marginal) of reduced areas corresponding to the periodic pattern of Figure 3.

Figure 4. To demonstrate c in Figure 4, yielding

these configurations,

we use (2.2) to compute b and

b-(a-c)=W/(2sa).

In the case where W>O initially,

(5.6)

the equations

(2.6) reduce to

&I_ K ---(2(A,-A,+ra2)2-(A,+A,-2A2)2], dt

2a2a2

(5.7) a% 3_ 6 ---(2(A3-A2+na2)2-(A,+Aj-2A2,2),

dt

2s2a2

where the maintenance requirements are ignored (u = 0) as in the simulations in [ 11. Subtracting the last equation from the first gives $=uZ(W+2/3),

(5.8)

where v=K/(7ra)2 and /3= 7ra2. Adding

the first and last of the equations

(5.9) (5.7) and subtracting

the second

D. J. GATES

30 gives ~~‘y(z’-W’+4pW) Equations (5.8) equations which A,‘s themselves combined with a

for

WaO.

(5.10)

and (5.10) are a pair of coupled, nonlinear differential completely describe the time development for W z 0. The can be obtained from the solutions to these equations third equation (5.11)

$(A,+2A2+A,)=4ra2,

with the known solution (1.8). For W

Equation

w2+4pw) conclusions

for

(5.8) is unchanged,

WGO.

may be immediately

while

(5.12) drawn

from

w+zp=g>o,

the

(5.13)

as long as plant 1 or 3 survives. Then (5.8) implies that Z is increasing, so that the spread A, -A, widens with 1. If 0~ W(O)<4/?, then from (5.10), W(O)>O. If W(O)>0 andA,>A,>A, at t=O, then Z> W, so that W(O)> 4apW(O). If Z(O)>2/3 and W(O)> -2p, then from (5.12) W(o)>y(w+2p)2>o.

(5.14)

Thus W(t) is initially increasing under a variety of conditions, demonstrates the observed increase in skewness after canopy closure. Eliminating Z from the equations gives W=16a/U++8a2W( W=8a(

W-4/3)(

W+2P)W-8a2(

W+Zp)

W+Zp)(

W+4p)

for

WaO,

which

(5.15)

for WC 0. (5.16)

These are equations of nonlinear oscillators with an additional force term (negative damping) proportional to the velocity. They can exhibit a variety of behavior, but the following are most relevant here: Suppose A, > A, > A, at t = 0. Then (a) W(t) increases indefinitely 0 and Z(0) > Z/3;

and exponentially

if W(0) > 0 or if W(0) <

GROWTH

AND COMPETITION

31

IN PLANTATIONS

2-

-30

5I

10 I

t

15 I

20,

25 I

30

FIG. 5. Variation of skewness SK( 1) and related function W(t) for the periodic of Figure 3 with initial conditions W(0) = - 1.28 and Z(0) = 1.5/3.

(b) if W(0) < 0 and Z(0) < 2/3 and 2/I-W(0) W(t) decreases initially to a minimum > -2 exponentially.

pattern

is sufficiently large, then and thenceforth increases

These properties follow from the nature of the conservative forces in (5.15) and (5.16), whose potential has maxima at W= -2p and W=4p. The former maximum causes a reversal of direction under the initial conditions in (b), but the latter is overtaken provided A, > A 2 > A, initially. An example of (b) is illustrated in Figure 5, where the skewness obtained from (5.5) is also shown. This agrees qualitatively with the skewness behavior after canopy closure reported in [l] and [6]. I am indebted to Mark Westcott for helpful discussions. REFERENCES

I 2 3

D. J. Gates, Competition and skewness in plantations, J. Theoret. Biol. to appear. E. D. Ford, Competition and stand structure in some even aged plant monocultures, J. Ecol. 63:311-333 (1975). D. J. Gates, A. .I. O’Connor, and M. Westcott, Partitioning the union of discs in plant competition models, Proc. Roy. Sot. London Ser. A 367:59-79 (1979).

32 4 5 6

D. J. GATES L. D. Landau and E. M. Lifshitz, Theor?, of Elusticiry, Pergamon, Oxford, 1975. M. Abramowitz and I. A. Stegun, Handbook of Mathematical Functions, Dover, New York, 1968. D. J. Gates, R. McMurtrie, and C. Borough. Negative skewness of DBH in Pinw rudiata plantations. Submitted for publication.

7

D. J. Gates competition

and M. Westcott, Negative skewness models, J. Math. Biol., to appear.

and negative

correlations

8

K. Yoda, T. Kira, and K. Hozumi, Intraspecific competition among Further analysis of the competitive interaction between adjacent

in spatial

higher plants. IX. individuals, Inst.

Polyrech. Osaka Ci{l; Unio. 8(D): 161- 178 (1957). 9 10

C. F. Cooper, Pattern in ponderosa M. G. R. Cannell, G. K. Njunga, Ross-Parker),

II

Variation

pine forests, Ecology 42:439-499 (1961). E. D. Ford, and R. Smith (appendix by H. M.

in yield among

competing

individuals

within

mixed genotype

stands of tea; a selection problem, J. Appl. Ecol. 49:969-985 (1977). E. D. Ford and P. J. Diggle, Competition for light in a plant monoculture a spatial

stochastic

process,

Ann. Botany, to appear.

modelled

as