---
title: "Population Dynamics"
author: "SBM"
date: "February 5, 2018"
output: pdf_document
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)

```

This week and next we are thinking about dynamical models for both population sizes and gene frequencies.  This material is very nicely covered by Otto and Day in chapters 3-5 for 1-d systems.  

To begin, let's consider population size and how it changes through time.  We will discuss both discrete and continuous time models, but to start out, let's just say that we are thinking about a change in population size that happens over a time increment of size $\Delta t$.  

In a closed population, we have the very general statement that the numbers of individuals in a given area at time $t+\Delta t$ must be given by the numbers that are there now, the plus the numbers that are born minus the numbers that die.  That is, 

$$N(t+\Delta) = N(t)+Births(t) - Deaths(t)$$  

This is totally general, absolutely correct, and at the same time completely useless.  In this formulation, we only know the next population size after we've determined births and deaths.  Not very predictive.  

In order to make it more useful we are going to need to determine models for births and deaths.  To do this, we are going to need to make some assumptions.  The simplest thing to assume is that the _per capita_ rates of natality and mortality are constant in time.  Let's call them $b$ and $d$.  The units for these rates are usually daughters per mother per unit time, so we need to multiply by the number of mothers and the time interval to get to numbers of daughters.  This gets us to 

$$N(t+\Delta t) = N(t)+b\Delta t N(t) - d\Delta t N(t)$$

or 

(@DiscreteExp) $$N(t+\Delta t) = [1+(b-d)\Delta t] N(t) $$

Now, setting the initial population size to $N(0)$ we can ask what the population size is some time in the future by iterating  (@DiscreteExp)  to find

(@DExp1) $$N(n\Delta t) = [1+(b-d)\Delta t]^n N(0) $$

During the life history section, we saw that population size grows without bounds if $R>1$, shrinks to 0 is $R<1$ and stays the same if $R=1$. 


# Moving to continuous time
Here's where our discrete time and continuous time ideas diverge.  
If we think that $\Delta t$ is the relevant time step and that not much of interest is happening between $t$ and $t+\Delta t$ then we can stick with (@Dexp1), perhaps setting $\Delta t$ to 1 generation and replacing $[1+(b-d)\Delta t]$ with the per-generation rate of increase $R$  so that (@Dexp1) becomes a model for the population size n generations into the future, i.e. 

(@Disc_n) $$N(n) = R^n N(0)$$.  


On the other hand, if we think that births and deaths can happen at any instant in time, it probably makes more sense to use continuous time.  To do this, we imagine that $\Delta t$ is really small and take the limit $\Delta t \to 0$.  We can either do this starting from (@DiscreteExp) or (@DExp1).  We will ultimately get the same answer either way, but it is probably instructive to do it both ways.  

Here's what happens if we start from (@DiscreteExp).  First, we subtract $N(t)$ from both sides and divide by $\Delta t$.  This gets us to

$$\frac{N(t+\Delta t)-N(t)}{\Delta t} = (b-d) N(t) $$
To simplify notation, let's replace $b-d$ with the instantaneous rate of per capita population growth, $r=b-d$.  Then taking the limit, and remembering the definition of the derivative, we get


(@cont1) $$\frac{dN}{dt} = (b-d) N(t) = r N(t) $$

The solution to this differential equation is given by 

(@contsol) $$N(t) = N(0)e^{rt} $$

which can be found by separating variables - i.e. rewriting (@cont1) as 

$$\frac{1}{N}dN = r dt $$
and then integrating both sides from 0 to t.  On the right side, we have $\int r dt$ which is just $rt$, assuming we started at time 0.  For the left side, recall that the integral of $x^{-1}dx$ is $\ln \left(x \right)$ (this is actually the definition of $\ln \left(x \right)$).  So we have 

$$ \int\limits_{N(0)}^{N(t)} \frac{1}{N}dN = \ln\left[N(t)\right]-\ln\left[N(0)\right] =  \int\limits_{0}^{t} rds = rt $$
__MATHEMATICAL ASIDE:__  Note that this expression opens the door to handling time-varying population growth.  If, for instance, population growth varied seasonally, with a period of $\tau$, then   

$$ \ln\left[N(\tau)\right]-\ln\left[N(0)\right] =  \int\limits_{0}^{\tau} r(s) ds = \bar r \tau$$
where the last piece on the right side comes from setting 

$$\bar r =\frac{1}{\tau}  \int\limits_{0}^{\tau} r(s) ds $$ 
which is just the average growth rate over one cycle.  Since, r is periodic, we will get exactly the same expression for $\bar r$ if we integrated over the interval $\left[\tau,2\tau\right]$.  This means that if we chose a discrete time step of size $\tau$, the dynamics would be given by $N(t+\tau)=e^{\bar r \tau}N(t)$ so (even though the growth rate is fluctuating continuously) the population would grow without bound if $\bar r > 0$.  

_Another way to take the limit_
The second way to get to a continuous time solution is to start from (@DExp1) and make use of this handy fact:


(@lim1) $$ \lim_{n \to \infty} (1+\frac{x}{n})^n = e^x $$

In order to apply this handy fact to (@DExp1) we need to think a bit about time.  We are no longer thinking about a number of steps of a fixed size, but instead thinking about what happens to the number of steps in a fixed interval of time as the step size becomes very small.  That is, if we think of the total time $t$ as fixed and also have $t=n\Delta t$, what has to happen to n as $\Delta t \to 0$? n has to go to infinity.  So rather we'll replace the $\Delta t$ bits in (@DExp1) by $\Delta t= t/n$ to get

$$N(t) = [1+(b-d)\frac{t}{n}]^n N(0) $$
Now we are ready to apply our handy fact (pretend $x = (b-d)t$) to find that


$$N(t) = \lim_{n \to \infty} [1+(b-d)\frac{t}{n}]^n N(0) =e^{(b-d)t}N(0)$$
which is the same as (@contsol).  By comparing (@consol) with (@Disc_n) we see that there is a clear relationship between $R$ and $r$. That is, $R=e^{rn}$ 

In both the discrete time and continuous time cases, the dynamical possibilities of this model are quite limited.  The population either grows or shrinks exponentially.  In continuous time, the population grows to infinity if $r>0$ (i.e. births exceed deaths), stays the same if $r=0$ (i.e. births and deaths balance exactly), or shrinks to 0 if $r<0$ (deaths exceed births).  In discrete time, the same things are possible: Increase ($R>1$), stay same ($R=1$), and shrink ($R<1$).  These are not particularly satisfying predictions since most living things dont do them...From a theoretical point of view, this tells us something important about both our model and real living things: birth and death rates cant be constant!


Now, there are a bunch of things embedded within the assumption of constant birth and death rates: no interactions with conspecifics, other species or the environment. Food is unlimited. All individuals are the same. There is no age structure, so offspring mature instantly, etc, etc.

Okay, that's alot of assumptions!  Let's try to relax a few, starting with the idea there are no interactions with conspecifics or limiting resources.

#Density dependence in continuous time
The most simple-minded way to proceed is to assert that $r$ must be some decreasing function of $N$.  It would be a little more rigorous to specify how the birth and death rates depend on $N$ separately, but the easiest way to get going is to say that $r(N) = r_{max}(1-N/K)$ so that $r$ is greatest at $N=0$, zero when $N=K$, and negative for $N>K$.  In continuous time, this means that 

(@cont_log) $$\frac{dN}{dt} = r_{max}(1-\frac{N}{K}) N $$

This is a parabola that has a max at $N=K/2$ and zeros at $N=0$ and $N=K$.  A plot with $r=2$ and $K=1$ would look like

```{r, echo=FALSE}
r=2
K=1
f<-function(x) r*x*(1-x/K)
curve(f, from=0, to=1.1,xlab="N", ylab="dN/dt",col="blue")
lines(c(0,1.1),c(0,0))
```

Let's think about what this implies for population growth.  If the population starts at some $N$ between 0 and $K$, growth is positive and the population approaches $K$.  If population size starts at some value greater than $K$ growth is negative and the population shrinks towards $K$.  

Both $N=0$ and $N=K$ are equilibria for this model; if we start there, we stay there forever.  But there's a big difference between these equilibria.  If we start at some point near $K$, the distance to $K$ will shrink through time.  This equilibrium is _stable_.  On the other hand, if we started at some population size arbitrarily close to $0$, the distance to $0$ would grow.  So, $N=0$ is an unstable equilibrium.

Here's a plot of the solution to the logistic (i.e. a plot of $N$ as a function of time) for a population starting near $0$ and another starting above $K$.

```{r,echo=FALSE}

