Chapter 3 Modeling with Rates of Change
Chapter 1 provided examples for modeling with rates of change, and Chapter 2 introduced the computational and visualization software R and RStudio, and how we can translate equations with rates of change to understand phenomena. The focus for this chapter will be on taking a contextual description and starting to develop differential equation models for them.
Oftentimes when we construct differential equations from a contextual description we bring our own understanding and knowledge to this situation. How you may write down the differential equation may be different from someone else - do not worry! This is the fun part of modeling: models can be considered testable hypotheses that can be refined when confronted with data. Let’s get started
3.1 Competing plant species and equilibrium solutions
Consider the following context to develop a mathematical model:
A newly introduced plant species is introduced to a region. It competes with another established species for nutrients (and is a better competitor). However, the growth rate of the new species is proportional to the difference between the current number of established species and the number of new species. You may assume that the number of established species is a constant E.
For this problem we will start by naming our variables. Let \(N\) represent number of new species and \(E\) the number of established species. We will break this down accordingly:
“the growth rate of the new species” describes the rate of change, or derivative, expressed as \(\displaystyle \frac{dN}{dt}\).
“is proportional to the difference between the current number of established species and the number of new species” means \(\displaystyle \alpha \cdot (E-N)\), where \(\alpha\) is the proportionality constant. Including this parameter helps to avoid assuming we have a 1:1 correspondence between the growth rate of the new species and the population difference.
“and is a better competitor” helps to explain why the term is \(\displaystyle \alpha \cdot (E-N)\) instead of \(\displaystyle \alpha \cdot (N-E)\). We know that the newly established species will start out in much smaller numbers than \(N\). But since it is a better competitor, we would expect its rate to increase initially. So \(\displaystyle \frac{dN}{dt}\) should be positive rather than negative.
Taking all these assumptions together, Equation (3.1) shows the differential equation to model this context:
\[\begin{equation} \frac{dN}{dt} = \alpha \cdot (E-N) \tag{3.1} \end{equation}\]
You may recognize that Equation (3.1) is similar to Equation (1.4) in Chapter 1 for the spread of Ebola. It is not surprising to have similar differential equations appear in different contexts. We will see throughout this book that it is more advantageous to learn techniques to analyze models qualitatively rather than memorize several different types of models and not see the connections between them.
An interesting solution to a differential equation is the steady state or equilibrium solution. Equilibrium solutions occur where the rates of change are zero. For Equation (3.1), this means that we are solving \(\displaystyle \frac{dN}{dt} = \alpha \cdot (E-N) = 0\). Granted, the expression \(\alpha \cdot (E-N)\) may look like alphabet soup, but it is helpful to remember that \(\alpha\) and \(E\) are both parameters; the steady state occurs when the expression \(E-N\) equals zero, or when \(N = E\). We may consider the new species \(N\) to be established when it reaches the same population level as \(E\). Identifying steady states in a model aids in understanding the behavior of any solutions for a differential equation. Chapters 5 and 6 dig deeper into steady states and their calculation.
3.2 The Law of Mass Action
Our next example focuses on how to generate a model that borrows concepts from modeling chemical reactions. For example let’s say you have a substrate A that reacts with enzyme B to form a product S. One common way to represent this process is with a reaction equation (Equation (3.2)):
\[\begin{equation} A+B \rightarrow S \tag{3.2} \end{equation}\]
Figure 3.1 is a schematic diagram of Equation (3.2):
FIGURE 3.1: Schematic diagram of a substrate-enzyme reaction.
One key quantity is the rate of formation for the product \(P\), which we express by Equation (3.3):
\[\begin{equation} \frac{dP}{dt}= kAB, \tag{3.3} \end{equation}\]
where \(k\) is the proportionality constant or the rate constant associated with the reaction. Notice how we express the interaction between \(A\) and \(B\) as a product - if either the substrate \(A\) or enzyme \(B\) is not present (i.e. \(A\) or \(B\) equals zero), then product \(P\) is not formed. Equation (3.3) is an example of the law of mass action.
Modeling interactions (whether between susceptible and infected individuals, enzymes and substrates, or predators and prey) with the law of mass action is always a good first assumption to understand the system, which can be subsequently refined. For example, if we consider that the substrate might decay, we can revise Figure 3.1 to Figure 3.2:
FIGURE 3.2: Revised schematic diagram of substrate-enzyme reaction with decay of the product \(P\).
In this instance the rate of change of \(P\) would then include a term \(dP\) (Equation (3.4):
\[\begin{equation} \frac{dP}{dt}= kAB - dP \tag{3.4} \end{equation}\]
3.3 Coupled differential equations: lynx and hares
Another example is a system of differential equations. The context is between the snowshoe hare and the Canadian lynx, shown in Figure 3.3. Figure 3.4 also displays a timeseries of the two populations overlaid. Notice how in Figure 3.4 both populations show regular periodic fluctuations. One plausible reason is that the lynx prey on the snowshoe hares, which causes the population to initially decline. Once the snowshoe hare population declines, then there is less food for the lynx to survive, so their population declines. The decline in the lynx population causes the hare population to increase, and the cycle repeats.11
FIGURE 3.3: Examples of lynx and hare - aren’t they beautiful?
FIGURE 3.4: Timeseries of the combined lynx and hare populations. Notice how the populations are coupled with each other.
In summary it is safe to say that the two populations are coupled to one another, yielding a coupled system of equations. But in order to understand how they are coupled together, first let’s consider the two populations separately.
To develop the mathematical model we will make some simplifying assumptions. The hares grow much more quickly than then lynx - in fact some hares have been known to reproduce several times a year. A reasonable assumption for large hare populations is that rate of change of the hares is proportional to the hare population. Based on this assumption Equation (3.5) describes the rate of change of the hare population, with \(H\) as the population of the hares:
\[\begin{equation} \frac{dH}{dt} = r H \tag{3.5} \end{equation}\]
Since the growth rate \(r\) is positive, so then the rate of change (\(H'\)) will be positive as well, and \(H\) will be increasing. A representative value for \(r\) is 0.5 year\(^{-1}\) (Mahaffy 2010; Brady and Butler 2021). You may be thinking that the units on \(r\) seem odd - (year\(^{-1}\)), but that unit on \(r\) makes the term \(rH\) dimensionally consistent to be a rate of change.
Let’s consider the lynx now. An approach is to assume their population declines exponentially, or changes at the rate proportional to the current population. Let’s consider \(L\) to be the lynx population, with the following differential equation (Equation (3.6)):
\[\begin{equation} \frac{dL}{dt} = -dL \tag{3.6} \end{equation}\]
We assume the death rate \(d\) in Equation (3.6) is positive, leading to a negative rate of change for the Lynx population (and a decreasing value for \(L\)). A typical value of \(d\) is 0.9 yr\(^{-1}\) (Mahaffy 2010; Brady and Butler 2021).
The next part to consider is how the lynx and hare interact. Since the hares are prey for the lynx, when the lynx hunt, the hare population decreases. We can represent the process of hunting with the following adjustment to our hare equation:
\[\begin{equation} \frac{dH}{dt} = r H - b HL \end{equation}\]
So the parameter \(b\) represents the hunting rate. Notice how we have the term \(HL\) for this interaction. This term injects a sense of realism: if the lynx are not present (\(L=0\)), then the hare population can’t decrease due to hunting. We model the interaction between the hares and the lynx with multiplication between the \(H\) and \(L\). A typical value for \(b\) is .024 lynx\(^{-1}\) year\(^{-1}\). It is okay if that unit seems a little odd to you - it should be! As before, if we multiply out the units on \(bHL\) we would get units of hares per year.
How does hunting affect the lynx population? One possibility is that it increases the lynx population:
\[\begin{equation} \frac{dL}{dt} =bHL -dL \end{equation}\]
Notice the symmetry between the rate of change for the hares and the lynx equations. In many cases this makes sense - if you subtract a rate from one population, then that rate should be added to the receiving population. You could also argue that there is some efficiency loss in converting the hares to lynx - not all of the hare is converted into lynx biomass. In this situation we sometimes like to adjust the hunting term for the lynx equation with another parameter \(e\), representing the efficiency that hares are converted into lynx:
\[\begin{equation} \frac{dL}{dt} =e \, bHL -dL \end{equation}\]
(sometimes people just make a new parameter \(c=e \, b\), but for now we will just leave it as is and set \(e=0.2\)). Equation (3.7) shows the coupled system of differential equations:
\[\begin{equation} \begin{split} \frac{dH}{dt} &= r H - b HL \\ \frac{dL}{dt} &=e\,bHL -dL \end{split} \tag{3.7} \end{equation}\]
The schematic diagram representing these interactions is shown in Figure 3.5:
FIGURE 3.5: Schematic diagram Lynx-Hare system.
Equation (3.7) is a classical model in mathematical biology and differential equations - it is called the predator-prey model, also known as the Lotka-Volterra model (Lotka 1920, 1926; Volterra 1926).
3.4 Functional responses
In several examples we have seen a rate of change proportional to the current population, as, for example, the rate of growth of the hare population is \(rH\). This is one example of what we would call a functional response. Another type of functional response assumes that the rate reaches a limiting value proportional to the population size, so \(\displaystyle \frac{dH}{dt} = \frac{rH}{1+arH}\). This is an example of a type II functional response. Finally, the type II response has also been generalized (a type III functional response) \(\displaystyle \frac{dH}{dt} = \frac{rH^{2}}{1+arH^{2}}\). Figure 3.6 shows all three functional responses together:
FIGURE 3.6: Comparison between examples of Type I - Type III functional responses. For a Type I functional response the rate grows proportional to population size H, whereas for Types II and III the rate reaches a saturating value.
Notice the limiting behavior in the Type II and Type III functional responses. These responses are commonly used in ecology and predator-prey dynamics and in problems of how animals search for food.
3.5 Exercises
Exercise 3.1 Consider the following types of functional responses:
\[\begin{equation} \begin{split} \mbox{ Type I: } \frac{dP}{dt} &= 0.1 P \\ \mbox{ Type II: } \frac{dP}{dt} &= \frac{0.1P}{1+.03P} \\ \mbox{ Type III: } \frac{dP}{dt} &= \frac{0.1P^{2}}{1+.05P^{2}} \end{split} \end{equation}\]
For each of the functional responses evaluate \(\displaystyle \lim_{P \rightarrow \infty} \frac{dP}{dt}\). Since these functional responses represent a rate of change of a population, what are some examples (hypothetical or actual) in which each of these responses would be appropriate?
Exercise 3.2 A population grows according to the equation:
\[\begin{equation} \frac{dP}{dt} = \frac{P}{1+.05P} -.1P = f(P) - g(P) \end{equation}\]
On the same axis, plot the equations \(f(P)\) and \(g(P)\). What are the two positive values of \(P\) where \(f(P)\) and \(g(P)\) intersect?
Next algebraically determine the two steady state values of \(P\), that is solve \(\displaystyle \frac{dP}{dt}=0\) for \(P\). (Hint: factor a \(P\) out of the expression \(\displaystyle f(P)-g(P)\).)
Does your algebraic solution match your graphical solutions?
Exercise 3.3 A population grows according to the equation:
\[\begin{equation} \frac{dP}{dt} = 2P - \frac{4P^{2}}{1+P^{2}} = r(P)-d(P) \end{equation}\]
- On the same axis, plot the equations \(r(P)\) and \(d(P)\). What are the two positive values of \(P\) where \(r(P)\) and \(d(P)\) intersect?
- Next algebraically determine the two steady state values of \(P\), that is solve \(\displaystyle \frac{dP}{dt}=0\) for \(P\). (Hint: factor a \(P\) out of the expression \(r(P)-d(P)\).)
- Does your algebraic solution match your graphical solutions?
Exercise 3.4 A population grows according to the equation:
\[\begin{equation} \frac{dP}{dt} = \frac{aP}{1+abP} - dP, \end{equation}\]
where \(a\), \(b\), and \(d\) are all positive parameters. Determine the two steady state values of \(P\), that is solve \(\displaystyle \frac{dP}{dt}=0\) for \(P\).
Exercise 3.5 A chemical reaction takes two chemicals \(X\) and \(Y\) to form a substrate \(Z\) through the law of mass action. However the substrate can also disassociate. The reaction schematic is the following:
\[\begin{equation} X + Y \rightleftharpoons Z, \end{equation}\]
where you may define the proportionality constant \(k_+\) as associated with the formation of the substrate \(Z\) and \(k_-\) the disassociation (\(Z\) decays back to \(X\) and \(Y\)).
Write down a differential equation that represents the rate of reaction \(\displaystyle \frac{dZ}{dt}\).
Exercise 3.6 (Inspired from Thornley and Johnson (1990) and Logan and Wolesensky (2009)) For each of the following exercises consider the following contextual situations modeling rates of change. For each problem you will need to:
- Name and describe all variables and parameters;
- Determine a differential equation representing the context;
- Write a brief one-two sentence explanation of why your differential equation models the situation at hand;
- Hand sketch a rough graph of what you think the solution is as a function of time, consistent with the context given.
- The rate of change of an animal’s body temperature is proportional to the difference in temperature between the environment and the current body temperature of the animal.
- A plant grows proportional to its current length \(L\). Assume this proportionality constant is \(\mu\), whose rate also decreases proportional to its current value. You will need to write down a system of two equations with variables \(L\) and \(\mu\).
- A patient undergoing chemotherapy receives an injection at rate \(I\). This injection decreases the rate that a tumor accumulates mass. Independent of the injection, the tumor accumulates mass at a rate proportional to the mass of the tumor.
- A cell with radius \(r\) assimilates nutrients at a rate proportional to its surface area, but uses nutrients proportional to its volume. Determine an equation that represents the rate of change of the radius.
- The rate that a cancer cell divides (increases in amount) is proportional to the number of healthy cells in its surrounding environment. You may assume that a healthy cell has mortality \(\delta_{H}\) and a cancer cell has mortality \(\delta_{C}\). Be sure to write down a system of differential equations for the population of cancer cells \(C\) and healthy cells \(H\).
- The rate that a virus is spread to the population is proportional to the probability that a person is sick (out of \(N\) total sick and healthy individuals).
FIGURE 3.7: Modeled reaction schemes representing the potential effect of a pesticide on water quality.
Exercise 3.7 (Inspired by Burnham and Anderson (2002)) You are tasked with the job of investigating the effect of a pesticide on water quality, in terms of its effects on the health of the plants and fish in the ecosystem. Different models can be created that investigate the effect of the pesticide. Different types of reaction schemes for this system are shown in Figure 3.7, where \(F\) represents the amount of pesticide in the fish, \(W\) the amount of pesticide in the water, and \(S\) the amount of pesticide in the soil. The prime (e.g. \(F'\), \(W'\), and \(S'\) represent other bound forms of the respective state). In all seven different models can be derived. For each of the model schematics, apply the Law of Mass Action to write down a system of differential equations.
References
.jpg){kind=link}
.jpg){kind=link}
There is a lot more nuance for reasons behind periodic fluctuations in these two populations, which includes more complicated food web interactions and climate variation. MacLulich (1937), Stenseth et al. (1997), and King and Schaffer (2001) are good places to dig into the complexity of this fascinating biological system. Image sources for Figure 3.3: Kilby (2012) and Frank, Jacob W. (2021). Image source for Figure 3.4: OpenStax (2016)↩︎