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,...)
}
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
M <- rbinom_raster(N,prob.seed)
pl(M)
P <- migrate_raster(N,kern=pollen.kernel,sigma=pollen.sigma,radius=max.pollen.dist)
pl(P)
S <- seed_production_raster(seeders=M,pollen=P,mating=mating.tensor,
fecundity=seed.fecundity)
pl(S)
SD <- migrate_raster(S,kern=seed.kernel,sigma=seed.sigma,radius=max.seed.dist,
normalize=TRUE)
pl(SD)
5. 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)
G <- rpois_raster(SD*K)
names(G) <- names(SD) # multiplication removes names
pl(G)
NN <- rbinom_raster(N,prob.survival) + G
pl(NN)
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.args <- list(
r0 = 0.01,
competition = migration(
kern="gaussian",
sigma=100,
radius=300,
normalize=1,
do.M=TRUE,
population=habitat
)
)
germination_fun <- function (N, carrying.capacity, ...) {
germination.args$r0 / ( 1 + migrate(rowSums(N),germination.args$competition)/carrying.capacity )
}
this.demography <- demography(
prob.seed = 0.2,
fecundity = 200,
prob.germination = germination_fun,
prob.survival = 0.9,
pollen.migration = migration(
kern = function (x) { exp(-sqrt(x)) },
sigma = 100,
radius = 1000,
normalize = NULL,
do.M=TRUE,
population=habitat
),
seed.migration = migration(
kern = "gaussian",
sigma = 20,
radius = 400,
normalize = 1,
do.M=TRUE,
population=habitat
),
genotypes = c("aa","aA","AA"),
mating = mating_tensor( c("aa","aA","AA") )
)
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) )
)
pl(pop)
for (k in 1:5) {
pop$N <- generation(pop,this.demography,carrying.capacity=values(pop$habitat)[pop$habitable])
pl(pop)
}