What controls the dynamics of species populations through time?

An essential question in population ecology, with clear fundamental and applied aspects. For instance, fundamentally, understanding what causes populations to fluctuate over time and across environmental gradients underlies much of population ecology. In an applied sense, the ability to forecast population dynamics and estimate population growth rates is essential to the management and conservation of species (e.g., think of fisheries).

We will start from absolute baseline assumptions about the ecology of populations, and eventually build up to incorporate a bit more realism. However, it is important to note that we will not include everything into these models, and the goal of this exercise is not to recapitulate the complexity of a natural system into a model. If we were to try to do this, we would maybe have a really well performing model for a single species in a single location. This course is about generalizing and linking natural systems to underlying theory.

A generalized model of population dynamics

\[\begin{equation} N_{t+1} = N_{t} + Births + Immigration - Deaths - Emigration \end{equation}\]

In the case of a closed population (no movement of individuals into or out of the patch), this reduces to

\[\begin{equation} N_{t+1} = N_{t} + Births - Deaths \end{equation}\]

This is to say, that populations change in terms of the numbers of individuals as a function of two things; birth and death. This is a discrete time model, as each timestep \(t\) is assumed to be on the generational time of the population of interest. Below, we will go into two of the simplest population dynamic models, but at least one which gives some interesting insight that is a bit unexpected. In doing so, we will introduce several terms, which I will try to put in bold.

Exponential growth

Discrete model

The simplest model of population dynamics is based on an exponential increase in population size given a positive growth rate. That is, the population at the next time point (\(N_{t+1}\)) is based on the population size at the current time (\(N_t\)) times the growth rate of the population (\(\lambda\)).

\[\begin{equation} N_{t+1} = \lambda * N_t \end{equation}\]

Here, \(\lambda\) is essentially \(R+1\), where \(R\) refers to the net per capita rate of growth. Why is this?

Because at each time step, we must consider the individuals that are still in the population \(N_t\) and those individuals do not die (or they might?), but instead give birth at some per capita rate \(R\). This means that the equation initially looks like:

\[\begin{equation} N_{t+1} = N_t + (R*N_t) \end{equation}\]

We then let \(\lambda = 1+R\) to reduce the equation. A nice feature of the exponential model is that it can be used to forecast population size after any number of generations. By example, let’s say we have an initial population size of 50 individuals at time \(t\). What is the population size after 3 generations when \(\lambda\) = 1.25

\[\begin{equation} N_{t+1} = \lambda * N_t = 1.25 * 50 = 62.5 \\ N_{t+2} = \lambda * N_t+1 = 1.25 * 62.25 = 78.125\\ N_{t+3} = \lambda * N_t+2 = 1.25 * 78.125 = 97.66 \\ \end{equation}\]

Or another way to look at it: \[\begin{equation} N_{t+1} = \lambda * N_t \\ N_{t+2} = \lambda * \lambda * N_t \\ N_{t+3} = \lambda * \lambda * \lambda * N_t \\ = 1.25*1.25*1.25*50 = 97.66 \\ or \\ N_{t+T} = \lambda^{T} * N_t \\ \end{equation}\]

Continuous model

A discrete model makes sense when this assumption matches the species biology. That is, if the species reproduces once per year, or if generation time can be bounded within some time window, then a discrete model might capture the relevant dynamics well. Let’s consider a system where we want that time window to be incredibly small.

\[\begin{equation} \frac{dN}{dt} = rN \end{equation}\]

where \(r\) is equal to \(b\) - \(d\) (), where \(b\) and \(d\) are per capita measures (births or deaths per individual per unit time). This \(r\) is the . When \(r < 0\), the population decreases towards 0. When \(r > 0\) the population increases exponentially (essentially geometrically, but in continuous time). This equation can be simplified back to discrete time (see Box 1.1 from the Case text), and we see the population size at time \(t\) (\(N_t\)) is

\[\begin{equation} N_{t} = N(0)e^{rt} \end{equation}\]

