Applying a chemical modelling technique to ecological problems

Ecological Modelling, 5 (1978) 259--268 © Elsevier Scientific Publishing Company, Amsterdam -- Printed in The Netherlands

259

APPLYING A CHEMICAL M O D E L L I N G TECHNIQUE TO ECOLOGICAL PROBLEMS G.A.M. KING

Physics and Engineering Laboratory, Department of Scientific and Industrial Research, Lower Hurt (New Zealand) (Received 10 May 1977)

ABSTRACT King, G.A.M., 1978. Applying a chemical modelling technique to ecological problems. Ecol. Modelling, 5: 259--268. A very simple computer program can be adapted to handle both complicated chemical and ecological problems, by using procedures which reduce different reaction types and flows of material all to the one form. The process of reduction discloses similarities between ecological situations which, at first sight, seem very different.

INTRODUCTION

A computer modelling technique, which was developed first for describing the chemistry of ozone in the stratosphere (King, 1975), has proved useful in biological studies (King, 1977 a and b). It is both simple and very general. The generality means that, in order to adapt it for ecological work, certain devices must be used. However, this turns out to be not a drawback b u t rather an advantage because it directs attention to inter-relations and possibilities which might otherwise be overlooked. The devices used for ecological modelling are extensions of those used to reduce complex chemical networks so that the computer program need handle only one kind of chemical equation. We introduce the program, therefore, by describing its application in chemistry -- the next section. The following section carries the ideas over into ecology. CHEMICAL EXAMPLE

Appendix I gives a print-out of the program in the BASIC language used b y the Hewlett--Packard 2100A computer. This particular program will handle up to 20 chemical species and up to 20 chemical equations. It can be changed readily to accommodate more, simply by altering the DIMENSION statement -- number 200. However, if the problem is much larger, it is probably best put on a more powerful computer to cut d o w n the running time.

260 All the chemical equations are first put in the bimolecular form: A+B-+ C+D.

(1)

Each time step in the program then reduces the concentrations of A and B, [A] and [B], by k. [A] • [B], where k is the rate coefficient for Eq. (1). [C] and [D] are augmented by the same amount. The following devices reduce all chemical equations to the form of (1): If a reaction is reversible, this fact is included by writing the reverse reaction as a separate equation. Thus: C+D-~A+B.

(2)

If a reaction is an addition, yielding only one product, a fictitious second product is added -- Z, in Eq. (3) A+B-+C+Z.

(3)

First order reactions, or decompositions, also use the fictitious reactant device to generate the appearance of a second order reaction. Thus; A-+ C+ D

(4)

becomes A+Z-~C+D.

(5)

For this circumstance, the concentration of Z must remain constant in order that the rate coefficient does not change during the progress of reaction. Z is therefore given unit concentration and is not included in the program loop which updates reactant concentrations-- numbers 690 to 710; Z is chemical species number 1, and the loop starts at 2. Combinations of the forms in Eqs. (1), (3) and (5) allow higher order reactions to be expressed as chains of bimolecular reactions. The fictitious reactant also allows us to remove material from the reacting system, A + Z-~ Z + Z

(6)

or to introduce material at rates determined by some transport process: Z+Z-+Z+A.

(7)

When all equations have been put in the bimolecular form, they are converted to numbers for entering into the program as data. This is done by giving each chemical species a number in a continuous sequence, number 1 being reserved for the fictitious species (see above). For example, an over-simplified version of ozone chemistry in the stratosphere is given by this set of equations (Chapman, 1930): 02 -+ O + O, O+O2+M-~O

photo-dissociation, 3+M,

(8a) (8b)

261 O + Os -* 02 + 0 2 .

(8c)

(8b) represents a " 3 - b o d y addition" mediated by some atmospheric molecule, M. We shall assume that M is present in large concentrations so that, although we would normally resolve (8b) into two bimolecular reactions, it is enough to use just O + 02 ~ Oa.

(9)

Number the chemical species:

Z

0

02

03,

1

2

3

4.

(I0)

The equations become:

3+I-> 2+2,

(lla)

2+3-~4+1,

(llb)

2+4-~ 3+3.

(11c)

