SCOPE 34 - Practitioner's Handbook on the Modelling of Dynamic Change in Ecosystems

5.1 THE ORIGINS AND DEVELOPMENT OF DYNAMIC MODELS

Today's dynamic models and sub-models of ecosystems have their origins in the classical analytical models of physics, chemistry and biology. These analytical models depend on mathematical procedures for finding exact solutions to differential and other equations. As the behaviour of many kinds of mathematical equations is well known, such analytical solutions provide considerable scope for the investigation of the response of mathematical expressions to changes in basic parameters or variables. They require little more than the use of pencil and paper and a knowledge of mathematics, although that knowledge may need to be considerable for anything but the simplest of problems. Analytical solutions to such problems have the enormous intellectual appeal which characterizes the aesthetics of mathematics, as well as the predictive power of scientific hypotheses. Indeed, Maynard-Smith (1974) makes a clear distinction between what he calls models, i.e. mathematical descriptions of ecological systems which are essentially theoretical, and simulations, which have a practical purpose and whose utility lies mainly in the description of particular cases.

Undoubtedly, analytical models have been, and remain, useful in some fields of ecology, especially in population biology and population dynamics, as well as in the development of the theoretical models which are the basis for many developments in ecological theory. A simple example of the use of analytical models is provided by the Volterra equations describing the interactions between a prey species and its predator. These equations were described by Volterra (1926), and Maynard-Smith (1974) gives the following analysis.

If we assume that the density of the prey species is represented by the variable x, and the density of the predator species by the variable y, then one model for the interactions between predator and prey is given by the equations:

This model makes a number of important assumptions, which may be summarized as follows.

1. The density of a species (the number of individuals per unit area) can be adequately represented by a single variable.

2. Changes in density can be described adequately by deterministic equations.

3. The effects of interactions within and between species are instantaneous.

4. In the absence of predation, the prey species would follow the logistic function, with an intrinsic rate of increase and a carrying capacity of a/b.

5. The rate at which prey are eaten as proportional to the product of the densities of predator and prey.

How does the system described by these equations behave? At any given time, t, the state of the system is described by the values of x and y, so that there is a corresponding point describing the system in the x and y phase plane. To each point we can also attach an arrow to indicate the direction in which a system at that point will move, and, if these arrows are joined up, they form trajectories describing how the system behaves. Formal analysis of the system in this way, derived from knowledge of the properties of differential equations, demonstrates that the system converges to a steady state with prey and predator both present only when the carrying capacity a/b is high enough to support the predator. Modification of the original equations enables various other influences on the relationship between prey and predator to be investigated, for example the effects of cover for the prey, the effect of a predator with a constant food intake, etc.

See Maynard-Smith, 1974, for a more complete account.

Despite their aesthetic attractions, analytical models are useful only under certain rather restricted conditions. First, the equations describing the ecological processes must be linear if more than a small number of these equations are to be solved simultaneously. Second, the solution of even a small number of non-linear equations may become difficult or even impossible under particular circumstances. Table 1 gives a simple decision table for the solution of mathematical equations by analytical methods, ranging from the trivial to the impossible. Finally, there are considerable difficulties in the use of analytical models for the description of ecological processes which are discontinuous, as are many of the processes of plant and animal reproduction. Analytical models have not, therefore, been found to be of particular value in the study of ecosystem dynamics and in the modelling of whole ecosystems.

As a result, and because of the greater interest in modelling complex ecological processes, increased attention has been given to simulation models which use simpler forms of mathematics in conjunction with the speed and power of modern electronic computers. By placing any equation set on a computer and exploring interactively the consequences of changing the parameters of the equations, or the values taken by the variables, it is possible to obtain increased understanding and improved prediction of the dynamics of ecological systems and their components. Such simulation models do not usually give an exact solution to an equation over time as does an analytical model, and they are frequently subject to errors which are introduced by the inexact nature of the techniques of obtaining a solution. However, it is usually possible to solve many equations nearly simultaneously, and it is also possible to include a wide variety of non-linear relationships with the equations. The remainder of this chapter describes the use of simulation models of dynamic systems.

Table 1 Mathematical Problems and their ease of solution by analytical methods

5.2 DEFINITIONS AND NOTATION

Some careful definition of terms and of mathematical notation is necessary for the understanding of dynamic models. Indeed, failure to conform to reasonable standards and definitions has done much to make dynamic modelling difficult to understand by the non-mathematician, and hence unpopular. In this section, we will therefore define some of the more important terms and their mathematical notation.

The variables which describe the state or condition of an ecological system are called state variables, and are usually shown as a capital letter qualified by a subscript e.g. Y_i, X_j, etc. A model which has more than one state variable is called multi-dimensional.

Rates of movement or energy between factors described by state variables of a system are called flow rates, and the symbol M(i, j) denotes the flow rate of material from, say, state variable Y_i to Y_j. The symbol M(i, i) represents that particular change in the state variable Y_i which is not a result of flows from other state variables. As an example, consider the simple predator-prey model in Figure 5. Y₁ in this figure is the population of some organism, and Y₂ is the population of another organism which acts as a predator on the first. Then M(1,2) is the rate at which the prey Y₁ are being consumed by predators Y₂. Similarly, M(1,1) is the net growth (births minus deaths) of the prey population, and M(2,2) is the net growth of the predator population in the absence of predation. There is no M(2,1) as no material moves from the predator to the prey.

