Introduction

An interesting population system to analyze is bacterial growth in a chemostat, or biological reactor.  This is a system used in many labs, where it is necessary to maintain a fresh culture of a desired bacteria.  The idea is to supply a culture, or reaction tank, with a constant supply of fresh nutrient through an inlet valve.  The  reactor is stirred constantly, and mixed effluent is tapped off at the same bulk rate as fluid is supplied from the nutrient tank.  A mathematical model for this system was proposed by Segel (1984).

Let N(t) be the relative population density of the bacteria in the reactor, and C(t) be the concentration of nutrient in the reactor.  The relative population density is a measure of the number of bacteria per unit volume in the chemostat.  From a practical point of view it would be impossible to count the total number of bacteria in the tank.  Instead, we will estimate the density.  Then an estimate of the total number of bacteria can be computed based on the volume of the tank.  The concentration of nutrient is relative to the influent nutrient concentration.  The influent concentration will be normalized to be one.  For example, if the concentration of nutrient in the chemostat is 0.98, this implies that the concentration level is 98% of the concentration of nutrient being pumped into the system.

One example of a set of modeling equations describing the temporal dynamics of the process is

(dn/dt)=f(n,c)=(Subscript[α, 2]) (c/(1+c))n-n,

(dc/dt)=g(n,c)=-n (c/(1+c))-c+(Subscript[α, 1]).

The negative linear term on the right-hand-side of the first equation represents the loss of bacteria due to material being drained out of the chemostat.  The negative terms on the right-hand-side of the second equation characterize the reduction in concentration of the nutrient due to consumption by the bacteria and due to material being drained off the tank.  The positive term in the first equation reflects the increase in the bacteria due to the nutrient being injected into the tank.  The positive constant, (Subscript[α, 2]), is the growth rate constant for the bacteria.  The positive constant, (Subscript[α, 1]), is the amount of nutrient being added to the tank.  We should also note that both N and C are reserved Mathematica terms so we will use lower case for the dependent variables.

What do we want to determine about the physical system from the mathematical model we are developing?

1. Is there a way to maintain a certain density of bacteria that is useful for lab experiments?

2. What parameter values will maintain a specified density of bacteria in the chemostat?

3. In these types of equilibrium states, can we guarantee stability of the process?  That is, if the system gets bumped or the concentration of the influent nutrient is perturbed slightly, will the equilibrium state be maintained?

The computations to answer these questions can be rather complicated to do by hand.  Conceptually, we already have some idea of what we need to do to analyze the model.  For equilibrium states we need the derivatives to be zero and for stability we need to look at the solution of the linearized, perturbed problem.  We can use Mathematica to do the calculations:

f[n_,c_]=(Subscript[α, 2])*c/(1+c)*n-n
g[n_,c_]=-n*c/(1+c)-c+(Subscript[α, 1])

Note that f and g are functions of the two variables n and c.  Next we want to find equilibrium solutions by setting f and g to zero and solve for n and c.

eqsols=Together[Solve[{f[n,c]==0,g[n,c]==0},{n,c}]]

The first steady state is the trivial case where there are no bacteria and we are just flushing nutrient through the system.  If this steady state solution is stable then we will have a hard time raising bacteria.  If it is unstable and the non-trivial steady state is stable then we can introduce (or start with) a small sample of the bacteria and over time the concentration of bacteria will grow to a much larger amount.  This second situation is the desired outcome.

There is more information that can be gleaned from the above results on steady states.  In order for the non-trivial steady state to be positive (i.e. make sense), we must have

(1/((Subscript[α, 2])-1))>0 ⟹ (Subscript[α, 2])>1.

From the n equation also note that

-(Subscript[α, 2])-(Subscript[α, 1]) (Subscript[α, 2])+(Subscript[α, 1]) (Subscript[α, 2]^2)>0 ⟹ (Subscript[α, 1])>(1/((Subscript[α, 2])-1))

