Lab 11: Stochastic Simulation of the Logistic Equation

In this lab, we will simulate the logistic model, allowing for the fact that parents can
have a variable number of offspring.  That is, we will incorporate demographic stochasticity.

In any particular generation, all the parents will have, on average, n[t] + r n[t] (1-n[t]/K) offspring
(from the logistic), or 1 + r (1-n[t]/K) offspring per parent.  For example, if r=0.5, then, on average, 
parents will have 1.5 kids/parent when the density of the population is low.  However,
any particular parent might have 0, 1, 2, ... kids.

We will let the probability that a parent has X kids be given by the Poisson distribution
with mean equal to the expected number of kids per parent (= mu).

(We will assume that this population is asexual or hermaphroditic, so that we don't
have the added complication of having to keep track of the numbers of males and females.)

We will be using some unusual graphics commands and asking Mathematica about discrete distributions,
so we first need to load some additional packages:

<<Graphics`Graphics`
<<Statistics`DiscreteDistributions`

Some preliminaries

After these packages have been loaded, Mathematica knows about the PoissonDistribution:

?PoissonDistribution

We'll often ask Mathematica to tell us the probability of getting X offspring based on the
PoissonDistribution. To do this, we will use the command, PDF:

?PDF

For example, let's say that we wanted to show a bar chart for the probability of having 0, 1, 2, ... kids,
if the offspring distribution is Poisson and if 1.5 is the mean number of kids:

?BarChart

BarChart[Table[{PDF[PoissonDistribution[1.5],x],x},{x,0,4}]]

QUESTION 1:  From the above graph, what is the chance of having 0 kids?

In our simulations, we will use the random number generator in Mathematica:

?Random

For example, the following command gives us a random number of offspring based on the 
Poisson Distribution with mean 1.5:

Random[PoissonDistribution[1.5]]

QUESTION 2:  Hit shift-enter on the above equation 10 times and count how many times you get 0.
Is this answer close to your answer to QUESTION 1?
[It would be even closer if you repeated this 100 times...and even closer if you did it 1000 times --
but you would get awfully bored!!]

Simulating the Logistic

The answer to QUESTION 2 tells us that if we started a population with only 1 individual 
where r=0.5 (1.5 kids per parent), there is a good chance that the population will go extinct
within one generation. 

Let's now develop a Mathematica program that simulates logistic growth over time
when the number of children that a parent has is distributed according to the Poisson.

The Poisson distribution is special in that the sum of n Poisson distributions with mean
mu is equal to a Poisson distribution with mean n*mu.  We will use this fact in the following
way.  At time t, there are n[t] parents that each have 1 + r (1-n[t]/K) kids on average.
If the number of kids per parent is Poisson distributed, then the number of kids produced
by the whole population will be Poisson distributed with mean n[t]*(1 + r (1-n[t]/K)).

Clear[num, t, r, k, n0, trial]
num[t_,r_,k_,n0_,trial_]:=num[t,r,k,n0,trial]=
    Random[PoissonDistribution[
	num[t-1,r,k,n0,trial]+r*num[t-1,r,k,n0,trial]*
	(1-num[t-1,r,k,n0,trial]/k)
   + 10^(-16) ]]
num[0,r_,k_,n0_,trial_]:=n0

There are two odd features about this simulation program.  One is the 10^(-16)
added to the basic logistic equation.  The reason I did this is that the Poisson Distribution
is not actually defined for mu=0, so if the population ever goes extinct then 
Mathematica doesn't know what to do.  By adding, 10^(-16) to the expected number
of offspring, I can avoid mu ever hitting zero exactly without influencing the dynamics
(since this number is so tiny).

The second is the variable "trial" that has been included.  The reason for this is that
you may want to run more than one trial with the same set of parameters.  To do this, just 
run the program once setting trial to one and then repeat setting trial to 2, 3, etc:

ListPlot[Table[num[t,0.5,100,10,1],{t,0,50}],
	PlotJoined->True,
	PlotRange->{{0,50},{0,150}}]

[The above graph sets r=0.5 K=100 and n0=10.]

QUESTION 3:  Rather than approaching carrying capacity and staying 
there forever, the population hovers about carrying capacity.
Roughly what range of population sizes is seen once the initial phase
of approaching the carrying capacity is over?

QUESTION 4:  With K=100, as in the above graph, extinction is fairly
unlikely.  Explore decreasing values of K to see what population size is
small enough to have a reasonably high risk of extinction within 50 generations.
Start these simulations with n0 equal to the value of the carrying capacity.
If you have time, quantify your results, saying for example, that the population
went extinct within 50 generations in x/10 trials with K=...

NOTE:  In this lab, you have simulated only one source of uncertainty in population 
growth.  That is, you have accounted for the fact that each parent does not always
have the number of offsparents expected.

In nature, there are additional sources of uncertainty/randomness in the growth of a
population, including variable environmental conditions over time and space and
variable sex ratios.  Furthermore, when the population gets very small, additional
effects could hamper growth, including inbreeding depression and difficulties finding
a mate.  Therefore, your answer to QUESTION 4 probably gives a very conservative 
picture of the risk of extinction for small populations.