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) <- sample( 100*c(1,2,5,NA), length(habitat), replace=TRUE )
plot(habitat,asp=1)

# plotting function for typical layers below
pl <- function (x,nr=1,zlim,...) {
    if (missing(zlim)) {
        if (inherits(x,"Raster")) {
            zlim  <- c(0,max(values(x),na.rm=TRUE))
        } else if (inherits(x,"population")) {
            zlim  <- c(0,max(x$N,na.rm=TRUE))
        }
    }
    plot(x,nr=nr,zlim=zlim,...)
}

A simulation step

Beginning with a population:

N <- rpois_raster( brick( habitat, habitat, habitat )/3 )
genotypes <- names(N) <- c("aa","aA","AA")
pl(N)

Here are the population parameters. To maintain population densities, we'll make the germination rate depend on the local population density. Concretely, let's use something motivated by the Beverton-Holt model, which without age structure has $n' = r_0/(1/n+1/M)$, and carrying capacity $(r_0-1)M$: if $f(x)$ is the nearby population density (possibly smoothed somehow), then the germination rate is $r_0/(1+f(x)/M(x))$.

# probability of seeding in a given year
prob.seed <- 0.2
# pollen dispersal kernel
pollen.kernel <- function (x) { exp(-sqrt(x)) }
pollen.sigma <- 100
max.pollen.dist <- 1000
# seed dispersal kernel
seed.kernel <- "gaussian"
seed.sigma <- 20
max.seed.dist <- 400
# number of seeds per individual
seed.fecundity <- 200
# genotype offspring tensor: 
#  [i,j,k] is the probability that genotypes i and j combine to make offspring genotype k
mating.tensor <- array( c(
        # aa offspring
        1,   1/2,   0, # aa
        1/2, 1/4,   0, # aA
        0,     0,   0, # AA
        # aA offspring
        0,   1/2,   1, # aa
        1/2, 1/2, 1/2, # aA
        1,   1/2,   0, # AA
        # AA offspring
        0,     0,   0, # aa
        0,   1/4, 1/2, # aA
        0,   1/2,   1  # AA
        ), dim=c(3,3,3) )
dimnames(mating.tensor) <- list( genotypes, genotypes, genotypes )
# raw germination rate
prob.germination <- 0.01
# competition kernel for germination
competition.kernel <- "gaussian"
competition.sigma <- 100
max.competition.dist <- 300
# sort of carrying capacity
carrying.capacity <- habitat
# probability of survival
prob.survival <- 0.9

Sample the number of seed-producing individuals:

M <- rbinom_raster(N,prob.seed)
pl(M)

Find the mean pollen flux:

P <- migrate_raster(N,kern=pollen.kernel,sigma=pollen.sigma,radius=max.pollen.dist)
pl(P)

Find the mean seed production field:

S <- seed_production_raster(seeders=M,pollen=P,mating=mating.tensor,
                     fecundity=seed.fecundity)
pl(S)

Disperse seeds:

SD <- migrate_raster(S,kern=seed.kernel,sigma=seed.sigma,radius=max.seed.dist,
                     normalize=TRUE)
pl(SD)

Find competition kernel: according to $r_0/(1+f(x)/M(x))$:

K <- migrate_raster(sum(N),kern=competition.kernel,sigma=competition.sigma,
                    radius=max.competition.dist,normalize=TRUE)
K <- prob.germination/(1+K/carrying.capacity)
pl(K)

Sample new individuals:

G <- rpois_raster(SD*K)
names(G) <- names(SD) # multiplication removes names
pl(G)

Kill off old individuals and add in new ones:

NN <- rbinom_raster(N,prob.survival) + G
pl(NN)

Demography set-up

Now, let's use tools from the package to set up the simulation we ran above. The function that determines current germination probabilities is left unspecified, to be passed in at runtime, because it is a raster, and we want to demographic setup to be independent of the raster being used.

germination_fun <- function (N, carrying.capacity, ...) {
    r0 / ( 1 + migrate(rowSums(N),competition)/carrying.capacity )
}

this.demography <- demography(
        prob.seed = 0.2,
        fecundity = 200,
        prob.germination = vital( germination_fun, 
                r0 = 0.01,  # one in 100 seeds will germinate at low densities
                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") )
    )

Demography application

Now, doing the above simulation step is just a matter of calling generation(). Importantly, note that we pass in carrying.capacity, since it was left unspecified in the function germination_fun above.

pop <- population( 
                  habitat = habitat,
                  genotypes = c("aa","aA","AA"),
                  N = cbind( aa=rpois_raster(habitat/4,only.values=TRUE),
                             aA=rpois_raster(habitat/2,only.values=TRUE),
                             AA=rpois_raster(habitat/4,only.values=TRUE) )
             )
this.demography <- setup_demography(this.demography,pop)

pl(pop)
for (k in 1:5) {
    pop$N <- generation(pop,this.demography,carrying.capacity=values(pop$habitat)[pop$habitable])
    pl(pop)
}