In petrelharp/landsim: Simulate Populations on Landscapes

fig.dim <- 8
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)

Set-up

Here is the habitat we will work with. Note that units are in meters, and the resolution of the raster is 100m.

pop <- make_population(
            habitat = random_habitat(),
            inaccessible.value = NA,
            uninhabitable.value = NA,
            genotypes = c("aa","aA","AA"),
            N = 0
        )
pop$N[,"aa"] <- rpois(nrow(pop$N),values(pop$habitat)[pop$habitable]/4)
pop$N[,"aA"] <- rpois(nrow(pop$N),values(pop$habitat)[pop$habitable]/2)
pop$N[,"AA"] <- rpois(nrow(pop$N),values(pop$habitat)[pop$habitable]/4)
plot(pop$habitat)

Here is the basic, default demography, modified to have population regulation through germination probabilities:

basic.migr <- migration(
                    kern = "gaussian",
                    sigma = 300,
                    radius = 1000,
                    normalize = 1
             )
basic.demog <- demography(
        prob.seed = 0.05,
        fecundity = 200,
        prob.germination = 0.4,
        prob.survival = 0.6,
        pollen.migration = basic.migr,
        seed.migration = basic.migr,
        genotypes = c("aa","aA","AA")
    )

demog <- basic.demog
demog$prob.germination <- vital(
                    function (N,...) {
                        out <- r0 / ( 1 + migrate(competition,x=rowSums(N))/K )
                        cbind( aa=out, aA=out, AA=out )
                    },
                    r0 = 0.4,
                    K = values(pop$habitat)[pop$habitable]/5,
                    competition = migration(
                                kern="gaussian",
                                sigma=200,
                                radius=400,
                                normalize=1
                        )
                )

We'll start things off by running the model to stationarity, as seen by the total numbers of the three genotypes:

demog <- setup_demography( demog, pop )
sim <- simulate_pop( pop, demog, times=seq(0,100,length.out=101),
                summaries=list( totals=function(N){colSums(N)} )
            )
matplot(sim$summaries[["totals"]],type='l',lty=1, log='y', ylab="number of individuals")
legend("bottomright",lty=1,col=1:3,legend=pop$genotypes)

pop$N[] <- sim$N[,,101]

Lineages

Now we have a population history and the demographic model that produced it. Let's start a bit earlier. We can ask the generation function to return everything, not just total numbers, namely:

1. seeders : number of individuals in each location that produce seeds
2. pollen : expected pollen flux (this is pollen.migration applied to N; since this is in expected numbers we can treat pollen as diploid, averaging over meiosis later)
3. seed.production : expected numbers of seeds of each genotype produced by mating between pollen and local genotypes (using the mating tensor)
4. seeds.dispersed : seed.production with seed.migration applied to it (still expected numbers)
5. germination : numbers of new individuals germinating in each location
6. death : number of previously present individuals that have died

The numbers of genotypes (N) in the next generation are survivors + germination, where survivors = N - death.

Lineages each are labeled by an allele, not a genotype; to keep things simple, we will follow lineages of a single allele; for the discussion below this will be allele $A$. Suppose there genotype proportions among the seed parents are $p_{AA}$, $p_{Aa}$, and $p_{aa}$; they are $q_{AA}$, $q_{Aa}$, and $q_{aa}$ among the pollen parents, the probability that parents $u$ and $v$ produce offspring genotype $g$ is $t_{u,v,g}$. Lineages in $AA$ offspring choose pollen or seed parent with equal probability (1/2 each). The numbers of $Aa$ offspring coming from (seed,pollen) parents $(u,v)$ is proportional to $p_u q_v t(u,v,Aa)$; an $A$ lineage in an $Aa$ offspring will choose a pair of parents $(u,v)$ with probability proportional to this. Let $A(u)$ be the number of $A$ alleles in genotype $u$. Assuming that meiosis is fair, given parents $(u,v)$ the $A$ lineage will choose parent $u$ with probability $A(u)/(A(u)+A(v))$.

Therefore, to sample lineage movement backwards through this we just need to:

1. Choose genotypes for the lineages to be in.
2. Stay in the same place with probability survivors/(survivors+germination).
3. Otherwise, do a backwards migration step with seed.migration weighted by seed.production, then
4. choose whether to follow a seed or pollen parent, and their genotype, by a backwards step through the mating tensor as described above, using seeders and pollen as weights.
5. If the seed parent is chosen, stay put.
6. If the pollen parent is chosen, do another backwards migration step with pollen.migration weighted by N.

These steps are implemented in lineage_generation; here we apply them starting with a sample of lineages chosen uniformly from the possible $A$ gametes to make sure that things make sense.

one.gen <- generation( pop, demog, return.everything=TRUE )

# current generation
now.N <- pop$N - one.gen$death + one.gen$germination