Where \(N(0)\) is the initial population size, \(r\) is the instantaneous rate of increase, and \(t\) is the number of timesteps. This looks quite similar to the discrete time case, except the growth rate is slightly different. This can also be used to project the expected population growth over time, where \(t\) can be any number greater than 1.

\[\begin{equation} r = ln(\lambda) \end{equation}\]

For more information on this relationship, particularly the relationship between \(R\) and \(r\), see page 9 of the Case text. Importantly, \(R\) is often referred to as ‘annual yield’ while \(r\) is a rate, often referred to as ‘annual rate’. Play with the continuous and discrete time models to explore how \(\lambda\), \(r\), and \(R\) are related.

\(\lambda <\) 1, \(r <\) 0: population decrease to 0

\(\lambda =\) 1, \(r =\) 0: population unchanging

\(\lambda >\) 1, \(r >\) 0: population increase to infinity

Equilibria

\(N = 0, r > 0\) (unstable)

\(N = 0 , r < 0\) (stable)

Assumptions of the exponential model:

The issue with this model is that there is nothing to stop it, so the time series of the population size quickly becomes exponential (as we’ll see in the coding demonstration).

Let’s explore population growth using population size in the US (as also done in the Case text). Here, we will use \(R\), which we have previously introduced. We will use a dataset that comes with the base installation on US population size from 1790 to 1970.

data(uspop)
par(mar=c(4,4,0.5,0.5))
plot(uspop, type='b', pch=16, col='dodgerblue',
  ylab='Population size')

These data consist of 19 evenly-spaced time points based on US census data. In the plot, we might look at this and characterize it as exponential growth. Let’s explore that assumption. To do so, we will visualize the data again, but this time in log space, as the slope of the line in natural log space allows us to estimate the rate of population increase (\(r\)).

data(uspop)
par(mar=c(4,4,0.5,0.5))
plot(log(uspop), type='b', pch=16, col='dodgerblue', 
  ylab='Population size')
mod <- lm(log(uspop) ~ time(uspop))
abline(mod)

summary(mod)
## 
## Call:
## lm(formula = log(uspop) ~ time(uspop))
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.33854 -0.17104  0.03883  0.20369  0.25081 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -3.774e+01  1.661e+00  -22.73 3.68e-14 ***
## time(uspop)  2.202e-02  8.829e-04   24.95 7.87e-15 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.2108 on 17 degrees of freedom
## Multiple R-squared:  0.9734, Adjusted R-squared:  0.9718 
## F-statistic: 622.3 on 1 and 17 DF,  p-value: 7.868e-15

The black line corresponds to the model fit, and the slope of the line is 0.022, meaning that there is an average increase in population size each census period of 2.2% of the population. It is a bit clearer in the Case text that the human population is not quite increasing exponentially in more recent years.

data(uspop)
par(mar=c(4,4,0.5,0.5))
plot(uspop, type='b', pch=16, col='dodgerblue', 
  ylab='Population size')
lines(as.vector(time(uspop)), exp(predict(mod)))

Incorporating variation in demographic rates

This lecture material will start on page 30 of the Case text, skipping over the role of individual movement and spatial variation in demographic rates. This is so that we can better focus on single populations and ignore space for now.

In the model above, we assume a per capita rate of increase in population size. That is, we assume that every individual in the population creates \(\lambda\) individuals in a time step \(t\). This is obviously an oversimplification, as individuals may vary in their demographic rates (typically demographic rates mean birth and death, but in this model we do not consider death, or implicitly assume that \(r=b-d\)). This individual-level variation in demographic rate is referred to as between-individual variation or demographic variability. By example, let us consider that we have a population of 10 individuals, each with their own \(\lambda_i\) value.

Further, we can consider that \(\lambda\) values may vary through time, referred to as temporal variation or environmental variation. We will explore both of these ideas independently and together, but it is important to consider how these forces are similar and different in their assumptions, structures, and outcomes.

We can draw these values from some distribution to assemble our population. Let’s consider them drawn from a normal distribution with mean 1 and variance 0.25.

# we start with 10 individuals
n = 10
# draw lambda values from a normal distribution
lambda <- rnorm(n, 1, 0.25)
hist(lambda)

mean(lambda)
## [1] 0.9951459