The flow rates M(i, j) will usually be affected by factors other than those described by the state variables. Model variables whose values are determined by management decisions are called decision variables. Such variables are controlled solely by human decisions. Flow rates may also be influenced by factors outside the system, but which nevertheless affect the behaviour of the system. These variables are called exogenous variables or driving variables. In dynamic ecological models, driving variables describing environmental conditions, such as temperature, precipitation and solar radiation, often have an important influence on flow rates like primary productivity, decomposition or nutrient uptake. It is not unusual for a flow rate to depend on several state variables, and on several exogeneous or decision variables. Rates that depend on several factors are thus represented mathematically by adding together or multiplying together effects of factors that are assumed to act independently of one another.

The total rate of change of a state variable, Y_i, is the sum of effects due to flow into Y_i from other variables less the flows from Y_i to the other variables.

where F_i is the total rate of change of a state variable Y_i; a_ji is a conversion ratio of material Y_i to material Y_j;

Figure 5. Simple example of state variable and flows in predator-prey model

and

M(j, k) are the flow rates from Y_j to Y_k

Rates of change are of fundamental importance in dynamic models because they are the basis of differential and difference equations used to predict the behaviour of ecological systems.

Some, or all, of the variables included in a model may be random variables, with each of which is associated a probability distribution that describes quantitatively the likelihood that the variable has a specified value, or falls within a specified range. Models which contain random variables are called stochastic models, and the output of such models represent probabilities rather than the fixed values which will be obtained from a deterministic model.

Finally, systems in which the variables change with time are called dynamic. Dynamic systems, and models of dynamic systems, may eventually reach a state at which the values of the variables do not change. A system, or model, in which all the state variables remain constant with time is said to be in a steady state or equilibrium. Long-term studies of systems may therefore concentrate on their steady-state behaviour rather than on transient fluctuations, and the calculation of the values of the state variables at which the system reaches a steady state is of particular importance.

5.3 DIFFERENTIAL AND DIFFERENCE EQUATIONS

Models of dynamic change incorporate rates of change, F_i with values of state variables, Y_i, so as to predict the altered states of those variables with time. If time is divided into discrete units such as days, generations or years, the value of the state variable, Y_i, at time t + Dt is the value of the state variable Y_i, at time t, plus the rate of change, P_i, multiplied by Dt, i.e. the number of time units that have elapsed. In mathematical terms, therefore:      

Y_i (t + Dt)=Y_i(t) + F_i(t)   Dt

Such an equation is called a difference equation, and can be used to predict changes at discrete points in time.

In contrast, equations that describe the dynamics of a system continuously through time are called differential equations, and are based on the derivative of Y_i(t). The derivative is a mathematical expression of the instantaneous rate of change of Y_i(t), and is denoted by the function d Y_i(t)/dt, or alternatively by Y.

Under defined mathematical conditions, this limit can be shown to exist and the function Y_i(t) is said to be differentiable at the point t.

The purpose of predictive equations in dynamic modelling is to determine the values of Y_i at each time t. An expression that gives these values of Y_i(t) is called the solution of the predictive equations. One of the principal advantages of differential equations is that a closed form or analytical solution may be possible from the mathematical expression of the relationships. For example, the exponential equation:

where Y(t₀) = c, has the closed form solution:

In general, however, closed form solutions are only attainable for models which are very simple, or which have some special form. Because of the complexity of most ecological systems, models of dynamic systems must be solved numerically, and the modern computer has greatly increased our ability to find numerical solutions to such problems.

The numerical solution of difference equations is found by recursive calculation of the values of the state variables, i.e. the values of the state variables are used to calculate the values at time (t + 1), which in turn are used to calculate the values at time (t + 2), etc. The process starts at t = t₀, for which initial values Y_i(t₀) are given. Comparison of recursive and analytical solutions shows that the state variables are consistently underestimated by the use of difference equations. This underestimation is caused by the false assumption that the growth rate remains constant for the whole of the time period of a single step, when, in fact, the rate changes continuously. The discrepancy between recursive and analytical solutions decreases as the time interval is made smaller. One of the principal problems of obtaining solutions by the use of difference equations, therefore, is the choice of an appropriate time interval. If it is too short, the length of time needed for the solution increases dramatically; if it is too long, the underestimate of the state variable distorts the solution.

It may be argued that an acceptable time interval is one which gives an answer within a specified range from the answer obtained by an analytical solution. However, this is not a particularly helpful criterion. It is not necessary to find an approximate answer if an analytical solution exists, and, if no analytical solution can be found, we have no standard of comparison for the approximation. An alternative approach, therefore, is to ensure than an acceptable time interval has been reached when halving its value does not change the relevant results of the simulation by more than a pre-set amount. This pre-set amount should not, in general, be smaller than is warranted by the accuracy of the parameters, initial values and tabulated functions. In this way, a reasonable compromise may be achieved between the apparent accuracy of the answer and the computational effort, which increases at least linearly with the decreasing magnitude of the time interval.

More sophisticated methods of integration of differential equations are, however, frequently desirable, especially where it is necessary to vary the size of the time interval during the simulation, as frequently occurs in the more complex models. A wide range of methods exist for calculating the solutions to differential equations. Many of these methods are based on a Taylor series expansion of the function, Y(t). From Taylor's theorem, any finite real-valued function, Y(t), for which a finite derivative dQ/dt exists over the interval t = s + t, can be represented by the expansion:

where m is any integer greater than zero and e(Dt^{m+ 1}) is an error term less than a constant time (Dt)^{m+ 1} Provided that t < 1, this error term becomes very small as m increases. If m = 1, then the expansion is reduced to:

One of the simplest techniques for solving differential equations numerically is based on the reduced Taylor expansion, and is due to Euler. The value of the state variable Q at time t + Dt is calculated from the equation:

where F(Y, t) is the rate of change. The error of this approximation, e(Dt²), will be small if t is less than 1. In essence, the solution becomes a difference equation that is used recursively to calculate the values of Y(t) for t₀ £ t £ T if Y(t) is differentiable within the limits t₀ < t < T.

Euler's method requires at least ( T- t₀ )/ Dt recursive calculations in the interval t₀ <t < T, and, if Dt is small, the number of calculations may be large. Even though the error at each single step is small, the accumulation of errors over many steps may lead to a considerable loss of accuracy. A method which has a smaller error than e(Dt²) at each step is known as the Runge-Kutta method. This uses a difference approximation of the form:

where C_i, a_i and b_i are chosen so that the expansion is equivalent to that of the Taylor series. When N = 4, the commonly used Runge-Kutta approximation is:

The approximation is a difference equation which can be used recursively, but which has an error at each step of the order (Dt)⁵. A computer subroutine for solving differential equations by the Runge-Kutta method is therefore an essential part of the ecologist's toolkit for the modelling dynamic systems. An alternative and rather simpler method of mid-point integration also gives quite acceptable accuracy, and a computer subroutine for this method is also generally useful.

BASIC algorithms for the Euler, Runge-Kutta and mid-point methods are given in the Appendix to this Handbook.

The individual steps in constructing and solving a dynamic model may be summarized as follows.

1. Define the state variables, exogenous variables and control variables of the system to be modelled.

2. Estimate, or determine experimentally, the flow rates M(j, i) between the state variables.

3. Calculate the rate of change of each variable as the sum of flows into the variable minus the flows from the variable, i.e.

4. Describe the dynamics of the system by substituting each rate of change, Fi, into a series of difference or differential equations.

5. Find a numerical solution to the predictive equations of the dynamic model. If only the equilibrium solution is of interest, it can be obtained by solving the equation F( Y) = 0.

6. Test the sensitivity of the model to small changes in the initial values of the state variables, and of the exogeneous and control variables. Test also the sensitivity of the model to small changes in the estimated flow rates.

7. When the model has been shown to be reasonably robust to small changes in its basic parameters, the effects of various management policies may be tested by trial and error simulation, or by optimization methods.

5.4 MODEL TYPES AND APPLICATIONS

5.4.1 Population models

One obvious application of dynamic models is the prediction of populations of organisms in response to various environmental effects. Such applications have a long history, and many of these are referred to in the case studies in Section 5.6. Here, therefore, we will concentrate on some simple, but important, examples.

The simple model for the growth of a single species is the exponential growth equation:

in which the rate of change F is equivalent to the number, mass or density of the population multiplied by a constant growth rate, r.

The exponential growth model has a closed form solution:

Y(t + 1) = Y(t) + rY(t)

Table 2 gives the computed values of the growth of a hypothetical population with a mass of 1 g at t₀, and a growth rate of 0.1 g/hour. The first column of this table gives the time from zero to 10. The second column gives the exact solution to the differential equation, while the second, third and fourth, respectively, give the Euler, mid-point and Runge-Kutta solutions. The Euler method quickly leads to quite serious underestimates of the mass of the population, while both the mid-points and Runge-Kutta methods give acceptable estimates.

A more realistic model of population growth of a single species is given by the logistic growth equation:

Table 2 Comparison of solutions of exponential growth equation y = e^0.1t

where the rate of change is equal to the number, mass or density of the population multiplied by a constant growth rate, r, and by a factor (K- Y)/K which approaches unity as Y approaches the value of K. Again, this differential equation has an exact or closed form solution:

where

The computed values of the growth of a hypothetical population with a mass of 0.45 g at t_o, a growth rate of 0.05 g/hour, and an upper value for K of 6 g are given in Table 3. Again, the Euler solution is less satisfactory as an approximation to the exact solution than either the mid-point or Runge-Kutta method.

For both the exponential and logistic growth models, exact or closed form solutions exist, but, for most population models, such solutions are either mathematically difficult or impossible. As a simple example of such a model,

Table 3 Comparison of solutions of logistic growth equation y = K/(l + ce^-0.05t)

consider the growth of a population of bacteria which is not limited by either space or food, but for which the relative growth rate is dependent on a fluctuating daily temperature (de Wit and Goudriaan, 1974). If the temperature is assumed to be prediced by the equation:

                             T= 20 + 10(sin(27p X/24)) °c

where X is the number of hours, it will fluctuate between 10° C and 30 ° C every 24 hours. If, further, the relative growth rate is assumed to be related to temperature by the quadratic equation:

                      G = 0.00554 + 0.01019T - 0.0001009T ²

then the exponential growth model has a variable relative growth rate:

Table 4 Comparison of solutions for temperature-dependent exponential growth of a bacterial colony

Table 4 summarises the predicted growth of a population of bacteria with a starting mass of 1 g under these assumptions, based on Euler, mid-point, and Runge-Kutta solutions. Here, the simple Euler solution gives a substantial underestimate of the growth of the bacterial colony after about 12 hours. The mid-point and Runge-Kutta solutions are roughly comparable, within the limits of the basic assumptions.