The program requests the number of chemical species -- in this case, four; the number of equations -- three; the equations themselves -- three sets of four numbers; the rate coefficients for the equations and the initial concentrations of the chemical species, as well as details of the o u t p u t desired. The rate coefficients for ( l l b ) and ( l l c ) come from laboratory data while that for ( l l a ) is the product of the intensity of incoming radiation and its absorption coefficient. (In a good modelling job, this would be an integral over wavelength rather than just a product.) The program always integrates in unit steps of time. One would not, therefore, enter the rates for this problem in conventional units with a time base of one second. Instead, the rate coefficients are scaled by a factor which alters their time base so t h a t each computation changes concentrations by a significant a m o u n t -- say, 1%. Too small steps means long machine running times, while too big steps means lowered accuracy because the program uses only a very simple integration procedure. (Earlier versions did use a more exact procedure, but the additional complication is seldom warranted for ecological work.) In many problems, one does not know which chemical species are going to be the interesting ones, and some may react so fast that the choice of unit time step causes t h e m to disappear from the model mixture. The program contains a safeguard against instability in this event -- numbers 720 to 750. It also presents various handling options, using statement 790, which allow flexibility during running. USE IN ECOLOGY Application to ecology depends on ignoring the differences between nutrients and living things, treating the latter as just complex chemical species.

262

However, it is certainly convenient to label the chemical types distinctively, so we shall label nutrients R, R1, R2, R3, etc., and populations of living things C, C1, C2, C3, etc. The simplest representation of a biological species recognizes that it takes up nutrients and reproduces: R+C-*C+C.

(12)

It also dies: C -~

(13)

Neither equation balances chemically, but the program doesn't mind. Number the chemical species: Z

R

C,

1

2

3.

(14)

Eqs. (12) and (13) become 2+3-,3+3,

(15a)

3+1-*1+1.

(15b)

Fig. 1 shows a typical run representing, say, batch culture of bacteria. It uses the rates, 0.015 for Eq. (15a) and 0.005 for (15b), the former being based on a time step of one minute which, when [R] = 1 at the start of the run, corresponds to a population doubling time of 46 min (46 × 0.015 = 0.69 ~ In 2). The "log phase" (or exponential increase) and the "phase of decline" (exponential decrease) are well portrayed. A more complicated equation set is needed to bring out the "lag phase" and the "stationary phase" characteristic of a bacterial culture.

r~ i 3F

0

500

1000 TIME

Fig. 1. Modelling of the batch culture of bacteria using the simplest representations of growth, reproduction and death.

263 Notice that the different chemical types in Eq. (15) are treated alike by the program, the living one being identified only by the "branching" process, (12) or (15a), which gives the potential for its self-reproductive increase. For many purposes, it is useful to distinguish the growing and reproductive stages in the life cycle of a biological species: R1 + C1 -* C2,

growth,

(16a)

R2 + C2 -* C1 + C1,

reproduction,

(16b)

C1 -* ,

decay,

(16c)

C2 -* ,

decay.

(16d)

Number the chemical species: Z

R1

R2

C1

C2,

1

2

3

4

5,

(17)

giving 2+4-*

5+1,

3+5-*4+4, 4+1-*

1+1,

(18)

5+1--1+1. This equation set can be interpreted in several ways. An interesting one for the cell biologist is the "synchronizing" of cells. C1 and C2 are taken to be distinct "states" of the biological species and the population can be synchronized in C1 by starving with respect to R1; it can be synchronized in C2 by cutting off R2. This raises the question of h o w nutrients are supplied. The "batch culture" situation was handled earlier (Eqs. (12)--(15)) just by specifying initial concentrations. The present " c h e m o s t a t " type of problem calls for further equations representing supply: -* R1,

(19a)

-* R2,

(19b)

inserted in the program as: 1+1-* 2+1,

(20a)

1+1-* 3+1.

(20b)

A quite different interpretation of Eqs. (16)--(20), of more interest to the animal ecologist, would say that C1 is a juvenile population while C2 is an adult population. It might even be useful to sub-divide the biological species further b y introducing states which distinguish between males and females whose interaction causes the expansion of the species. The role of sex in stabilizing population density then becomes manifest.

264 Let us elaborate on the chemical condition for a self-reproducing species. Consider a biological species which uses four nutrients, with four corresponding "states": R1 + C1 -* C2,

(21a}

R2 + C2 -~ C3 + C1,

(21b)

R3 + C3 -* C4,

(21c)

R4 + C4 -* C1.

(21d)

The conspicuous feature of Eq. (21) is that it is cyclic in the C's. In their selfreproductive aspect, living things are, chemically speaking, auto-catalytic (Hinshelwood, 1946) and an essential c o m p o n e n t of this is the catalytic, or cyclic, character of their chemistry. Earlier, in Eqs. (12) and (16b), the selfreproductive feature was built in with a branching reaction which yielded two "daughters". Eq. (21) has been arranged to show that the daughters need not appear in the same reaction -- C1 is consumed in (21a) and produced in both (21b) and (21d). But there must be a branching reaction in the cyclic set -- in this case (21b) -- wherein one of the cycling intermediates is consumed and two are produced. This kind of sub-division can be applied to species which have complicated life-cycles. When a single biological species occupies two (or more} regions with limited exchange between them, it is treated as two (or more) species. Exchange is then treated as transformation between species. Thus, if a biological species is labelled C8 in one region, and C9 in another, these equations would appear:

C8 ~ C9,

(22a)

C9 -+ C8.

(22b)

The ecologist is much concerned with the circulation of materials in the environment. If he were able to write down the set of equations for all materials in a closed environment (suffering neither loss nor gain), this set would be cyclic in all materials. Certain sub-sets of this enormous set would be characterized by three features: (1) They are cyclic with respect to a set of chemical species. (2) Moreover, these cyclic sub-sets contain branching reactions so that they form self-reproducing systems. (3) And they are also "packaged" in physical units -- the biological individuals -- see King (1977a and b) for a suggested explanation. Eqs. (16a) and (16b) stand as a very simple example. However, the long term average population of such a species is steady in spite of the auto-catalytic growth because the decay reactions, e.g. (16c) and (16d), feed the materials back into the over-all set for the environment. Now, the over-all set also contains sub-sets showing the first of the above characteristics and probably the second but not the third. These are the networks of associations between biological species in the environment. It may be useful to appreciate the

265 similarities between the chemistry of the environment and the chemistry of the cell. The discussion has become rather abstract, so let us end with an example familiar to all theoretical ecologists -- the predator--prey scheme (Goel et al., 1971}. We will give it the extra touch of different rates of predation on juvenile and adult prey. R

+ C1 -* C2

/

R + C2-* C1 + C1

prey

C1 + C3 -~ C3 + C3,

predation on juveniles,

C2 + C3 ~ C3 + C3,

predation on adults,

(23)

C3 -* . Following tradition, we keep R constant and it is therefore convenient to use the fictitious reactant for it. Notice that the predators suffer death from old age, while adult prey die on reproducing. R

C1

C2

C3

1

2

3

4

(24) Eq. (23) becomes 1 + 2-* 3 + 1,

rate = 0.01,

1+3-* 2+2,

0.01,

2 + 4 -* 4 + 4,

0.005,

3 + 4 -~ 4 + 4,

0.0025,

4 + 1 -~ 1 + 1,

0.005.

(25)

Eq. (25) lists the rates on the right. The first two rates correspond to a time step of 1 day if the prey has the potential of doubling its population in about 2 months. Again following tradition, we have ignored many necessary scaling factors such as the inefficient use of food. We could include this particular factor by supposing that each unit of prey represents, perhaps, 30 individuals. Fig. 2 shows o u t p u t as total populations of prey (. . . . ) and of predators ( ......... ) together with the juvenile/adult ratio in the prey. This last measure started at a value of 2/3 but quickly shifted into an oscillation, driven by the predators, around the equilibrium value which is 1.281 for the particular rate coefficients used. Naturally, the shift in juvenile/adult ratio is accompanied by a small change in the amplitude of predator and prey oscillations. Numerical modelling gives some " f e e l " for the behaviour of non-linear systems. It is valuable to complement it with analysis, even though this may not be taken very far. Thus, the set of chemical equations, (23), generates the

266

3

. |

FiI~EY~

|

i

~ ~". 0

.

1

/

.

.

--'

I

','

~"~PREDATOR

/

\ " ',

X

",

!

'"

J '\

'"

',

"\\

~

/

: \

"",,,./,/

/"

l ....................... 500 1000 1500

2000

'

", |

//~

.OU~TS

TIME

"

'\ 2500

1

1

Fig. 2. T h e i n t e r a c t i o n b e t w e e n p r e y ( . . . . ) and p r e d a t o r s ( ......... ) with d i f f e r e n t rates o f p r e d a t i o n o n juvenile and adult prey. T h e j u v e n i l e / a d u l t ratio is s h o w n ( ).

following differential equations:

CI=--kl.R.CI+2.k2-R-C2--k3.C1.C3, C2 = k l • R • C1 - - k 2 • R • C 2 - - k 4

• C2 • C3,

(26)

C3 = k3 • C1 • C3 + k4 • C2 • C 3 - - k s • C3, w h e r e t h e k ' s are t h e r a t e c o n s t a n t s f o r t h e five c h e m i c a l e q u a t i o n s a n d conc e n t r a t i o n b r a c k e t s have b e e n o m i t t e d . T h e e q u i l i b r i u m values are C I * = 0.719,

C2" = 0.562,