In the Case text, he refers to the universe distribution of individual contributions to \(\lambda\) (i.e, \(\lambda_i\) values). That is, if we could characterize the distribution of values for all individuals at a given time, we could then downscale to smaller subpopulations by taking a random sample from this known distribution. Here, our known distribution is the normal distribution with mean of 1 and variance of 0.25 (\(Normal(1,0.25)\)).

We will revist this in a minute, but let’s first think about temporal variation. In the simplest case, we can imagine that some times are better for a population than others (e.g., maybe some climatic differences affect reproduction). This would create a situation in which we draw a value of \(\lambda\) not for every individual, but for every time \(t\). What underlying assumption does this make?

Temporal variation in demographic rates

We saw the exponential growth model before. Let us now revisit this model by incorporating temporal variation in \(\lambda\). Here, we define a \(\lambda_t\) by sampling from the universe distribution, which we previously defined as being normal with mean 1 and variance 0.25.

expoGrowth <- function(n, lam){
  n*lam
}

expoDynamics <- function(n,lambda, steps=100){
  if(length(lambda) < steps){ 
    lambda <- rep(lambda, steps)
  }
  ret <- c()
  ret[1] <- n
  for(i in 1:steps){
    ret[i+1] <- expoGrowth(ret[i], lambda[i])
  }
  return(ret)
}
n0 <- 20
t <- 50

# this is the model as we previously introduced it, with a constant lambda over time
plot(1:(t+1), expoDynamics(n0, 1, steps=t),
  type='l', las=1, 
  xlab='Time', 
  ylab='Population size', 
  col='dodgerblue',
  ylim=c(0,50))

# now we incorporate temporal variation in lambda
lambdas <- rnorm(t, 1, 0.25)
lines(1:(t+1), expoDynamics(n0, lambdas, steps=t),
  col=adjustcolor(1,0.5))

# and no we do this for many times
sapply(1:100, function(x){
  lambdas <- rnorm(t, 1, 0.25)
  lines(1:(t+1), expoDynamics(n0, lambdas, steps=t),
    col=adjustcolor(1,0.5))
})

legend('topleft', bty='n', 
  c("Constant lambda", "Time varying"),
  pch=16, col=c('dodgerblue', 'black'))

Notice that the population declines to extinction in some cases, even though the mean of all the \(\lambda\) values at each time should have a mean that is similar to the true mean \(\mu\). Let’s play around with the universe distribution to clearly demonstrate how this distribution constrains the population dynamics we observe. Below, we still assume that the \(\lambda\) values are drawn from a normal distribution with mean 1, but now we will change the variance to reduce the spread in those values, from 0.25 to 0.05.

n0 <- 20
t <- 50

# this is the model as we previously introduced it, with a constant lambda over time
plot(1:(t+1), expoDynamics(n0, 1, steps=t),
  type='l', las=1, 
  xlab='Time', 
  ylab='Population size', 
  col='dodgerblue',
  ylim=c(0,50))

# now we incorporate temporal variation in lambda
lambdas <- rnorm(t, 1, 0.05)
lines(1:(t+1), expoDynamics(n0, lambdas, steps=t),
  col=adjustcolor(1,0.5))

# and no we do this for many times
sapply(1:100, function(x){
  lambdas <- rnorm(t, 1, 0.05)
  lines(1:(t+1), expoDynamics(n0, lambdas, steps=t),
    col=adjustcolor(1,0.5))
})

legend('topleft', bty='n', 
  c("Constant lambda", "Time varying"),
  pch=16, col=c('dodgerblue', 'black'))

And we see about what we’d expect, right?

The difference between arithmetic and geometric means

Population growth is a geometric process, leading two populations with the same arithmetic mean in temporal \(\lambda\) values to have markedly different dynamics. Let’s explore the example of a constant environment, corresponding to a constant \(\lambda\) to a time-varying \(\lambda_t\).

n0 <- 20
t <- 50

lambdas <- rnorm(t, 1, 0.25)
lambdaStatic <- rep(mean(lambdas), length(lambdas))