f<-function(x,N0) N0*exp(2*x)/(1+N0*(exp(2*x)-1)/1)

#par(mfrow=c(1,2))
curve(f(x,.01), from =0, to = 10,ylim=c(0,2),col="blue")

curve(f(x,1.9), from =0, to = 10,add=TRUE,col="red")

```

It seems like growth up to $K$ is sigmoid (starts slow, increases quickly, then slows down as $N$ approaches $K$), but the decay to $K$ doesn't have an inflection.  Can you figure out why?


If we are thinking more generally about the possible behaviors this kind of model can produce,  it is somewhat better than the exponential. The population no longer goes to $0$ or $\infty$.  Now it goes to some postive population size and stays there. But how much does this result depend on the particular shape (i.e. parabola) we wrote down for $dN/dt$?  The answer, if we are interested in the long-run behavior of the model, is 'not much'.  Let's think about this some more.  What if the plot of growth v. population size looked like any of the things below?  _{I left the r code visible here so that you can see the functions corresponding to each of the graphs}_

```{r}

f1<-function(x) 6*x*(1-x)^2*sign(1-x)
f2<-function(x) 6*x^2*(1-x)
f3<-function(x) 1.5*x^.2*abs(1-x)^.5*sign(1-x)
f4<-function(x) 5*x^2*(1-x)*(1+.3*sin(30*x))

