The Montagues, the Capulets, and the Zombies

A socio-historical study in population dynamics (or not)

Marcus McGuff, A.Mus.D.
Austin Community College

Introduction

Systems of differential equations can be used to model many mechanical systems.  Some researchers have also tried to use these same methods to examine other less mechanistic systems, such as modeling population growth and interactions in biology.  A little contemplation will probably convince you that this is really one of those Hard (Subscript[Problems, TM]) that can give scientists fits:  There are so many different variables involved, many of which are difficult to measure (assuming that we even know what all of them are), that finding an accurate model for population growth seems moderately hopeless.  This doesn't mean that attempting to model these systems is a waste of time; we can start with fairly simple models and then work our way up to more comprehensive ones.  Along the way, we can capture some basic properties of population dynamics by examining these models.  This lab will take you through some basic ideas to model population using some very simple equations and systems using several "species" from historical fiction.

Enter the Montagues - Some basic equations to model a population

The Montague species is descended from 13th Century Italian roots.  The simplest model for their population might be that the rate at which their population increases is proportional to their size at any time, or mathematically: 

(dx/dt)=a x 

(where a>0 is some constant growth rate).

StreamPlot[{1,0.1y},{x,0,10},{y,0,10},Axes->True,AxesOrigin->{0,0},Frame->False,PlotLabel->"Exponential Growth",AxesLabel->{"t","x"}]

This exponential growth model captures the basic idea that the more Montagues there are, the more Montagues there will be; they start off growing relatively slowly, but grow faster over time as their numbers increase.  However, it makes the idealistic assumption that there will always be food for any number of Montagues to eat.  In most ecosystems, there is a definite upper limit to the size of a population that the environment can support (due to limited food supplies, build up of waste materials, overcrowding, etc.), called the carrying capacity of the environment.  So perhaps it might be better to consider a growth rate that is proportional to how far away the population is from the maximum carrying capacity:

(dx/dt)=a (1-b x)  

(where a>0 is some constant growth rate and b>0 is a constant that gives an equilibrium population value of (1/b))

StreamPlot[{1,1-0.2y},{x,0,10},{y,0,10},Axes->True,AxesOrigin->{0,0},Frame->False,PlotLabel->"Limited Growth",AxesLabel->{"t","x"}]

This limited growth model captures the idea of an upper limit for the population (why is that limit (1/b)?), so it better captures the likely long-term behavior of the Montagues.  However, this method loses some of the useful short-term features of the exponential growth model used earlier (limited growth curves always start out growing or shrinking the fastest at the very beginning).  So, we can combine those behaviors together by multiplying the two models together:

(dx/dt)=a x(1-b x)

(where a>0 and b>0 are as above).

StreamPlot[{1,y(1-0.2y)},{x,0,10},{y,0,10},Axes->True,AxesOrigin->{0,0},Frame->False,PlotLabel->"Logistic Growth",AxesLabel->{"t","x"}]

Notice that this has the short-term growth characteristics of exponential growth (as long as you start below the carrying capacity) and the long-term behavior of limited growth.  This model is called the logistic growth model.

Of course, none of these models is particularly accurate, since they completely ignore the impact of other species on the Montagues (and the impact of the Montagues on those species...).  However, we will use variations on these models to build up models for the behavior of interdependent species.

Night of the living zombies - An exploration of the Predator-Prey system

Sadly, a feral pack of zombies has just moved in near the Montagues.  As everyone knows, zombies eat people, and Montagues are particularly tasty.  What isn't as widely known is that zombie's must eat people to survive (well, you get the idea).  So, let's start with a fairly simple model; we will make the following assumptions about the Montagues (x(t)≡ the number of Montagues at time t, y(t)≡ the number of zombies at time t in some reasonable units - 1000's, maybe):