All of the above models can be solved relatively easily by the use of the simple computer algorithms DIFFEU, DIFFEQ and DIFFER which are given in the Appendix. The only changes that need to be made in the programs are the substitution of the required function in the DEF statement at the beginning of the program. Similarly, one or more of the algorithms can be incorporated into a larger program to simulate the growth of a population.

More complex population models may require competition between species, or even between plants of the same species, to be taken into account. For example, the number of plants of each species is a mixed stand of an agricultural or forestry crop is determined at the time of planting. The competition between each crop species can then be expressed in terms of 'relative space', a dimensionless variable which characterises the effects of crowding on available root and foliage space, nutrients, sunlight and associated factors. The actual production of dry matter can then be obtained from the product of 'relative space' and the maximum possible yield for dense monocultures. All three quantities are therefore functions of time.

To simulate, for example, the growth and competition of barley and oats planted as a mixed stand, the differential equations for the relative space r_b and r_o for barley and oats respectively are as follows:

dr_b = G(r_b - r_b² - r_br_o)
dr_o = G_o(r_o - r_o² - r_br_o)

where G_b and G_o are relative growth rates of barley and oats as empirical functions of time in the absence of competition summarized in Table 5.

Table 5 Empirical growth rates for barley and oats

Table 6 Computed proportions of growing space for oats and barley growing in competition

The calculated proportions of the growing space occupied by the two species when sown at an initial spacing of 0.02 cm are given in Table 6. Because the barley maintains its relatively fast growth rate, a disproportionate share of the available space is occupied by barley at an early stage when the two species are grown together. By the time the oats begin to grow, there is insufficient space to accommodate them.

5.4.2 Predator-prey models

While population models, with or without competition or the constraints on resources to support growth, represent one of the main areas of application for dynamic differential or difference equation models, one of the most popular areas for such models has always been that of the predator-prey relationship. In this simple review of predator-prey relationships, we will first consider predator-prey systems without an age structure and then extend the models to include systems which have an age structure.

5.4.2.1 Predator-prey species without an age structure

Dynamic models of predator-prey interaction are usually based on several simplifying assumptions, including the following.

1. The density of species, usually expressed as the number of individuals per unit area, can be adequately represented by a single variable. Differences of age, sex, genotype, etc., are all assumed to be capable of being ignored for the purpose of this model.

2. Changes in the density of the organisms are assumed to be adequately described by deterministic equations, i.e. there are no obvious random or probabilistic events.

3. The effects of the interactions between the organisms are assumed to be instantaneous, i.e. there is no lag between the consumption of one organism by another and the response of the consuming organism to the increase of available energy.

If these relationships are true, or approximately true, then the predator-prey relationships may be capable of being represented by differential or difference equations.

5.4.2.2 Volterra equations

Volterra (1926) described the interactions between a prey species, x, and its predator, y, by the following equations:

We have already considered the analytical solution to these equations in Section 5.1, and suggested that the system converges to a steady state with prey and predator both present only when the carrying capacity at b is high enough to support the predator (Maynard-Smith, 1974). However, it will be interesting to explore this simple representation of a predator-prey relationship by the same methods described for population models. If we can show that these methods give the same results as the analytical solutions, where these exist, we may have more confidence in the use of the numerical methods in situations where the analytical solutions do not exist.

Analytical solution of these equations showed that both predator and prey numbers oscillate with decreasing amplitude, the predator oscillations lagging in phase behind the prey. If y = 0, the prey species is limited only by the predator, and the prey increases exponentially in the absence of the predator . The oscillations are then of constant amplitude, depending on the initial conditions - a system started close to its steady state will have small-amplitude oscillations, and one started far from its steady state will have large-amplitude oscillations. Such a system is called 'conservative' because there is a quality which is conserved during the motion, as energy is conserved in simple harmonic motion.

The term bx², expressing the inhibiting effect of a species on its own growth, is referred to as a 'damping' term. In ecology, the main factor reducing oscillations is the presence of self-inhibiting effects. The Volterra equations, without damping, are written as:

From these equations can then be derived the equation: 

e ln(x) -c' x + a ln(y) = constant

which represents a family of closed curves in which each member of the family corresponds to a different value of the constant. Choice of a starting point, i.e. of initial values of x and y, determines the value of the constant. Any population will continue indefinitely to follow the cycle on which it starts. When the cycles are of very small amplitude, they may be approximated by ellipses. For simple models of the kind illustrated by the Volterra equations, it is possible to derive the equation of such an ellipse, and some interesting properties can be inferred from it.

Four properties are true of small cycles (in the neighbourhood of the equilibrium point).

1. The sizes of both the prey and the predator populations vary sinusoidally with period T = 2p / ae. The period is the same for both species and depends only on the parameters a and e.

2. The two populations are always a quarter-cycle out of phase. Thus, the prey population begins to decrease from its maximum size at the instant when the predator population is attaining its fastest growth. One quarter-cycle later, when the prey population is declining most rapidly, the predator population attains its maximum size and starts to decline.

3. The amplitudes of the oscillations, unlike the period of the oscillations, depend on the initial sizes of the populations.

4. The mean size of the prey population over a period T is e/c', and its identical to the equilibrium size of the prey population. Similarly, the mean size of the predator population is a/c.

Consider a difference equation model of a predator and prey equation whose rate of change is defined by the equations:

F₁ = r₁ Y₁ -a₁ Y₁ Y₂
F₂=r₂Y₂+a₂Y₁Y₂

Then:

Y₁(t + 1) = Y₁(t) + r₁ Y₁(t) - a₁ Y₁(t) Y₂(t)
Y₂(t + 1) = Y₂(t) + r₂QY₂(t) - a₂ Y₁ (t) Y₂(t) 

If r₁=0.5, a₁=a₂=0.001, and r₂= -1,and Y₂= -1, and Y₁(0)=1000 and Y₂(0) = 400, recursive calculation of these difference equations gives the results shown in Table 7.

Table 7

The properties of dynamic models of this kind and readily explored by computers, especially when provided with graphics terminals. Figure 6 gives a simple difference equation solution, using the Euler approximation, derived from the BASIC algorithm in the Appendix. In this program, the solution is given for the equations:

x = 4x - 2xy
y = 3y + xy

with the starting values of x = 3,  y = 0.5; x = 3,  y = 1; x = 3,  y = 1.5; and x = 3,  y = 2.

The results in Figure 6 show three roughly ellipsoid trajectories and a single point, this point corresponding to the starting value of x = 3, y = 2 and representing equilibrium. If the two populations start at this point, no change ever occurs. For any other starting values, the populations of both prey and predator circle about the equilibrium point in a counter-clockwise direction.

Figure 6. Trajectories of predator-prey model for four different starting points

A modification of this algorithm shows the effect of the damping term -bx², by changing statement 120 of the program to X = X + (4 * X - 0.25 * X * X - 2 * X * Y) * 0.01. If the starting values are set to x = 3, y = 0.5, the sizes of the two populations quickly lead to the equilibrium point of the solution (Figure 7).

Figure 7. Trajectories of predator-prey model when damping term is included

Use of computer models and difference equations in this way, even with relatively inefficient methods such as the Euler approximation, can therefore be helpful in showing the behaviour of systems of equations. They are particularly useful when, in contrast to the simple model of the Volterra equations, analytical solutions cannot be readily derived. A general-purpose algorithm for plotting the trajectory of predator-prey populations is given as the program PREDAT in the Appendix.

Leslie (1948) gave an alternative formulation of the predator-prey equations: 

The equation for the prey is the same as Volterra's equation with damping. The equation for the predator resembles the logistic equation, but the second term has been modified to allow for the density of the prey. If x/y is large (many prey per predator), the predator increases exponentially. If x/y = f/e, the predator is at equilibrium. If x/y is less than f/e, the predator decreases in numbers. The equations lead to rapidly damped oscillations.

Volterra's equations are usually to be preferred for two reasons.

1. In the Volterra equations, whether the predator increases or decreases depends only on the density of the prey, while, for the Leslie equations, it depends on the number of prey per predator .

