From Population Dynamics to Partial Differential Equations

Michael Kerckhove

Differential equation models for population dynamics are now standard fare in single-variable calculus. Building on these ordinary differential equation (ODE) models provides the opportunity for a meaningful and intuitive introduction to partial differential equations (PDEs). This article illustrates PDE models for location-dependent carrying capacities, migrations, and the dispersion of a population. The PDE models themselves are built from the logistic equation with location-dependent parameters, the traveling wave equation, and the diffusion equation. The approach presented here is suitable for a single lecture, reading assignment, and exercise set in multivariable calculus. Interactive examples accompany the text and the article is designed for use as a CDF document in which some of the input can remain hidden.

Introduction: Ordinary Differential Equations and Population Dynamics

The ordinary differential equations with which students are most familiar are the equations for exponential and logistic population growth (see [], for example). Historically, Thomas Malthus initiated the mathematical treatment of population dynamics []. His investigations into the consequences of exponential growth in human populations coupled with linear growth in agricultural production led him to a pessimistic view of the eventual fate of humankind, involving misery and starvation.

Pierre-Francois Verhulst and others considered how to model factors that might retard exponential growth, among which the roles of disease and famine had been suggested already by Malthus. In his 1838 paper [], Verhulst transformed the discussion into a mathematical model by writing

y'=m y-φ(y).

Here m y corresponds to exponential growth and the term -φ(y) corresponds to the effect of obstacles that slow the growth rate and that depend on the population size y. A simple form for the obstacle term is φ(y)=k (y^2), and setting m=k A yields the factored form

y'=k y(A-y),

the form we use in this article. Verhulst solved this equation; showed how it led, in the long run, to a bounded population size; and used actual population data to calculate the size of the limiting population for Belgium according to the model. He termed the graph of the solution the logistic curve: the word logistic has its origins in the Greek “logistikē”—the art of calculating—and Verhulst was the first to be able to calculate the effects of obstacles on the growth of human populations.

The logistic equation is especially useful for introducing ideas involving the qualitative behavior of solutions, in particular the notion of a stable equilibrium. We now write the relevant ODE as

y'=k y(A-y); k,A>0,y(0)=(Subscript[y, 0])≥0.

In the logistic model, y(t) is the population size at time t. Using the terminology of population dynamics, the parameters in this model are the growth constant, k, which controls the rate of population growth near y=0; the carrying capacity, A, which limits the eventual size of the population; and the initial population size (Subscript[y, 0]). DSolve gives an explicit formula for the solution to the corresponding initial value problem.

odeLogistic=D[y[t],t]== k y[t] (A-y[t]);
yLogistic=y[t]/.First[Quiet@DSolve[{odeLogistic, y[0]== y0},y,t]]

The usual textbook formula for the solution divides both numerator and denominator by (e^(A k t)) to obtain y(t)=((A (Subscript[y, o]))/((Subscript[y, 0])+(A-(Subscript[y, 0]))(e^(-A k t)))). From this simplified formula it is easy to check that y(0)=(Subscript[y, 0]) and that (lim_(t→∞))y(t)=A for (Subscript[y, 0])≠0.

