Evaluating predator prey models
Duncan Golicher
23/07/2021
Introduction
This document provides a guide for your project. It is not citeable material and is written in an informal style to provide advice on how to address the project. Do not quote directly from it, but do feel free to rewrite the ideas expressed in your own words following the style of a scientific paper.
Introduction
The response to the pandemic of 2020/21 was heavily influenced by the use of ecological modelling in order to predict the dynamics of a complex system. Much of the modelling was derived from work conducted by a team of modellers at Imperial college(Ferguson et al. 2006). This has heightened public awareness of both the strengths and potential weaknesses of ecological modelling. A particular challenge arises when models are used to predict the result of some change to a complex system in the future. As data on the future result of the intervention are clearly not available, the model cannot be validated directly. In many respects expecting any model to predict the future is asking far too much. Ecological models provide frameworks for thought. If used carefully they can be useful for improving understanding of a complex system. If used without understanding they may produce misleading results. This directed IRP provides an opportunity to look at ecological modelling in some depth. There is a very rich literature on the subject to be explored with many differing viewpoints and many critiques of models.
The IRP is based on the use of an adapted classical ecological model to predict the impact of introducing pine martens in order to control grey squirrels. Modelling provides a good source of material for a project as, unlike field work, models can be relied upon to produce “data.” The thought that goes into parameterising a model can be at least as useful as the model results themselves, and often much more so. All models are wrong, but some are useful. They are only useful if their limitations are well understood.
The Lotka Volterra predator prey model
The classic “text book” model that links the dynamics of a prey population to that of a predator population was developed by Lotka and Volterra in 1927 (Berryman 1992). It is extremely simple.
\(\frac {dN}{dt} = rN-\alpha NP\)
\(\frac {dP}{dt} = \alpha \beta NP-\gamma P\)
There are many explanations of the derivation of these equations in text books and on the internet. They can be written in several different ways, but the basic idea is always the same. I’ll explain it in simple terms here.
\(\frac {dN}{dt}\) and \(\frac {dP}{dt}\) are the rates of change of the prey and predator populations with respect to time. If the model time step in a year the equations represent the increase or decrease of the numbers of prey and predators over a year. This is a dynamic model, so the rates change as the model runs as a result of changes in both populations simultaneously. It is this aspect of the model that makes the dynamic interesting and hard to predict without running the model as a numerical simulation. You can’t do the calculations in your head!
The \(r\) parameter represents the intrinsic rate of growth of the prey population \(N\) in the absence of predation. If there were no predators then the prey population would increase exponentially. However it doesn’t do this, as prey are consumed by predators in the model.
The \(\alpha\) parameter controls mortality due to hunting. The loss of prey is modelled by multiplying \(\alpha\) by the number of prey and the number of predators. To understand this, think about the case when there is just one predator and 100 prey. If the predator attack and kill 10% of the prey (i.e. ten prey) then \(\alpha\) would be 0.1. The problem with this in terms of realism is that if there were 1000 prey then a single predator would be assumed to catch 100 of them. If \(\alpha\) is set to be 1 then this results in one single predator wiping out the whole prey population in one time step! So behavioural aspects are not being included in the model. This can be adjusted for through including a functional response model, at the cost of making the model framework more complex and difficult to understand as a whole.
Once the prey are eaten they can be used to increase the predator population. The very simple model assumes that there is a conversion factor \(\beta\) that determines how many new prey are added to the population.
The predator population dies off through an exponential decay process modelled as \(\gamma P\). So if half the predator populations dies each year \(\gamma\) would equal 0.5
Although the text book treatments of these equations can give the impression of great sophistication and lead to beautiful mathematical treatments, in reality the model is very unstable and unrealistic. It is clearly glossly oversimplified and mainly of theoretical interest.
There is no limitation to the growth of the prey population apart from predation. In other words food availability and other similar limiting factors are not included.
Mortality of prey only takes place through predation.
Predators only consume one form of prey. So in the absence of all prey they would become extinct.
The mathematical model allows fractions of individuals. So prey populations can fall below one in the model.
The \(\alpha\) and \(\beta\) parameters lead to a model in which there are no biological contraints to the attack rate, nor to the conversion of prey to predators.
If the simple model is implemented in R using parameters as stated above it behaves like this (code to run the model is embedded in this document).
These extreme cycles can be modified by changing the parameters, but the instability and lack of realism of the model means that it is not a suitable representation of any real life system. It is a mathematical toy. The classic example of rather similar cycles found in nature is the Canadian snowshoe hare and Lynx data provided by trapping records. Although this has frequently been used to justify the model closer inspection of the data led to some suspicions. In fact the cycles didn’t match the predictions at all as the peak in prey numbers occurred after the peak in predators. This led one author to suggest (rather tongue in cheek) that hares must eat Lynx! (Gilpin 1973).
The Lotka Volterra model with prey carrying capacity
One element of the model that can be improved very easily is to assume that the prey population is controlled by some other environmental limitations. To do that we need to know the prey carrying capacity (K).
\(\frac {dN}{dt} = r_1N (1-\frac N K) -\alpha NP\)
\(\frac {dP}{dt} = \alpha \beta NP-\gamma P\)
Now if we run the model we get a much more boring result from the mathematical perspective. When \(\alpha \beta NP = -\gamma P\) then the predator population will stabilise. At some point the prey population will be limited by both the carrying capacity and the off take from the predator population.
The cyclical behaviour changes into one of long term equilibria between predator and prey. However, although less interesting mathematically, this is now much more useful if we are thinking about how changing an ecosystem through introducing (or reintroducing) a predator may lead to an impact though the establishment of a new equilibrium prey population which is below the carrying capacity established by food availability. This may provide opportunities for the food to be exploited by another species and may stabilise boom and bust cycles in the prey.
Alternative prey
How can we take into consideration the possibility that the predator itself does not depend only on the prey species in order to survive? One simple way is to allow the predator to give birth to young without consuming the prey. Additional prey are a bonus. Consuming the prey species does add some predators, but the predator is not dependent on them. This can be added to the model by allowing the predators an intrinsic exponential increase. Providing this is below the intrinsic death rate some prey will still be required in order to allow the population to increase, but the number will be less.
\(\frac {dN}{dt} = r_1N (1-\frac N K) -\alpha NP\)
\(\frac {dP}{dt} = r_2P +\alpha \beta NP-\gamma P\)
An extension to this is to add an additional carrying capacity for the predator species.
Now even though the two systems could be uncoupled if predation levels are low, high levels of predation can lead to cycles again. The predator population can survive both with and without prey. The effect of the predator population on depressing the prey will depend on how much the birth of predators is increased through the consumption of the additional prey.
\(\frac {dN}{dt} = r_1N (1-\frac N K_1) -\alpha NP -\gamma_1N\)
\(\frac {dP}{dt} = r_2P (1-\frac P K_2) +\alpha \beta NP -\gamma_2 P\)
This is quite a complex model that can show a wide range of behaviour again. However this model still leaves out many aspects of the real system.
- The model doesn’t capture switching between prey types as one form of prey increases in abundance relative to another.
- The model does not include a term for a density dependent functional response that includes the concept of handling time which imposes limits to consumption at high prey densities.
- More complex elements such as change in behaviour and “behavioural refuges” (Křivan 2013) are overlooked.
Functional response
Some of the weaknesses of the model with respect to the omission of behaviour can be overcome by including a term for the functional response. This refers to the intrinsic limitation on the number of prey that any single predator can consume. The basic idea is that consumption of prey saturates at some upper limit. Without this the model uses the unrealistic assumption that consumption is always proportional to the number of prey multiplied by the number of predators.
\(N_c = \frac{A N} {F + N}\)
In this case the number of prey consumed is \(N_c\). The maximum number of prey that one predator can consume in a model time step is \(A\) and the number of prey at which the consumption reaches half this asymptotic level is F.
So when this is added to the equations we get the following.
\(\frac {dN}{dt} = r_1N (1-\frac N K_1) -\alpha \frac{A N} {F + N}P -\gamma_1N\)
\(\frac {dP}{dt} = r_2P (1-\frac P K_2) +\alpha \beta \frac{A N} {F + N}P -\gamma_2 P\)
This model can still produce wave like behaviour with some combinations of parameter, but in general the model leads to more stable co-existence, as unrealistc booms in the predator population through overconsumption are avoided.
This in essence is the model you will use for the project, expressed in mathematical terms. You can use these equations in your methods section in order to give more formality to the work. However if you don’t feel comfortable with these dry formal mathematics, don’t worry. The mathematical formulation remains very simple in essence. To develop a more intuitive model and extend it so that parameters are given more direct meaning we will use the mathematics as a framework, but we will implement a model based on them in a more expressive form using compartment flow modelling.
Estimating intrinsic rate of growth using a Leslie matrix
One of the challenges when estimating parameters for a model is finding values for the intrinsic rate of growth of the population. It might be tempting to assume that the rate of growth is simply the birth rate per individual. However in an animal population this would be misleading as many of the individuals added to the population through birth do not survive to reproduce. So this would overestimate the rate of growth, even without taking into account the effect of the carrying capacity. In order to model age structured populations demographers often use a Leslie matrix. The same concept could also be modelled as
## [,1] [,2] [,3] [,4]
## [1,] 0.0 3.0 2.0 1.0
## [2,] 0.5 0.0 0.0 0.0
## [3,] 0.0 0.8 0.0 0.0
## [4,] 0.0 0.0 0.5 0.2## [1] 1.467244Compartment flow modelling
The sets of differential equations that form the Lotka Volterra model, and its modifications are an example of compartment flow models. Modern software can implement the equations in a more intuitive and flexible form using a graphical interface. The basic concept involves considering populations as stocks represented by compartments. Changes to the value held in the compartments takes place through flows in and out. The equations are the flows. The software can produce more flexible models through implementing additional rules to modify flows and can also link to real data. More complex models can be built up in this manner that can be difficult to write down as well defined mathematical equations, even though they are still based around formal mathematics. A free online tool is available for implementing compartment flow models here. https://insightmaker.com/
The site includes a set of tutorials explaining how to build models from first principles. https://insightmaker.com/tutorials
It is worth working through these tutorials in order to learn how to build your own models. However for the purposes of the IRP you are provided with a working model that can be modified if required in order to include more detail. Modifying models is challenging, as the behaviour can be difficult to predict and parameters hard to estimate.
https://insightmaker.com/insight/CNapV0sS7UfmHaN4NrSfZ/IRP-Martens-Squirrels
Challenges
There are many challenges involved when applying models to real life systems. Models may appear to be mathematically sophisticated but they still omit many key aspects of the real system. You may well conclude that models are not suitable tools at all. This would be a reasonable conclusion, but to justify it you do need to gain a good understanding of the model first. One benefit of modelling is that it helps to point out parameter of a system that need to be measured, or at least estimated, in order to manage it effectively. The carrying capacity is such a parameter. Most populations rarely exist at their theoretical carrying capacity. In the presence of predators a prey population will ususually be depressed below the level it could reach in their absence. There are many benefits for management in depressing populations in order to prevent them exhausting their food supply. In some situations predators can even elliminate an undesirable prey species completely, but this does require an alternative food source for the predator if maintenance of the predator is desirable.
Squirrels and pine martens
The system we will look at potentially involves at least three key species, and in reality would include additional ones. The issue involves the potential of reintroduction of pine martens as a means of depressing or even eliminating grey squirrels in order to allow populations of red squirrels to be maintained. A paper by Jones et al. (2017) provides a good template for the study. In this paper squirrel pox is also included as a key factor involved in the interaction between grey squirrel populations and red squirrels. Grey squirrels outcompete red squirrels in many habitats but may also harbour squirrelpox, which is typically an asymptomatic infection which is harmless to grey squirrels but can produce pathological disease in red squirrels.
This emphasises the complexity of ecosystems and the potential for overlooking important drivers when only considering trophic interactions. You may want to consider these additional elements in your discussion in order to critically appraise the strengths and weaknesses of modelling. However, even if squirrel pox is not considered explicitly in a model it could be assumed that effective control of grey squirrels prior to the reintroduction of red squirrels could produce circumstances in which red squirrel populations may develop without being impacted by the disease. So we will consider a system that consists of grey squirrels, pine martens and additional food sources that allow pine marten populations to develop which prey on grey squirrels. The impact of pine martens on red squirrels can be excluded from the model for the purposes of the current exercise. Again, this is an element that can be added to the work in order to enrich and deepen the discussion.
Carrying capacity
The paper by Jones et al. (2017) provides some estimates of squirrel carrying capacity that may be useful. For the purposes of the exercise you should select an area of the new forest (see map tool added into this document for measurement) and use the estimates from Jones to get an approximate carrying capacity for grey squirrels. This only needs to be a “ball park” figure, as sensitivity analysis can be used to evaluate how the model responds to changes in assumptions.
Pine marten populations
As pine martens are not present in most of the UK in order to estimate the impact of martens on squirrels you will have to look at a range of literature, much of it derived from studies in Europe.
Sensitivity analysis
When models are used to make predictions that could have implications in real life it is vitally important to run sensitivity analysis and conduct a critical appraisal of the models. Model do not need to include all factors of a complex system. This would be too difficult to achieve. However factors that influence the system in ways that are difficult to quantify must add to uncertainty about the models performance.
In the case of the pine marten squirrel model the carrying capacity for the squirrel population includes a multitude of complex factors. The upper bound for the number of squirrels in an area is not simply linked to food supply. Squirrel behaviour also plays a role as does the availability of suitable habitat that squirrels can use. The birth and death rates of squirrels can be estimated within some bounds, but these will vary from year to year. The proportion of the pine marten’s diet that consists of squirrels will also vary greatly. The more elements you consider the more difficult it becomes to trust the model. The way to deal with this is not to use single parameters at all. Instead the model is run using a range of values for all the parameters that are uncertain. Running the model thousands of times produces a range of outcomes. A plausible set of assumptions and parameter estimates may lead to a percentage of the models predicting the extinction of the squirrels and some models predicting extinction of the martens, with a range of other outcomes possible. Although this means that the results of the model cannot be used to suggest any outcome with certainty, the model output is now more useful, as at least it may rule out some implausible scenarios.
(Park, Rabinovich, and Lek 2007)(Rose 1983)