C3" = 1.123.

(27)

T h e linearized version o f (26) yields eigenvalues, --2.83×10 -2

and

--6.5×10 -5-+4.52×10

-ai,

suggesting a p e r i o d i c s o l u t i o n d e c a y i n g slowly t o w a r d s t h e e q u i l i b r i u m values. ACKNOWLEDGEMENT

M y t h a n k s t o t h e j o u r n a l referee, w h o suggested t h e discussion o n Eq. (26), a n d t o L. F r a d k i n w h o c h e c k e d t h e n o n - l i n e a r characteristics o f this e q u a t i o n set. REFERENCES Chapman, S., 1930. A theory of upper-atmospheric ozone. Mere. R. Meteorol. Soc., 111: 103--125. Goel, N.S., Maitra, S.C. and Montroll, E.W., 1971. O n the Volterra and other non-linear models of interacting populations. Rev. Mod. Phys., 43: 231--276. Hinshelwood, C.M., 1946. Chemical Kinetics of the Bacterial Cell. Oxford University Press, London, 274 pp. King, G.A.M., 1975. Halomethanes and the Ozone Layer -- A Survey. P.E.L. Report 480, Dep. Sci. Ind. Res., Lower Hutt, 16 pp. King, G.A.M., 1977 a. Symbiosis and the origin of life.Origins Life, 8: 39--53. King, G.A.M., 1977 b. Symbiosis and the evolution of prokaryotes, BioSystems, 9: 35--42.

267 APPENDIX

Program listing CHEMEQ 100 110 120 130 140 150 160 170 180 190 200 210 220 230 240 250 260 270 280 290 300 310 320 330 340 350 360 370 380 390 400 410 420 430 440 450 460 470 480 490 500 510 520 530 540 550 560

REM Chemical equations -- CHEMEQ R E M P r e p a r e list of c h e m i c a l species R E M First species is a b l a n k " b o x " R E M Prepare list of e q u a t i o n s REM Convert equations to number form R E M List rate c o e f f i c i e n t s R E M List initial c o n c e n t r a t i o n s for species R E M First species has u n i t c o n c e n t r a t i o n R E M Species a n d e q u a t i o n s d i m e n s i o n e d t o 20 REM DIM E [ 2 0 , 4 ] , R [ 2 0 ] , V [ 2 0 ] , W [ 2 0 ] P r i n t " e n t e r n u m b e r o f c h e m i c a l species I n p u t N1 Print "enter number of equations I n p u t I1 P r i n t " e n t e r e q u a t i o n s - - 4 integers ForI=ltoI1 Input E[I,1], E[I,2], Eli,31, E[I,4] Next I Print P r i n t " e n t e r rates ForI=ltoI1 Input R[I] Next I P r i n t " e n t e r initial c o n c e n t r a t i o n s ForN=ltoN1 Input V[N] Next N T9 = 0 P r i n t " s p e c i f y w h i c h species p r i n t e d - - c h o o s e 4 Input P[1], P[2],P[3],P[4] P r i n t " s p e c i f y p r i n t interval a n d limit I n p u t T1, T2 Print Print P[1], P[2], P[3], P[4], "T" REM T = T9 ForJ=lto4 K = P[J] X[J] = V[K] Next J Print X[1], X[2], X[3], X[4], T If T > = T 2 t h e n 7 8 0 REM T---INT(T+I'I) ForN=ltoN1 W [ N ] -- 0

- - e n t r y PT 1 " -- e n t r y PT 2 " - - e n t r y PT 3 "

- - e n t r y PT 4 "

-- e n t r y PT 5 "

- - e n t r y PT 6 " - - e n t r y PT 7 "

268 570 580 590 600 610 620 630 640 650 660 670 680 690 700 710 720 730 740 750 760 770 780 790 800 810 820 830

Next N ForI=ltoI1 A-- E[I,I] B = E[I,2] C = E[I,3] D = Eli,4] Z = R[I]*V[A]*V[B] W[A] = W[A]-- Z W[B] = W[B] - - Z W[C] = W[C] + Z W[D] = W [ D ] + Z Next I ForN=2toN1 V [ N ] = V [ N ] + W[N] Next N ForN=2toN1 If V [ N ] > 0 t h e n 750 V [ N ] -- 0 Next N If I N T ( T / T 1 ) * T 1 = T t h e n 470 GOTO 540 T9 = T Print " c h o o s e e n t r y PT for n e x t run - - use 8 t o e n d s e t " Input H GOTO H o f 210, 230, 250, 300, 340, 3 9 0 , 4 1 0 , 8 3 0 GOTO 790 End