to insure that n > 0. Now we are faced with a more difficult question:  is this a stable steady state?  To answer this we liniarize the general solution about the steady state and examine solutions to the linear equations as follows.

Stability for the System of Ordinary Differential Equations

First some notation.  Let u be a vector whose components are n and c.  That is

u=(({{n}, {c}}))

Let F(u) be the vector valued function whose components are the right-hand-side of the system of differential equations.  That is

F(u)=(({{f(n,c)}, {g(n,c)}})).

Then the system of differential equations can be represented in matrix form as

(d/dt)(({{n}, {c}}))=(({{f(n,c)}, {g(n,c)}})),

or the general system is

(du/dt)=F(u).

Now suppose that we have found an equilibrium solution, (Subscript[u, 0]), of the general system.  This means that

(d/dt)(Subscript[u, 0])=0=F((Subscript[u, 0])),

where 0 is the zero vector.  Now if we introduce a perturbation to the equilibrium solution of the form ϵ v, where v is a vector function and ϵ≪ 1, then the perturbed vector would be

u=(Subscript[u, 0])+ϵ v.

Substituting into the the left-hand-side of the general system gives

(d/dt)u=(d/dt)((Subscript[u, 0])+ϵ v)=ϵ (dv/dt),

since (d/dt)(Subscript[u, 0])=0.

Next we need to deal with the right-hand-side of the general solution.  In the one dimensional analysis of the Logistic equation we used Taylor series expansion in the parameter  about the equilibrium solution.  An analogous procedure will be used in the case of systems of ODE's.  In particular, we let the parameter  vary and we expand F((Subscript[u, 0])+ϵ v) about the point (Subscript[u, 0]).  The formal Taylor series can be found in a calculus book as

F((Subscript[u, 0])+ϵ v)=F((Subscript[u, 0]))+ϵ (({{(df/dn), (df/dc)}, {(dg/dn), (dg/dc)}})) v+𝔒((ϵ^2)).

Again, the term 𝔒((ϵ^2)) represents all terms in the Taylor series that are multiplied by at least a factor of (ϵ^2).  Note that the matrix

J =(({{(df/dn), (df/dc)}, {(dg/dn), (dg/dc)}}))

is called the Jacobian matrix for the funtion F(u), and is the analogy of the derivative in the usual one dimensional Taylor series.  At this point we will not focus on the derivative expansion above.  We will assume that it is well defined and proceed to the next step.  The entries in the matrix are all evaluated at the equilibrium solution (Subscript[u, 0]).

Using a shorthand notation for the Taylor expansion in terms of the Jacobian we can now write the general equation for the perturbed system as follows:

ϵ(dv/dt)=F((Subscript[u, 0]))+ϵ J((Subscript[u, 0])) v+𝔒((ϵ^2)).

Barring notation this is virtually the same as the result we got in the one dimensional analysis.  If we note that F((Subscript[u, 0]))=0, divide through by ϵ and then drop any remaining term in ϵ (since ϵ ≪ 1), then we get,

(dv/dt)=J((Subscript[u, 0])) v.

Since the Jacobian has fixed entries, this is a constant coefficient, linear system of ODE's.  Before we explore the issue of stability, let's take a look at how one goes about solving (in closed form) a system like this.

Solving Constant Coefficient Linear Systems of ODE's

For this analysis we will restrict ourselves to 2 by 2 systems.  Suppose that the 2 by 2 matrix J has two eigenvalues, (Subscript[λ, 1]), (Subscript[λ, 2]), and two corresponding eigenvectors, (Subscript[μ, 1]), (Subscript[μ, 2]).  If the iegenvectors are distinct, then they can be used as a basis and any two dimentional vector can be written as a linear combination of the eigenvectors.  Thus, we can write

v(t)=(Subscript[v, 1])(t) (Subscript[μ, 1])+(Subscript[v, 2])(t) (Subscript[μ, 2]),

for any perturbation vector v.  In terms of the linearized system

(dv/dt)=J v

we have