plot(1:(t+1), expoDynamics(n0, lambdas, steps=t),
  type='l', las=1, 
  xlab='Time', 
  ylab='Population size', lwd=5,
  col=adjustcolor('dodgerblue', 0.5),
  ylim=c(0,50))

lines(1:(t+1), expoDynamics(n0, lambdaStatic, steps=t),
  col=adjustcolor(1,0.5), lwd=5)

# arithmetic mean
mean(lambdas)
## [1] 0.9921282
mean(lambdaStatic)
## [1] 0.9921282
# geometric mean
prod(lambdas)**(1/length(lambdas))
## [1] 0.9611235
prod(lambdaStatic)**(1/length(lambdaStatic))
## [1] 0.9921282

Individual-level demographic variability

The above dynamics assume that population growth rate is the same for all individuals in the population \(i\), but varies through time \(t\). Here, we consider the case where each individual contributes differentially, with individual \(\lambda\) values drawn from the universe distribution. There are numerous ways to think about incorporating this. For instance, the same individual could have the same population growth rate over time, but we have no way to track that individuals identity without an individual-based model. For simplicity, we will consider the overall population growth rate at time \(t\) to be the mean value of a draw from the universe distribution for every individual in the population.

By example, we had previously defined the universe distribution as normal with mean 1 and variance 0.25. So if we start a population with 20 individuals, then at the first time point, we draw 20 random samples from this distribution.

hist(rnorm(20, 1, 0.25))

We do this at every time step, conditional upon the population size \(N_t\).

expoGrowth <- function(n, lam){
  n*lam
}

expoDynamics2 <- function(n, lambda, variance=0.1, steps=100){
  ret <- c()
  ret[1] <- n
  for(i in 1:steps){
    ret[i+1] <- expoGrowth(ret[i], mean(rnorm(ret[i], lambda, variance)))
  }
  return(ret)
}
n0 <- 20
t <- 50

# this is the model as we previously introduced it, with a constant lambda over time
plot(1:(t+1), expoDynamics2(n0, 1, variance=0, steps=t),
  type='l', las=1, 
  xlab='Time', 
  ylab='Population size', 
  col='dodgerblue',
  ylim=c(0,50))

# now we incorporate individual variation in lambda
lines(1:(t+1), expoDynamics2(n0, 1, variance=0.25, steps=t),
  col=adjustcolor(1,0.5))

# and no we do this for many times
sapply(1:100, function(x){
  lines(1:(t+1), expoDynamics2(n0, 1, variance=0.25, steps=t),
    col=adjustcolor(1,0.5))
})

legend('topleft', bty='n', 
  c("Constant lambda", "Time varying"),
  pch=16, col=c('dodgerblue', 'black'))

There is a key difference between these two types of variation; demographic variability and environmental variability. It is important to note that in the case of demographic variability we get closer to the true value of the mean of the universe distribution the more samples we take from it. That is, large populations will tend to be more like the deterministic model than smaller populations. This is to say that demographic variability, as we have incorporated it, is density-dependent.

n0 <- 20
t <- 50

layout(matrix(1:2, ncol=2))
par(mar=c(4,4,0.5,0.5))
# this is the model as we previously introduced it, with a constant lambda over time
plot(1:(t+1), expoDynamics2(n0, 1, variance=0, steps=t),
  type='l', las=1, 
  xlab='Time', 
  ylab='Population size', 
  col='dodgerblue',
  ylim=c(0,50))

# now we incorporate individual variation in lambda
lines(1:(t+1), expoDynamics2(n0, 1, variance=0.25, steps=t),
  col=adjustcolor(1,0.5))

# and no we do this for many times
sapply(1:100, function(x){
  lines(1:(t+1), expoDynamics2(n0, 1, variance=0.25, steps=t),
    col=adjustcolor(1,0.5))
})



n0 <- 2000
t <- 50
par(mar=c(4,4,0.5,0.5))
# this is the model as we previously introduced it, with a constant lambda over time
plot(1:(t+1), expoDynamics2(n0, 1, variance=0, steps=t),
  type='l', las=1, 
  xlab='Time', 
  ylab='Population size', 
  col='dodgerblue',
  ylim=c(1800,2200))