2. The Volterra equations relate the rate of increase of the predators to the rate (c' xy) at which the prey is being eaten. In the Leslie equations, there is no relationship between the rate at which a predator eats and the rate at which it reproduces.

5.4.2.3 The effect of protection for the prey

Suppose that some number, x_r, of the prey can be protected so as to make them inaccessible to the predator. Volterra's equations (without damping) then becomes:

If the number of prey which are protected is a constant fraction of the total, so that x_r = kx, the equations can be written as:

Such a change clearly does alter the nature of the equilibrium. If, however, the number of the prey which are protected is constant, so that x_r = k, then it can be readily shown that the effects of the protection are stabilizing, changing a 'conservative' system into one with convergent oscillation.

5.4.2.4 Predator with constant food intake

For some predators, a more reasonable assumption that an individual predator takes prey at a constant rate is determined by the food requirements of the predator and not by the density of the prey. The equations then become:

However, the equation for the predator is unrealistic, as it provides for an exponential increase if c' > e, or an exponential decrease if c' < e. If we assume that a predator has a constant food intake, independent of prey density, the predator cannot be limited by its prey. The predator must, therefore, be self-limited in some way, with a damping term in the equation for y, i.e.:

Such a system has two stationary points when x and y are non-zero. One of these points is a stable non-oscillatory equilibrium and the other is an unstable equilibrium. Which of the two stationary points the system will reach depends on the initial conditions.

The model is unrealistic for any initial conditions which are not close to the stable equilibrium point, because the assumption that every predator takes prey at a constant rate, regardless of prey density, must ultimately break down when the prey become rare.

5.4.2.5 General case - Rozenzweig-MacArthur model

A more general case is given by the equations :

where f(x) is the rate of change of x in the absence of predators;
ø(x, y) is the rate of predation;
k is the conversion efficiency of prey into predator; and
e is the mortality rate of the predator.

Analysis of these equations for continuously reproducing systems gives the following conclusions.

1. The interaction of predators and prey, when the predator is limited only by the supply of prey, leads to regular fluctuations in the numbers of both predators and prey.

2. If the numbers of prey are limited by resources, and not by the predator, the oscillations will tend to be damped out.

3. If the predator is limited by some other factor than the prey, the oscillations will also tend to be damped out.

4. If some form of protection makes a constant number of prey unavailable to the predators, the oscillations will again tend to be damped out.

5. Fluctuations are likely to increase in amplitude, perhaps leading to the extinction of one or both species, if the predator is able to maintain itself when the prey density falls below the level necessary for the survival of the predator.

5.4.2.6 Stochastic simulation of predator-prey populations

A stochastic version of the purely deterministic Volterra equations:

can also be simulated, assuming that these equations represent the next outcome of a large number of random events.

Suppose that a is the average birth rate for the prey, and that all prey deaths are due to attack by the predators. Similarly, we suppose that c is the average death rate of the predator, and that the birth of a predator always coincides with the death of a prey animal, e.g. if the predators are parasitic wasps that lay their eggs in the larvae of host insects so that development of each parasite necessarily involves destruction of host larvae or pupae. With this limitation, the number of possible events is either three of four, the case with three possible events arising when every death of a prey animal is accompanied by the birth of a predator, so that b = d = p (say). The events, and their associated probabilities, are summarized below, where c is a constant.

If some prey deaths do not coincide with predator births, and only a proportion of predator attacks are successful, then d = q b, where q is the fixed proportion. The event x®  x - 1 has two possible outcomes: either y®  y + 1 with probability proportional to q bxy = dxy, or y® y with probability proportional to (1- q )bxy = (b -d)xy. The case q = 1 is the special case above.

The risk of extinction of one or the other of the two species is clearly greater when the trajectory of the combined populations passes close to either of the axes. When a population reaches a point in its trajectory such that one of the species has only a few members, the probability that species will fail to recover is high. If the predator dies out, the prey population will grow without restriction until some density-dependent or regulating function comes into operation. If the prey dies, the predator must also die, unless it can transfer its dependence to some other species.

5.4.2.7 Predator-prey systems with age structure

So far, we have considered predator-prey systems with no defined age structure, thus keeping the mathematical representation of the systems as simple as possible. For most ecological systems, however, the effects of breeding seasons and age structure cannot be ignored. These effects often cause oscillations of large amplitude to occur in the behaviour of the systems through the operation of the phenomenon which is known to engineers, as 'feedback'. As a general rule, if the duration of the delay in a feedback loop is longer than the 'natural period' of a system, oscillations of large amplitude will result. The 'natural period' for the system is readily defined as follows.

If the growth of a population in the absence of a feedback loop governing its regulation is represented by the differential equation:

then the natural period is 1/r.

5.4.2.8 Delayed regulation of systems

Delay in the regulation of ecosystems may arise from one or more of three causes.

1. Development time. Changes in the environment of a species may produce an immediate change in (say) the rate at which eggs are produced by mature females. However, the corresponding change in the number of adults will be delayed by the length of time that it takes for an egg to develop into an adult, i.e.:

   dx/dt = f(Xt -to)

   where Xt- T is the adult breeding population at some time t₀ in the past.

2. Discrete breeding seasons. If a species breeds only at a particular time, e.g. only in the spring of each year, some delay is necessarily introduced into the processes which regulate the ecosystem, even where the individual organisms usually survive in breeding again in successive years. Where the individuals of a species live for many years, and produce relatively few young each year, the delay of one year is likely to be short compared to the natural period for the species, so that any oscillations in numbers are convergent. Where, however, the adults breeding in one season survive only rarely to breed in the next year, the equation becomes:

   xn+1 = øxn

   where x_n is the population size in year n. Even where the differential equation has a stable equilibrium, the corresponding difference equation will show divergent oscillations.

3. Delayed response by limiting factors. The numbers of a species may oscillate if there is a delay in the response to the factors limiting the population size of the species, even if the species itself responds quickly to the environmental conditions.

Maynard-Smith (1974) gives some simple examples of predator-prey systems incorporating such delays.

5.4.2.9 Matrix representation of population growth with age-dependent birth and death rates

Where a population can be divided naturally into a series of discrete age classes, a matrix representation of the population has been suggested by Leslie (1945, 1948). At any given time t, the population can be represented by a column vector: 

In this vector, n_it is the number of females of age i at time t .

The population structure at the next time period t + 1 is then given by equation:

In this equation, F₀, F₁, F_m-1, F_m are the fecundities of females of different ages, i.e. F_i is the number of daughters born to a female of age i which survive to the next time interval, and so contribute to the number n_0,t+1. Similarly, Pi is the probability that a female of age i will survive to age i + 1, i.e. n_i+1,t+1 = P_in_it.

The matrix notation is particularly convenient for the purposes of calculation and for the analysis of the basic properties of ecological systems (Pielou, 1969; Williamson, 1972). Usher (1972) gives an extensive account of the developments of the Leslie matrix model for a wide range of applications, including predator-prey systems.

5.4.3 Energy flow and nutrient cycles

In ecosystem research, we are often concerned with such processes as the cycling of nutrients and the flow of energy. Dynamic models can be well adapted so as to deal with such processes, representing the input of energy or nutrients into the ecosystem, and the transfer of energy or nutrients within the ecosystem. The natural compartmentation of the ecosystem into its species composition or into its trophic levels facilitates the formation of equations describing the distribution of energy or nutrients in the ecosystem at any given time (Phillipson, 1966). Losses from the ecosystem may be specified explicitly, or they may be assumed to be the difference between input and the sum of the output from, the storage in, anyone compartment. Initially, it is simpler to assume that the method of creating compartments for an ecosystem is through trophic levels.

Figure 8 gives a simple example of the flow of phosphorus in a three- compartment system taken from Smith (1970). The parameters of this model are defined as:

Figure 8. Representation of the phosphorus cycle in three compartment ecosystem

x_i = the amount of phosphorus in the i^th compartment at any specified time; 
a_i = the rate of inflow of phosphorus into the i^th compartment;
f_ij = the rate of flow of phosphorus from the i^th to the j^th compartment;
f_ij = the proportion of phosphorus that is stored in the i^th compartment during a given period of time.
 z = the rate of outflow of phosphorus from the i^th compartment.

In matrix notation, the equations for this model are given in Figure 9. The two vectors in this equation denote X_t,i of energy or nutrients in the ith compartment at time t. The matrix elements of the form f_ij(i ¹ j) denote the transitions between the ith and jth compartments. The elements in the leading diagonal are composed of two factors, namely the energy or nutrient not transferring between compartments (f_i) and the input to the ith compartment, which is, in this model, independent of the amount of energy or nutrients already in that compartment, and is therefore represented as a_i/x_t,i. For the parameters summarized in Table 8, it can be shown that this ecosystem is in a steady state, with no increase in the total amount of phosphate.

Usher (1972) summarizes methods of analysis appropriate to matrix models of energy and nutrient cycles, and gives an example of the energy flow in a four-compartment ecosystem. Alternatively, the algorithm in the Appendix can be used to explore the changes in the state variables.

Figure 9. Matrix representation of phosphorus cycle

Table 8 Parameter values for three-compartment ecosystem

5.5 COMPUTATIONAL SOLUTIONS OF DYNAMIC MODELS

Scientists with a sound knowledge of computer programming can quite easily write efficient algorithms for the computational solution of dynamic models. With the development of high-level computing languages like FORTRAN or BASIC, it is a relatively simple matter to write down the difference equations describing the effects of flows on the system variables, and subroutines can be used recursively in order to obtain the solution of dynamic models, and also to explore their consequences for the management of natural systems. Some of these simple algorithms are given in this Handbook.

However, as interest has grown in the use of dynamic models in a wide range of applications, special-purpose languages have been developed to simplify the task of programming the computers used in deriving the solutions. One of the earliest of these languages was CSMP, standing for Continuous System Modeling Program. CSMP is a superset of the general-purpose language FORTRAN, and, as a result, is generally well maintained and supported by computer centres. The practical advantages of the language, and of other similar languages, is that the rate equations can be written down in a relatively simple way, and without any regard to the order in which they must be calculated. Given the necessary parameters for the model and the time limits within which the solutions are required, the translation of the model into the underlying FORTRAN is performed automatically. An added advantage of CSMP, in particular, is that it contains a wide choice of approximate integrating procedures which can be selected by the modeller, depending on the nature and purpose of the model. This opportunity to choose can be a most attractive characteristic in modelling biological systems, with their frequently complex and interacting continuous relationships. Like most languages of this type, it incorporates automatic time-keeping routines and sophisticated output facilities, including the production of graphs and diagrams. Anyone wanting to use CSMP for the modelling of dynamic systems should consult their local computing centre. The books by de Wit and Goudriaan (1974), Dent and Blackie (1979) and Brockington (1979) all contain examples of CSMP programs.

A perhaps even more widely used integrating simulation language is DYNAMO (an acronym for DYNAmic MOdels), which is available for most mainframe and microcomputers. Although DYNAMO is similar in appearance to CSMP, it has a more limited choice of integrating procedures. In contrast, DYNAMO provides excellent facilities for the exploration of models. This exploration is achieved by examining a whole series of specific simulations of the model with variations in parameter values, exogeneous data, model structure, etc., and provides an insight into overall system performance. Again, anyone wishing to use DYNAMO should consult their local computer centre. An excellent introduction to computer simulation through the use of DYNAMO is given by Roberts et al. (1983).

With the rapid development of microprocessors, improved special languages for the modelling of dynamic processes can be expected to be developed, exploiting, especially, the high-resolution graphics that have become almost commonplace on such machines. On the other hand, the interactive facilities provided by microprocessors for both model development and model exploration make it even more desirable than before that scientists learn to program computers for themselves. New general- and special-purpose languages are constantly being devised which will make it possible for dynamic processes to be modelled and explored.

5.6 CASE STUDIES IN DYNAMIC MODELLING

As dynamic modelling is by far the most popular technique for modelling change in ecosystems, there are numerous case studies-certainly too many to enumerate in full. For the interested reader, however, there are some useful collections, and notably the following.

Hall and Day (1977) have written introductory chapters on modelling theory and procedures, followed by case studies using models to analyse natural systems, assess environmental impact, and design optimal (or, at least, better) interactions between man and nature. Four of the introductory chapters are of particular value. They begin with definitions of a system and systems analysis and then go through some basic steps generally used to develop conceptual, diagrammatic, mathematical and computer models. The presentation reflects the authors' belief that the mechanics of building a model are relatively simple, particularly when compared with the complexity of determining the important attributes and relations within the ecosystem. Hall et al. (1977) describe a symbolic language intended to assist in the conveyance of complicated and voluminous ecological information and in the analysis of complex systems. The energy flow language is based on a series of modules, but represents both system processes and mathematical functions, connected by lines representing transfer particles of energy, materials, or information. Overton (1977) emphasizes the structure and strategic aspects of the model building process, with little attention to tactics. The tactical matters of solving equations and choosing the best form for an explicit relationship are technical problems, but the availability of suitable tactics limits the possible strategies, and the selection of a strategy constrains the tactics so that the two aspects are closely interrelated. Finally, Shoemaker (1977a, b) describes the fundamental mathematical steps in the construction and solutions of models of an ecological system intended to understand and predict the behaviour of the system as a whole by quantitatively describing interactions between parts of the system. The discussion of mathematical models is limited to those models designed to use data as well as ecological principles to predict the behaviour of specific ecological systems.

Of the remaining papers describing particular case studies, one of the most important is that by Botkin (1977) which describes a computer model of forest growth, designed to simulate the growth of individual trees on small forest plots. Experience with the model suggests that is reproduces the major dynamic characteristics of a forest community, and that it can be used to construct simulated experiments about the response of a forest to perturbations and manipulation. These experiments given insight into the importance of species interactions to the success and survival of anyone species, and to the persistence of the entire forest community. Further details of this model are given by Botkin et al. (1972a, b).

Jeffers (1972) contains the proceedings of a symposium on the use of mathematical models in ecology and related sciences. The volume contains some important papers, ranging from the philosophical discussion of various aspects of mathematical modelling in empirical science by Skellam (1972) to the detailed description of developments of the Leslie matrix model by Usher (1972) with its extensive bibliography, already quoted extensively in this Handbook. Goodall (1972) discusses the building of dynamic models of ecosystems for predictive purposes in terms of the identification of state variables, the processes involved in their changes and factors influencing their rate of change. Conway and Murdie (1972) illustrate the use of population models as an aid to pest control at several different levels, with examples drawn from models of component processes, including reproduction and predator-prey interactions, and population models of mosquitoes and red bollworm. Other papers emphasize the importance of parameter sensitivity (Plinston, 1972), stochastic model fitting by evolutionary operation (Ross, 1972) and the use of problem-orientated packages in the formulation of ecological models (Radford, 1972), in addition to a wide variety of applications to both large-scale and small-scale ecological problems.

Perhaps the most extensive set of case studies, however, occurs in the four volumes entitled Systems Analysis and Simulation in Ecology, edited by B. C. Patten. Patten's primer for ecological modelling and simulation with analogue and digital computers remains one of the best introductions to dynamic modelling, and gives a sufficient treatment of programming elements to permit ecological models of no small significance (systems of coupled differential equations) to be implemented effectively for simulation or systems analysis studies on both digital and analogue computers (Patten, 1971). Other chapters in the same volume, for example O'Neill (1971), apply these techniques to actual prediction of ecological phenomena, with varying success.

The second volume of the series (Patten, 1972) contains chapters describing applications of selected dynamic analysis procedures, in addition to a section on ecological systems theory, and an extensive final section on applications of social relevance. The theoretical discussions include, for the first time, such topics as the steady-state equilibrium in simple non-linear food webs, the structural properties of food webs, and the concept of niche pattern. The third volume (Patten, 1975) is broadly divided between a discussion of ecosystem modelling in the US International Biological Programme, with chapters on grassland, eastern deciduous forest and desert, tundra and coniferous forest biomes, and models of freshwater-estuarine ecosystems. The case studies contain considerably more detail than those in earlier volumes, consistent with the general development of the dynamic modelling techniques.

The fourth volume of the series (Patten, 1976) is divided into two main sections, the first dealing with models of estuarine marine ecosystems, terresterial ecosystems, and human ecosystems, and the second dealing with the special problems and theory of ecosystem modelling. The four models together provide a succinct review of the development of dynamic modelling of ecosystems over the most active period of this development, and the papers contained in these volumes provide an extensive bibliography of ecosystem modelling. The final volume, in particular, amply demonstrates the current scope of systems ecology from its past and present emphasis on the parts and mechanisms in simulation modelling, and its movement towards systems analysis and novel, more formal consideration of ecosystem theory. The seventeen chapters make clear that, although the systems approach is young in ecology, it has substantially enriched the science, both methodologically and conceptually.

One of the largest collections of mathematical models of air pollution, water pollution, ecology and other environmental topics is the proceedings volume of a conference on environmental modelling simulation held in Cincinnati, Ohio, in Apri11976. Edited by W. R. Ott (1976), this volume contains 164 papers, covering a very wide range of environmental concerns. Many of these papers are too short to give adequate detail of the formulation and use of models in the environmental sciences, but, collectively, they provide a very extensive bibliography to guide further reading and investigation within highly specialized fields of application.

Jorgensen (1979a) also contains the proceedings of a conference, this time on the state-of-the-art in ecology modelling, held in Copenhagen in 1978. The first part of this book contains invited state-of-the-art reviews of river models (Hahn and Eppler, 1979), graphical methods (Maguire, 1979), microcosms (Lassiter, 1979), predator-prey (Dubois, 1979), water quality and irrigation (Lassiter et al., 1979; Skogerboe et al., 1979; Jorgensen, 1979b). The second part contains 29 original papers presented at the conference. Again, such a volume provides a review of present knowledge, and an introduction to the relevant specialist literature.

Shugart and O'Neill (1979) also provide a series of benchmark papers, i.e. studies which are thought to have had the greatest influence on later research. The authors have selected some early papers that had a significant influence, as well as more recent papers that represent innovative approaches. They have also included several papers that contain information or points of view that are important for the understanding of systems ecology.

Innis and O'Neill (1979) address some of the problems of constructing and analysing ecosystem models, and present the papers delivered on this topic at the Satellite Programme in Statistical Ecology held at various locations in 1976-1978.

Finally, for the reader interested in further exploration of the systems ecology literature, and particularly of models available only as internal reports, the bibliography of O'Neill et al. (1977) contains over 900 references to the modelling literature.