((Subscript[dv, 1])/dt) (Subscript[μ, 1])+((Subscript[dv, 2])/dt) (Subscript[μ, 2])=J((Subscript[v, 1])(Subscript[μ, 1])+(Subscript[v, 2])(Subscript[μ, 2]))=(Subscript[λ, 1])(Subscript[v, 1])(Subscript[μ, 1])+(Subscript[λ, 2])(Subscript[v, 2])(Subscript[μ, 2]),

where we have used results from linear algebra about eigenvalues and eigenvectors.  Since the eigenvectors are linearly independent, we can "decouple" the equations along the directions of the eigenvectors.  That is, we get two, one dimensional linear ODE's

((Subscript[dv, 1])/dt)=(Subscript[λ, 1])(Subscript[v, 1]),

((Subscript[dv, 2])/dt)=(Subscript[λ, 2])(Subscript[v, 2]).

We can solve each of these equations separately and get

(Subscript[v, 1])=(Subscript[c, 1])(e^((Subscript[λ, 1])t)) and (Subscript[v, 2])=(Subscript[c, 2])(e^((Subscript[λ, 2])t)),

where (Subscript[c, 1]) and (Subscript[c, 2]) are arbitrary constants.  Then the solution to the linearized system is

v(t)=(Subscript[c, 1])(e^((Subscript[λ, 1])t))(Subscript[μ, 1])+(Subscript[c, 2])(e^((Subscript[λ, 2])t))(Subscript[μ, 2]).

Stability of the System of ODE's

Now that we have the solution to the linearized problem we can investigate the stability of the chemostat problem.  In order to guarantee stability we must have v(t) shrink to the zero vector.  This will only occur if both eigenvalues of J are negative (or at least have negative real part).

At this point we can proceed in one of two ways.  The first is to compute the eigenvalues directly from the matrix J.  The other is to use some properties of matrices and the characteristic polynomials of 2 by 2 matrices to find the relevant information about the eigenvalues.  Consider a general 2 by 2 matrix

A={{a,b},{c,d}};
MatrixForm[A]

The characteristic polynomial for this matrix, p(λ)=Det(A-λ I), can be found in Mathematica by using the CharacteristicPolynomial command as follows:

p=CharacteristicPolynomial[A,λ]

Note that when simplified, the characteristic polynomial is

p(λ)=(λ^2)-(a+d)λ+(a d-b c).

The eigenvalues are the roots of this polynomial.  If we new the eigenvalues then we could write p(λ) as

p(λ)=(λ-(Subscript[λ, 1]))(λ-(Subscript[λ, 2])).

Multiplying out this expression gives

p(λ)=(λ^2)-((Subscript[λ, 1])+(Subscript[λ, 2]))λ+(Subscript[λ, 1])(Subscript[λ, 2]).

Compare this with the characteristic polynomial above gives two conditions on the eigenvalues:

a+b=(Subscript[λ, 1])+(Subscript[λ, 2]), and  a d-b c=(Subscript[λ, 1])(Subscript[λ, 2]).

Now notice the following Mathematica commands:

Tr[A]

Det[A]

For the 2 by 2 matrix A the Tr and Det commands five us the left-hand-sides of the above expressions.  Thus, both eigenvalues will be negative if the determinate is positive and the trace is negative.

Stability of the Chemostat

In the chemostat problem the jacobiam matrix is

jacm[n_,c_]=D[{f[n,c],g[n,c]},{{n,c}}]

MatrixForm[jacm[n,c]]

The determinate and trace of this matrix are

tracejacm=Simplify[Tr[jacm[n,c]]]

detjacm=Simplify[Det[jacm[n,c]]]

We need to evaluate these expressions at the equilibrium solutions as follows

traceATeq=Simplify[tracejacm/.eqsols]

detATeq=Simplify[detjacm/.eqsols]