The built-in Mathematica function Manipulate lets you vary the parameters k, A, and (Subscript[y, 0]) interactively to see the effect of varying those parameters graphically. Equilibrium solutions for the model (values of y for which y'=0) occur at y=0 and y=A. The carrying capacity A is indicated by the dashed horizontal line (slope zero) in the figure, while the solution curve is purple. The initial value (Subscript[y, 0]) is the height of the curve at its intersection with the y axis. How would you describe the effect of varying the growth constant k?

Manipulate[
Plot[{carryingCapacity,yLogistic/.{k->growthConstant,A->carryingCapacity,y0->initialPopulation}},
 {t,0,600},PlotRange ->{-.3, 2},AxesOrigin ->{0, 0},AxesLabel->{t,y},
 Ticks->{{{150," "},{300," "},{450," "},{600," "}},{{.6," "},{1.2," "},{1.8," "}}},
 PlotStyle ->{Dashed,Thick},PlotLabel ->"Logistic Growth",Epilog->{Text[Row[{Style["y",Italic]," = ",Style["A",Italic],"\nstable"}],{30,carryingCapacity+.15}],
Text[Row[{Style["y",Italic]," = 0\n  unstable"}],{30,-.15}]}
],
 {{growthConstant,0.01,"growth constant k"},0.01,0.2,.001,Appearance->"Labeled"},
 {{carryingCapacity,1,"carrying capacity A"},0.5,1.5,.001,Appearance->"Labeled"},
 {{initialPopulation,0.1,"initial population (Subscript[y, 0])"},0.01,1.9,.001,Appearance->"Labeled"},SaveDefinitions->True]

The equilibrium solution y=A is stable, since whenever (Subscript[y, 0])>0, (lim_(t→∞))y(t)=A. (However, y(t)=A for some time t only if (Subscript[y, 0])=A, and hence y(t)=A for all times t.) Since solution curves starting near y=0 do not satisfy (lim_(t->∞))y(t)=0, this equilibrium solution is not stable. This information can be captured in a “phase diagram,” shown in the next section. The diagram visually condenses the claims just made: solution curves diverge away from y=0 and converge toward y=A. Use the sliders above to verify these claims, based on the plots of the various solution curves. Since this article is primarily concerned with the shapes of solution curves for the differential equations models we study, graphs do not include units on the axes.

A Key Concept: The Time Evolution of a Spatial Profile Curve

Now consider the case of a population that is spread out along a narrow strip of land or sea, which we will think of as a line. An initial study of the population might result in a function (Subscript[u, 0])(x) giving the size of the population at location x along the strip. But population size is a dynamic quantity, so we should really employ a function of two variables, time and location, in our attempt to model its evolution. Thus, the population size at time t and location x is u(t,x) with u(0,x)=(Subscript[u, 0])(x) recording the observed initial population distribution.

We model the time evolution of the population using the partial differential equation

((∂u)/(∂t))=k u ( A(x)-u); k,A(x)>0, u(0,x)=(Subscript[u, 0])(x)≥0.

This is just a modification of the logistic equation in which the carrying capacity A now depends on location A(x). This allows the new model to account for environmental conditions that vary from place to place. The function u=u(t,x) is now the unknown function that solves the PDE. In this setting, the role of the initial value, formerly played by the number (Subscript[y, 0]), is being played by the function u(0,x)=(Subscript[u, 0])(x). Essentially, what we have is a separate initial value problem at each location.

We will call the graph of u(0,x) the initial profile curve. Think of u(t,x) for a given value of t as the profile curve for the population at time t. Its graph shows how the total population is distributed along the line. The left-hand side of our PDE gives the time rate of change of the function u(t,x), and we will interpret the right-hand side as the rule governing the time evolution of the initial profile curve u(0,x). This idea, the time evolution of a profile curve, is crucial to what follows.

Here is an example that illustrates this point of view. The location-dependent carrying capacity is given by the function A(x)=300+100cos(((2π)/50)x), and the growth constant is k=0.04. Plots of A(x) and (Subscript[u, 0])(x) are shown below.

GraphicsRow[{
Plot[3+Cos[((2π)/50)x],{x,0,100},PlotRange->{0,4.5},
 Ticks->{{{25," "},{50," "},{75," "},{100," "}},{{1.5," "},{3," "},{4.5," "}}},
 PlotLabel->Row[{"Carrying Capacity ",HoldForm[A[x]]}]],
Plot[((3+Cos[((2π x)/50)])/(1+Exp[3+Cos[((2π x)/50)]])),{x,0,100},PlotRange->{0,4.5},
 Ticks->{{{25," "},{50," "},{75," "},{100," "}},{{1.5," "},{3," "},{4.5," "}}},
 PlotLabel->"Initial Population (Subscript[u, 0])(x)"]
},ImageSize->{400,150}]

We might predict that the initial population will evolve toward the carrying capacity. To see whether this prediction is accurate, we use Manipulate to show the evolving profile curve as t increases. Use the “Reset and Play” button to show the time evolution.

Module[{u},uLogistic=u[t,x]/.First[Quiet@DSolve[{D[u[t,x],t]==k u[t,x](A-u[t,x]) (*logistc PDE *),u[0,x]==u0},u[t,x],{t,x}]];

Manipulate[Plot[{3+Cos[((2π x)/50)] (* carrying capacity *),
 uLogistic/.{k->0.04,A->3+Cos[((2π x)/50)],u0->((3+Cos[((2π x)/50)])/(1+Exp[3+Cos[((2π x)/50)]])),t->time}},{x,0,100},
 AxesLabel->{x,Style["u",Italic]},PlotRange->{0,4.5},PlotStyle->{Dashed,Automatic},
 Ticks->{{{25," "},{50," "},{75," "},{100," "}},{{1.5," "},{3," "},{4.5," "}}}],
{time,0,100},SaveDefinitions->True,ControlType->Trigger]
]

Here is a thought experiment: suppose we were to change the initial population curve, leaving the carrying capacity at each location x the same; would you expect the new initial profile curve to evolve to A(x)? Can you explain why this is so?

Another thought experiment: what would be the effect of a larger value of k on the time it takes for the initial population to approach equilibrium?

Once again, we can summarize the investigation with a “phase diagram.”

Module[{Α,plt,offset,arrows,lbl},
Α[x_]:=3+Cos[((2π)/50)x];
plt=Plot[Α[x],{x,0,100}];
offset=.2;
arrows=Graphics[
{Table[{Arrowheads[.025],Arrow[{{x,Α[x]+4offset},{x,Α[x]+offset}}]},{x,5,95,18}],
Table[{Arrowheads[.025],Arrow[{{x,Α[x]-4offset},{x,Α[x]-offset}}]},{x,5,95,18}],
{Gray,Table[{Arrowheads[.025],Arrow[{{x,0+offset},{x,0+4offset}}]},{x,5,95,18}]},
{Gray,Table[{Arrowheads[.025],Arrow[{{x,0-offset},{x,0-4offset}}]},{x,5,95,18}]}}];
lbl=Graphics[Text[HoldForm[A[x]],{105,4}]];
Show[{plt,arrows,lbl},AxesLabel->{x,u},PlotRange->{{0,110},{-1,5}},AxesOrigin->{0,0},
 Ticks->{{{25," "},{50," "},{75," "},{100," "}},{{1.5," "},{3," "},{4.5," "}}},
 PlotLabel->Row[{"Stability of Equilibria for ",HoldForm[D[u,t]==k u (A[x]-u)]}]]
]

As before, the diagram indicates the equilibrium solutions, which are the curves along which ((∂u)/(∂t)) is zero at every point. These are graphs of the functions u(t,x)=0 and u(t,x)=A(x). The arrows indicate that if (Subscript[u, 0])(x)>0 at every location x, then (lim_(t→∞))u(x,t)=A(x), so that A(x) can be termed a stable equilibrium for this PDE.

A final note: another variation of logistic growth is given by ((∂u)/(∂t))=k(x) u(A-u), where the carrying capacity is constant, but the growth constant k, that is, the speed with which the initial population approaches the equilibrium population, depends on location.

Migrating Populations

In the previous example, the time evolution of the profile curve was governed by a right-hand side that involved only the values of the function u(t,x). The next PDE we consider has the time evolution at location x governed by the slope of the profile curve at x. The simplest equation of this type has the time rate of change of u proportional to its slope, ((∂u)/(∂t))=-c((∂u)/(∂x)), with u(0,x)=(Subscript[u, 0])(x).

We can interpret the solution to this PDE in terms of a migratory population with the parameter c controlling the speed of the migration as shown in the output from Manipulate displayed below (use the “speed” buttons to vary c). Data for the migration of gray whales along the California coastline would provide a real-life example of the phenomenon we are modeling. In the graphic, it is best to use the “Reset” button before changing the speed.

Module[{U},UMigration=U[t,x]/.First[DSolve[{D[U[t,x],t]==-c D[U[t,x],x](* migration PDE *),U[0,x]==PDF[NormalDistribution[25,6],x]},
U[t,x],{t,x}]];
Manipulate[Plot[UMigration /.{c ->speed,t->time},{x,0,100},AxesLabel->{x,u},
 Ticks->{{{25," "},{50," "},{75," "},{100," "}},{{.02," "},{.04," "},{.06," "}}},
 PlotLabel->Row[{"Migrating Population\n",HoldForm[D[u,t] = -c D[u,x]]}]],
 {{speed,10,"speed c"},{10, 50}},{{time,0,"time t"},0,1},
 ControlType->{SetterBar,Trigger},SaveDefinitions->True]
]

DSolve gives the following formula for the solution to the PDE ((∂u)/(∂t))=-c((∂u)/(∂x)), which we might call the migration equation. As usual, the initial population profile is u(0,x)=(Subscript[u, 0])(x).

Module[{u},
pdeMigration=D[u[t,x],t]==-c D[u[t,x],x];
uMigrationSoln[c_,u0_][t_,x_]:=u[t,x]/.First[DSolve[{pdeMigration,u[0,x]==u0[x]},u[t,x],{t,x}]];
Row[{"u[t,x] = ",uMigrationSoln[c,u0][t,x]}]
]

You may have seen in previous courses how changing f(x) to f(x-a) changes the graph. Consider the y value (Subscript[y, 0])=f(0) on the graph of f(x); this same y value occurs on the graph of f(x-a), but now it occurs at x=a, since (Subscript[y, 0])=f(0)=f(a-a). Likewise, you see that (Subscript[y, 1])=f((Subscript[x, 1]))=f(((Subscript[x, 1])+a)-a), so that the same y value, (Subscript[y, 1]), that occurred at (Subscript[x, 1]) along the graph of f(x) will occur along the graph of f(x-a) at input value (Subscript[x, 1])+a. Since this relationship holds for every (Subscript[x, 1]), we conclude that the graph of f(x-a) is a horizontal translation of the graph of f(x) to the right by a units. In our setting, (Subscript[u, 0])(x-c t) gives a graph that is the graph of (Subscript[u, 0])(x) displaced to the right by c t units. Note how displacement = speed × time.

An application of the chain rule shows that (Subscript[u, 0])(x-c t) satisfies the migration equation for any initial population distribution (Subscript[u, 0]). This model also applies to the familiar stadium wave and is commonly referred to as the traveling wave equation. Jean le Rond d’Alembert first introduced PDEs for waves and methods for solving them in his study of vibrating strings []. In fluid dynamics, ((∂u)/(∂t))=-c((∂u)/(∂x)) is a simplified version of the transport equation that models how dissolved solids travel along with a current.

Seasonal Migration

The migration of gray whales is seasonal, ranging from the southern Baja peninsula near the Tropic of Cancer to the Chukchi Sea north of the Arctic Circle []. To model the periodic nature of the whales’ migration, we can translate the initial graph using a periodic function to obtain u(t,x)=(Subscript[u, 0])(x-A sin(c t)). Now u(t,x) satisfies the periodic migration equation

((∂u)/(∂t))=-A c cos(c t)((∂u)/(∂x));u(0,x)=(Subscript[u, 0])(x)≥0,

with (Subscript[u, 0])(x) as the initial profile curve. You can see how the factor -A c cos(c t) comes from the chain rule. The parameter c still governs the speed of the migration; can you predict the feature of the graph that is related to the parameter A?

Module[{U},USeasonal=U[t,x]/.First[DSolve[{D[U[t,x],t]==-A c Cos[c t]D[U[t,x](* seasonal PDE *),x](* seasonal PDE *),U[0,x]==PDF[NormalDistribution[50,6],x]},
 U[t,x],{t,x}]];
Manipulate[Plot[USeasonal/.{A->amplitude,c->speed,t->time},{x,0,100},
 AxesLabel->{x,u},PlotLabel->Row[{"Seasonal Migration\n",
HoldForm@D[u,t]," = ",-A c Cos[c t],HoldForm@D[u,x]}],
 Ticks->{{{25," "},{50," "},{75," "},{100," "}},{{.02," "},{.04," "},{.06," "}}}], 
 {{amplitude,10,"A"},{10,30}},{{speed,10,"speed c"},{ 10,20}},
 {{time,0,"time t"},0,2(π/10)},ControlType->{SetterBar,SetterBar,Trigger},SaveDefinitions->True]
]

Dispersion

In response to overcrowding in one location, a population may disperse over a wider area. A model is provided by making the time evolution of the population graph proportional to the concavity of the graph. To see why this works, recall that a local maximum is characterized via the second derivative test as occurring where the tangent line is horizontal and the graph is concave down (negative second derivative). Thus, the model stipulates that near a local maximum, the future population will decrease. We write the dispersion equation as

((∂u)/(∂t))= d(((∂^2)u)/(∂(x^2))); d>0,

and refer to the proportionality constant d as the dispersion coefficient. The graph of the so-called fundamental solution to this equation,

u(t,x)=(1/(Sqrt[4 π d t]))exp(((-(x^2))/(4 d t))),

is a bell-shaped curve whose standard deviation increases as time moves forward, starting from some time (Subscript[t, 0])>0.

This verifies that this function solves the dispersion equation.

dispersionEquation=D[u,t]==d D[u,x,x];
u[t_,x_]:=(1/(Sqrt[4π d t]))Exp[((-(x^2))/(4d t))];
dispersionEquation/.{u->u[t,x]}

The larger the dispersion coefficient, the faster the dispersion takes place; use the buttons for the dispersion coefficient and the plus and minus buttons for time to show that the graph with dispersion coefficient 5 at time 20 is identical to the graph with dispersion coefficient 20 at time 5.

uDispersion[d_][t_,x_]:=(20/(Sqrt[4π d t]))Exp[((-(((x-50))^2))/(4d t))];
Manipulate[Plot[uDispersion[d][t,x],{x,0,100},PlotRange->{0,2.75},AxesLabel->{x,u}, Ticks->{{{25," "},{50," "},{75," "},{100," "}},{{.8," "},{1.6," "},{2.4," "}}}], 
 {{d, 5, "dispersion coefficient d"},{ 5, 20}}, 
 {{t,1,"time t"},1,100,Appearance->"Open",
AppearanceElements->{"InputField","StepLeftButton","PlayPauseButton","StepRightButton"}},
ControlType->{SetterBar,Automatic},SaveDefinitions->True]

In these graphs, the total population is represented by the area under the curve and does not change. As time moves forward, the dispersion process leads to a population that is, for practical purposes, uniformly distributed. In chemistry, our dispersion equation is called the diffusion equation: the concentration of dye in a liquid diffuses toward a uniform concentration over time. The origins of this equation are in the study of how heat dissipates from a source. Joseph Fourier began circulating results on the mathematics of heat conduction in 1807, when the heat equation was born. He published a treatise on the subject in 1822 [].

Combined Models

For a population that disperses as it migrates, we can combine the two previous models via the advection-dispersion equation

((∂u)/(∂t))=d(((∂^2)u)/(∂(x^2)))-c((∂u)/(∂x)); d>0.

As you might expect, a fundamental solution involves replacing x by x-c t:

u(x,t)=(1/(Sqrt[4π d t]))exp(((-(((x-c t))^2))/(4 d t))).

Use the “Reset and Play” button to see how the initial data evolves.

uAdvectionDispersion[c_,d_][t_,x_]:=(20/(Sqrt[4π d t]))Exp[((-(((x-20-c t))^2))/(4d t))];
Manipulate[Plot[uAdvectionDispersion[5,10][t,x],{x,0,100},PlotRange->{0, 2},
 AxesLabel->{x,u},
 Ticks->{{{25," "},{50," "},{75," "},{100," "}},{{.6," "},{1.2," "},{1.8," "}}}],
 {{t,1,"time t"},1,10},ControlType->Trigger,SaveDefinitions->True]

A more interesting combined model can be employed in the study of invasive species. The figure below shows data related to the spread of gypsy moths in Wisconsin over a ten-year period from 1996 to 2006 [].

The darkest regions on the maps represent 11 or more moths per trapping area, which we will take as the statewide carrying capacity. We let x represent the distance from the Mississippi River on the western edge of the state, with the total distance to Lake Michigan roughly 180 miles. We combine dispersion with logistic growth to obtain the PDE

((∂u)/(∂t))=d(((∂^2)u)/(∂(x^2)))+k u (A-u); d,k,A>0.

This equation is solved numerically (think of applying Euler’s method to the entire profile curve), and the time evolution of the initial profile curve is plotted below.

Module[{U,f},
f[t_,x_,AA_,kk_]:=AA 30/(30+(AA-30)Exp[-AA kk (x-144)]);
solnFK[kk_,AA_,dd_]:=NDSolve[{D[U[t,x],t]==k U[t,x](A-U[t,x])+d D[U[t,x],x,x](* Fisher–Kolmogorov PDE *)/.{k->kk,A->AA,d->dd},
 U[0,x]==f[0,x,AA,kk],U[t,0]==f[t,0,AA,kk],U[t,180]==f[t,180,AA,kk]},U,{x,0,180},{t,1,10}];
UFisherKolmogorov=U/.First[solnFK[.01,100,8]];
Manipulate[Plot[UFisherKolmogorov[t,x],{x,0,180},PlotRange->{0,100},
 AxesLabel->{x,u},Ticks->{{{45," "},{90," "},{135," "},{180,"180"}},{{33," "},{66," "},{100,"11"}}}],
 {{t,1,"time t"},1,10},ControlType->Trigger,SaveDefinitions->True]]

You can see that the model does not fit the data precisely: the “edges” of the gypsy moth population shown on the two maps are not as uniformly sharp as the model would predict. Biological systems are complicated, particularly because an element of randomness may come into play. Still, the concept of a model that includes dispersion and logistic growth, a variant of the Fisher–Kolmogorov equation, seems to capture most of the essential features of the gypsy moth data.

Finally, note that the total gypsy moth population is not constant, as it would be under a dispersion-only model. Here the growth term k u (A-u) on the right-hand side of the equation overcomes dispersion, and as time moves forward, the population spreads across the whole state, while its size approaches the carrying capacity A at each location.

Summary

Using DSolve and Manipulate, we have illustrated three basic partial differential equations and interpreted the equations via a profile curve (a function of x at a specific time t) that evolves with time, according to a rule that depends on function value, slope, or concavity. In terms of population dynamics, these equations were used to model location-dependent carrying capacity, migration, and dispersion. Combinations of the basic models were used to capture more complicated phenomena in a conceptually simple way that is suitable for students in a multivariate calculus class and that takes advantage of their previous experience with ODEs. A beginning textbook for PDEs is []. Further applications to population dynamics appear in [].

References

J. Stewart, Calculus: Concepts and Contexts, 4th ed., Belmont, CA: Brooks/Cole, 2009.

T. Malthus. An Essay on the Principle of Population, London: J. Johnson, 1798.

P.-F. Verhulst, “Notice sur la loi que la population poursuit dans son accroissement,” Correspondance mathématique et physique, 10, 1838 pp. 113–121. books.google.com/books?hl=fr&id=8GsEAAAAYAAJ&jtp=113# v=onepage&q&f=false.

J. d’Alembert, “Recherches sur la courbe que forme une corde tendue mise en vibration,” Histoire de l’Académie royale des sciences et belles-lettres de Berlin, 3, 1747 pp. 214–219. www.math.jussieu.fr/~daubin/cours/Textes/Dalembert-HAB-1746-TFA.pdf.

Ocean Institute. “Whales of the World.” (Sep 2, 2011) www.ocean-institute.org/visitor/gray_whale.html.

J. Fourier, Théorie analytique de la chaleur, Paris: F. Didot, 1822.

P. Tobin and L. Blackburn, “Long-Distance Dispersal of the Gypsy Moth (Lepidoptera: Lymantriidae) Facilitated Its Initial Invasion of Wisconsin,” Environmental Entomology, 37(1), 2008 pp. 87–93. Figure used with permission of the authors. www.treesearch.fs.fed.us/pubs/14617.

W. Strauss, Partial Differential Equations: An Introduction, 2nd ed., Hoboken, NJ: John Wiley & Sons, 2009.

E. Holmes, M. Lewis, J. Banks, and R. Veit, “Partial Differential Equations in Ecology: Spatial Interactions and Population Dynamics,” Ecology, 75(1), 1994 pp. 17–29. www.esajournals.org/doi/abs/10.2307/1939378.