If left on their own, the Montagues population will grow exponentially (so, we're assuming an unlimited carrying capacity for Montagues at this stage): 

(dx/dt)=a x  

where a>0 would be the standard growth rate (natural birth rate - natural death rate) for the Montagues.

When zombies meet Montagues, this is not good for those Montagues.  This can be modeled by an "interaction term" of x y.  Since the Montagues get eaten, this term clearly needs to be negative, which diminishes the rate of growth of the Montagues: 

(dx/dt)=a x-b x y  

where a,b>0 and b controls the rate of "lethality" of these interactions in some sense.

For the zombies, let's make the following assumptions:

If left on their own, the zombies will just "die" away to nothing exponentially (well, decay away or whatever it is that zombies do...):

(dy/dt)=-p y

where p>0 is the rate of decay for the zombies.

When zombies meet Montagues, this is very good for the zombies.  We can introduce an "interaction term" of x y for the zombies as well, which increases the rate of growth of the zombie population:

(dy/dt)=-p y+q x y

where p,q>0 and q controls the rate at which zombies are created when they meet.

This system of equations (the second one in each case) is called the Lotka-Volterra model (after two early 20th Century mathematicians who developed them to study population growth).  Clearly, they are an over-simplification of this system, but let's examine what predictions they make about the respective populations of Montagues and zombies.

Questions:

Explain why x y makes sense as an "interaction term".  Why might this make more sense than x+y or just x or y?  (There could, of course, be far more complicated models for this, like (x^2)y or y(Sqrt[x]), but this one is a good simple first approximation.)

From the general form of these equations (see above), find all equilibrium points in terms of the parameters a,b,p,q.

For each of those equilibrium points, find the eigenvalues (again, in terms of a,b,p,q) for the local linearization (Jacobian) at those points.  Clearly state what those eigenvalues tell you about the solutions near those points.  If that behavior depends on the value of the parameters, show this clearly.  (You may assume that a,b,p,q>0, so you don't need to consider negative values of these or equilibrium points outside of the first quadrant and positive axes.)

In particular, what sort of long-term coexistence is possible for the Montagues and the zombies?  How will their populations relate to each other over time?

Use the Phase Field widget at the end of this lab to examine the impact of  the values of the parameters on the shapes of the solution curves.  In particular, briefly explain the impact of the following situations on the solutions and why those impacts do or do not make sense physically:  (For simplicity, stick with all parameters between 0 and 1.)

Montague growth rate increases?

Zombie decay rate increases?

Interactions become more lethal to the Montagues?

Interactions become more beneficial to the zombies?

The nightmare is over...  or is it?

Finally, the zombies have been out-sourced overseas, leaving room in the community for a new species, the Capulets, to move into the ecosystem.  Unlike the zombies, the Capulets don't hunt the Montagues (well, you do have to watch them closely at parties and balls...), but they are competing with the Montagues for resources (food, water, money, jobs, publicity).  So, let's take a look at one possible way to model this new system.  Let x(t)≡ the number of Montagues at time t and y(t)≡the number of Capulets at time t (in some reasonable units - 1000's, maybe).  Let's treat both groups exactly the same way; so, for each group, let's assume the following:

By themselves, let's assume they grow logistically:

(dx/dt)=x(a-b x)

where a,b>0 and a corresponds to the growth rate and the carrying capacity for the Montagues alone (for example) would be (a/b).  (I'm using the parameters a little differently here than earlier; it doesn't really matter whether you factor the a our or not.)

Competition from the Capulets serves to reduce the carrying capacity for the Montagues in proportion to the size of the Capulet population:

(dx/dt)=x(a-b x-c y)

where a,b,c>0 and c reflects the relative impact of the Capulets on the Montague growth rate.

Another way to look at this would be to multiply through:

(dx/dt)=a x-b (x^2)-c x y

where: 

a x - is the usual exponential growth term

b (x^2) - is the "brake" imposed on the population by the environmental limitations

c x y - is an "interaction term" as it was with the zombies, but here, the negative impact isn't due to being eaten (the Capulets don't do those sorts of things), but due to resources being used up by the Capulets

Of course, there is a similar equation for the Capulets:

(dy/dt)=y(p-q y-r x)

where p,q,r>0.

(Notice that if you allow for some of the constants to be 0 or negative, you can get the Lotka-Volterra equations from these as well; so while we're using this to model a very different situation, the equations we are using are really just generalizations of the previous set.)

So, there's good news and bad news here:  This is a more sophisticated model, which shows a wider range of interactions than our earlier system, but it's also more complicated to analyze what's going on (there are 6 parameters to deal with, just for a start).

Questions:

Explain why (x^2) makes sense as a "limiting term" for the population.  Why might this make more sense than x by itself?  (There could, of course, be far more complicated models for this, but this one is a good simple first approximation.)

Find the equilibrium points for this system.  There are 3 of these that always exist, regardless of the parameter values (assuming all parameter values are >0) plus one that sometimes exists, depending on the parameters.  Find a general formula for the 3 that always exists, in terms of the parameters and briefly explain which physical parameters impact the location of each one.  Explain geometrically what determines whether there is a fourth equilibrium point.  (I suggest you turn on the "nullclines" option in the Phase Field widget and look at where the green and orange nullclines intersect as you vary the parameters.)  If you were assigned to use Mathematica commands to do your lab, find a formula for that fourth equilibrium point and explain what is required mathematically for it to exist.

Find general formulas for the eigenvalues for the local linearization at the first 3 equilibrium points (the ones that always exist). For each of them, clearly show (with examples) what this tells you about the behavior of the system at those points.  (Hint:  For some of them, that behavior will change, depending on the values of the parameters, so pay attention to when each one is positive, negative, or complex.)

The fourth equilibrium point might be considered the most socially "hopeful" solution.  (Why?)  You can totally analyze this system algebraically (especially if you know how to use Mathematica to do the grunt-work), but honestly, the best way to see what's going on here is to examine how the eigenvalues change at the equilibrium points as you move the nullclines.  The nullclines that really matter here are the ones that are out in the 1st quadrant.  Consider all the different ways these two nullclines could be related (i.e., intersecting, not intersecting, X nullclines steeper, Y nullclines steeper, and combinations of these) and, for each case, show an example and clearly classify the eigenvalues of the local linearization at each equilibrium point.  (Hint:  you can move the nullclines without changing their slope by varying a and p.)  For each of these configurations, clearly explain what the long-term behavior of the system would be.

Given your work in the previous question, under what conditions can the Montagues and the Capulets live in peaceful coexistence?  (You do not have to give restrictions on the variables, but rather explain what must be true about the nullclines for this to happen.)

It has been suggested that the Montagues and the Capulets can live together in harmony if b q>c r (assuming there is a fourth equilibrium point and it is in the first quadrant).  Does this agree with your conclusion from above?  (Show why or why not...)  Explain what this inequality implies about the physical situation.  (Hint:  Try to figure out the meaning of each side in terms of what the different parameters do.)

Extra Credit - Happily ever after...

So, the zombies have heard about the happy new community with the Montagues and the Capulets and, deciding that there's no place like home, they move back into the neighborhood.  Can these three groups live together in -- well, can they live together?  Here are three possible ways we could adjust our model to incorporate the zombies into the population, the easy way, the hard way, and the insane way:

One way to incorporate the zombies would be to not worry about how many of them there were, but to just focus on their impact on the other two groups.  The simplest way to do this would be to include a "harvesting" term in the Montague and Capulet equations.  The simplest of these would be a "constant effort" harvesting term, which means simply modifying the Montague equation, for example to be:

(dx/dt)=a x-b (x^2)-c x y-d x

where d is rate of attrition due to zombies ("constant effort" refers to a fixed percentage of the population).  Of course, the Capulets would need a similar modification.  Without doing a whole lot of computation, what is the impact of this type of "harvesting" on the populations?  How does it change your conclusions from the previous section about mutual coexistence (if at all)?  (Hint:  If you have to do much work here, you probably aren't thinking about this the right way...)

A more obvious way to do this might be to simply assume the zombies harvest Montagues and Capulets at a constant rate each year.  So, the Montague equation would look like:

(dx/dt)=a x-b (x^2)-c x y-d

(and similarly for the Capulets).  This is easy to understand, but it has some very extreme consequences.  Take a look at these equations in the Phase Field widget and see what happens to the nullclines and critical points (hint:  keep your harvesting constants quite small or there's no hope).  Explain why equilibrium points seem so much harder to come by here (give some examples) and what sorts of long-term prospects there are for the two groups.  (You do not have to give an exhaustive treatment of all possible outcomes; you are primarily looking to understand why such a simple change has such severe consequences.)

Of course, the best thing to do would be to set up a 3-dimensional system to model the populations for all three.  Using the models developed earlier in the lab, write out a system of 3 equations (x ≡ Montagues, y ≡ Capulets, z ≡ zombies) to model the complete system.  Assume that the Montagues and Capulets behave in a logistic manner if they were alone and the zombies behave and interact with the other two species in the same way we did earlier (so, combine the two models together).   Good news:  You do not need to try to solve this system or even find its equilibrium points; a real examination of this system might make for a good research paper, but it's far beyond the scope of this lab.

The phase field widget