xm<-1.5
par(mfrow=c(1,4))
curve(f1,from=0, to=xm,col="blue",ylim=c(-.25,1),xlab="N",ylab="dN/dt")
lines(c(0,xm),c(0,0))
curve(f2,from=0, to=xm,col="blue",ylim=c(-.25,1),xlab="N",ylab="dN/dt")
lines(c(0,xm),c(0,0))
curve(f3,from=0, to=xm,col="blue",ylim=c(-.25,1),xlab="N",ylab="dN/dt")
lines(c(0,xm),c(0,0))
curve(f4,from=0, to=xm,col="blue",ylim=c(-.25,1),xlab="N",ylab="dN/dt")
lines(c(0,xm),c(0,0))

```

The solutions for these are below _{again, I left the code visible so you could see how to use the differential equation solver in r - you first need to install the package deSolve}_

```{r}
library(deSolve)
f1<-function(t,y,parms) list(6*y*(1-y)^2*sign(1-y))
f2<-function(t,y,parms) list(6*y^2*(1-y))
f3<-function(t,y,parms) list(1.5*y^.2*abs(1-y)^.5*sign(1-y))
f4<-function(t,y,parms) list(5*y^2*(1-y)*(1+.3*sin(30*y)))


times <- seq(from = 0, to =5, by = 0.001)
out1 <- ode(y = .1, times = times, f1, parms = NULL)
out2 <- ode(y = .1, times = times, f2, parms = NULL)
out3 <- ode(y = .1, times = times, f3, parms = NULL)
out4 <- ode(y = .1, times = times, f4, parms = NULL)


par(mfrow=c(1,4))
plot(out1[,1],out1[,2],type="line",col="blue",ylab="N",xlab="t")
plot(out2[,1],out2[,2],type="line",col="blue",ylab="N",xlab="t")
plot(out3[,1],out3[,2],type="line",col="blue",ylab="N",xlab="t")
plot(out4[,1],out4[,2],type="line",col="blue",ylab="N",xlab="t")