# now we incorporate individual variation in lambda
lines(1:(t+1), expoDynamics2(n0, 1, variance=0.25, steps=t),
  col=adjustcolor(1,0.5))

# and no we do this for many times
sapply(1:100, function(x){
  lines(1:(t+1), expoDynamics2(n0, 1, variance=0.25, steps=t),
    col=adjustcolor(1,0.5))
})

The above plots may look similar, but consider the relative change in population sizes. In the left panel, we have upper values of around 50 individuals, more than double than the starting population size. In the right panel, it appears to be a lot of change, but it is around 10% change in either direction.

We will not cover the last section of this chapter, which goes into the combination of temporal variation and demographic variability. If we have time in the semester, we will revisit some of these concepts, juxtaposing them under the umbrella of demographic stochasticity (density-dependent processes introducing randomness into population models through probabilistic processes of whole-numbered individuals) and environmental stochasticity (density-independent processes introducing randomness into population models through probabilistic processes dependent on environmental or temporal variation in demographic rates). These are incredibly related to the concepts introduced above, though arguably there are many ways to introduce these two forces into population models other than the ones described above.

Movement

Depending on what we cover on Tuesday, fill this in to cover the material at the beginning of Chapter 2 on movement/diffusion/random walks, excluding some of the partial diff eq bits.

walkTime <- function(p, lower, upper) {
  vals <- c(0)
  current <- 0
  while (current > lower & current < upper) {
    current <- current + ifelse(runif(1) < p, 1, -1)
    vals <- c(vals, current)
  }
  return(vals)
}

walkSpace <- function(seed=1, size=100, times=100, maxStepSize=1){
  set.seed(seed)
  landscape <- matrix(0,nrow=size, ncol=size)
  landscape[round(size/2,0), round(size/2,0)] <- 1
  for(i in 1:times){
    vert <- sample(-maxStepSize:maxStepSize, 1)
    horiz <- sample(-maxStepSize:maxStepSize, 1)
    loc <- which(landscape==i, arr.ind=TRUE)
    landscape[loc[1]+horiz, loc[2]+vert] <- i+1
  }
  return(landscape)
}


image(walkSpace(times=1))

image(walkSpace(times=100))

image(walkSpace(times=1000))

sims <- lapply(1:500, walkSpace, size=100, times=100)
sims2 <- Reduce('+', lapply(sims, function(x){x>0}))

plot(sims2[,50], type='b', pch=16, ylab='Population size', xlab='Space')

plot(sims2[50,], type='b', pch=16, ylab='Population size', xlab='Space')

Movement plus geometric growth

Above, we only consider spatial movement without reproduction. But what if we treat each cell as a potential ‘population’. That is, birth/death processes may occur at each timestep prior to dispersal, dependent on the presence of individuals in that cell.

What assumptions does this approach make?

# 1 dimensional landscapes 
getSpatial <- function(nVec, lambdaVec, times=100, p=0.2){
  ret <- matrix(0, ncol=length(nVec), nrow=times+1)
  ret[1,] <- nVec
  
  move <- function(vals){
    ret <- rep(0, length(vals))
    disp <- t(sapply(vals, function(x){
      movers <- sample(-1:1, x, prob=c(p/2, 1-p, p/2), replace=TRUE)
      return(c(sum(movers==-1), sum(movers==0), sum(movers==1)))
    }))
    
    vals <- disp[,2]
    # move left
    vals[1:(length(vals)-1)] <- vals[1:(length(vals)-1)] + disp[2:length(vals),1]
    # move right 
    vals[2:(length(vals))] <- vals[2:(length(vals))] + disp[1:(length(vals)-1),3]
    return(vals)
  }

  for(i in 1:times){
    ret[i+1,] <- ret[i,] * lambdaVec 
    ret[i+1,] <- move(ret[i+1, ])
  }
  return(ret)
}


n <- c(0,0,0,10,0,0,0)
lambdas <- rep(1, length(n))

test <- getSpatial(n, lambdas)

