Introduction

The interaction between mathematics and theorectical biology has been long and rich.  Some of the earliest differential equations involved descriptions of how populations grow and change in time.  Perhaps more interestingly, theoretical biology has challenged mathematicians with novel nonlinear equations and nonlinear phenomena.  Even now, some of the most interesting nonlinear equations which mathematicians study are derived from models of population biology.

In this chapter, you will be introduced to a range of population models, from a simple exponential model of bacteria growth to a complicated nonlinear set of equations describing the population in a chemostat, or biological reactor.  Early on we will discover how much simpler it is to use symbolic aids to take care of mechanical manipulations.

Malthusian Models - Exponential Growth

In 1798 Malthus applied a simple differential model to the growth of the human species and predicted that, barring natural disaster and disease, the human population would grow exponentially, use all available resources, and then die of starvation.  This model failed to include limitations due to disease and other population density dependent processes.  To begin the modeling work, let P(t) denote the number of individuals in a population at any given time, t.  We will assume that the model measures the population continuously over time.  We must realize that this model is at best an approximation of the real population.  In reality, births and deaths occur discretely in time.  However, if we model a large population over a long period of time, the process is continuous.

The Malthusiam model states that a population grows in proportion to the number of individuals in the population.  Mathematically that is,

(dP/dt)=k P(t)

where the proportionality constant, k, is called the growth rate constant or net reproduction rate.  To complete this first order ODE model we need initial conditions.  Thus, let t = 0 be a reference time at which we know that the population has (Subscript[P, 0]) individuals.  So the full model of population growths (as of now) is

(dP/dt)=k P(t),  P(0)=(Subscript[P, 0]).

This is one of the first types of problems one encounters in a basic course in ODE's.  One method to solve such a problem is the method of separation of variables.  Using the algebra of differentials the ODE becomes

dP=k P dt

Next, divide both sides by P to get

(dP/P)=k dt

Integrate both side

ln(P)=k t +(Subscript[C, 1]).

Exponentiating both sides gives

P(t)=(e^(k t+(Subscript[C, 1])))=(e^(Subscript[C, 1]))(e^(k t))=C (e^(k t)).

This is the general solution (or family of solutions) to the ODE.  In order to find the particular solution that corresponds to the initial conditions, we substitute the initial conditions into the general solution and solve for C.

P(0)=(Subscript[P, 0])⟹C (e^(k·0))=(Subscript[P, 0])

So C=(Subscript[P, 0]), and the solution to the initial value problem is

P(t)=(Subscript[P, 0])(e^kt).

In Mathematica we would arrive at this solution by doing the following:

sols=DSolve[{P'[t]==k*P[t],P[0]==P0},P[t],t]

We could graph this solution (after plugging in values for the parameter k and the initial condition P0) by issueing the following command:

Plot[P[t]/.sols/.{k->0.1,P0->1.0},{t,0,20},PlotRange->{0,8},AxesLabel->{"t","P(t)"}]

If we want to see how the growth rate effects the solution, we could graph the solution for a variety of values of k.  In what follows we use the Table command to generate several solutions in the same graph.

Plot[P[t]/.sols/.Table[{k->0.1*(2*j-1),P0->1.0},{j,6}],{t,0,20},PlotRange->{0,20},AxesLabel->{"t","P(t)"}]

Here we used k values of .1, .3, .5, .7, .9, and 1.2.  As k increases, so does the rate at which the population grows (hence we call k the growth rate parameter).

The Logistic Model

Clearly it is not realistic to expect ongoing exponential growth.  As a population grows, it consumes resources.  That in turn influences how rapidly it can continue to grow.  One model, first proposed by Verhulst in 1838, that reflects this consumption of resources is the logistic model.  Again, let P(t) denote the number of individuals in some population and let C(t) denote the amount of some limiting resource, such as arable land in the case of human populations.  As the population increases the amount of land needed to support the population must also increase.  So, if there is a limited amount of land then the population will also be limited.

Mathematically we must write equations that relate the rate of change in population and the use of the limiting resource.  To do this we assume that the growth rate of the population varies jointly with the current population and the unused amount of the limiting resource available.  That is

(dP/dt)=a C P.

The other relationship is that of how the unused amount of the limiting resource changes when the population changes.  In this case as the population increases the amount of unused resources will decrease.  If we assume that these rates are proportional then

(dC/dt)=-b((dP )/dt),

with b > 0.  The initial conditions that go along with this system of ODE's would be

P(0)=(Subscript[P, 0]), and C(0)=(Subscript[C, 0]).

In this case the system of ODE's can easily be reduced to a system of one equation.  We do this by integrating the second ODE to get

C(t)=-b P(t)+(Subscript[C, 1]).

Plugging in the initial conditions for C gives us

(Subscript[C, 0])=-b (Subscript[P, 0])+(Subscript[C, 1]).

This means that (Subscript[C, 1])=(Subscript[C, 0])+b (Subscript[P, 0]) and

C(t)=(Subscript[C, 0])+b (Subscript[P, 0])-b P(t).

Next substitute the expression for C into the first ODE to reduce the problem to

(dP/dt)=a P ((Subscript[C, 1])-b P), P(0)=(Subscript[P, 0]).

Note that we keep the (Subscript[C, 1]) to shorten the notation, but that (Subscript[C, 1]) is a positive constant.  Let's see if Mathematica will solve this initial value problem for us:

Clear[P,P0]
Logistic=DSolve[{P'[t]==a*P[t]*(C1-b*P[t]),P[0]==P0},P[t],t]

This looks somewhat complicated, so we will assume that the initial conditions are both unitary (that is, 1 unit of population and 1unit of available resources).

C1=1;
P0=1;

Then the solution to the Logistic Model looks like

Logistic

There are still two parameters that effect the shape of the solution.  To see how the parameter a effects the solution we can look at several graphs holding b fixed (at say 0.1) and varying a.

p1=Plot[P[t]/.Logistic/.{a->0.1,b->0.1},{t,0,50},PlotStyle->{Red}];
p2=Plot[P[t]/.Logistic/.{a->0.2,b->0.1},{t,0,50},PlotStyle->{Green}];
p3=Plot[P[t]/.Logistic/.{a->0.3,b->0.1},{t,0,50},PlotStyle->{Blue}];
Show[p1,p2,p3,PlotRange->{0,10}]

Note that in all of the above cases the solution approaches the value 10 as time increases.  The parameter a effects how fast the solution approaches it's limiting value.  To see how b effects the solution, graph it holding a fixed at .2 and vary b.

p4=Plot[P[t]/.Logistic/.{a->0.2,b->0.05},{t,0,50},PlotStyle->{Red}];
p5=Plot[P[t]/.Logistic/.{a->0.2,b->0.2},{t,0,50},PlotStyle->{Black}];
p6=Plot[P[t]/.Logistic/.{a->0.2,b->0.3},{t,0,50},PlotStyle->{Blue}];
Show[p2,p4,p5,p6,PlotRange->{0,20}]

In this case the limiting value of the solution is effected by changing b.  Now lets look at changing the initial population.  Let's use a = .2 and b = .1, and graph initial populations from 1 to 15.

Clear[P0]
Plot[P[t]/.Logistic/.Table[{a->.2,b->.1,P0->j},{j,15}],{t,0,30},PlotRange->{0,16}]

It appears that starting with a positive initial population leads to a population of 10 (for the values of the other parameters we have used).  Let's see why.

Equilibrium Solutions

An equilibrium solution to a differential equation is a constant that solves the equation.  If (Subscript[P, eq]) is an equilibrium solution then for this solution we have

((Subscript[dP, eq])/dt)=0

That means the right hand side of the differential equation is zero at an equilibrium solution.  That is

a (Subscript[P, eq])((Subscript[C, 1])-b (Subscript[P, eq]))=0,

where (Subscript[P, eq]) is the equilibrium solution.  In this case we have two solutions to the quadratic equation in (Subscript[P, eq]), namely

(Subscript[P, eq])=0, or (Subscript[P, eq])=((Subscript[C, 1])/b).

The first equilibrium solution is the trivial solution.  If you start with zero population in a closed system (no migration), then you always have zero population.  The second equilibrium solution is called the carrying capacity for the system.  From the graphics we had above the behavior of these two equilibrium solution is quite different.

In particular, the stability of these two solutions is different.  Stability in this case means: if we start with initial conditions "close" to the equilibrium solution then does the solution stay close to the equilibrium solution?  For the trivial solution the answer to this question is no, so we say that (Subscript[P, eq])=0 is an unstable equilibrium solution.  For the other equilibrium solution, small changes in initial condition do not change the solution much, so this equilibrium solution is called stable.  Mathematically we "prove" this as follows.

Consider the differential equation

(dP/dt)=f(P),

where f is a differentiable function of P and that f((Subscript[P, eq]))=0 by our definition of equilibrium.  To analyze the stability of an equilibrium solution we add a "small" perturbation to the equilibrium solution and see what happens to the solution of the differential equation.  Formally, we start by setting

P(t)=(Subscript[P, eq])+ϵ p(t),

and take ϵ to be small: 0 < ϵ ≪ 1.  Next substitute the perturbed function into the differential equation.  This will give an equation of the form

(d/dt)((Subscript[P, eq])+ϵ p(t))=f((Subscript[P, eq])+ϵ p(t)).

The linearity of the derivative and (Subscript[P, eq]) a constant implies that

ϵ(d/dt)(p(t))=f((Subscript[P, eq])+ϵ p(t)).

Expand the right-hand-side in a Taylor series centered at (Subscript[P, eq]) to get

ϵ(dp/dt)=f((Subscript[P, eq]))+f'((Subscript[P, eq])) ϵ p+(1/2)f''((Subscript[P, eq])) (((ϵ p))^2)+⋯

The first term on the right is zero.  If we divide by ϵ we get

(dp/dt)=f'((Subscript[P, eq])) p+𝔒(ϵ).

The "big-O" notation indicates that the remaining terms are of order ϵ.  If we think of ϵ as a small number then we have approximately

(dp/dt)=α p, where α=f'((Subscript[P, eq])).

this has solution p(t)=p(0) (e^(αt)) and there are three possible types of behaior that can occur.

1. For α < 0, the solution, p(t), will decay to zero.

2. For α = 0, the solution remains constant, p(t) = p(0).

3. For α > 0, the solution will grow rapidly away from zero.

These behaviors tell us that the perturbed solution will in (1) decay to the equilibrium solution (the equilibrium is stable), (2) stay approximately the same distance from the equilibrium (could be stable or unstable, but we can not tell by this analysis), or in (3) move away from the equilibrium solution (unstable).  Thus, to check the stability of an equilibrium we must look at the derivative of the right hand side evaluated at the equilibrium.  Let's check the Logistic model.

First define the right hand side of the Logistic model as a function of P

Clear[f,a,b,P,C1]
f[P_]=a*P*(C1-b*P)

Then the derivative of the right-hand-side function is

df[P_]=f'[P]

Evaluating this at the zero equilibrium gives us

df[0]

Considering that all paramenters in the model are positive, this indicates that the zero equilibrium is unstable.  Evaluating the positive equilibrium gives us

df[C1/b]

Again, because all parameters in the model are positive, we have that the linear stability analysis gives that the positive equilibrium is stable.

Linear Stability in a General Example

Consider the following example

(dx/dt)=(((x^2)(1-x))^3).

If we wish to explore the equilibrium solutions using linear stability analysis, we would do the following:

Clear[f]
f[x_]=x^2*(1-x)^3

Let's look at a graph of this function to get an idea of the equilibrium solutions for the ODE.

Plot[f[x],{x,-1,2},PlotRange->{-0.2,0.2},AxesLabel->{"x","f"}]

Based on the graph, which can easily be confirmed by some simple algebra, there are two equilibrium solutions: namely x = 0 and x = 1.  If we try the linear stability analysis on these solutions we get

df[x_]=f'[x]

df[0]

df[1]

In both cases the linear stability analysis fails to provide any information about the stability of the steady-state (equilibrium) solutions.  However, we can get information about stability from the graph of f.  Namely, when x < 0, then f(x) > 0 so x(t) is increasing.  When 0 < x < 1, then f(x) is still possitive so x(t) is increasing.  Finally for x > 1, then f(x) is negative so x(t) is decreasing.

Recall from the theory of differential equations that solutions through a point are unique (provided the right-hand-side is "nice").  Another way of saying that is that two different solutions can not intersect.  That is, in a state space graph, if a solution starts below an equilibrium solution, then it must stay below the wquilibrium solution.  If a solutions starts above an equilibrium solution, then it must stay above the equilibrium solution.

So, in this example, if we have IC's x(0) < 0, then the solution increases towards 0, but can not cross the t-axis.  If we have IC's between 0 and 1, then the solution will increase toward 1, but can not cross x = 1.  If we start with IC's greater than 1 then the solution will decrease towards 1, but will not cross x = 1.  Put that all together to conclude that x = 1 is a stable steady-state solution and x = 0 is a semi-stable steady-state solution (unstable).

Next, graph several solutions to verify our analysis.

sols=Table[NDSolve[{x'[t]==f[x[t]],x[0]==0.2*k},x[t],{t,0,20}],{k,-10,15}];

Plot[x[t]/.sols,{t,0,20},PlotRange->{-2,3},AxesLabel->{"t","x"}]

In conclusion, if solutions are increasing for values below the steady-state and decreasing for values above the steady-state, then the steady-state is stable.  Otherwise, the steady-state is not stable, but could be semi-stable (or stable on one side).

Exercise 5

1. For the general system (dx/dt)=f(x) we have learned to start a qualitative analysis by determining the locations of the equilibrium solutions and the stability of the equilibrium solutions.  For each of the cases below, determine the equilibrium solutions and their stability characteristics.  In some cases, there will be infinitely many equilibrium.  In these cases, characterize the first few equilibrium solutions closest to the origin, x = 0.  Some Hints:

1. It may not be algebraically possible to compute solutions.  So, you might have to consider settling for approximations.  In these cases consider using Solve and/or NSolve commands to find the equilibrium.

2. To analyze the stability you can compute the derivative of the right-hand-side function using appropriate Mathematica commands.

3. Recall that if the derivative is zero at an equilibrium solution we cannot tell about the stability of the equilibrium using linear stability analysis.  That is, we would say the stability is still indeterminate.  See problem 2 and 3.

a.

(dx/dt)=(x^3)(x-1)(x+1).

b.

(dx/dt)=sin(x).

c.

(dx/dt)=(e^x)-cos(x).

d.

(dx/dt)=tan(x)-x.

e.

(dx/dt)=x tan(x).

2. For this problem we will consider the differential equation (dx/dt)=(x^2)(1-(x^2)).  Perform the following steps toward analyzing the stability of equilibrium solutions.

a. Graph the right-hand-side function as a means of visually identifying the locations of the equilibrium solutions (i.e. right-hand-side equals zero).

b. Use the derivative test to ascertain the stability of the two non-trivial equilibrium solutions.  In addition, show that the test does not tell us anything aout the stability of the equilibrium solution x = 0.

c. Now use the graph of the function to determine the stability of the solution x = 0.  To understand this, consider what happens when (dx/dt)>0, and (dx/dt)<0.  In the first case x must be increasing and in the second case, x must be decreasing.  If x is increasing to the left of an equilibrium solution and x is decreasing to the right of the same equilibrium solution, then the equilibrium solution must be stable.  If the opposite occurs, the equilibrium must be unstable.

3. Repeat problem 2 for the system defined in problem 1a.