```

They all do more or less the same thing; rise quickly to an asymptote. The solutions in the first and third panels are less 'sigmoidal' than the other two because the growth rate near zero is very large.  But maybe the most interesting thing about these plots is how similar the second and third panels are, given the big difference in the models for population growth.  

Ultimately, what I hope to get across from this is that the long run behavior is not determined by the details of the function we write down.  What matters is that there is an unstable equilibrium at $0$ and a single stable equilibrium at $N=K$ and that growth for $N\in[0,K]$ is positive and for $N>K$ it's negative.  

_The logistic as a first-order approximation_ 
To try and formalize the idea that growth near equilibrium is approximately logistic, no matter what the model is, let's say that we have a generic per-capita growth function $r(N)$

$$ \frac{dN}{dt}=Nr(N) $$
that has a single zero at $N=K$, i.e. $r(K)=0$.  If we use a linear approximation to $r(N)$ near $K$, we would have $r(N) \approx r(K)+r'(K)(N-K)$.  Since $r(K)=0$, this is really, just $r(N) \approx r'(K)(N-K)$  Now if we factored out $-K$ from the term in parentheses, we would have $r(N) \approx -r'(K)K(1-N/K)$ and our population model would be approximately

$$\frac{dN}{dt} \approx -r'(K)KN(1-N/K) $$

Which kind of looks like the logistic, if we make the analogy $r_{max} \to -r'(K)K$  
(If you are wondering what the heck that negative sign is going there, hang on a minute, it ought to become clear in the next few paragraphs.  In any event, I wouldnt spend too much time trying to intuit what this means).  The point is just that ANY ode is approximately logistic in the vicinity of the stable equilibrium $K$.


Let's take a step back and try to think about what exactly makes $N=K$ stable and $N=0$ unstable.  In each case, the population growth rate goes from positive to negative at the stable point and from negative to positive at the unstable point.  That is, the slope is  negative for a stable euqilibrium and positive for an unstable equilibrium.  Let's see if we can make this a bit more formal.

##Stability in continuous time
Let's say we have a generic population growth model, $dN/dt = f(N)$ with an equilibrium at $N=n^*$, i.e. $f(n^*)=0$.  [Note the difference between this and the 'generic' model we had in the previous section. Here we've left the whole model unspecified, rather than the per capita growth rate.  The only reason to do this is to make the next steps a little simpler. You would get to the same conclusion using $dN/dt = Nr(N)$ but with an extra step or two].  

Now, let's think about what happens to the population near, but not at $n^*$.  In fact, let's introduce $\delta = N-n^*$ to represent the distance to the equilibrium.  Then we'll ask whether this distance grows or shrinks when it starts out close to $0$.  To do this, we'll rewrite our generic model in terms of $\delta$.  Since $n^*$ is constant, $d\delta/dt=dN/dt$ so we can just plug in to get $d\delta/dt = f(N) = f(n^*+\delta)$  Since $\delta$ starts out small, we can use a linear approximation for $f$ and get

$$\frac{d\delta}{dt} \approx f(n^*)+f'(n^*)\delta$$
where we've used the deriviative of f, $f'(n^*)$ for the slope at $n^*$.  Since $n^*$ is an equilibrium, $f(n^*)=0$ and we are left with

$$\frac{d\delta}{dt} \approx f'(n^*)\delta$$
If we look back at (@cont1),  this looks alot like our model for exponential growth.  From this we know that $\delta$ will growth if $f'(n^*)>0$ and shrink of $f'(n^*)<0$.  So, stability is determined by the sign of the derivative at the equilibrium.  Not the value - it can be small or large - just the sign.    

Let's try to apply this to the logistic model.  Here we have $f(N) = r_{max}N(1-N/K)$. The derivative of $f$ is $f'(N)=r_{max}(1-2N/K)$.  Next we evaluate this at our two euqilibria. Taking $n^*=0$ first, we have $f'(0) = r_{max}$ which is always positive, no matter what we make $r$ or $K$ (provided that $r>0$, of course!).  Next, let's evaluate stability for $n^*=K$. Plugging into the derivative, we get $f'(K)=r_{max}(1-2K/K)=-r_{max}$ which is always negative.  So the equilibrium at $n^*=K$ is always stable.  

Let's consider where we've ended up - we set out to allow population growth to decrease with population size.  The simplest way to do this is to say that the per capita growth rate decreases linearly with population size and this lead us to the logistic model.  We decided that the only thing this model can do is converge monotonically to the equilibrium. We then showed that this behavior is much more general; it will happen in any model that has a single stable equilibrium.  Last, we zoomed in on the idea of 'stability' and found that it is determined soleley by the sign of the derivative at the equilibrium. 

# Discrete time

Now, we're going to move on to discrete time dynamics with density dependence.  Recall that the density-independent models had the same behavior (exponential growth / decay) regardless of whether time was discrete or continuous. This is no longer true for models with density dependence.  Ultimately this is because in continous time, the change in population size at any instant is infinitessimally small. Conequently, in order to go from the initial population size $N(0)$ up to $K$, we need to move through every value in between.   In discrete time, however, the changes in population size are finite.  This allows us to skip over intervening population sizes which means that it may be possible to leapfrog over the equilibrium. Because of this, the dynamical possibilities of 1-d models in discrete time is quite alot broader than that for their continuous counterparts: In addition to monotonic approaches to equilibrium, they can also exhibit damped oscillations, stable limit cycles, and chaos.

The logistic model in discrete time is given by

(@disclog) $$N(t+\Delta t) = N(t) + r\Delta t N(t)\left[1-\frac{N(t)}{K}\right] $$

To make life easier we are going to assume that time is scaled such that $\Delta t = 1$, so that $N(t+1) = N(t)+rN(t)[1-N(t)/K ]$.

In discrete time, the equilibrium occurs where $N(t+1) = N(t) = n^*$. In a graph, this is where the logistic map (this is another term applied to discrete time models) crosses the $1:1$ line.

```{r,echo=FALSE}
f<-function(x) (x+2*x*(1-x))*(x<(1+2)/2)
xm=2
curve(f,from=0,to=xm,col="blue",xlab="N(t)",ylab="N(t+1)")
lines(c(0,xm),c(0,0),col="black")
lines(c(0,xm),c(0,xm),col="red")
points(1,0,col="red")

