1. Continuous Population Models for Single Species The increasing study of realistic mathematical models in ecology (basically the study of the relation between species and their environment) is a reflection of their use in helping to understand the dynamic processes involved in such areas as predator-prey and competition interactions, renewable resource management, evolution of pesticide resistant strains, ecological control of pests, multi-species societies, plant-herbivore systems and so on. The continually expanding list of applications is extensive. There are also interesting and useful applications of single species models in the biomedical sciences: in Section 1.5 we discuss two practical examples of these which arise in physiology. Here, and in the following three chapters, we shall consider some deterministic models. The book edited by May (1981) gives an overview of theoretical ecology from a variety of different aspects; experts in diverse fields review their areas. The book by Nisbet and Gurney (1982) is a comprehensive account of mathematical modelling in population dynamics: a good elementary introduction is given in the textbook by Edelstein-Keshet (1988). 1.1 Continuous Growth Models Single species models are of relevance to laboratory studies in particular but, in the real world, can reflect a telescoping of effects which influence the population dynamics. Let N(t) be the population of the species at time t, then the rate of change dN -t— = births — deaths 4- roioratioti - (11) dt " ° x ' is a conservation equation for the population. The form of the various terms on the right hand side of (1.1) necessitates modelling the situation that we are concerned with. The simplest model has no migration and the birth and death terms are proportional to N. That is ^ = bN~~dN N{t) = No e^11* dt where 6, d are positive constants and the initial population N(Q) — No. Thus if 6 > d the population grows exponentially while if 6 < d it dies out. This 2 1. Continuous Population Models for Single Species approach, due to Malthus in 1798 but actually suggested earlier by Euler, is pretty unrealistic. However if we consider the past and predicted growth estimates for the total world population from the 17th to 21st centuries it is perhaps less unrealistic as seen in the following table. Date Mid 17th Early 19th 1918-27 1960 1974 1987 1999 2010 2022 Century Century Population in billions 0.5 1 2 3 4 5 6 7 8 In the long run of course there must be some adjustment to such exponential growth. Verhulst in 1836 proposed that a self-limiting process should operate when a population becomes too large. He suggested dN dt = rN(l - N/K) , (1.2) where r and K are positive constants. This is called logistic growth in a population. In this model the per capita birth rate is r(l — N/K), that is, it is dependent on JV. The constant K is the carrying capacity of the environment, which is usually determined by the available sustaining resources. There are two steady states or equilibrium states for (1.2), namely JV = 0 and N ~ K, that is where dN/di = 0. JV = 0 is unstable since linearization about it (that is JV2 is neglected compared with JV) gives dN/dt rJV, and so N grows exponentially from any initial value. The other equilibrium N = K is stable: linearization about it (that is (JV — K)2 is neglected compared with |JV - K\) gives d(N - K)/dt « -r(N - K) and so JV -> K as t -+ oo. The carrying capacity K determines the size of the stable steady state population while r is a measure of the rate at which it is reached, that is, it is a measure of the dynamics: we could incorporate it in the time by a transformation from t to rt. Thus 1/r is a representative time scale of the response of the model to any change in the population. N(i) K No K/2 0 Kg. 1.1. Logistic population growth. Note the qualitative difference for the two cases JVb < K/2 and K > Na > K/2. 1.1 Continuous Growth Models 3 If N(0) = No the solution of (1.2) is "rc-gt^Sff-ilT* M t~iM- (L3) and is illustrated in Fig. 1.1. Prom (1.2), if No < K, N{t) simply increases monotonically to K while if No > K it decreases monotonically to K. In the former case there is a qualitative difference depending on whether No > K/2 or Nq < K/2: with Nq < K/2 the form has a typical sigmoid character, which is commonly observed. In the case where Nq > K this would imply that the per capita birth rate is negative! Of course all it is really saying is that in (1.1) the births plus immigration is less than the deaths plus emigration. The point about (1.2) is that it is more like a metaphor for a class of population models with density dependent regulatory mechanisms - a kind of compensating effect of overcrowding - and must not be taken too literally as the equation governing the population dynamics. It is a particularly convenient form to take when seeking qualitative dynamic behaviour in populations in which N = 0 is an unstable steady state and N(t) tends to a finite positive stable steady state. The logistic form will occur in a variety of different contexts throughout the book. In general if we consider a population to be governed by f-/W, 0. This is clear from linearizing about N* by writing n(t) & N{t) - N*, K*)|<1 and (1.4) becomes dn m = f{N* + n) » f(N*) + nf'(N*) + ... , which to first order in n(i} gives dn dt nf'(N*) n(t) K exp[/'(JV*)i] . (1.5) So n grows or decays according as f(N*) > 0 or f'(N*) < 0. The time scale of the response of the population to a disturbance is of the order of l/J/'(jV*)|: it is the time to change the initial disturbance by a factor e. There may be several equilibrium, or steady state populations N* which are solutions of f(N) = 0: it depends on the system f(N) models. Graphically plotting f(N) against N immediately gives the equilibria. The gradient f'(N*) Unstable 0 n Fig. 1.2. Population dynamics model dN/dt = /(n) with several steady states. The gradient f'(N) at the steady state, that is where f(N) = 0, determines the linear stability. at each steady state then determines its linear stability. Such steady states may, however, be unstable to finite disturbances. Suppose, for example, that f(N) is as illustrated in Fig. 1.2. The gradients f'{N) at N = 0,N% are positive so these equilibria are unstable while those at N = JVi,iV3 are stable to small perturbations: the arrows symbolically indicate stability or instability. If, for example, we now perturb the population from its equilibrium N\ so that N is in the range N2 < N < Nz then N -+ N% rather than returning to Ni. A similar perturbation from N3 to a value in the range 0 < N < N2 would result in N(t) —> N\. Qualitatively there is a threshold perturbation below which the steady states are always stable, and this threshold depends on the full nonlinear form of f{N). For JVi, for example, the necessary threshold perturbation is Ni-Ni. 1.2 Insect Outbreak Model: Spruce Budworm A practical model which exhibits two positive linearly stable steady state populations is that for the spruce budworm which can, with ferocious efficiency, defoliate the balsam fir: it is a major problem in Canada. Ludwig et al. (1978) considered the budworm population dynamics to be governed by the equation P{N) 0 n Fig. 1.3. Typical functional form of the predation in the spruce budworm model: note the sigmoid character. The population value JVe is an approximate threshold value. For N < Nc predation is small, while for N > Nc it is "switched on". 1.2 Insect Outbreak Model: Spruce Budworm 5 Here r_ is the linear birth rate of the budworm and Kg is the carrying capacity which is related to the density of foliage available on the trees. The p(iV)-term represents predation, generally by birds: the qualitative form of it is important and is illustrated in Fig. 1.3. Predation usually saturates for large enough N. There is an approximate threshold value Nc, below which the predation is small, while above it the predation is close to its saturation value: such a functional form is like a switch with Nc being the critical switch value. For small population densities N, the birds tend to seek food elsewhere and so the predation term p(N) drops more rapidly, as N —> 0, than a linear rate proportional to N, To be specific we take the form for p(N) suggested by Ludwig et al. (1978) namely BN2/(A2 + N2) where A and B are positive constants, and the dynamics of N(t) is then governed by dN ( N\ BN2 „ c. „=ra^^_„j _______ (1.6) This equation has four parameters, rg, Kg, B and A, with A and Kg having the same dimensions as N, rg has dimension (time)-1 and B has the dimensions of AT(time)"1. A is a measure of the threshold where the predation is 'switched on', that is Nc in Fig. 1.3. If A is small the 'threshold' is small, but the effect is dramatic. Before analysing the model it is essential, or rather obligatory, to express it in nondimensionol terms. This has several advantages. For example, the units used in the analysis are then unimportant and the adjectives small and large have a definite relative meaning. It also always reduces the number of relevant parameters to dimensionless groupings which determine the dynamics. A pedagogical article with several practical examples by Segel (1972) discusses the necessity and advantages for nondimensionalisation and scaling in general. Here we introduce nondimensional quantities by N ArS Kb Bt " = T r = ^' q=~A> T ~ ~A (L7) which on substituting into (1.6) becomes ^ = rH1_^"iT^^/(u;r'?)' (L8) where / is defined by this equation. Note that it has only two parameters r and q, which are pure numbers, as also is u of course. Now, for example, if u < 1 it means simply that N ■< A. In real terms it means that predation is negligible in this population range. In any model there are usually several different nondimensionalisations possible. The dimensionless groupings to choose depends on the aspects you want to investigate. The reasons for the particular form (1.7) will become clear below. The steady states are solutions of u 2 f{u;r,q) = 0 =*- ru 1+u2' (1.9) Clearly u = 0 is one solution with other solutions, if they exist, satisfying Pig. 1.4a,b. Equilibrium states for the spruce budworm population model (1.8). The positive equilibria are given by the intersections of the straight line r(l — u/q) and u/(l + u2). With the solid straight line in (a) there are 3 steady states with f(u;r,q) typically as in (b). Although we know the analytical solutions of a cubic (Appendix 2), they are often clumsy to use because of their algebraic complexity; this is one of these cases. It is convenient here to determine the existence of solutions of (1.10) graphically as shown in Fig. 1.4 (a). We have plotted the straight line, the left of (1.10), and the function on the right of (1.10); the intersections give the solutions. The actual expressions are not important here. What is important, however, is the existence of one, three, or again, one solution as r increases for a fixed q, as in Fig. 1.4 (a), or as also happens for a fixed r and a varying q. When r is in the appropriate range, which depends on q, there are three equilibria with a typical corresponding f(u; r,q) as shown in Fig. 1.4 (b). The nondimensional groupings which leave the two parameters appearing only in the straight line part of Fig. 1.4 is particularly helpful and was the motivation for the form introduced in (1.7). By inspection u = 0, u = ui are linearly unstable, since df/du > 0 at u — 0,w2> while ui and u3 are stable steady states, since at these df/du < 0. There is a domain in the r,q parameter space where three roots of (1.10) exist. This is shown in Fig. 1.5: the analytical derivation of the boundary curves is left as an exercise (Exercise 1). (1.10) 1.3 Delay Models One of the deficiencies of single population models like (1.4) is that the birth rate is considered to act instantaneously whereas there may be a time delay to take account of the time to reach maturity, the finite gestation period and so on. We can incorporate such delays by considering delay differential equation models of the form ^~p-=f(m,x(t-T))} ci.li) where T > 0, the delay, is a parameter. One such model, which has been used, is an extension of the logistic growth model (1.2), namely the differential delay VUU14U1V1I M dt = rN{t) 1 - N{t-T)~ K (1.12) where r, K and T are positive constants. This says that the regulatory effect depends on the population at an earlier time, t — T, rather than that at t. This equation is itself a model for a delay effect which should really be an average over past populations and which results in an integro-differential equation. Thus a more accurate model than (1.12) is, for example, the convolution type model dN dt = rN{t) 1 - K'1 f J—c w(t - s)N(s) ds (1.13) where w(t) is a weighting factor which says how much emphasis should be given to the size of the population at earlier times to determine the present effect Kg. 1.7. Typical weighting function w(t) for an integrated delay effect on growth limitation for the delay model represented by (1.13). on resource availability. Practically w(i) will tend to zero for large negative and positive t and will probably have a maximum at some representative time T. Typically w(t) is as illustrated in Fig. 1.7. If w(t) is sharper in the sense that the region around T is narrower or larger then in the limit we can think of w(t) as approximating the Dirac function 6(T — t), where j —I 6(T-t)f{t)dt = f{T). Equation (1.13) in this case then reduces to (1.12) j—i 8{t-T~ s)N{s) ds=N(t- T) The character of the solutions of (1.12), and the type of boundary conditions required are quite different to those of (1.2). Even with the seemingly innocuous equation (1.12) the solutions in general have to be found numerically. Note that to compute the solution for t > 0 we require N(t) for all —T < t < 0. We can however get some qualitative impression of the kind of solutions of (1.12) which are possible, by the following heuristic reasoning. Refer now to Fig. 1.8 and suppose that for some t = ti, N(ti) — K and that for some time t < *i, N(t — T) 0, dN(t)/dt > 0 and so N(t) at ti is still increasing. When < = ii+T, N(t-T) = N(ti) = K and so dN/dt = 0. For *j + T < t < i2, N(t ~T) > K and so dN/di < 0 and N(t) decreases until t = ii + T since then dN/di = 0 again because N(h + T-T) = N(t2) = K. There is therefore the possibility of oscillatory behaviour. For example, with the simple linear delay equation dN ir „, m. , . 7r£ -— = -—N(t - T) ^ N(t) = A cos — , dt 2t k ' w 2T which is periodic in time and which can be easily checked is a solution. 1.6 Harvesting a Single Natural Population It is clearly necessary to develop an ecologically acceptable strategy for harvesting any renewable resource be it animals, fish, plants or whatever. We also usually want the maximum sustainable yield with the minimum effort. The inclusion of economic factors in population models of renewable resources is increasing and these introduce important constraints: see for example the book by Clark (1976a). The collection of papers edited by Vincent and Skowronski (1981) specifically deals with renewable resource management. The review article by Plant and Mangel (1987) is concerned with insect pest management. The model we describe here is a simple logistic one with the inclusion of a harvesting contribution: it was discussed by Beddington and May (1977). Although It is a particularly simple one it brings out several interesting and important points which more sophisticated models must also take into account. Rotenberg (1987) also considered the logistic model with harvesting, with a view to making the model more quantitative. He extinction. Most species have a growth rate, depending on the population, which more or less maintains a constant population equal to the environment's carrying capacity K. That is the growth and death rates are about equal. Harvesting the species affects the mortality rate and, if it is not excessive, the population adjusts and settles down to a new equilibrium state < K. The modelling problem is how to maximize the sustained yield by determining the population growth dynamics so as to fix the harvesting rate which keeps the population at its maximum growth rate. We discuss here a basic model which consists of the logistic population model (1.2) in which the mortality rate is enhanced, as a result of harvesting, by a term 1.6 Harvesting a Single Natural Population 25 linearly proportional to JV, namely, TW (l ~ |^ m = /(JV). (1.41) Here r, K and J3 are positive constants and EN is trie harvesting yield per unit time with E a measure of the effort expended. K and r are the natural carrying capacity and the linear per capita growth rate respectively. The new non-zero steady state from (1-41) is Nk{E) = K[l--)>Q if E < r (1.42) which gives a yield Y(E) = ENh{E) = EK (l - ~) . (1.43) Clearly if the harvesting effort is sufficiently large so that it is greater than the linear growth rate when the population is low the species will die out. That is if E > r the only realistic steady state is JV = 0. If E < r (which was possibly not the case, for example, with whaling in the 1970's) the maximum sustained yield and the new harvesting steady state are, from (1.43) and (1.42), YM=Y(E)]E=r/2 = —, Nk]YM = ~. . (1.44) Does an analysis of the dynamic behavior tell us anything different from the naive, and often used, steady state analysis just given here? Fig. 1.15 illustrates the growth rate f{N) in (1.41) as a function of N for various efforts E. Linearising (1.41) about Nh{E) gives d{NjtNk) « tWJBW* - m ft®- r)(N - Nh) , (1.45) which shows linear stability if E < r: arrows indicate stability or instability in Fig. 1.15. We can consider the dynamic aspects of the process by determining the time scale of the recovery after harvesting. If E = 0 then, from (1.41), the recovery time T — 0(1/r), namely the time scale of the reproductive growth. This is the order of magnitude of the recovery time of N to its carrying capacity K after a small perturbation from K since, for N(t) — K small and iVfc(O) = K, (1.45) shows d^N~ K) ~ -r(N - K) => N{t)-K K/2, the equilibrium population for the maximum yield Ym- The recovery time ratio Tr(E)/Tr(0) from (1-47) is then approximately 1. So increasing 22, and hence the yield, we are on branch £+. As E increases further, Nk(E) decreases towards K/2, the value for the maximum sustained yield Ym and we reach the point A 1.6 Harvesting a Single Natural Population 27 in Fig. 1.16 (a) when Nh(E) = if/2. As E is increased further, Nh(E) < K/2 and the recovery time is further increased but with a decreasing yield: we are now on the L_ branch. We can now see what an optimal harvesting strategy could be, at least with this deterministic point of view. An effort E should be made which keeps the equilibrium population density N^{E) > K/2, but as close as possible to K/2, the value for the maximum sustained yield. The closer to K/2 however the more delicate the situation becomes since we might inadvertantly move onto branch i_ in Fig. 1.16 (a). At this stage, when N^E) is close to K/2 a stochastic analysis should be carried out as has been done by Beddington and May (1977). Stochastic elements of course reduce the predictability of the catch. In fact they decrease the average yield for a given effort. (a) (b) Fig. 1.16a,b. (a) Ratio of the recovery times as a function of the yield for the logistic growth model, with yield proportional to the population: equation (1-41). (b) Graphical method for determining the steady state yield Y for the harvested logistic model (1.41). As an alternative harvesting resource strategy suppose we harvest with a constant yield Yq as our goal. The model equation is then dN ( N\ Fig. 1.17 (a) shows the graphical way of determining the steady states as Yq varies. It is trivial to find the equilibria analytically of course, but often the behavioral traits as a parameter varies are more obvious from a figure, such as here. If 0 < Yo < = Ym> the maximum sustainable yield here, there are two positive steady states Ni(Yq), Ni(Yq) > iVi(Fo) which from Fig. 1.17 are respectively unstable and stable. As Yq —> rK/4, the maximum sustainable yield of the previous model, this model is even more sensitive to fluctuations since if a perturbation from N2 takes N to a value N < iVj the mechanism then drives N to zero: see Fig. 1.17 (b). Not only that, N —y 0 in a finite time since for 28 1. Continuous Population Models for Single Species Fig. 1.17a,b. (a) Equilibrium states for the logistic growth harvested with a constant yield Yq: equation (1.49). (b) Growth rate f(N) in (1.49) as the yield Y0 increases. small enough N, (1-49) becomes dN/dt rj —Yq and so for any starting Nq at to, For easy comparison with the constant effort model we evaluate the equivalent recovery time ratio Tr(Yq)/Tr(0). The recovery time Tr(Yq) is only relevant to the stable equilibrium N2{Yq) which from (1.49) is N2{Y0) = K 1 + Kr 1/2' F0 < rüT/4 . The linearized form of (1.49) is then d(N-N2(Y0)) dt (N - N2) df ON -(N-N2)r 1- rK 1/2 Thus Tr(Yq) Ym) 1/2- rK (1.50) which shows that Tr(Yo)/Tr(0) —* co as Fo —> Ym- This model is thus a much more sensitive one and, as a harvesting strategy, is not very good. The main conclusion from this modelling exercise is that a constant effort rather than a constant yield harvesting strategy is less potentially catastrophic. It calls into question, even with this simple model, the fishing laws, for example, which regulate catches. A more realistic model, on the lines described here, should take into account the economic costs of harvesting and other factors. This implies a feedback mechanism which can be a stabilizing factor: see Clark (1976a). With the unpredictability of the real world it is probably essential to include feedback. Nevertheless such simple models can pose highly relevant ecological and long term financial questions which have to be considered in any more realistic and more sophisticated model.