colz <- terrain.colors(20)
plot(test[1,], 
  pch=16, type='b', 
  ylab='pop size', xlab='spatial position')

for(i in 1:20){
  lines(test[i,], col=colz[i])
}

This results in the center of the species range having larger and larger population sizes, and the population will spread out in space in a predictable manner. That is, the spread will occur with constant velocity, under the assumption that \(\lambda\) is the same for each cell in the spatial landscape. This also means that the radius of populations occupied area will grow linearly with time. That is, there is a scaling between the time it takes to grow \(n\) times as large as it starts with (e.g., for area to 4 times as large as it currently is, it needs twice as long, if it needs to grow 9 times as large, it needs 3 times as long).

This is influenced by a number of things, with one of the most obvious things that species rarely move a single cell away.

What happens when we increase the distance that any individual may travel?

What happens if growth rate is larger?

Both of these things influence the radius of the population spread.

\[ Radius^{2}(t) = 4D_t ln(N(t)) \]

But \(ln(N(t))\) is a function of geometric growth. That is,

\[ ln(N(t)) = rt + ln(N_0) \]

where \(r\) is the growth rate, and \(N_0\) is the initial population size. \(ln(N_0)\) will be a very small number relative to \(rt\), so previous authors have argued that it can simply be ignored for simplicity at getting at spatial spread.

\[ radius^{2}(t) = 4Dt * rt\]

Simplifying this and putting in terms of units of distance/time leads to \(2rD^{0.5}\), removing time from the equation (i.e., the rate of spread of the population is independent of time, only depending on displacement (\(D\)) and population growth rate (\(r\))). The rate also does not depend on the initial population size, but note that this is the asymptotic rate of spread. We can go back to the case where there is no population growth, and this makes the initial population size matter a bit more.

\[ Radius^{2}(t) = 4D_t ln(N(0)) \]

This suggests that area occupied should increase linearly with time, as the radius of the spread increases as a function of the square root of time. Specifically, the radius increases at rate

\[ \dfrac{2\sqrt{Dln(N_0)}}{\sqrt{t}}\]

This suggests that initial population size does influence spread, and that the rate of spread (in terms of the rate of change in the radius) depends strongly on time, suggesting that the rate of change declines over time.

Spatial variation

Population growth rates may vary in space, such that some cells are just better for the population in terms of population growth rate \(\lambda\). How do we incorporate this, and what are the implications?

## let's explore a bigger landscape (where edge effects are not pronounced)
n <- rep(0, 500)
n[250] <- 100
lambdas <- rep(1, length(n))

test <- getSpatial(n, lambdas)

colz <- terrain.colors(20)
plot(test[1,], 
  pch=16, type='b', 
  ylab='pop size', xlab='spatial position')

seqz <- round(seq(1, nrow(test), length.out=20),0)
for(i in 1:20){
  lines(test[seqz[i],], col=colz[i])
}

let’s allow them to spread for longer

n <- rep(0, 500)
n[250] <- 100
lambdas <- rep(1, length(n))

test <- getSpatial(n, lambdas, times=1000)

colz <- terrain.colors(20)
plot(test[1,], 
  pch=16, type='b', 
  ylab='pop size', xlab='spatial position')

seqz <- round(seq(1, nrow(test), length.out=20),0)
for(i in 1:20){
  lines(test[seqz[i],], col=colz[i])
}

## let's incorporate variation in lambda across space 
n <- rep(0, 25)
n[round(length(n)/2,0)] <- 50
lambdas <- rep(1, length(n))
lambdas[12:17] <- 1.25

test <- getSpatial(n, lambdas, times=20)

colz <- terrain.colors(20)
plot(test[1,], 
  pch=16, type='l', ylim=c(0, 600), col=colz[1],
  ylab='pop size', xlab='spatial position')

seqz <- round(seq(1, nrow(test), length.out=20),0)
for(i in 1:20){
  lines(test[seqz[i],], col=colz[i])
}

What’s going on here?

What happens if lambda is only greater than 1 in a very small window?

n <- rep(0, 30)
n[round(length(n)/2,0)] <- 50
lambdas <- rep(0, length(n))
lambdas[15] <- 1.25

