How Mathematics Describes Life
The circle of life is something we have all heard of from somewhere, but we don't usually try to calculate it. For some time we have been working on analyzing a predator-prey model to better understand how mathematics can describe life, in particular the interaction between two different species (a predator and a prey). The model we are analyzing is called the Holling-Tanner model, and it can't be solved analytically. The Holling-Tanner model is a very common model in population dynamics because it is a simple descriptor of how predators and prey interact. The model is a system of differential equations that is controlled by constants of how the species interact. This model is a generalized math concept, so it can describe any predator prey interaction given the right constants; from lions and zebras to white blood cells and diseases.The model is a system of two differential equations with solution properties described by three parameters which are characteristics of the interacting species and their environment. The model is not specific to any particular set of species and so can describe predator-prey species ranging from lions and zebras to white blood cells and infections.
One thing all these systems have in common is they all have critical points. A critical point is a value for both populations that keeps both populations constant. It is a pair of equilibrium populations for the predator and the prey; this is important because for such a point the right hand sides for both differential equations are equal to zero. For this model there are two critical points, a predator free critical point and a coexistence critical point. Most of the analysis we did is on the coexistence critical point because the predator free critical point is always unstable and frankly less interesting than the coexistence critical point.
What we did is consider two regimes for the differential equations, large $B$ and small $B$. $B$, $A$, and $C$ are constants in the differential equations that control the nature of the solution where $B$ measures the responsiveness of the predators population to change; as $B$ gets larger the predators become more responsive. $A$ represents predation of the prey, and $C$ represents the satiation point of the prey population. For the large $B$ case we were able to approximate the system of differential equations by a single scalar equation. For the small $B$ case we were able to predict the limit cycle, as well as look at the type of solution all regimes of $B$ give us. The limit cycle is what I referred to earlier as the circle of life; it is a process of the predator and prey populations growing and shrinking periodically. This model has a limit cycle in the regime of small $B$, that we solved for numerically. With some assumptions to reduce the differential equations, we were able to create a system of equations and unknowns to predict the behavior of the limit cycle for small
Why
One fundamental problem in population dynamics is the predator-prey
problem. In its simplest form there are only two species - the prey
and the predator. The predators survive by devouring the prey. When
the prey population is large, the predators will increase because
there is more for them to eat. Conversely, when the prey population is
small the predator population will decrease because they do not have
enough to eat. In many instances the decrease of the predator population
causes the prey to recover and then subsequently the predator population
will increase again. The result is a cyclic behavior of the predators
and the prey with the predators lagging the prey. While this problem
has been stated in terms of population dynamics, similar behavior
can occur at the cellular level, for example white blood cells acting
as predators and some external intrusion, i.e., an infection, acting
as the prey.
Let $P(t)$ and $V(t)$ represent the populations of the prey and predators,
respectively, as functions of time $t$. Many different predator-prey
models have been proposed and studied, but one very important model
is called the Holling-Tanner model. In its most general form this model
is a coupled system of two differential equations,
$$
\frac{dP}{dt}=P(\alpha-\beta P-A\frac{V}{P+C}) \,, \quad \frac{dV}{dt}=BV(1-\gamma \frac{V}{P}) \, ,
$$
where $\alpha$, $\beta$, $\gamma$, $A$, $B$ and $C$ are
positive constants. This model can be simplified by appropriate
nondimensionalization of $P$, $V$ and $t$ resulting in a set of equations
where $\alpha$, $\beta$ and $\gamma$ are set to 1,
$$
\frac{dP}{dt}=P(1- P-A\frac{V}{P+C}) \,, \quad \frac{dV}{dt}=Bv(1- \frac{V}{P}) \, .
$$
In order to understand what these equations mean, consider first the $P$
equation.
The term in parenthesis is the net birth rate per unit population.
For small values of $P$ and $V$ this term is just one and is just the
natural birth rate of the species. The next term $-P$ represents
a crowding effect. For any species, if the population is too large
then intraspecies competition for a scarce resource decreases
the birth rate (hence the negative sign). The last term is predation
(again with a minus sign) and represents the decrease in the prey
population when they interact with predators. The denominator
in this term represents a satiation effect and simply says that
there are limits to how much the predators can consume, i.e., for large
$P$ predation (per unit population) decreases.
Next consider the $V$ equation. Again, the term in parenthesis
is the net birth rate per unit population for the predators and
the equation simply states that when the predator population is greater
than the prey population $V>P$ there is not enough prey for the predators
to consume and their population decreases (term in parenthesis is negative), whereas if $V< P$ the predator population will grow (term in parenthesis is positive). This system is a highly nonlinear system of differential equations and in general there is no analytic solution - it can only be solved numerically. However, it can be analyzed mathematically by looking for critical (equilibrium) populations where the right hand side of each equation is zero. If you find such a point then the populations are in equilibrium, since their time derivative is zero.
Note that the only solutions of interest are when $P(V)\ge 0 )$ since $P$ and $V$ represent populations. Negative values of $P$ or $V$ describe non-physical solutions. It is easy to see that one such critical point is $(P,V)=(1,0)$. This is called the predator-free state since $V=0$. Simple math analysis shows that the predator-free state is always unstable, small perturbations will lead to a solution which moves away from the predator-free state ($V$ grows). A second critical point $(P_{*},V_{*})$ with $V_{*}=P_{*}>0$
also exists. The right
hand side of the second equation is zero when $V=P$ so that we must
have $V_{*}=P_{*}$. We can then solve for $P_{*}$ from a simple
algebraic manipulation of the right hand side of the $P$ equation.
This critical point is called the coexistence critical point
since the two species exist together in equilibrium. For a large
range of parameter values the coexistence critical point is
unstable, and the solution will evolve to a periodic solution
$P(t),V(t)$ representing cyclic growth and decline of the
two species. This solution is called a limit cycle because
it is the limit (as $t\rightarrow \infty$) of general initial conditions
and because it is periodic (cyclic).
The limit cycle can be studied in two ways, (i) compute and plot
$P(t)$ and $V(t)$ and (ii) consider the curve in the $P-V$ plane
which this periodic solution traces out (by eliminating $t$).
This last approach is called a phase plane analysis and is useful
because you can directly see the relationship between $P$ and $V$.
In general, the limit cycle and solutions in general can only be
attained through numerical computation. However, in this project
I show how you can mathematically analyze the system in the two
regimes (i) $B$ small and (ii) $B$ large. In the first case I show
how you can get a simple analytic approximation to the limit cycle.
In the second case I show how you can reduce the system of equations
from two equations to just one scalar
How
To look at the stability of the critical points of these linear, autonomous differential equations we used two methods. The first is to solve for the direction field numerically. The second is to make a Jacobian matrix and find the eigenvalues. The equation that describes them is the following: $$ \frac{dP}{dt} = P\left(1-P-\frac{AV}{P+C}\right) $$ $$ \frac{dV}{dt} = BV(1-\frac{V}{P}) $$ These two equation describe the growth rate of the two populations with respect to time and each other. When these two equations equal zero is when they are in equilibrium (have constant population); that is what the critical point describes. There are two critical points, when $(V,P) = (0,1)$ and $(V,P) = (P_c,P_c)$. The $P_c$ we are looking at is positive and real; the value describes what population the two species can coexist. To calculate the critical point we can do the following. $$\frac{dP}{dt} = P\left(1-P-\frac{AV}{P+C}\right) = 0~~~~~~~ \frac{P(1-P)(P+C)-AP^2}{P+C} = 0$$ We set $V = P$ because that is where the critical point is, and now we can simply solve for the numerator using the quadratic formula; $\gamma = A + C$. $$P_c = \frac{1-\gamma}{2} + \frac{\sqrt{(1-\gamma)^2 + 4C}}{2}$$ To confirm what kind of critical points we found we checked the Jacobian matrix of the differential equations and then found their eigenvalues. Depending on the eigenvalues we can determine the type of critical point. The Jacobian matrix of two differential equations can be described as the following: $$\frac{dP}{dt} = P\left(1-P-\frac{AV}{P+C}\right) = f$$ $$\frac{dV}{dt} = BV(1-\frac{V}{P}) = g$$ $$J = \begin{bmatrix} \frac{\partial f}{\partial P} & \frac{\partial f}{\partial V}\\ \frac{\partial g}{\partial P} & \frac{\partial g}{\partial V} \end{bmatrix} $$ We can then use this information to categorize the solutions for all regimes of all the constants. To approximate these differential equations into a single scalar equation we started with the rate of change of the predators. $$ \frac{dV}{dt} = BV(1-\frac{V}{P}) $$ If we assume that $B$ is large then we can bring it to the other side of the equation, and we can assume at large enough $B$: $$ \frac{1}{B} \frac{dV}{dt}=V(1-\frac{V}{P})\approx 0 $$ Since we are not investigating the predator free state we can remove that as a solution. $$(1-\frac{V}{P})\approx 0$$ What is left is the assumtion $V\approx P$ and we can then have a single scalar equation that describes the rate of change of the two populations. $$\frac{dP}{dt} = P\left(1-P-\frac{AP}{P+C}\right) $$
To understand how we predicted the shape of the limit cycle for small $B$, the first step is to reformulate the system so that the small number $B$ is in front of one of the derivatives. To do this make the change of variables $\tau=Bt $, so that we now get the system of equations $$ B\frac{dP}{d \tau}=P(1-P-A\frac{V}{P+C})\, , \quad \frac{dV}{d \tau}=V(1-\frac{V}{P}) \, . $$ Now consider the direction field for this system.
The plot above is for $B = 10^{-5}$. In the plot above the halo red dot represents the critical point, and the solid red dot represents where the numerical solution stopped solving for the limit cycle.
Since for $P$ we now
have $B$ in the denominator
$(\dfrac{dP}{d \tau}=\dfrac{P}{B}(1-P-A\dfrac{V}{P+C})\, )$,
we expect the direction field to be nearly horizontal ($dP$ predominates)
unless either $P$ or $(1-P-A\dfrac{V}{P+C})$ is close to zero.
The equation $1-P-A\dfrac{V}{P+C}=0$ describes a
downward facing parabola. Thus, portions of
the limit cycle where $V$ really changes (not horizontal in the $P-V$ plane)
have to be either near this parabola or have $P$ small.
Consider first the parabola. Along the right branch of the parabola
$dV$ is positive so that the limit cycle will include a curve very
close (not exact of course) to the right branch of the parabola.
When the maximum of the parabola is reached (call it $V_{max}$) the
limit cycle can no longer follow the parabola since the
parabola turns down
but the direction field still points upward ($dV$ still positive).
The only possibility is that the limit cycle is now nearly horizontal until
it reaches a small value of $P$ ($P_{turn}$)
when it will turn down. We approximate
this turning point by looking at a point where the direction field
points in the $45^{\circ}$ direction and $V=V_{max}$
($dP=dV$ and $V=V_{max}$ - one equation and one unknown for $P_{turn}$).
Next consider how $V$ changes when $P$ is small.
Since $dP$ will now be small, we expect this part
of the limit cycle to be nearly vertical ($dV$ predominates).
For small $P$ we can
get an approximate equation for $\dfrac{dP}{dV}$ which we can solve.
$$
\frac{dP}{dV} = -\frac{P(1-\frac{A}{C}V)}{B\frac{V^2}{P}}
$$
To solve we will seperate the variables and integrate and set $\frac{dP}{dV} = -1$. This will give us two equations for $P$ as a function of $V$ which will allow us to solve for $V_{small}$. This will also allow us to find $P_{turn 2}$. The limit cycle
will follow this nearly vertical line until it turns back again at
$P_{turn2}$, and $V_{small}$. When it turns again it will be nearly
horizontal until it intersects the parabola
at $P=P_{intersect}$. We can then calculate $P_{intersect}$ by setting the equation for the parabola equal to $V_{min}$. These steps allow us to construct a good approximation to the
limit cycle for small
Results
To showcase the results we will display many plots. The first of which is a region plot that characterizes the solutions to the Holling-Tanner model.
.png)
.png)
It is well known that a critical point can be classified
into several different possibilities. These possibilities
are called nodes, spirals and saddles. Nodes and
spirals can also be either stable or unstable.
Limit cycles occur when the coexistence critical point
is an unstable spiral. Since there are three parameters
(A, B, and C) in order to visualize regions where
the critical point has different behavior we reduced the
visualization to two dimensions. We did this by freezing
the value of the critical point (for the figures above the
critical point was (0.1,0.1)). The value of the critical
point does not depend on B, only on A. For the
left most figure we considered different values
of A and B, computing what C would have to be in order
to keep the critical point fixed. For the right most figure
we did the same but giving C and computing A.
If this resulted in negative parameters we did not try
to classify the critical point, rather we just considered
the case unphysical. The region of the unstable spiral
is where cyclic behavior occurs and there is where limit cycles occur.
In the plot below$ P(t)$ is the population of prey, $V(t)$ is the population
of predators and in blue we plot $P$ obtained from the
single scalar equation that I derived for the large $B$ case.
The first plot is for B=2 and the solution for the system almost
immediately latches on to the solution to the single
equation. The second plot is for small B and while
asymptotically the solution to the system approaches that
of the scalar equation. The lower plot is the same curve except for $B=0.001$. This displays how the solution degradates as $B$ gets smaller. We infact took an interest in how much the solution degredates as a function $B$, so we plotted the relative error of the solutions as a function of $B$.
For the small $B$ regime we were able to analytically predict the limit cycle without solving the system of differential equations. The plot below shows the two solutions overlayed on top of eachother.
The plot above shows the two solutions in the $V-P$ plane to better illistrate the limit cycle. This result is very exiciting because it gives us an analytical method to approximate the limit cycle. The plot above has the critical point $(0.1,0.1)$, and $B = 10^{-5}$. Notice that the approximation is not perfect; at the top and bottom horizontal regions the approximation is slightly off of the numerical solution. The more we raise the value of $B$ the larger these errors become.
The plot above has the parameter $B = 10^{-3}$. The error might become larger, but is still a reasonable approximation. This method can be used to better understand the system of equations, and the limit cycle. Hopefully this might even inspire new ideas for other models.
Acknowledgements
This approach was suggested by Applied Mathematics Professors Vladimir A. Volpert and Alvin Bayliss. I also want to acknowledge help from Applied Math graduate student Eric A. Autry.
