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 randomly generated landscape; the white areas (which are NA) are inaccessible areas. Note that units are in meters.

diam <- 1e4
habitat <- raster(xmn=-diam, xmx=diam, ymn=-diam, ymx=diam, 
      resolution=100,
      crs="+proj=utm +zone=11 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs")
values(habitat) <- pmin(20,(2+rcauchy(ncell(habitat))))
habitat <- 20*migrate_raster( habitat, kern="gaussian", sigma=300, radius=1500 )
values(habitat)[values(habitat)<0] <- NA
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)
        names(x) <- colnames(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

demog <- demography(
        prob.seed = 0.2,
        fecundity = 100,
        prob.germination = germination_fun <- vital(
                function (N, carrying.capacity, ...) {
                    out <- r0 / ( 1 + migrate(rowSums(N),competition)/carrying.capacity )
                    return( cbind( aa=out, aA=s*out, AA=s^2*out ) )
                },
                r0 = 0.01,  # one in ten 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 = 30,
                            radius = 500,
                            normalize = NULL
                     ),
        seed.migration = migration(
                            kern = "gaussian",
                            sigma = 50,
                            radius = 400,
                            normalize = 1
                     ),
        genotypes = 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*40/total.habitat,only.values=TRUE),
                             AA=0 )
             )
demog <- setup_demography( demog, pop )

pl(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 is the number of offspring per individual at low density, minus one:

base.r <- intrinsic_growth( pop, demog, 
                         carrying.capacity=values(pop$habitat)[pop$habitable] )
pl(pop$habitat, v=base.r-1)

To check for stability, here's the same thing, at density 500:

base.r <- intrinsic_growth( pop, demog, density=500,
                         carrying.capacity=values(pop$habitat)[pop$habitable] )
pl(pop$habitat, v=base.r-1)

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,700,length.out=71)
sim <- simulate_pop( pop, demog, 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",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"]]))