```
From the figure, we see that again, $N=K$ and $N=0$ are the equilibria. Note that the parabola doesnt hit the horizontal axis at $N=K$ in this figure because we are plotting next population size, not the change in size.  A little algebra shows that it hits the axis at $N=K(1+r)/r$. Without modification, this model says that if $N(t)>K(1+r)/r$, then $N(t+1)$ would be negative.  This doesnt make much sense and is sometimes interpretted as saying that the population goes extinct.  Alternately, we could constrain the $r$ so that this is not possible.  Let's think about that for a second.  Over values of $N \in [0,K(1+r)/r]$ the maximum population size next time occurs where $dN(t+1)/dN(t) = 0$.  Since $dN(t+1)/dN(t)=1+r(1-2N(t)/K)$ the maximum occurs when $N(t) = K(1+r)/(2r)$.  Plugging this into (@disclog) and re-arranging, we find that $N(t+1)$ will always stay in the viable range if $r \leq 3$.

# Stability in discrete time
To evaluate stability in discrete time, we use the same reasoning that we did before.  We imagine that we are starting the population near an equilibrium and ask whether the distance to the equilibrium grows or shrinks through time.  Let's start with a generic discrete time model, say $N(t+1)=f[N(t)]$ and that the equilibirum, $n^*$ satisfies, $n^*=f[n^*]$.  Now we'll say that we're starting the population close to equilibrium so $N(t) = n^* + \delta(t)$ where $\delta$ is small.  Plugging this into our generic model, for both $N(t)$ and $N(t+1)$, we get $n^*+\delta(t+1) = f[n^*+\delta(t)]$.  Since $\delta$ is small, we can use a linear approximation for $f$ to get

$$n^*+\delta(t+1) \approx f[n^*]+f'[n^*]\delta(t)]$$
Since $n^*=f[n^*]$, this reduces to

$$\delta(t+1) \approx f'[n^*]\delta(t)]$$
So, just like population sizes in (@DExp1), the distance to the equilibrium grows (i.e. $n^*$ is unstable) if $|f'[n^*]|>1$ and shrinks (stable) if $|f'[n^*]|<1$.

You might wonder why the absolute value is here.  Let's think about that for a minute. Let's say that $f'[n^*]=-1/2$. After 1 step, the deviation from equilibrium is $\delta(1) \approx -1/2 \delta(0)$. So it is half as big, but opposite in sign to the previous step. After two steps, we have  $\delta(2) \approx 1/4 \delta(0)$, so it's back to the side it started from, but only a quarter as big.

When we think about this more generally, we see that if $|f'[n^*]|$ is negative, we could write the long-run case as $\delta(t) \approx f'[n^*]^t \delta(0)=(-1)^t|f'[n^*]|^t \delta(0)$.  That is, the negative sign means that we  swap sides each step, but the part that determines whether or not we get close to equilibrium is the  magnitude of $f'[n^*]$.

Let's apply this to the logistic map. In this case, $f[N] = N + r N\left[1-N/K\right]$ and the derivative is $f'[N] = 1+r-2rN/K$.  Evaluating the stability of the equilibrium at $0$, we get $f'[0]=1+r$. So $n^*=0$ is always unstable (Unless $r=0$ in which case nothing happens ever!).  What about the equilibrium at $n^*=K$?  Here we find that

(@eq) $$f'[K] = 1-r$$

Now this is sort of interesting...Our previous thinking has led us to say that $r$ is between $0$ and $3$.  So there are three intervals of interest here.  For $r$ in $[0,1]$, the derivative at equilibrium is positive and less than one, so deviations shrink monotonically.  For $r$ in $[1,2]$ the abolute value is still less than one, so the equilibrium is stable, but since the sign is negative, the approach to equilibrium will be oscilatory.  And for $r \in [2,3]$ the equilibrium is unstable.  But, curiously, this doesnt mean that the trajectory will head off to infinity.  Actually the first thing that happens as $r$ becomes just slightly bigger than 2 is that the oscillations back and forth across the equilibrium stop decaying to zero and instead become a stable limit cycle with period 2.  Here are few graphs illustrating this range of behaviors.  Going from left to right, we have $r = {0.75, 1.75, 2.25}$

```{r}
#parameters
r<-c(.75, 1.75, 2.25)
K=1

