In petrelharp/landsim: Simulate Populations on Landscapes

fig.dim <- 3
knitr::opts_chunk$set(fig.width=3*fig.dim,fig.height=fig.dim,fig.align='center')
library(Matrix)
library(raster)
library(landsim)
set.seed(42)

Here's a simple layer we'll work with. Note that units are in meters.

habitat <- raster(xmn=-1000, xmx=1000, ymn=-1000, ymx=1000, 
      resolution=100,
      crs="+proj=utm +zone=11 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs")
values(habitat) <- exp(2+rnorm(ncell(habitat)))
habitat <- migrate_raster( habitat, kern="gaussian", sigma=100, radius=1500 )
plot(habitat)

The total carrying capacity is r sum(values(habitat),na.rm=TRUE).

# plotting function for typical layers below
pl <- function (x,v,nr=1,zlim, ...) {
    if (!missing(v)) {
        if (nlayers(x)<NCOL(v)) {
            x <- do.call(stack,list(x)[rep(1,NCOL(v))])
        }
        values(x)[!is.na(values(x))] <- as.numeric(v)
    }
    if (missing(zlim)) {
        if (inherits(x,"Raster")) {
            zlim  <- range(0,values(x),finite=TRUE)
        } else if (inherits(x,"population")) {
            zlim  <- range(0,as.numeric(x$N),finite=TRUE)
        }
    }
    plot(x,nr=nr,zlim=zlim,...)
}

Genotype-dependent recruitment

Let's suppose there's strong, additive selection: heterozygous Aa alleles are 10% more likely to germinate than aa alleles, and AA alleles are 20% more likely to germinate. The germination_fun gets called with the current state of the population as an argument (N), whose result is multiplied by the total seed flux, per genotype and per location, to get the mean number of newly recruited individuals.

Here's the demographic set-up

germination_fun <- function (N, carrying.capacity, ...) {
    out <- r0 / ( 1 + migrate(rowSums(N),competition)/carrying.capacity )
    return( cbind( aa=out, aA=s*out, AA=s^2*out ) )
}

this.demography <- demography(
        prob.seed = 0.2,
        fecundity = 100,
        prob.germination = vital( germination_fun, 
                r0 = 0.01,  # one in 100 seeds will germinate at low densities
                s = 1.5,    # multiplicative selective benefit of the A allele
                competition = migration(
                                     kern="gaussian",
                                     sigma=100,
                                     radius=300,
                                     normalize=1
                                 )
                ),
        prob.survival = 0.9,
        pollen.migration = migration(
                            kern = function (x) { exp(-sqrt(x)) },
                            sigma = 100,
                            radius = 1000,
                            normalize = NULL
                     ),
        seed.migration = migration(
                            kern = "gaussian",
                            sigma = 20,
                            radius = 400,
                            normalize = 1
                     ),
        genotypes = c("aa","aA","AA"),
        mating = mating_tensor( c("aa","aA","AA") )
    )

We'll start the population with only a few of the advantageous alleles:

total.habitat <- sum(values(habitat),na.rm=TRUE)
pop <- population( 
                  habitat = habitat,
                  genotypes = c("aa","aA","AA"),
                  N = cbind( aa=rpois_raster(habitat,only.values=TRUE),
                             aA=rpois_raster(habitat*100/total.habitat,only.values=TRUE),
                             AA=0 )
             )
this.demography <- setup_demography( this.demography, pop )

plot(pop)

Let's check this looks to allow a stable population. The mean number of offspring per individual the probability of seeding multiplied by the fecundity, multiplied by the probability of germination; this must be greater than the probability of death, calculated at low density. Here this is:

exp.pop <- pop
exp.pop$N[] <- 1
vitals <- generation(exp.pop,this.demography,
                     carrying.capacity=values(pop$habitat)[pop$habitable],
                     expected=TRUE,return.everything=TRUE)
for (k in seq_along(vitals)) {
    pl(exp.pop$habitat,v=vitals[[k]],main=c("",paste(names(vitals)[k],"density 1"),""))
}

Just for comparison, here's the number of germinating seeds per capita at a much higher density

exp.pop$N[] <- 500
vitals <- generation(exp.pop,this.demography,
                     carrying.capacity=values(pop$habitat)[pop$habitable],
                     expected=TRUE,return.everything=TRUE)
pl(exp.pop$habitat,v=vitals$germination/exp.pop$N,main=c("","per capita germination at density 500","") )

Now, we can run the simulation. Note that we have to pass carrying.capacity in. We'll record the full state of the population at several times, and record the total number of each genotype in every generation.

plot.times <- seq(0,200,length.out=51)
sim <- simulate_pop( pop, this.demography, times=plot.times, 
                carrying.capacity=values(pop$habitat)[pop$habitable],
                summaries=list( totals=function(N){colSums(N)} )
            )

Here is a video of the simulation:

for (k in seq_along(plot.times)) {
    pl(pop$habitat,v=sim$N[,,k],main=c("",sprintf("t=%d",floor(plot.times[k]))),zlim=range(sim$N,finite=TRUE))
}

And, here are the total numbers of the three genotypes:

matplot( sim$summaries[["totals"]], type='l', xlab='time', ylab='numbers', lty=1 )
legend("topright",lty=1,col=1:3,legend=colnames(sim$summaries[["totals"]]))