Here the first values are for the trivial equilibrium and the second values are for the positive equilibrium.  Recall that (Subscript[α, 2])>1 and (Subscript[α, 1])(Subscript[α, 2])-(Subscript[α, 1])-1>0 in order to have a positive equilibrium.  Thus the determinant of the jacobian is positive.  That means that the eigenvalues of the jacobiam have the same sign, they are either both positive or both negative.  However, the trace of the jacobian, which is the sum of the eigenvalues, is clearly negative given the restrictions on the parameters.  Thus both eigenvalues have negative real parts and the nontrivial equilibrium is stable.

Solutions to the Chemostat Problem

Now let's solve the chemostat problem.  First pick values for (Subscript[α, 1]) and (Subscript[α, 2]).  Ordinarily, these values would depend on the types of bacteria and nutrient used, but we will just assign arbitrary values (keeping in mind the restrictions from above).

(Subscript[α, 1])=3
(Subscript[α, 2])=2

Next let's build a "function" that allows us to solve the system with different initial conditions.

Clear[n,c]
chemo[n0_,c0_]:=NDSolve[
{n'[t]==f[n[t],c[t]],c'[t]==g[n[t],c[t]],
n[0]==n0,c[0]==c0},
{n,c},{t,0,20}
]

Notice the " := ".  This is a delayed assignment.  Mathematica does not attempt to execute the NDSolve command until we call the function chemo.  We do this because the initial conditions are parameters of the function and Mathematica would not be able to perform the NDSolve command until we specify numeric values for the initial conditions.  Lets test this function with small initial conditions:

sols1=chemo[0.1,0.1]

ParametricPlot[Evaluate[{c[t],n[t]}/.sols1],{t,0,20},PlotRange->{{0,4},{0,4}}]

Recall that this is a phase space plot.  For the values of the parameters and initial conditions used we see that first the concentration of the nutrient increases with little change in the concentration of the bacteria.  When the concentration of nutrient reaches 2-2.5, the the concentration of bacteria starts to increase.  Then as the number of bacteria increase, the concentration of nutrient decreases.  The system reaches an equilibrium at about (1,4).

If we wanted to graph the solutions using several initial conditions, we could do the following.

sols=Table[chemo[.5*j,.5*j],{j,1,16}];

Show[
ListPlot[{Table[{.5*j,.5*j},{j,1,16}]},PlotMarkers->"*",PlotStyle->Red],
ParametricPlot[Evaluate[{c[t],n[t]}/.sols],{t,0,20},
PlotStyle->Blue],
ListPlot[{{1,4}},PlotMarkers->"o",PlotStyle->Black],
PlotRange->{{0,8},{0,12}}
]

The nonzero equilibrium for this system is plotted as a black circle at (1,4), the initial conditions are graphed as red stars.

Exercise 6

Beddington and May (1982) proposes a model to describe the interaction between baleen whales and krill.  Krill are the primary food resource for baleen whales which means the whales are the predators in the system and the krill are the prey in the system.  Let's let y(t) denote the population density of baleen whales and x(t) denote the population density of krill.  The model proposed was defined by the following system of ordinary differential equations.

(dx/dt)=r x (1-(x/K))-((a x y)/(x+d)),

(dy/dt)=s y (1-(y/(b x))).

a, b, d, r, s, and K are parameters of the system that depend on the characteristics of the whales and krill.  r and s are reproductions rates.  K is a carrying capacity for the krill.  b is a relative carrying capacity for the whales given x number of krill.  a and d are constants that describe the predation of the whales on the krill.  x and y are densities which might have values of tons per square mile of ocean?

Answer the following questions:

1. Determine all potential steady states (equilibrium solutions) of the system using Mathematica.  Discuss the solutions in terms of the real system.

2. Determine the stability of all steady states for the system.  Be sure to use the Simplify command to reduce your answers.

3. Discuss the qualitative behavior of the system, using phase space plots where appropriate.

Note:  Parts of this problem will require a little extra thought.  You will have to make assumptions about parameter values (maybe in order to guarantee stability).  Try different parameter values and see if you can generate different dynamics for the solutions.  The steps you will follow should be similar to the ones used in the chemostat problem.  Good luck and tally-ho!