#logistic map
f<-function(x,r) {x+r*x*(1-x/K)}


T=30 #number of generations
x=array(data=NA,T) # for storing iterations


x[1]=.011  #initial value


#make some plots
par(mfrow=c(1,3))

for (i in 1:3){

#loop to iterate map
for (t in 1:(T-1)){x[t+1]<-f(x[t],r[i])}

#right plot-resulting time series
plot(c(1:T),x,xlab="Time",ylab="N", pch=c(16), col="red")
lines(c(1:T),x,col="black")

}
```




A complete analysis of this model is beyond the scope of this course, but let's just say that as we continue to increase $r$ things get even more interesting.  The points in the period two cycle split, becoming a period 4 cycle, then an 8-cycle, and so on.  Around $r=2.57$ the trajectory becomes chaotic; it never repeats and two trajectories that start close to one another diverge rapidly.  Let's try to make some graphs to look at this.  From left to right, $r={2.5, 2.545, and 3}$ corresponding to a stable 4-cycle, 8-cycle, and chaos.

```{r, echo=FALSE}
#parameters
r<-c(2.5, 2.545,3)
K=1

#logistic map
f<-function(x,r) {x+r*x*(1-x/K)}


T=30 #number of generations
x=array(data=NA,T) # for storing iterations


x[1]=.72  #initial value


#make some plots
par(mfrow=c(1,3))

for (i in 1:3){

#loop to iterate map
for (t in 1:(T-1)){x[t+1]<-f(x[t],r[i])}
x[T]
#right plot-resulting time series
plot(c(1:T),x,xlab="Time",ylab="N", pch=c(16), col="red")
lines(c(1:T),x,col="black")

}
```


So to sum up the main ideas, we found that without density dependence, populations growth without bound or shrink to zero in both discrete and continuous time.  If we impose density dependence in continuous time, any model we write down such that there is a single stable equilibrium that is greater than zero, will exhibit monotonic growth to this equilibrium and stay there. Stability is determined by the sign of the derivative at the equilibrium.  Oscillations and chaos are impossible in 1-d.  At least 2-d are needed for cycles and 3+ for chaos.  This contrasts with the discrete time case, where it is possible to have limit cycles and chaos in just 1-d. Stability is determined by whether or not the absolute value of the derivative is less than 1 at the equilibrium.  