# choose lineage locations: coordinates are (genotype, location index)
mean.n.lineages <- 1000
n.lineages <- rpois_matrix( sweep(now.N,2,c(0,1,2),"*")*mean.n.lineages/sum(now.N) )
lineages <- arrayInd( rep( seq_along(n.lineages), n.lineages ), .dim=dim(now.N) )
colnames(lineages) <- c("location","genotype")

parents <- lineage_generation( lineages, N=pop$N, gen=one.gen, num.alleles=c(0,1,2), demog=demog )

now.locs <- xyFromCell(pop$habitat,which(pop$habitable)[lineages[,1]])
parent.locs <- xyFromCell(pop$habitat,which(pop$habitable)[parents[,1]])

Here are single-generation steps for a bunch of lineages:

plot(pop$habitat)
suppressWarnings( arrows( x0=now.locs[,'x'], x1=parent.locs[,'x'], 
           y0=now.locs[,'y'], y1=parent.locs[,'y'], 
           length=0.05 ) )

Now to do this across many generations we need to keep track of everything across those generations. Here are a few lineages for 100 generations:

many.gens <- simulate_pop( pop, demog, times=1:100, tinit=1, return.everything=TRUE )
mean.n.lineages <- 20
now.N <- many.gens$N[,,length(many.gens$times)]
n.lineages <- rpois_matrix( sweep(now.N,2,c(0,1,2),"*")*mean.n.lineages/sum(now.N) )
lineages <- arrayInd( rep( seq_along(n.lineages), n.lineages ), .dim=dim(now.N) )
colnames(lineages) <- c("location","genotype")

lin.hist <- simulate_lineages( lineages, many.gens, num.alleles=c(0,1,2), demog=demog )

plot_lineages( lin.hist, pop )

Note that these lineages are all of type A, but are jumping back and forth between Aa and AA genotypes; we're not looking at that.

Lineages in a traveling wave

Just for fun, let's do the same thing backwards in a traveling wave. Let's suppose there's strong, additive selection: heterozygous Aa and AA alleles are 20% more likely to germinate than aa alleles. 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$prob.germination <- vital(
                    function (N,...) {
                        out <- r0 / ( 1 + migrate(competition,x=rowSums(N))/K )
                        cbind( aa=out, aA=1.2*out, AA=1.2*out )
                    },
                    r0 = 0.4,
                    K = values(pop$habitat)[pop$habitable]/5,
                    competition = migration(
                                kern="gaussian",
                                sigma=200,
                                radius=400,
                                normalize=1
                        )
                )

We'll start the population with only a few of the advantageous alleles, by first running it to stationarity with only aa individuals; then adding a few aAs.

demog <- setup_demography( demog, pop )
pop$N[,c("aA","AA")] <- 0
sim <- simulate_pop( pop, demog, times=seq(0,100),
                summaries=list( totals=function(N){colSums(N)} )
            )
matplot(sim$summaries[["totals"]],type='l',lty=1, log='y', ylab="number of individuals", main='warmup')
legend("bottomright",lty=1,col=1:3,legend=pop$genotypes)

start.N <- sim$N[,,dim(sim$N)[3]]

# try until establishment
for (ntries in 1:200) {
    pop$N[] <- start.N
    wave.start <- sample.int(nrow(pop$N),1)
    pop$N[wave.start,c("aa","aA")] <- pop$N[wave.start,c("aA","aa")]  
    sim <- simulate_pop( pop, demog, times=seq(0,600),
                 summaries = list( totals=function(N){colSums(N)} ),
                 stop.fun=function(N){ sum(N)==0 },
                 return.everything=TRUE )
    total.aA <- sim$summaries[[1]][nrow(sim$summaries[[1]]),"aA"]
    if ( total.aA>0 ) { break }
}

matplot(sim$summaries[["totals"]],type='l',lty=1, log='y', ylab="number of individuals", main='sweep')
legend("bottomright",lty=1,col=1:3,legend=pop$genotypes)

Here are a few lineages for those 600 generations; the origin of the wave is shown with a black dot.

mean.n.lineages <- 20
now.N <- sim$N[,,length(sim$times)]
nA <- sweep(now.N,2,c(0,1,2),"*")
lineages <- arrayInd( sample( seq_along(now.N), size=mean.n.lineages, prob=nA, replace=TRUE ), .dim=dim(now.N) )
colnames(lineages) <- c("location","genotype")

lin.hist <- simulate_lineages( lineages, sim, num.alleles=c(0,1,2), demog=demog )

plot_lineages( lin.hist, pop )
points( xyFromCell(pop$habitat, which(pop$habitable)[wave.start]), cex=2, pch=20 )

petrelharp/landsim documentation built on May 25, 2019, 2:53 a.m.