test <- getSpatial(n, lambdas, times=20)

colz <- terrain.colors(20)
plot(test[1,], 
  pch=16, type='l', ylim=c(0, 600), col=colz[1],
  ylab='pop size', xlab='spatial position')

seqz <- round(seq(1, nrow(test), length.out=20),0)
for(i in 1:20){
  lines(test[seqz[i],], col=colz[i])
}

What’s going on here?

What happens if we increase growth rate to be really high?

n <- rep(0, 30)
n[round(length(n)/2,0)] <- 50
lambdas <- rep(0, length(n))
lambdas[15] <- 2

test <- getSpatial(n, lambdas, times=20)

colz <- terrain.colors(20)
plot((test[1,]+1), ylim=c(3,400000),
  pch=16, type='l', log='y', col=colz[1],
  ylab='pop size', xlab='spatial position')

seqz <- round(seq(1, nrow(test), length.out=20),0)
for(i in 1:20){
  lines((test[seqz[i],]+1), col=colz[i])
}

what’s going on here?

How much variation is there in spatial spread?

n <- rep(0, 30)
n[round(length(n)/2,0)] <- 50
lambdas <- rep(0, length(n))
lambdas[10:20] <- 1.05

test <- list()

for(i in 1:100){
  test[[i]] <- getSpatial(n, lambdas, times=100)
  print(i)
}
## [1] 1
## [1] 2
## [1] 3
## [1] 4
## [1] 5
## [1] 6
## [1] 7
## [1] 8
## [1] 9
## [1] 10
## [1] 11
## [1] 12
## [1] 13
## [1] 14
## [1] 15
## [1] 16
## [1] 17
## [1] 18
## [1] 19
## [1] 20
## [1] 21
## [1] 22
## [1] 23
## [1] 24
## [1] 25
## [1] 26
## [1] 27
## [1] 28
## [1] 29
## [1] 30
## [1] 31
## [1] 32
## [1] 33
## [1] 34
## [1] 35
## [1] 36
## [1] 37
## [1] 38
## [1] 39
## [1] 40
## [1] 41
## [1] 42
## [1] 43
## [1] 44
## [1] 45
## [1] 46
## [1] 47
## [1] 48
## [1] 49
## [1] 50
## [1] 51
## [1] 52
## [1] 53
## [1] 54
## [1] 55
## [1] 56
## [1] 57
## [1] 58
## [1] 59
## [1] 60
## [1] 61
## [1] 62
## [1] 63
## [1] 64
## [1] 65
## [1] 66
## [1] 67
## [1] 68
## [1] 69
## [1] 70
## [1] 71
## [1] 72
## [1] 73
## [1] 74
## [1] 75
## [1] 76
## [1] 77
## [1] 78
## [1] 79
## [1] 80
## [1] 81
## [1] 82
## [1] 83
## [1] 84
## [1] 85
## [1] 86
## [1] 87
## [1] 88
## [1] 89
## [1] 90
## [1] 91
## [1] 92
## [1] 93
## [1] 94
## [1] 95
## [1] 96
## [1] 97
## [1] 98
## [1] 99
## [1] 100
layout(matrix(c(1:3), ncol=3))
# at 10 time steps
plot((test[[1]][10,]), ylim=c(1,50),
  pch=16, type='l', col=colz[1],
  ylab='pop size', xlab='spatial position')
for(i in 1:length(test)){
  lines(test[[i]][10,], col=colz[i])
}
title('10 time steps', line=-1)

# at 50 time steps
plot((test[[1]][50,]), ylim=c(1,50),
  pch=16, type='l', col=colz[1],
  ylab='pop size', xlab='spatial position')
for(i in 1:length(test)){
  lines(test[[i]][50,], col=colz[i])
}
title('50 time steps', line=-1)


# at 100 time steps
plot((test[[1]][100,]), ylim=c(1,50),
  pch=16, type='l', col=colz[1],
  ylab='pop size', xlab='spatial position')
for(i in 1:length(test)){
  lines(test[[i]][100,], col=colz[i])
}
title('100 time steps', line=-1)