Effects of population- and seed bank noise on neutral evolution and efficacy of natural selection
Lukas Heinrich Johannes Müller Aurélien Tellier Daniel Zivković
q-bio.PE 92D10, 39A50, 60H10







TUMGarching]Lukas Heinrich 
TUMGarching,ICB]Johannes Müllermycorrespondingauthor 
TUMWeihen]Aurélien Tellier 
TUMWeihen]Daniel Živković 


[TUMGarching]Center for Mathematics, Technische Universität München, 85748 Garching, Germany
[TUMWeihen]Section of Population Genetics, Center of Life and Food Sciences Weihenstephan, Technische Universität München, 85354 Freising, Germany
[ICB]Institute for Computational Biology, Helmholtz Center Munich, 85764 Neuherberg, Germany


[mycorrespondingauthor]Corresponding author








Population genetics models typically consider a fixed population size and 
a unique selection coefficient. However, population dynamics inherently 
generate noise in numbers of individuals and selection acts on various 
components of the individuals' fitness. In plant species with seed banks, 
the size of both the above- and below-ground compartments present noise 
depending on seed production and the state of the seed bank. We investigate 
if this noise has consequences on 1) the rate of genetic drift, 
and 2) the efficacy of selection. We consider four variants of two-allele 
Moran-type models defined by combinations of presence and absence of noise 
in above-ground and seed bank compartments. Time scale analysis and dimension 
reduction methods allow us to reduce the corresponding Fokker-Planck equation 
to a one-dimensional diffusion approximation of a Moran model. 
We first show that if 
the above-ground noise classically affects the rate of genetic drift, 
below-ground noise reduces the diversity storage effect of the seed bank. 
Second, we consider that selection can act on four different components 
of the plant fitness: plant or seed death rate, seed production or seed 
germination. Our striking result is that the efficacy of selection for seed 
death rate or germination rate is reduced by seed bank noise, whereas 
selection occurring on plant death rate or seed production is not affected. 
We derive the expected site-frequency spectrum reflecting this heterogeneity 
in selection efficacy between genes underpinning different plant fitness 
components. Our results highlight the importance to consider the effect 
of ecological noise to predict the impact of seed banks on neutral and selective 
evolution.




Diffusion Moran model seed bank selection site-frequency spectrum 
[2010] 92D10, 39A50,  60H10








§ INTRODUCTION


Genetic drift and natural selection are prominent forces shaping the amount of genetic 
diversity in populations. In diploid dioecious organisms, natural selection can be 
decomposed in different components: 1) viability selection as the differential 
survival of the genotypes from zygotes to adults, 2) fecundity (or fertility) 
selection as the differential zygote production, 3) sexual selection as the differential 
success of the genotypes at mating, and 4) gametic selection as the distorted segregation 
in heterozygotes <cit.>.


In effect, population genetic models with 
discrete generations or with Malthusian parameter ignoring age-structure lump these 
components into one unique selection parameter. Experimental or genomic population 
studies describing thus changes in allele frequencies often fail to describe 
and dissect the respective effects of these selective modes. Several theoretical 
studies on fertility <cit.> 

or on sexual selection <cit.> 

as well as experimental work on animal <cit.>


and plant <cit.>

populations have attempted 
to disentangle the respective influence of these selection components.



In age-structured populations, however, genetic drift and selection 
can act differently than predicted by models without age structure 
(overview in the book by ).

The first type of age-structured models are simply obtained 
by individuals' life span and reproduction overlapping several 
generations. Here, genetic drift acts equally on all individuals 
of all age classes at any generation. The magnitude of genetic drift 
is defined by the population size which can be fixed, or fluctuating 
following a logistic dynamic size constrained by the population carrying 
capacity. As an extreme type of overlapping generation model, the Moran 
model exhibits a rate of random genetic drift all but the same as 
the Wright-Fisher (WF) model up to a rescaling of the population 
size (see, e.g., the book by ). Meanwhile, in age-structured 
populations, selection for fecundity and 
viability can show different outcomes, such as time to allele fixation and/or maintenance 
of alleles, compared to discrete models. This occurs if fecundity or longevity act 
at different ages of the structured population and under a logistic population size 
dynamic <cit.>. 


An interesting question 
arising from the current increasing availability of genomic data is whether 
in age-structured populations selection for fecundity can be disentangled from that 
of viability using population genomics statistics (such as the site-frequency 
spectrum, SFS).



A second type of age-structured model is obtained when considering that individuals may 
remain as dormant/quiescent structures spanning several generations. Quiescence 
in reproductive structure is in fact wide-spread, such that seeds or eggs can be persistent 
states that allow to buffer a variable environment <cit.>. 

The time at which offspring germinates or hedges can be variable, such that 
only some of the offspring live in detrimental periods, and most likely at least some in 
a beneficial environment. Bacterial spores or lysogenic states of temperate phages are also 
examples of such a bet-hedging strategy. Seedbanks represent thus a storage of genetic 
diversity decreasing the probability of population extinction <cit.>, 

diminishing the effect of genetic drift <cit.>, 

slowing down the 
action of natural selection <cit.> 


and favouring balancing selection <cit.>.

The strength of the seed bank effect clearly depends on the organism under consideration. 
Dormant seeds and diapausing eggs have a somehow short live span compared to 
the average coalescent time, 
and we call these seed banks “weak”. Dormant states of bacteria, however, can last 
a longer time (many generations) even than the average coalescent time. These seed banks 
are called “strong”, and modelled in  a similar way as weak selection: the time scale 
of the quiescent state is scaled by the inverse of the population size. In particular 
<cit.> investigated in a series of papers strong 
seed banks, and find mathematically appealing results as deviations from the Kingman 
coalescent.


In the present paper, we focus on weak seed banks aiming at applications to plant 
or invertebrate species. In a seminal paper <cit.> 
investigated the effect of weak seed banks on the coalescent. 
The key parameter here is G that denotes the average number of plant generations 
a seed rests in the soil. 
<cit.>  
have obtained 
the Kingman n-coalescent rescaled by (1+G)^-2 (the nomenclature in the paper 
of <cit.> is somewhat different, as they consider non-overlapping 
generations in a discrete time setting, while we formulate
the result already for a time continuous model with overlapping generations). 
It was subsequently shown 
that 
1+G can be estimated using polymorphism data and information on the census size 
of populations <cit.>. 

Furthermore, neglecting seed banks may yield 
distorted results for the inference of past demography using for example the 
SFS <cit.>. Interestingly, the effect 
of (weak) selection is enhanced by the slow-down of the time scale due to seed 
banks <cit.>. The effect of genetic drift 
and 
weak natural selection on allele frequencies can be computed in 
a diffusion framework 
in a Moran model with deterministic seed bank  <cit.>. While 
the diffusion term, defining genetic drift is scaled by (1+G)^-2, matching 
the backward coalescent result of <cit.>, the coefficient 
of natural selection in the drift term is multiplied only by (1+G)^-1. 
In biological 
terms, this means that the strength of selection, as defined by a unique 
selective 
coefficient, is enhanced by the seed bank compared to the effect 
of genetic drift, 
even though the time to reach fixation for an allele is 
increased <cit.>.
We investigate two additions to this current body of theoretical literature on weak seed 
banks. 
First, we compute the effect of realistic models relaxing the hypotheses of a fixed size 
for the population above ground, and of deterministically large for the seed bank 
compartment. By doing so, we generate noise in the population size above or below ground. 
This extends the classic Moran or Wright-Fisher models, as the importance of the noise 
effect on population dynamics on the rate of genetic drift or selection is being 
recognized <cit.>. However, this approach is faced with an additional 
difficulty: if two alleles are present in a population of a constant size, 
it is sufficient to keep track of the number of individuals for one allele only. 
In the case of a fluctuating population size with logistic dynamics constrained 
by the carrying capacity, the state space becomes essentially two-dimensional. 
This more general case can still be examined via a time scale analysis and a dimension 
reduction by singular perturbation approaches. These methods are applicable 
in the case where the population size becomes 
large <cit.>. Though in those papers 
(as in the present one) the arguments are used in a formal way, the validity 
of this approach is proven <cit.>. This latter analysis reveals 
that the dynamics are well described by appropriate scaled 
diffusion approximation of a Moran models.
Second, we dissect the fitness of plants into four components, which can be possibly 
affected by genetic drift occurring above ground and in the seed bank. We compute 
classic population genetics results for neutral and selected alleles and derive 
the expected SFS for the alleles under different fitness components. 
The analysis 
reveals the effect of noise in the above-ground population (plants) as well as 
in the below-ground population (seeds). In particular, we find that genetic drift 
and selection are differently influenced by above and below-ground noise. Additionally, 
the selection coefficient of alleles involved in seed death rate or germination rate 
are reduced by seed bank noise, but not by above-ground noise, while the two others 
are not affected (plant death rate and seed production). Below-ground noise also 
reduces the seed bank storage effect of neutral genetic diversity.



§ METHODS

The aim of the present study is to investigate 
the effect of above-ground/below-ground population noise on evolution 
in presence of seed banks. We extend the results obtained 
in <cit.>, 
where a seedbank model with deterministic seedbank, fixed 
above ground population size and weak selection in seed survival 
and seed production 
has been considered. In the present paper, we consider more sources 
of variability and selection also in plant death rate and seed germination, 
and compare the results of four different models 
(Figure <ref>):  fluctuations/no fluctuations in the total
above ground population, intrinsic fluctuations/no fluctuations in the seed bank.

Noise in population size can be located either 
above ground (plant population), or below ground (seed bank or seed population). 
Traditionally, the plant above-ground population is assumed
to consist of N individuals, N being fixed. A noisy substitute of this 
assumption is a logistic (fluctuating) population 
model <cit.>. 



For seeds, we find in 
the literature
either the assumption of a finite, fixed number <cit.> of seeds,
or, a deterministic infinite seed density <cit.>. 
In the latter, the assumption is that the number of seeds per plant 
is large enough, 
such that intrinsic noise is negligible. We suggest, here, an alternative 
model, 
in which each plant produces single seeds at time points that are 
distributed according 
to a Poisson process. The seeds also die (or loose their ability
to germinate) after a random time. 
Throughout the paper, we only consider Markovian processes, that is, 
all waiting times are exponentially distributed. 
The “fluctuating  seedbank” assumption assumes 
stochastically varying seed bank size. We therefore consider four models 
defined 
by all possible combinations of fixed or logistic above-ground population 
and deterministic or fluctuating seed bank. 




Notation: 
X_1,t denotes the number of allele-A plants. In a model with fixed above-ground 
population size N, the number of allele-B plants reads X_2,t=N-X_1,t; 
in the logistic fluctuating version, we set up a stochastic process for X_1,t 
and X_2,t. In that case, N does denote the carrying capacity of the population 
(i.e., the maximal possible size). The average population size (conditioned on non-extinction) 
is κ N for some κ∈(0,1). 
Y_t (Z_t) always refers to the amount of allele A (allele B) seeds in the 
bank. For the deterministic seed bank, Y_t (Z_t) are real 
numbers that follow (conditioned on X_1,t, X_2,t) an ordinary differential 
equation (ODE). In the case of stochastically fluctuating seed banks, 
Y_t (Z_t) are non-negative integers, that follow a stochastic 
birth-death process.

We allow for weak natural selection. 
The rates for allele B individuals slightly differ from those for allele 
A individuals on a scale of 1/N, as it is usual for weak 
effects. If σ_i>0, allele B has a disadvantage in comparison 
with allele A in the respective process (where σ_1
addresses the death of a plant, σ_2 the death of seeds, 
σ_3 the production of seeds, and σ_4 the germination of seeds). Of course, the signs 
of σ_i can be chosen in an arbitrary way to consider 
a genotype B that has an advantage above genotype A (σ_i<0), 
or for example a situation where B has a disadvantage 
above ground (σ_1, σ_3>0) and an 
advantage below ground (σ_2, σ_4<0). 
The parameters of our models are summarized in Table <ref>.








meaning     symbol (Allele A)     symbol (Allele B)

death rate of  plants     ζ     ζ (1+σ_1/N)

death rate of seeds     μ     μ (1+σ_2/N)

production rate of seeds     β     β (1-σ_3/N)

germination rate of seeds (log. pop. only)     γ    γ (1-σ_4/N)




tableParameters of the models.

















 §.§ Fixed population size and deterministic seed bank


This model has been developed and discussed before <cit.>. The total 
above-ground population has size N, and the transitions for the 
allele-A plant population X_1,t∈{0,…,N} given the 
seed densities Y_t and Z_t∈_+ are presented in Table <ref>. 
E.g., an A-plant dies at rate ζ X_1,t. 
In the standard Moran model, it is instantaneously 
replaced by a B-individual with 
probability 1-X_t,1/N. In our case, the seeds determine the probability 
for a B-individual, where this probability is given by Z_t/(Y_t+Z_t). In the same way we obtain 
the rate at which a B-individual dies (the rate is ζ (1+σ_1/N) (N-X_1,t)) and is replaced by 
an A-individual (the probability is Y_t/(Y_t+Z_t)). 






event     offset    rate

death of A, birth of B   
X_1,t↦ X_1,t-1    
ζ X_1,t  Z_t/(Y_t+Z_t)

death of B, birth of A   
X_1,t↦ X_1,t+1    
ζ (1+σ_1/N) (N-X_1,t)  Y_t/(Y_t+Z_t)




tablePossible transitions and their rates.




In this model we assume that the number of seeds a plant produces 
is large (basically infinitely large), such that the seed density 
in the soil follows a deterministic process, given the history of 
the above ground population. The dynamics of seeds follows a Davis' piecewise 
deterministic process <cit.>, where plants produce seeds 
at rate β (resp. β(1-σ_3/N)) and seeds die at rate 
μ (resp. μ(1+σ_4/N))

    Ẏ_t = β X_1,t-μ Y_t,  Ż_t = β(1-σ_3/N) (N-X_1,t)-μ(1+σ_2/N) Z_t.




 §.§ Logistic population dynamics and deterministic seed bank


For the logistic model, we do not couple death and birth of a plant as 
it is usually done in Moran-type models to keep the total population 
size constant. We generalize the logistic dynamics as investigated, e.g., 
by <cit.>, or <cit.> 
for the situation at hand and separate birth and death events. If 
the seed densities Y_t,Z_t∈_+ are given, the transitions 
for X_1,t, X_2,t∈{0,…,N}, X_1,t+X_2,t≤ N, 
read as summarized in Table <ref>.






event     offset     rate 

death of A     X_1,t→ X_1,t-1     ζ X_1,t 

death of B     X_2,t→ X_2,t-1     (1+σ_1/N)ζ X_2,t 

birth of A     X_1,t→ X_1,t+1     γ(1-(X_1,t+X_2,t)/N)Y_t 

birth of B     X_2,t→ X_2,t+1     (1-σ_4/N)γ(1-(X_1,t+X_2,t)/N)Z_t




tablePossible transitions and their rates.




Conditioned on X_1,t and X_2,t, the dynamics of seeds are 
again deterministic, 

    Ẏ = β X_1 - μ Y,    Ż = (1-σ_3/N)β X_2 - (1+σ_2/N)μ Z.





 §.§ Fixed population size and fluctuating seed bank


For the above-ground population we return to a fixed population size, 
s.t. it is sufficient to follow X_1,t as  X_2,t=N-X_1,t. 
In the present model we address the noise in the number of seeds, Y_t,Z_t∈_0. 
The seeds follow a stochastic birth-death process, where the death rate 
is kept constant, and the birth rate is proportional to the number of 
corresponding above-ground plants. We obtain the transitions summarized in Table <ref>.







event     offset     rate 

death of A, birth of A     (X_1,t,Y_t)→ (X_1,t,Y_t-1)     ζ X_1,t Y_t/(Y_t+Z_t) 

death of A, birth of B     (X_1,t,Z_t)→ (X_1,t-1,Z_t-1)     ζ X_1,t Z_t/(Y_t+Z_t) 

death of B, birth of A     (X_1,t,Y_t)→ (X_1,t+1,Y_t-1)     (1+σ_1/N)ζ (N-X_1,t) Y_t/(Y_t+Z_t) 

death of B, birth of B     (X_1,t,Z_t)→ (X_1,t,Z_t-1)       (1+σ_1/N)ζ (N-X_1,t) Z_t/(Y_t+Z_t) 

birth of A-seed     Y_t → Y_t+1     β X_1,t 

death of A-seed     Y_t → Y_t-1     μ Y_t 

birth of B-seed     Z_t → Z_t+1     (1-σ_3/N)β (N-X_1,t) 

death of B-seed     Z_t → Z_t-1     (1+σ_2/N)μ Z_t 




tablePossible transitions and their rates.







 §.§ Logistic population dynamics and fluctuating seed bank



The last model incorporates logistic growth and a stochastically fluctuating 
seed bank. This model is an obvious combination of the last two models.






event     offset     rate

death of A     X_1,t→ X_1,t-1     ζ X_1,t 

death of B     X_2,t→ X_2,t-1     (1+σ_1/N)ζ X_2,t 

birth of A     (X_1,t,Y_t)→ (X_1,t+1,Y_t-1)     γ(1-(X_1,t+X_2,t)/N)Y_t

birth of B     (X_2,t,Z_t)→ (X_2,t+1,Z_t-1)     (1-σ_4/N)γ(1-(X_1,t+X_2,t)/N)Z_t

birth of A-seed     Y_t → Y_t+1     β X_1,t 

death of A-seed     Y_t → Y_t-1     μ Y_t 

birth of B-seed     Z_t → Z_t+1     (1-σ_3/N)β X_2,t 

death of B-seed     Z_t → Z_t-1     (1+σ_2/N)μ Z_t




tablePossible transitions and their rates.






 §.§ Strategy for the analysis of the models


The aim of the analysis is the reduction of the four models to 
a one-dimensional
diffusion approximation of a Moran model representing the fraction of allele-A individuals within 
the population. As the details of the analysis are tedious, we present them in detail in the appendix and only outline the basic idea in the present section.

The strategy of the analysis differs for the first model 
(fixed population size, deterministic seed bank) and the other three models. 
The reason is that there is one single stochastic state 
variable X_1,t in the first model, so that if we know the history of X_1,t, 
the state of the seed bank is known. No dimension 
reduction method is thus required. We basically can use the approach of <cit.> 
to derive the diffusion approximation of a Moran model. However, we outline in appendix <ref> 
a heuristic argument
based on a small-delay approximation, as this route seems to provide 
an appealing 
short-cut (though to our knowledge for this approximation, 
that is used in theoretical 
physics, no hard convergence theorem is available): if ecological 
time t is not considered but rather the evolutionary 
time τ=t/N (population size N large), the delay of 
a weak seed bank is small. In this case, 
the solution can be expanded w.r.t. the delay. 
As a result, the seed bank can be removed from 
the stochastic process and replaced by appropriately 
rescaled parameters. Since for the time beeing, the short delay approximation is only a heuristic approach, we formulate also a proof for the result based on time scale arguments (as explained next) in appendix <ref>.




The three other models have two, three or four stochastic state variables. 
Methods of dimension reduction are required to obtain a one-dimensional 
diffusion approximation of a 
Moran model. The key insight here is that any realization saddles fast 
on a one-dimensional manifold. If we consider, e.g., the deterministic version 
of the logistic population dynamics with stochastically fluctuating seed bank, 
we find, according to arguments by, e.g., <cit.>, for the deterministic
limit as N→∞ 
(with x_i(t)=X_t,i/N, y(t)=Y_t/N, z(t)=Z_t/N) 

    ẋ_̇1̇   = γ(1-x_1-x_2)y -ζ x_1 
    ẋ_̇2̇   = γ(1-x_1-x_2)z -ζ x_2 
    ẏ   = β x_1 - μ y
    ż   = β x_2 - μ z

It is straightforward  to show that a line of stable equilibria, 
the so-called coexistence line, exists:

Let ϑ := (β-ζ)/μ>0, κ := (γϑ-ζ)/(γϑ) ∈ [0,1]. Then, there is a line of stationary points in [0,κ]^2×ℝ^2_+ given by

    (x_1,x_2,y,z)=(x,κ-x,ϑ x,ϑ(κ-x)),       x∈[0,κ].

The line of stationary points is transversally stable 
(locally and globally in the positive cone).

It turns out that the stochastic process rapidly approaches  
this line of equilibria, and performs a random walk close to 
it (see Figure <ref>). The analysis reveals that 
the distribution on a transversal cut is just a normal distribution with 
a variance of O(1/N). Along the line of stationary points, however, 
the realizations will move 
according to a one-dimensional diffusion approximation of a Moran process. 
This approximative process is a combination 
of one component of the full process parallel to this line, 
and a second component that results from an interaction between a component 
perpendicular to this line with the deterministic vector field directed 
towards this line. 
In order to reveal this structure, we first 
use a large population limit (Kramers-Moyal expansion)
to obtain a Fokker-Planck/Kolmogorov 
forward equation for the full process. In the second step we apply singular 
perturbation methods as described, e.g., in <cit.> or <cit.> to perform the dimension 
reduction to the one-dimensional Moran model.




§ RESULTS




 §.§ Timescales for different seed bank models


For all of our models, the resulting one-dimensional Fokker-Planck equation 
assumes the form of a diffusion approximation of a Moran-model with weak selection, 

    ∂_τu = -σ a∂_{(1-)u} +1/2b∂_^2{(1-)u}

where a and b describe the speed of selection and genetic drift, respectively. The term 
σ represents selective coefficients and is a generic parameter for notation including σ_i, i=1,…,4. In order to formulate 
the results for a and b, let us introduce three composite parameters: 
G=ζ/μ is the number of plant generations a seed survives on average; 
Y=β/ζ is the average number of seeds produced by a plant; 
κ already defined in Proposition <ref> is the average 
fraction of the 
above-ground population size in the logistic model, 
in comparison with the maximal possible 
population size N. Furthermore, for deterministic and fluctuating seed banks, we respectively denote (1+G)^-1 and (1+(1-1/Y)G)^-1 as 𝒢, which can be seen as the number of plant generations that seeds survive on average corrected by the size of the seed bank. Using these abbreviations, the parameters a, b 
and σ for our four models are summarized in Table <ref>.






h]llll
model/scale term     a(selection)     b(genetic drift)     σ (selection coeff.)

fix. pop., det. s.b.      ζ 𝒢^-1     2ζ 𝒢^-2
    σ_1+σ_2+σ_3

fix. pop., fluct. s.b.     ζ 𝒢^-1     2ζ 𝒢^-2 
   σ_1+(1-1/Y)σ_2+σ_3

log. pop., det. s.b.      ζ 𝒢^-1     2ζ 𝒢^-2 κ^-1 
   σ_1+σ_2+σ_3+σ_4

log. pop., fluct. s.b.      ζ 𝒢^-1     2ζ 𝒢^-2 κ^-1 
    σ_1+σ_2+σ_3+(1-1/Y)σ_4




tableDrift, diffusion, and selection coefficients for the different population/seed bank models.







The basic seed bank model (fixed population size, deterministic seed bank) 
and the standard Moran model without seed bank 
can be used as reference models. 
A seed bank slows down the time scale of selection as well as that of genetic drift <cit.>, 
where selection is less affected (by a factor of (1+G)^-1) than genetic drift (by a factor of (1+G)^-2). 
We find that fluctuations in the above-ground population and in the seeds have different effects.Fluctuations in the seed number reduce the storage effect of seed banks. The additional noise yields a  
reduction of the effective time a seed spends within the seed bank, and thus increases the rate of genetic drift. For Y→∞ 
(noise in seed bank tends to zero), we obtain the result for the deterministic 
seed bank, for Y→ 1 (noise is maximized), the model converges towards 
the standard Moran model without seed bank. Note that Y 
is not the average number of seeds per plant directly measured but the effective number of seeds per plants. 
For example, a certain fraction of seeds might be getting lost due to other environmental reasons (abiotic or biotic factors) than their intrinsic mortality. 
These seeds do not contribute to the bank. The noise in the above-ground population only affects genetic drift and does not appear in the selection term. 
This result reflects that 
the actual competition between alleles A and B only happens above ground. 
Nonlinear terms in the transition rates only appear in the birth term of the plants. 
By increasing solely genetic drift, the above-ground noise can counteract the amplification of selection by seed banks.
The scaling of selection by 𝒢^-1 and that 
of  genetic drift by 𝒢^-2 is somehow expected. 
All mutations are affected in the same 
way by the above-ground noise. Our result concerning the lumped 
selection coefficient σ, however, is unexpected: 
Mutations for some fitness components (mortality of seeds, σ_2, 
and germination ability, σ_4) show reduced selection while this is 
not the case for selective coefficients of other fitness components 
(σ_1 and  σ_3). 
This means that if the number of seeds per plant is not too large, 
beneficial mutations in the mortality of seeds (σ_2) have 
a reduced chance to reach fixation compared with 
a beneficial mutation for the production of seeds (σ_3). 




 §.§ Site-frequency spectrum (SFS)


The SFS is a commonly used statistic for the analysis of genomewide distributed SNPs. It is defined as the distribution of the number of times a mutation is observed in a population or a sample of n sequences conditional on segregation. 
Herein, this distribution is taken over numerous unlinked sites and mutations occur only on previously monomorphic ones <cit.> at rate θ/2 per N generations.  Each mutant allele that arises from the wildtype at such an independent site is assumed to marginally follow the diffusion model specified in (<ref>) so that in particular all mutants have equivalent selective effects among sites. Mutations are allowed to occur in plants and seeds, but the following results can be easily adapted to the scenario, where mutations may only arise in plants. The proportion of sites at equilibrium, where the mutant frequency is in (y, y + dy), is routinely obtained as 
(e.g., )

    f̂(y)=θ/b  y (1-y)1-exp{-2 a/b σ (1-y)}/1-exp{-2 a/b σ},

where a, b and σ≠0 are given in Table <ref>. The sample SFS 
at equilibrium can be immediately obtained from (<ref>) via binomial 
sampling as 

    f̂_n,k=θn/b  k(n-k)1-_1F_1(k;n;2 a/b σ)e^-2 a/b σ/1-e^-2 a/b σ,

where _1F_1 denotes the confluent hypergeometric function of the first 
kind <cit.>. The neutral versions 
of (<ref>) and (<ref>) are respectively given 
by f̂(y)=θ/(b y) and f̂_n,k=θ/(b k). 

Fig. <ref> indicates the striking effects of the seed bank noise 
for the model with logistic population dynamics and fluctuating seed banks. 
When the noise is low (Y high) the number of segregating sites is expected 
to be high, and mutations involved in the selection coefficients σ_1 
and σ_4 lead to similar SFS. The SFS show the typical U-shape expected 
under positive pervasive selection (left panel). If, however, Y becomes small, the number 
of segregating sites decreases, and the selection coefficient  σ_1 shows 
a U-shaped SFS while mutations under σ_4 do not (right panel).






§ DISCUSSION


We considered four models to investigate the effect of combined noise 
in the above-ground population and below-ground seed bank. In all
four cases time scale arguments allow us a reduction to 
a diffusion approximation of a Moran model. Our results extend 
the findings that a seed bank without noise yields a change of 
time scale in selection and genetic
drift that amplifies the effect of weak selection <cit.>. 
The first main result of this present work is that there is no direct 
interaction between the noise above ground and below ground. 
The above-ground noise increases the effect of genetic drift 
compared to a fixed above-ground population size, but does 
not affect selection. One can propose a new definition of 
the effective population size N_e =  κ 𝒢^2 N/2 
describing the change in genetic drift due to this noise 
(according to <cit.>. 
This allows us to redefine the evolutionary time scale 
T = ζ t/N_e. 
As a result of this procedure the parameter κ, 
which represents the reduction of the average 
total population size and indicates the increase of the above-ground 
noise by the logistic population dynamics, appears as a factor 
in front of the selection 
term in the Moran model. Note that in this notation
we find that the terms for below-ground 
(𝒢^2) and above-ground (κ) noise 
are multiplied indicating their independence. 

The second result is that the below-ground noise affects the scaling of time 
for both the selection and the 
genetic drift term. We introduce the concept of the mean effective number of 
seeds per plant in a similar definition to the effective population 
size <cit.>. 
If this average seed number tends to infinity, we recover the effect of 
a deterministic seed bank (no noise), and if this number tends to one, 
the noise of the seed bank accelerates the time scale such that 
the seed bank has no effect at all. The magnitude of noise in the seed bank is 
thus tuned between the two extreme cases of “no seed bank” (minimal effect) 
and “deterministic seed bank” (maximal effect). 


The third result of the present study is 
the insight that below-ground noise may affect the four plants' fitness components 
of viability and fecundity differentially. If the 
below-ground noise is large, the selective effect of 
mutations involved in seed death and germination may even be cancelled out. 
Other mutations involved in fitness traits of the above-ground population, such 
as plant death or seed production are not affected. In other words, while finite 
size and noise above ground do not affect selection, noise and finite size of the 
seed bank do change the selection coefficients.
In biological terms we interpret this result as follows. 
Above ground, the fate of an allele under positive selection 
is classically determined by the strength of genetic drift, which 
depends on the population size (the diffusion term in the Fokker-Planck equation) 
compared to the strength of selection which depends on the selection 
coefficient (the drift term in the forward model). Under a deterministic 
seed bank (as in the present models and in <cit.>), selection 
is efficient because it occurs on plants, when they are above ground, with a probability (1+G)^-1 and genetic drift occurs on a coalescent scale of (1+G)^-2. 
Any change in the allele frequencies above ground translates directly into the deterministic 
seed bank, just with a small time delay. However, when the seed bank has 
a fluctuating finite size, the strength of selection on seed fitness (the 
coefficients σ_2 and σ_4) is decreased by the noise in the 
seed compartment. This occurs because selection and genetic drift in the bank 
occur at every generation and not only when seeds germinate.

This observation implies that competition experiments, which measure and 
compare the effect of mutations and/or allow to determine selection coefficients 
on the ecological time scale can hardly be used to extrapolate to the evolutionary 
time scale. Indeed, even if two mutations seem to have an equivalent value for 
the plant fitness in a competition experiment, the presence of seed bank noise may 
lead to different evolutionary outcomes.

We finally discuss two ways of using genome polymorphism data to estimate seed bank 
parameters as well as the selection coefficients. 
We can attempt to infer the seed bank parameters based on neutral genetic diversity under the idealized conditions that the sample SFS reached an equilibrium as f̂_n,k (i.e., there is no recent demographic impact) and that we can measure or estimate the population mutation rate θ and the death rate of plants ζ. It turns out that 𝒢 (the number of plant generations a seed survives on average) and κ (the average above-ground 
population size) can then be identified, so that it is possible to disentangle the effects 
of the seed bank from those of the above-ground noise on neutral evolution by utilizing f̂_n,k. However, 
we cannot identify the below-ground noise, as G 
and Y only appear in the composite parameter 𝒢. Note that if we only have information about the relative 
sample SFS f̂_n,k/∑_if̂_n,i, the multiplicative 
constant in f_n,· cancels out, and therefore only κ 𝒢 and the combined effect of above- and below-ground noise can be estimated. 
Extending the work by <cit.>, 

we suggest that the number of plant 
generations a seed survives on average can be estimated from the absolute SFS using for example a Bayesian inference method with 
priors on the census size of the above-ground population and the death rate of plants.

We may also aim to infer the selection coefficients underpinning the various 
plants' fitness components, namely the fecundity and viability of plants or seeds. 
With the abundance of gene expression, molecular and Gene Ontology data, it becomes feasible 
to group genes by categories of function or pathways, for example to know the genes involved in seed 
germination or seed integrity, viability and seed dormancy (e.g., ). These functional groups of genes actually underlie the different plants' fitness components investigated in this study. 

Using genome-wide polymorphism data of several individuals, the SFS for each functional groups of genes can be computed and used to draw inference of selection, for example as the distribution of selective effects (e.g., the method by ). Our prediction is thus that the SFS would reflect the differential selection on these 
fitness components and can be observed over the genes involved in the different functional groups. The limitation in current data lies so far on the functional side, as more gene expression study are needed to assigne genes to functional network and to different plant fitness components.
As an extension, our expected SFS shows that the behavior of the seed fitness coefficients can be affected 
by the below-ground noise. So, we predict that populations with a small sized seed bank 
should exhibit less selection signatures on genes related 
to seed fitness compared to populations with a larger seed bank compartment. The same analysis as above can be conducted, but now comparing the SFS and inferred selection for the different functional groups of genes across several populations with know ecological set-ups which define the seed bank size. We further suggest that such procedure can be applied to disentangle selection for fecundity from that for viability in age-structured populations, if functional groups of genes can be assigned to these traits (for example using gene expression at different life stages and ages).

Our results differ thus from those of classic age-structured populations by the overlap 
of generations, as the seed bank can present its own rate of genetic drift. Moreover, 
selection acts differently above ground and below ground on the different plants' fitness 
components, which may allow us to disentangle their effect on the overall selection coefficient.



§ ACKNOWLEDGEMENTS

This research is supported in part by Deutsche Forschungsgemeinschaft grants TE 809/1 (AT), STE 325/14 from the Priority Program 1590 (DZ), 
and MU2339/2-2 from the Priority Program 1617 (JM). 




§ REFERENCES

elsarticle-harv















§ ANALYSIS OF THE MODELS



We present the conceptions for the analysis of the model of a
fixed population size and a deterministic seed bank as well as for 
logistic population dynamics and a deterministic seed bank in detail. Note that the computations for the remaining two models  
are similar but even more extensive, so that we do not present them in full length but only mention the results of the 
main steps. The computations for the dimension reduction have been checked using 
the computer algebra package MAXIMA <cit.> (see the available supplementary files).



 §.§ Fixed population size and deterministic seed bank


To keep the demonstration short, we present a nice  
but heuristic argument using the 
idea of a short delay approximation <cit.>.
To our knowledge no rigorous 
approximation theorem is available, therefore we 
also give an alternative approach, based on time scale
analysis, at the end of the appendix (appendix <ref>). Both approaches yield
identical results. 

Using the variation-of constant formula, we find 

    Y_t = β∫_0^∞ e^-μ sX_t-s  ds,
        
    Z_t = β(1-σ_3/N)∫_0^∞ e^-μ(1+σ_2/N)s(N-X_t-s)  ds.

Let x_t=X_1,t/N and ε^2=1/N. Then, 

    X_1,t→ X_1,t+1  ε^-2(1+σ_1ε^2) ζ (1-x_t) ∫_0^∞ e^-μ sx_t-s  ds/∫_0^∞ e^-μ sx_t-s  ds+(1-σ_3 ε^2)∫_0^∞ e^-μ(1+σ_2 ε^2)s(1-x_t-s)  ds,


    X_1,t→ X_1,t-1  ε^-2 ζ x_t (1-σ_3^2)  ∫_0^∞ e^-μ(1+σ_2 ε^2)(t-s)(1-x_t-s)  ds/∫_0^∞ e^-μ sx_t-s  ds+(1-σ_3 ε^2)∫_0^∞ e^-μ(1+σ_2 ε^2)s(1-x_t-s)  ds.

With standard arguments, we obtain a stochastic 
delay differential equation (SDDE) at the evolutionary 
time scale τ = t ε^2 as

    dx_τ   =   (ε^-4(1+σ_1ε^2) ζ (1-x_τ) ∫_0^∞ e^-μ s/^2x_τ-s  ds/^-2∫_0^∞ e^-μ s/^2x_τ-s  ds+^-2(1-σ_3 ε^2)∫_0^∞ e^-μ(1+σ_2 ε^2) s/^2(1-x_τ-s)  ds
          -
    ε^-4 ζ x_τ ∫_0^∞ e^-μ(1+σ_2 ε^2)s/^2(1-x_τ-s)  ds/^-2∫_0^∞ e^-μ s/^2x_τ-s  ds+^-2(1-σ_3 ε^2)∫_0^∞ e^-μ(1+σ_2 ε^2)s/^2(1-x_τ-s)  ds)  dτ
       +   (^-2 (1+σ_1ε^2)
     ζ (1-x_τ) ∫_0^∞ e^-μ s
    /^2x_τ-s  ds/^-2∫_0^∞ 
    e^-μ s/^2x_τ-s  ds+^-2
    (1-σ_3 ε^2)∫_0^∞ 
    e^-μ(1+σ_2 ε^2)s/^2
    (1-x_τ-s)  ds
          +
    ^-2 ζ x_τ ∫_0^∞ e^-μ(1+σ_2 ε^2)s/^2(1-x_τ-s)  ds/^-2∫_0^∞ e^-μ s/^2x_τ-s  ds+^-2(1-σ_3 ε^2)∫_0^∞ e^-μ(1+σ_2 ε^2)s/^2(1-x_τ-s)  ds)^1/2 dW_τ.

We aim a small delay approximation. Therefore, we note that for a function Φ(t), which is sufficiently smooth and bounded, we have (for μ>0)

    ^-2 μ ∫_0^∞ e^-μ(τ-s)/^2Φ(τ-s)  ds
       =   ^-2∫_0^∞μ e^-μ s/^2(Φ(τ)-sΦ'(τ)+ O(s^2))  ds
       =   Φ(τ)-^2μ^-1Φ'(τ)+ O(^4).

Thus, at a formal level, 

    lim_→ 0ζ(^-4 μ ∫_0^∞ e^-μ s/^2x_τ-s  ds-^-2 x_τ)dτ
    =
    -ζ/μ  dx_τ
    =
    -G  dx_τ,

where G=ζ/μ. 
Note that this equation has to be interpreted in 
terms of the Euler-Maruyama approximation of an SDDE, 
where differential quotients are replaced by difference quotients. 
We again emphasize that this approach is only 
mend to be formal, as to our knowledge no rigorous approximation theorems 
for small delay approximations in the context of SDDE are available. The result is 
consistent with the time scale analysis of the 
present model (appendix <ref>); we conjecture that the short delay approach yields a valid approximation under suited (rather general) conditions. 
If we add 

    - ζ(^-4 μ ∫_0^∞ e^-μ s/^2x_τ-s  ds- ^-2 x_τ)dτ/^-2μ ∫_0^∞ e^-μ s/^2x_τ-s  ds+(1-σ_3 ε^2) ^-2 μ ∫_0^∞ e^-μ(1+σ_2 ε^2)s/^2(1-x_τ-s)  ds
 
to both sides of the SDDE and let → 0, we obtain

    (1+G)dx_τ   =   ζ(σ_1+σ_2+σ_3)  x_τ(1-x_τ)  dτ
    +
    (2ζ x_τ(1-x_τ))^1/2  dW_τ.

This equation yields the desired result with 
σ=σ_1+σ_2+σ_3.











 §.§ Logistic population dynamics and deterministic seed bank


Let p_i,j(k,l,t)=ℙ(X_1,t=i,X_2,t=j,Y_t∈(k,k+dk),Z_t∈(l,l+dl)) be the 
joint probability density of the resulting stochastic process (with discrete i,j and continuous k,l). The corresponding master equation reads:

    ṗ_i,j(k,l,t) + ∇[[                     β i-μ k; β (1-σ_3/N) j-μ (1+σ_2/N) l ]p_i,j(k,l,t)] 
    
    =    -[ζ (i + (1+σ_1/N)j) + γN-i-j/N(k+(1-σ_4/N)l]p_i,j(k,l,t) 
       + ζ(i+1)p_i+1,j(k,l,t) + (1+σ_1/N)ζ(j+1)p_i,j+1(k,l,t) 
       + γ(N-i-j+1)k/N p_i-1,j(k,l,t) + (1-σ_4/N)γ(N-i-j+1)l/N p_i,j-1(k,l,t),

where the operator ∇ acts with respect to continuous state space variables k,l. 
Standard arguments yield the Fokker-Plank-approximation 
for large populations, where  

    _̣tu(x_1,x_2,y,z,t) 
       =∂_x_1{[ζ x_1-γ(1-x_1-x_2)y]u(x_1,x_2,y,z,t)}
         +∂_x_2{[(1+σ_1/N)ζ x_2-(1-σ_4/N)γ(1-x_1-x_2)z]u(x_1,x_2,y,z,t)}
         +_̣y{[(μ y -β x_1)]u(x_1,x_2,y,z,t)}
         +_̣z{[(1+σ_2/N)μ z -(1-σ_3/N)β x_2]u(x_1,x_2,y,z,t)}
         +1/2N_̣x_1^2{[ζ x_1+γ(1-x_1-x_2)y]u(x_1,x_2,y,z,t)}
         +1/2N_̣x_2^2{[(1+σ_1/N)ζ x_2+(1-σ_4/N)γ(1-x_1-x_2)z]u(x_1,x_2,y,z,t)}.

Note that the second order noise terms are solely 
due to x_1 and x_2; no noise is added by the 
seed bank variables y and z. 



  §.§.§ Deterministic model

The corresponding deterministic model 
(drift terms only)  yields the ODEs

    ẋ_̇1̇   = γ(1-x_1-x_2)y -ζ x_1, 
    ẋ_̇2̇   = γ(1-x_1-x_2)z -ζ x_2, 
    ẏ   = β x_1 - μ y,
    ż   = β x_2 - μ z.

The lifetime reproductive success (or basic reproduction number) of a plant reads 
R_0=β γ/(μζ).
If R_0>1, the plant population can persist.
Under this condition, there is a line of stationary solutions:

Assume R_0>0. 
Let ϑ := β/μ, κ := (γϑ-ζ)/γϑ. 
 Then, κ∈[0,1], and there is a line of stationary points in [0,κ]^2×ℝ^2_+ given by

    (x_1,x_2,y,z)=(κ x,κ (1-x),ϑ κ x,ϑ κ (1-x)),       x∈[0,1].

The line of stationary points is transversally stable (locally and globally).

Proof: 
It is straightforward to check that the line indicated above 
consists indeed of stationary points with non-negative values. The local stability is a consequence of Hartman-Grobman 
and the analysis of the Jacobian (due to the block structure 
of this matrix the eigenvalues can be stated explicitly). 
Note that one eigenvalue necessarily is zero with an eigenvector pointing in the direction of the line of stationary points. Now, to show global stability, we prove that the system will approach the equilibrium line from any starting point. We first respectively denote P=x_1+x_2 and S=y+z as the total plant and seed populations and consider the resulting reduced system

    Ṗ = γ(1-P)S -ζ P,   Ṡ = β P - μ S.

Here we use that we only consider weak selection: all 
selection effects tend to zero for N→∞, and hence the ODE describes the neutral case. The divergence 
of this system is negative, and thus the combination of the 
theorems of Bendixon-Dulac and Pointcaré-Bendixon imply 
that trajectories (P,S) tend to stationary points. 
This observation yields the desired global stability.


Note that in equilibrium, x_1+x_2=κ. That is, κ  N 
represents the average above-ground population size of the model 
conditioned on non-extinction.




  §.§.§ Dimension reduction by time scale analysis

The computations in this section follow 
closely the calculations in the paper of <cit.> to perform a dimension reduction
by a time scale analysis. First, new local variables for the boundary layer around the equilibrium line are defined as

    x_1   =   κ  +ε/2 , 
        x_2 = κ (1-) +ε/2 , 
        y    =   κ ϑ +ε/2 (++ϑ ), 
        z =  κ ϑ (1-) +ε/2 (-+ϑ ),

where

    =x_1-x_2+κ/2κ,    = ε^-1(x_1+x_2-κ),    = ε^-1(y+z-ϑ(x_1+x_2)),    = ε^-1(y-z-ϑ(x_1-x_2),

and ε^2=1/N. 
For the transformed density ρ(,,,,t;ε)=u(x_1,x_2,y,z,t), we find 


    _̣tρ =      ε_̣{[-ζ/2ϑκ-γϑ/2(1-2)]ρ} + ε^2_̣{[γ/2κ+ζ/2(-σ_1-σ_4)(1-)]ρ}
       + _̣{[-ζ/ϑ+γϑκ]ρ} + ε_̣{[γ+γϑ^2+(σ_4+σ_1)ζκ(1-)]ρ} + ε^2_̣{…}
       + _̣{[μ+ζ-γϑ^2κ]ρ}
       + ε_̣{[-γϑ-γϑ^2^2+((-σ_1-σ_4)ζϑ
         +(σ_3+σ_2)β)κ(1-)]ρ} + ε^2_̣^2{[…]}
       + _̣{[μ+ζ+γϑ^2κ(1-2)]ρ}
       + ε_̣{[-γϑ+((σ_4+σ_1)ζϑ+(-σ_2-σ_3)β)κ(1-)]ρ} + ε^2_̣^2{[…]}
       + ε^2_̣^2{[ζ/4κ]ρ} + ε_̣_̣{[-ζ(1-2)]ρ} + ε^2_̣_̣{[…]}
       + ε_̣_̣{[ζϑ(1-2)]ρ} + ε^2_̣_̣{…} + ε_̣_̣{[-ζϑ]ρ} + ε^2_̣_̣{[…]}
       + _̣^2{[ζκ]ρ} + ε_̣^2{…} + _̣^2{[ζϑ^2κ]ρ} + ε_̣^2{[…]}
       + _̣^2{[ζϑ^2κ]ρ} + ε_̣^2{[…]ρ}
       + _̣_̣{[-2ζϑκ]ρ} + ε_̣_̣{…} + _̣_̣{[2ζϑκ(1-2)]ρ} + ε_̣_̣{[…]}
       + _̣_̣{[-2ζϑ^2κ(1-2)]ρ} + ε_̣_̣{[…]ρ},


where terms that will not be needed later on are represented by “…”. The reparametrized Fokker-Planck equation in operator notation is



    _̣tρ(,,,,t;ε) = (Ł+ε + ε^2L^(2))ρ(,,,,t) + 𝒪(ε^3),


with linear differential operators



    L^(0)ρ =    _̣{[-ζ/ϑ+γϑκ]ρ}+_̣{[(μ+ζ)-γϑ^2κ]ρ}
       +_̣{[(μ+ζ)+γϑ^2κ(1-2)]ρ}
       +_̣^2{ζκρ} +_̣^2{ζϑ^2κρ} + _̣^2{ζϑ^2κρ}
       + _̣_̣{-2ζϑκρ} + _̣_̣{2ζϑκ(1-2)ρ} + _̣_̣{-2ζϑ^2κ(1-2)ρ},
    
    
    
    
    
    L^(1)ρ=    _̣{[-ζ/2ϑκ-γϑ/2(1-2)]ρ}+_̣{[γ+γϑ^2+(σ_4+σ_1)ζκ(1-)]}
       +_̣{[-γϑ-γϑ^2^2-(σ_1+σ_4)ζϑκ(1-)+(σ_3+σ_2)βκ(1-)]ρ}
       +_̣{[-γϑ+(σ_4+σ_1)ζϑκ(1-)-(σ_2+σ_3)βκ(1-)]ρ}
       +_̣_̣{-ζ(1-2)ρ} + _̣_̣{ζϑ(1-2)ρ}+ _̣_̣{-ζϑρ}
       +_̣^2{…} +_̣^2{…} + _̣^2{…} +_̣_̣{…} +_̣_̣{…} + _̣_̣{…},
    
    
    
    
    
    L^(2)ρ =     _̣{[γ/2κ-(σ_1+σ_4)ζ/2(1-)]ρ}+_̣^2{ζ/4κρ}+_̣{…}+_̣{…}+_̣{…}.



We employ a time scale separation and focus on a solution evolving on the slow time τ = ε^2t=t/N using the Ansatz


    ρ(,,,,t) = ρ^(0)(,,,,ε^2t) + ερ^(1)(,,,,ε^2t) + ε^2ρ^(2)(,,,,ε^2t)+𝒪(ε^3).


Plugging this into (<ref>) and comparing same order terms, we have


    L^(0)ρ^(0)=0,   L^(0)ρ^(1)=-L^(1)ρ^(0),   L^(0)ρ^(2)=_̣τρ^(0)-L^(1)ρ^(1)-L^(2)ρ^(0),


which indicates that ρ^(0) can be written as


    ρ^(0)(,,,,τ) = f(,τ)ρ̂(,,;),

with a time independent normal distribution ρ̂(,,;) in ,, that satisfies L^(0)ρ̂=0. 
The function f(,τ) modifies ρ̂(,,;) 
and represents the time evolution of ρ^(0) along the line of stationary points. Integrating the last equation of (<ref>) from -∞ to ∞ with respect to ,, - the left hand side is a total derivative w.r.t. variables of integration and becomes zero - yields the evolution equation

    _̣τf = ∫ρ^(1)d(,,) + ∫ L^(2)d̊(,,).

We start by computing the second integral. Since all terms that are full derivatives w.r.t. ,  or  vanish upon integration, we have

    ∫ L^(2)ρ^(0) d(,,) 
       = _̣{∫[γ/2κ+(-σ_1-σ_4)ζ/2(1-)]ρ^(0)d(,,)} + _̣^2{∫ζ/4κ ρ^(0)d(,,)}
       = _̣{∫[γ/2κ]ρ^(0)d(,,)} + _̣{(-σ_1-σ_4)ζ/2(1-)f} + _̣^2{ζ/4κf}.

We take a closer look at the first integral. 
We aim to write =(L^(0))^+h^+-g() 
for suitable h^+, g, where (L^(0))^+ is the adjoint of L^(0). 
If we can do so, then

    ∫ρ^(0)  d(,,)
    = ∫ h^+L^(0)ρ^(0)  d(,,)-g()f(,τ)
    =-g()f(,τ)
    .

To identify h^+ and g, we reduce the 
problem to linear algebra. 
Define the finite-dimensional vector space (for k∈ℕ)

    H_k:={P:ℝ^4→ℝ, P polynomial homogenous of degree k (w.r.t. variables ,,)},

so that, e.g., γ/2κ ∈ H_2, while γ/2κ/;+c∉ H_2 for c∈ℝ∖{0}. Examining (Ł)^+, we find

    (Ł)^+h^+ = 0,    for h^+∈ H_0. 
        (Ł)^+h^+ = h∈ H_1,    for  h^+∈ H_1. 
        (Ł)^+h^+ = h+g∈ H_2 ⊕ H_0,    for  h^+∈ H_2.

In particular the last observation (Ł)^+H_2→ H_2 ⊕ H_0,
h^+↦(h,g) 
allows to define an operator 
M:H_2→ H_2, h^+↦ h. 
To simplify notation, we identify vectors w.r.t. a given basis 
in H_k and the corresponding polynomials. 


 W.r.t. the canonical basis 
(^2,^2,^2,,,), the operator M has the 
representation

    M = [        -2γϑκ            0            0        γϑ^2κ  -γϑ^2κ(1-2)            0;            0      -2(μ+ζ)            0          ζ/ϑ            0            0;            0            0      -2(μ+ζ)            0            0            0;         2ζ/ϑ       2γϑ^2κ            0   -(μ+ζ+γϑκ)            0  -γϑ^2κ(1-2);            0            0 -2γϑ^2κ(1-2)            0   -(μ+ζ+γϑκ)        γϑ^2κ;            0            0            0            0          ζ/ϑ      -2(μ+ζ) ],

and is invertible, if 
all parameters are positive. 
Let h=(a_1(),…,a_6())^T∈ H_2. 
Then, (L^(0))^+(h) = M h + g(), 
where g()∈ H_0 is uniquely defined by g()=Ĝ()  h and
Ĝ() denoting the row-vector

    Ĝ() = (2ζκ, 2ζϑ^2κ, 2ζϑ^2κ, 
    - 2ζϑκ, 2ζϑκ(1-2), - 2ζϑ^2κ(1-2)).
 
For h̃∈ H_2, we find (L^(0))^+(M^-1h̃) = h̃ + g() 
with g()=Ĝ()  M^-1h̃, so that

    ∫h̃ ρ(,,,,τ) d(,,)=-g()  f(,τ).
 

Proof:  It is straightforward to obtain the representation of M 
and g by applying Ł to the elements of the basis given above. Using, e.g., 
Gauß-elimination, we find that M is invertible if the parameters 
are positive. In order to obtain the i''th component of Ĝ(), 
consider the i'th entry of the canonical basis b_i, compute (L^(0))^+b_i,
and identify the component that is in H_0. E.g., for i=1 we find 
b_1=^2, and 
(L^(0))^+^2 
= -2 [(-ζ/ϑ)+γϑκ]+2ζκ, 
where -2 [(-ζ/ϑ)+γϑκ]∈ H_2 and 
2ζκ∈ H_0. 
Hence, (Ĝ())_1=2ζκ.
The equation (L^(0))^+(M^-1h) = h + g() with g()=Ĝ()  M^-1 h  
implies

    ∫ h ρ(,,,,τ) d(,,) 
    = 
    ∫ [(L^(0))^+(M^-1h)- g()] ρ(,,,,τ) d(,,) 
    = 0 - g()  f(,τ).




Solving the system for h(,,,)=γ/2κ, we obtain

    g() = -γζϑ/2(μ+γϑ)(1-2)

by using the computer algebra package MAXIMA <cit.>. In summary, we have

    ∫ L^(2)ρ^(0)d(,,) 
    = _̣{(
    γζϑ/2(μ+γϑ)
    +(-σ_1-σ_4)ζ/2)
    (1-2)f(,τ)} + _̣^2{ζ/4κf(,τ)}.

Now, we turn to the computation of the first integral in (<ref>). With

    h_0^+ := ζ/2ϑκ(μ+ζ)[+(1-2)+(β+ζϑ)/ζ(1-2)],

we obtain
(Ł)^+h_0^+ = -ζ/(2ϑκ)-γϑ/2(1-2).
Remember that full derivatives w.r.t. ,  or  vanish upon integration. So,

    ∫ρ^(1)d(,,) 
       = _̣∫[-ζ/2ϑκ-γϑ/2(1-2)]ρ^(1)d(,,) = _̣∫[(Ł)^+h_0^+]ρ^(1)d(,,) 
       =_̣∫ h_0^+ [Łρ^(1)]d(,,) = _̣∫-h_0^+[ρ^(0)]d(,,).

To handle ∫-h_0^+[ρ^(0)]d(,,), 
we use partial integration to move all derivatives in L^(1) 
w.r.t. ,, from ρ^(0) to h_0^+ (note that 
h_0^+ is linear in these variables, s.t. all second derivatives 
in these variables vanish). Let us consider one of these terms 
occurring in ∫-h_0^+[ρ^(0)]d(,,):

    ∫[+(1-2)+(β+ζϑ)/ζ(1-2)] ∂_{…}  d(,,) 
       =
    ∫∂_[+(1-2)+(β+ζϑ)/ζ(1-2)] {…}  d(,,) 
    = 0.

 It is not possible to use the same 
procedure for , as we do not integrate w.r.t. . Here
we use the product rule, e.g.,

    ∫[+(1-2)+(β+ζϑ)/ζ(1-2)] _̣{[-ζ/2ϑκ-γϑ/2(1-2)]ρ}  d(,,) 
       = 
    _̣x ∫[+(1-2)+(β+ζϑ)/ζ(1-2)] {[-ζ/2ϑκ-γϑ/2(1-2)]ρ}  d(,,) 
       -
    ∫_̣[+(1-2)+(β+ζϑ)/ζ(1-2)] {[-ζ/2ϑκ-γϑ/2(1-2)]ρ}  d(,,).

In this way, we obtain 

    ∫ρ^(1)d(,,) 
        = ζ/2ϑκ(μ+ζ)_̣^2∫[ [t]    γϑ^2(μ+ζ)/2ζ(1-2)^2^2+ζ/2ϑκ^2+γϑ/2(1-2)^2
       +μ+γϑ/2κ(1-2)+ζ/2ϑκ(1-2)]ρ^(0)d(,,) =.5pt-.9pta
          +ζ/2ϑκ(μ+ζ)_̣∫[γϑ^22μ+ζ/ζ(1-2)^2+γϑμ+ζ/ζ(1-2)+μ+ζ-γϑκ/κ+ζ/ϑκ]ρ^(0)d(,,)=.5pt-.9ptb
          +_̣∫[(σ_4+σ_1)ζ/2(1-)-(σ_4+σ_1+σ_3+σ_2)ζμ/μ+ζ(1-)]ρ^(0)d(,,)
    
    +ζ/2ϑκ(μ+ζ)_̣∫ -β(1-2)_̣{(1-2)ρ^(0)}-_̣{ζϑρ^(0)}d(,,)
    =.5pt-.9ptc.

We will treat all three terms separately starting with .5pt-.9ptc.

    .5pt-.9ptc   = _̣{[(σ_4+σ_1)ζ/2(1-)-(σ_4+σ_1+σ_3+σ_2)ζμ/μ+ζ(1-)]f}
         -ζ/2ϑκ(μ+ζ)_̣^2{[β(1-2)^2+ζϑ]f}-ζ/2ϑκ(μ+ζ)_̣{2β(1-2)f}.

To compute .5pt-.9ptb, we proceed as in the computations for ∫ L^(2)ρ^(0)d(,,), that is, we solve

    M[ a_1(); a_2(); a_3(); a_4(); a_5(); a_6() ]
    = [ γϑ^22μ+ζ/ζ(1-2);               0;               0;    γϑμ+ζ/ζ(1-2);       μ+ζ-γϑκ/κ;            ζ/ϑκ ]_=̂h,

and find

    .5pt-.9ptb = ζ/2ϑκ(μ+ζ)_̣{-g()f(,τ)} = ζ/2ϑκ(μ+ζ)_̣{[2β-γϑ^2κ(μ+ζ)/μ+γϑ](1-2)f(,τ)}.


For .5pt-.9pta, this can be done similarly, where the system to be solved is


    M[ a_1(); a_2(); a_3(); a_4(); a_5(); a_6() ]
    = [ γϑ^2(μ+ζ)/2ζ(1-2)^2;                   0;               ζ/2ϑκ;         γϑ/2(1-2)^2;        μ+γϑ/2κ(1-2);          ζ/2ϑκ(1-2) ]_=̂h,


and as a result, we have

    .5pt-.9pta   = ζ/2ϑκ(μ+ζ)_̣^2{-g()f(,τ)}
       = ζ/2ϑκ(μ+ζ)_̣^2{[β+ζϑ/2(1-2)^2-ζ^2ϑ/2(μ+ζ)[(1-2)^2-1]]f(,τ)}.

With 
∫ρ^(1)d(,,) = .5pt-.9pta +.5pt-.9ptb + .5pt-.9ptc 
and 

    _̣τf = ∫ρ^(1)d(,,)+∫ L^(2)ρ^(0)d(,,),

we obtain

    _̣τf = ζ/(1+G)^2_̣x^2{(1-)f} + ζ(-σ_1-σ_4-σ_2-σ_3)/1+G_̣x{(1-)f}.


















 §.§ Fixed population size and fluctuating seed bank


We only present the main steps for this and the following model,
as the computations are lengthy and resemble that of the last model.






Let p_i,j,k(t)=ℙ(X_t=i,Y_t=j,Z_t=k) be the probability distribution of the resulting stochastic process. The corresponding master equation reads:


    ṗ_i,j,k(t) 
    =    -[ζ N + σ_1/Nζ (N-i) + β N - σ_3/Nβ (N-i) + μ (j+k) + σ_2/Nμ k]p_i,j,k(t) 
       + [ζi(j+1)/j+k+1]p_i,j+1,k(t) 
     + [ζ(i+1)(k+1)/j+k+1]p_i+1,j,k+1(t) 
       + [(1+σ_1/N)ζ(N-i+1)(j+1)/j+k+1]p_i-1,j+1,k(t) 
       + [(1+σ_1/N)ζ(N-i)(k+1)/j+k+1]p_i,j,k+1(t) + [β i]p_i,j-1,k(t) + [μ (j+1)]p_i,j+1,k(t) 
       + [(1-σ_3/N)β (N-i)]p_i,j,k-1(t) + [(1+σ_2/N)μ (k+1)]p_i,j,k+1(t).


We now obtain the Fokker-Planck equation for the approximating diffusion process as described before. The first step is to transform the system to a quasi-continuous state space by scaling with N^-1, i.e., defining x=i/N, y=j/N, z=k/N, h=1/N 
and the quasi-continuous density u(x,y,z,t)=p_i,j,k(t). The resulting PDE reads


    _̣tu(x,y,z,t) 
    =     _̣x[(ζ x-ζy/y+z-σ_1/Nζ(1-x)y/y+z)u(x,y,z,t)] 
       + _̣y[(ζy/y+z+(μ y -β x) +σ_1/Nζ(1-x)y/y+z)u(x,y,z,t)] 
       + _̣z[(ζz/y+z+(μ z -β (1-x)) +σ_1/Nζ(1-x)z/y+z + σ_3/Nβ(1-x) +σ_2/Nμ z )u(x,y,z,t)] 
       + 1/2N_̣x^2[(ζxz/y+z+(1+σ_1/N)ζ(1-x)y/y+z)u(x,y,z,t)]
       + 1/2N_̣y^2[(ζxy/y+z+(1+σ_1/N)ζ(1-x)y/y+z+β x + μ y)u(x,y,z,t)]
       + 1/2N_̣z^2[(ζxz/y+z+(1+σ_1/N)ζ(1-x)z/y+z+(1-σ_3/N)β (1-x) + (1+σ_2/N)μ z)u(x,y,z,t)]
       + 1/N_̣x_̣y[-(1+σ_1/N)ζ(1-x)y/y+zu(x,y,z,t)]
    + 1/N_̣x_̣z[ζxz/y+zu(x,y,z,t)].



For N→∞, we obtain a deterministic
model, governed by the ODEs

    ẋ   = -ζ x + ζy/y+z = -ζ xz/y+z + ζ(1-x)y/y+z, 
    ẏ   = β x - μ y - ζy/y+z, 
    ż   = β(1-x) - μ z - ζz/y+z,

which is a neutral competition model since the selection terms vanish in the limit.

Let ϑ = (β-ζ)/μ>0. Then, there is a line of stationary points for (<ref>) in [0,1]×ℝ^2_+ given by

    (x,y,z)=(x,ϑ x,ϑ(1-x)),      x∈[0,1].

The line of stationary points is transversally stable. The eigenvectors perpendicular to the line of stationary points (together with the eigenvalues) read

    X_1 = [   0;   x; 1-x ],  λ_1 = -μ,
        
    X_2 = [  ζ; -β;  β ],  λ_2 = -ζ-β/ϑ.


The proof of this proposition is straight forward, along the lines
of the proof of Proposition <ref>. 
To formulate the inner solution, we introduce local coordinates


    =x,    = ε^-1(y+z-ϑ),
       = ε^-1(y-z-ϑ(2x-1),  ρ(,,,t)=u(x,y,z,t),


where ϑ=(β-ζ)/μ again and ε^2=1/N. Alternatively formulated, we have


    x = ,   y = ε/2(+)+ϑ x,
       z = ε/2(-)+ϑ(1-x).


 can be thought of as measuring the deviation of the total amount of seeds from its deterministic value ϑ and  as measuring the deviation of the allele ratio in seeds from the allele ratio in plants. Both are scaled by ε^-1 so we expect them to be of order 𝒪(1). By transforming derivatives, we have

    _̣x    = /x̣_̣ + /x̣_̣ + /x̣_̣
    = _̣ - 2ε^-1ϑ_̣,  _̣y = ε^-1(_̣+_̣),  _̣z 
    = ε^-1(_̣-_̣).


We now approximate _̣tρ(,,,t) by transforming all terms on the r.h.s. of (<ref>) and ignoring terms of 𝒪(ε^3). We find

    _̣tρ(,,,t)=(L^(0)+ε L^(1)+ε^2 L^(2))ρ(,,,t)+𝒪(ε^3),


with linear differential operators L^(0), L^(1) and L^(2) that take the form

    L^(0)ρ =   _̣{μρ}+_̣{[μ+ζ(1+1/ϑ)(+(1-2))]ρ}
       +_̣^2{βρ} -_̣_̣{2β(1-2)ρ} + _̣^2{[β+4ζϑ(1+ϑ)(1-)]ρ},
    
    L^(1)ρ=   -_̣{ζ/2ϑ(+(1-2))ρ}
        +_̣{[σ_1ζ(1-)+σ_3β(1-)+σ_2μϑ(1-)]ρ}
       -_̣{[ζ/ϑ(1+1/ϑ)(+(1-2)) -σ_1ζ (2(ϑ+1)-1) (1-)+σ_3β(1-)+σ_2μϑ(1-)]ρ}
       +_̣^2{μ/2ρ} 
     +_̣^2{[ζ(1+ϑ)(+(1-2))-4ζ(ϑ+1)(1-)+μ/2- 2ζ(-)(1+ϑ) ]ρ}
       -_̣_̣{2(2ϑ+1)ζ(1-)ρ} 
     +_̣_̣{(μ+ζ(1+1/ϑ)(+(1-2)))ρ},
    
    L^(2)ρ =   _̣{[ζ/2ϑ^2(+(1-2))
                      -σ_1ζ(1-)]ρ}+_̣^2{ζ(1-)ρ}+_̣{…}+_̣{…}.







Before we proceed, we define the functions

    h_0(,;)   =     -[+(1-2)]/μ+ζ(1+1/ϑ),
    
    h_1(,;)   =     -[(+(1-2))]/2μ+ζ(1+1/ϑ),
    
    h_2(,;)   =     ((1-2))^2ζ(1+1/ϑ)/2(2μ+ζ(1+1/ϑ))(μ+ζ(1+1/ϑ)
       - (1-2)μ/(2μ+ζ(1+1/ϑ))(μ+ζ(1+1/ϑ) - ^2/2(μ+ζ(1+1/ϑ)), 
    
    g()   =     4((1-))β+ζϑ(1+ϑ)/μ+ζ(1+1/ϑ),

and note that 


    (L^(0))^+h_0(,;)    = (+(1-2)),
    
    (L^(0))^+h_1(,;)    = (+(1-2)),
    
    (L^(0))^+h_2(,;) + g()    = (+(1-2)).

That is, h_0 and h_1 are eigenfunctions of (L^(0))^+. 
Observing the system on slow time τ=ε^2t=t/N, 
we see rapid dynamics for τ close to zero, τ∈𝒪(1/N) 
to be precise. This is the new boundary layer. For τ∈𝒪(1), 
only slow drift effects along the equilibrium line should remain. 
Being interested in the outer solution that develops on the slow time scale, 
we make the Ansatz 


    ρ(,,,t;ε) = ρ^(0)(,,,ε^2t) + ερ^(1)(,,,ε^2t) + ε^2ρ^(2)(,,,ε^2t)+𝒪(ε^3),


assuming _̣tρ=𝒪(ε^2) to focus on the outer solution and neglect the boundary layer. We plug this Ansatz into (<ref>), compare same order terms on both sides and obtain


    L^(0)ρ^(0)=0,   L^(0)ρ^(1)=-L^(1)ρ^(0),   L^(0)ρ^(2)=_̣τρ^(0)-L^(1)ρ^(1)-L^(2)ρ^(0).


We have ∫(+(1-2))ρ d(,)
=∫(L^(0))^+ h^+_1ρ d(,)=0, and hence

    ∫ L^(2)ρ^(0) d(,)    =   _̣∫[ζ/2ϑ(+(1-2))-σ_1ζ(1-)]ρ^(0) d(,)
    +_̣^2∫ζ(1-)ρ^(0) d(,)
       =    
    -σ_1 ζ _̣ ((1-) f(τ,))
    + ζ_̣^2(ζ(1-)f(τ,) ).

Using the same procedure as in Appendix <ref>, 
consider (<ref>)-(<ref>), to handle 
∫ L^1ρ^1 d(,,), and by applying (<ref>) - (<ref>), we find 

    ∫ L^1ρ^1 d(,,) = T_a+T_b+T_c,

where

    T_a    =   _̣-ζ^2/2ϑ^2((ϑ+1)ζ+μϑ)∫(+(1-2))ρ^(0) d(,,) = 0,
    
    T_b    =   _̣(
     -ζ(1-) (-ζσ_1ϑ - ζσ_2 + β(σ_2+σ_3))/(ϑ+1)ζ+μϑ∫ρ^(0) d(,,)),
    
     T_c    =   ^̣2_( 
     ζ^2/4ϑ^2((1+1/ϑ) ζ+μ)∫ (+(1-2)) + (+(1-2))  (1-2)) ρ^(0) d(,,))
           - ^̣2_( 
     ζ^2/ϑ ((1+1/ϑ) ζ+μ)∫ (2ϑ+1)(1-) ρ^(0) d(,,))
       =   ^̣2_( 
     -ζ^2 [(ϑ+1)^2ζ+μϑ(1+2ϑ)+β)] (1-) /(ϑ((1+1/ϑ)ζ+μ))^2∫ ρ^(0) d(,,)).

 With

    _̣τf = ∫ρ^(1)d(,,)+∫ L^(2)ρ^(0)d(,,),

we find

    _̣τf(,τ)  
        =    - ζ[σ_1+(1-1/Y)σ_2+σ_3]/(1+(1-1/Y)G)_̣{(1-)f(,τ)}
          + ζ/(1+(1-1/Y)G)^2_̣^2{(1-)
    f(,τ)}.












 §.§ Logistic population dynamics and fluctuating seed bank





Let p_i,j,k,l(t)=P(X_1,t=i,X_2,t=j,Y_t=k,Z_t=l) be the probability distribution of the resulting stochastic process. The corresponding master equation reads:


    ṗ_i,j,k,l(t)    = -[ [t]    ζ (i + (1+σ_1/N)j) + γN-i-j/N(k+(1-σ_4/N)l) 
       + μ (k + (1+σ_2/N)l)
     + β (i + (1-σ_3/N)j)]p_i,j(k,l,t) 
         + ζ(i-1)p_i-1,j,k,l(t)        
          + (1--σ_1/N)ζ(j-1)p_i,j-1,k,l(t) 
         + γ (1-(i+1)/N-j/N) (k+1) p_i-1,j,k+1,l(t) 
         + (1-σ_4/N)γ (1-i/N-(j+1)/N) (l+1) p_i,j-1,k,l+1(t) 
         + β i p_i,j,k-1,l(t) + (1-σ_3/N)β j p_i,j,k,l-1(t)
         + μ(k+1)p_i,j,k+1,l(t) 
          + (1+σ_2/N) μ(l+1)p_i,j,k,l+1(t).


We transform the system to the quasi-continuous state space with rescaled parameters x_1=i/N, x_2=j/N, y=k/N, z=l/N, h=1/N and quasi-continuous density u(x_1,x_2,y,z,t)=p_i,j,k,l(t). After expanding about (x_1,x_2,y,z) in terms up to order 𝒪(ε^2), we have a 4-dimensional Fokker-Planck equation again:


    _̣tu    =_̣x_1{[ζ x_1-γ(1-x_1-x_2)y]u}
         +_̣x_2{[(1+σ_1/N)ζ x_2-(1-σ_4/N)γ(1-x_1-x_2)z]u}
         +_̣y{[μ y -β x_1
                                    + γ(1-x_1-x_2)y]u}
         +_̣z{[(1+σ_2/N)μ z -(1-σ_3/N)β x_2 
                                    + (1-σ_4/N)γ(1-x_1-x_2)z]u}
         +1/2N_̣x_1^2{[ζ x_1+γ(1-x_1-x_2)y]u}
         +1/2N_̣x_2^2{[(1+σ_1/N)ζ x_2
                                    +(1-σ_4/N)γ(1-x_1-x_2)z]u}
         +1/2N_̣y^2{[μ y+β x_1
                                               + γ(1-x_1-x_2)y]u}
         +1/2N_̣z^2{[ (1+σ_2/N)μ z 
                                    + (1-σ_3/N)β x_2
                                    + (1-σ_4/N)γ(1-x_1-x_2)z
                                     ]u}
          -1/N_̣x1_̣y{[γ(1-x_1-x_2)y
                                     ]u}
           -1/N_̣x2_̣z{[(1-σ_4/N)γ(1-x_1-x_2)z
                                     ]u}.


The limiting deterministic system for N→∞ has dynamics

    ẋ_̇1̇   = γ(1-x_1-x_2)y -ζ x_1, 
    ẋ_̇2̇   = γ(1-x_1-x_2)z -ζ x_2, 
    ẏ   = β x_1 - μ y - γ(1-x_1-x_2)y,
    ż   = β x_2 - μ z - γ(1-x_1-x_2)z.

As before, we can show that a line of stable equilibria exists:


Let ϑ := (β-ζ)/μ>0, 
κ := (γϑ-ζ)/γϑ∈ [0,1]. 
Then, there is a line of stationary points 
in [0,κ]^2×ℝ^2_+ given by

    (x_1,x_2,y,z)=(x, κ-x, ϑ x, ϑ(κ-x)),       x∈[0,κ].


The line of stationary points is transversally stable (locally and globally).


Since the proof is similar to the proof of Proposition <ref>, 
it is omitted here.



  §.§.§ Perturbation approximation

As before, new local variables for the boundary layer around the equilibrium line are defined:

    x_1   =   κ  +1/2ε , 
       x_2 = κ (1-) +1/2ε , 
       y    =   κ ϑ +1/2ε (++ϑ ), 
       z =  κ ϑ (1-) +1/2ε (-+ϑ ),

where

    =   x_1-x_2+κ/2κ,    = ε^-1(x_1+x_2-κ), 
         =   ε^-1(y+z-ϑ(x_1+x_2)),    = ε^-1(y-z-ϑ(x_1-x_2)),

and ε^2=1/N. 
The considerations for the global behavior are the same as for model 2. We choose the same local variables for the boundary layer around the equilibrium line and transform (<ref>) up to terms of 𝒪(ε^3) and higher. 
The resulting reparametrized Fokker-Planck equation in the local variables is

    _̣tρ(,,,,t) = (Ł+ε + ε^2L^(0))ρ(,,,,t) + 𝒪(ε^2),

with linear differential operators


    L^(0)ρ   =     _̣{[-ζ/ϑ 
                        + γϑκ ] ρ}
       + _̣{[(β+ζϑ)/ϑ     
                     - γϑ (1-ϑ)κ] ρ}
       + _̣{[ (β+ζϑ)/ϑ     
             + γϑ (1+ϑ)κ (1-2)] ρ}
       + ^̣2_v{ζ κ  ρ}
       + ^̣2_{[ κ (β + ζϑ (1+ϑ))] ρ}
       + ^̣2_{[κ (β + ζϑ  (1+ϑ))] ρ}
       - _̣ _̣{ζ (1+2ϑ)κρ}
     +_̣ _̣{ζ(1+2ϑ)κ 
                         (1-2)ρ}
       -_̣ _̣{ 
        [2 (1+ϑ)ϑζ +2 β] κ (1-2) 
      ρ},
    
    
    
    
    
    L^(1)ρ=      
      - _̣{[
      ζ/2ϑκ  +γϑ/2  (1-2)]ρ}
     + _̣{[
        γ (ϑ +) + (σ_4+σ_1)ζκ(1-)
        ]ρ}
       +_̣{[ 
        - γ(1+ϑ)(ϑ+) 
         + (σ_2+σ_3)βκ(1-) 
         - (σ_4 (1+1/ϑ)+σ_1)(1-)ζϑκ]ρ}
       +_̣{[ 
           -  γ(1+ϑ) 
           - (σ_2+σ_3)βκ(1-) 
           + (σ_4 (1+1/ϑ)+σ_1)(1-)ζϑκ]ρ}
        - _̣_̣{ζ (1-2)ρ}
      + _̣_̣{ζ(1-2)(ϑ+1/2)ρ}
      - _̣_̣{ζ(ϑ +1/2)ρ}
    
       +_̣^2{…} +_̣^2{…} + _̣^2{…} +_̣_̣{…} +_̣_̣{…} + _̣_̣{…},     
    
    
    
    
    
    L^(2)ρ =      _̣{[γ/2κ-(σ_1+σ_4)ζ/2(1-)]ρ}+_̣^2{ζ/4κρ}+_̣{…}+_̣{…}+_̣{…}.


To handle the integrals below, we define an operator 
M:H_2→ H_0 resp. Ĝ:H_2→ H_0 in  a 
similar way as above (Proposition <ref>). 
Recall that we use 
(^2,^2,^2,,,) as the basis in H_2.
The operator M has the 
representation

    M = [           -2γϑκ               0               0        γϑ(1+ϑ)κ  -γϑ(1+ϑ)κ(1-2)               0;               0      -2(β+ζϑ)/ϑ               0             ζ/ϑ               0               0;               0               0      -2(β+ζϑ)/ϑ               0               0               0;            2ζ/ϑ       2γϑ(1+ϑ)κ               0   -γϑκ-(β+ζϑ)/ϑ               0  -γϑ(1+ϑ)κ(1-2);               0               0 -2γϑ(1+ϑ)κ(1-2)               0   -γϑκ-(β+ζϑ)/ϑ        γϑ(1+ϑ)κ;               0               0               0               0             ζ/ϑ      -2(β+ζϑ)/ϑ ].



If we define g via g()=Ĝ()  h, we find

    Ĝ() = κ   (
      2ζ,   
          2(ϑ(ϑ+1)ζ+β),   
         (2(ϑ(ϑ+1)ζ+β),   
       - (1+2ϑ)ζ,   
         (1+2ϑ)ζ (1-2),   
       -2(ϑ(1+ϑ)ζ+β)(1-2)
    ).



As before, we employ time scale separation and focus on a solution evolving on the slow time τ = ε^2t=t/N, using the Ansatz


    ρ(,,,,t) = ρ^(0)(,,,,ε^2t) + ερ^(1)(,,,,ε^2t) + ε^2ρ^(2)(,,,,ε^2t)+𝒪(ε^2).

If we define h_0^+ by 

    h_0^+ 
    =
    ζ/2κ(β+ζϑ)[+(1-2)+β+ζϑ/ζ (1-2) ],

we find

    (L^(0))^+h_0^+ =
    ζ/2ϑκ  +γϑ/2  (1-2).

By now, we have all ingredients together to go along the same route as in Appendix <ref>. 
We start with ∫ L^(2)ρ^(0)  d(,,). 
If we use 

    γ/2 κ 
    =
    (L^(0))^+
    M^-1(
    γ/2 κ )
    + 
    γζ(β+ζϑ-βϑ)/2μ(β+γϑ^2)(1-2),

we obtain 

    ∫ L^(2)ρ^(0)  d(,,) =
    _̣{[γζ(β+ζϑ-βϑ)/2μ(β+γϑ^2)(1-2)
          - (σ_1+σ_4)ζ(1-)/2] f }
    + 
    ^̣2_{ζ/4κ  f }.

Next, we turn to ∫ L^(1)ρ^(1)  d(,,)
=-∫ h_0^+ L^(1)ρ^(0)  d(,,)
with

    h_0^+ 
    = 
    ζ(+(1-2)+(β+ζϑ)(1-2)/ζ)/(2κ(β+ζϑ)).

We integrate by parts w.r.t. all derivatives 
_̣, _̣, _̣, and respectively move the derivatives 
_̣ in front of the integral by means of the chain rule, so that we obtain 

    -∫ h_0^+ L^(1)ρ^(0)  d(,,)
    = _̣(T_a+T_b+T_c+T_d,1)+^̣2_T_d,2,

with

    T_a    =   ∫ρ^(0) ζ/2κ(β+ζϑ) (
             γϑ(β+ζϑ)/2ζ  (1-2)^2^2
           + ζ^2/2ϑκ
           + γϑ/2   (1-2)^2  
                         + β+γϑ^2/2ϑκ  (1-2)
           + ζ/2ϑκ (1-2))
      d(,,),
    
    T_b    =   ∫ζρ^(0)/2κ(β+ζϑ)(
       γϑ(2β+ζϑ-ζ)/ζ  (1-2)^2
     + γϑ (μ+ζ)/ζ  (1-2)
                   + β+ζϑ-γϑ(1+ϑ)κ/ϑκ
     + ζ/ϑκ )  d(,,),
    
    T_c    =   ζ ∫(
               1/2 (σ_1+σ_4) (1-)
                -β (σ_1+σ_4(1-ζ/β) 
                     -σ_2-σ_3)/β+ζϑ  
                               (1-)
               )ρ^(0)  d(,,),
    
    T_d,1   =   ∫(
    ζ/2κ(β+ζϑ)
    (ζ-2β) (1-2)ρ^(0))  d(,,),
    
    T_d,2   =    - ∫(
    
    -ζ^2  (2(-1)-ϑ)
                      -ζβ   (4(1-)-1)
                     /2κ(β+ζϑ)) ρ^(0)   d(,,).

In particular, the integrals in T_a and T_b can be transformed 
using the operators M and Ĝ, 

    T_a    =   ζ/2κ(β+ζϑ)   ∫(
    β+ζϑ/2 (1-2)^2
             -
    2 ζ(β+ζϑ(1+ϑ))/2(β+ζϑ)   (1-)
            )
    ) ρ^(0)   d(,,),
    
    T_b    =   ζ/2κ(β+ζϑ)   ∫(
          2β-ζ
           +γκ (β+ζϑ-βϑ) (β+ζϑ)/μ(β+γϑ^2)))
                             (1-2)
     ρ^(0)   d(,,).

Recall that in the present variant of the model 
ϑ=(β-ζ)/μ, 
G=ζ/μ and Y=β/ζ.
With

    _̣τf = ∫ρ^(1)d(,,)+∫ L^(2)ρ^(0)d(,,),

we find

    _̣τf = ζ/(1+(1-1/Y)G)^2_̣x^2{(1-)f} + ζ(-σ_1-(1-1/Y)σ_4-σ_2-σ_3)/1+(1-1/Y)G_̣x{(1-)f}.





 §.§ Fixed population size and deterministic seedbank - Singular perturbation approach



Let p_i(k,l,t)=ℙ(X_t=i, Y_t∈(k,k+dk), Z_t∈(l,l+dl)). The master equations 
are given by 

    ṗ_i(k,l,t)	+∇[ ([                       β i-μ l; β (1-σ_3/N)(N-i)-μ(1+σ_2/N) k ])p_i(k,l,t)   ]
       = - [ 
    ζ i l/k+l+ζ(1+σ_1/N)(N-i)k/k+l] p_i(k,l,t)
       	+ [ζ(i+1)l/k+l] p_i+1(k,l,t)
    + [ζ (1+σ_1/N)(N-i+1)k/k+l] p_i-1(k,l,t)

The usual expansion yields the corresponding Fokker-Planck equation 
(x=i/N, y=k/N, z=l/N)

    ∂_t u    =
    ∂_x{[ 
    ζxz/y+z-ζ(1+σ_1/N)(1-x)y/y+z]u}
       	-
    ∇_y,z{[[                       β x-μ y; β (1-σ_3/N)(1-x)-μ(1+σ_2/N) z ]]u}
        +	1/2 N∂_x^2{[
    ζxz/y+z+ζ(1+σ_1/N)(1-x)y/y+z]u}




  §.§.§ Deterministic model


The drift term of the Fokker-Planck equation 
define the ODE model

    ẋ   =    -ζxz/x+z+ζ(1-x)y/y+z
    ẏ   =   β x -μ y
    ż   =   β(1-x)-μ z

with the line of stationary points

    (x,y,z) = (y,ϑ x, ϑ(1-x)),    ϑ = β/μ,   x∈[0,1].
 




  §.§.§ Dimension reduction by time scale analysis

We introduce
new coordinates, 

    x = x̃,     
    	y = ϑ+1/2  (+),   
    	z = ϑ (1-) + 1/2  (-),  ρ(t,x̃,ỹ,z̃)=u(t,x,y,z)

where, as before, ϑ= β/μ and ^2=1/N.
With

    =  x,      = ^-1(y+z-ϑ),    = ^-1(y-z+ϑ (1-2x)),

we obtain

    ∂_x 
    = ∂_-2^-1ϑ∂_,    ∂_y 
    = ^-1(∂_ +∂_),   ∂_z
    = ^-1(∂_ -∂_).

We transform the Fokker-Planck equation, neglecting terms of O(^3). 
For ρ(,,,t;) we obtain ∂_tρ=L^(0)ρ+ L^(1)ρ+^2L^(2)ρ with 

    L^(0)ρ   =    ∂_[(μ )ρ] 
    + ∂_[(μ+ζ(+(1-2)))ρ] 
    + ∂_^2[(4ζϑ^2(1-)))ρ]


    L^(1)ρ =    
    - ∂_[( ζ(+(1-2))/(2ϑ))ρ] 
    	- ∂_[( - σ_3β(1-ß) - σ_2μϑ(1-))ρ]  
       	- ∂_[(  ζ(1/ϑ)(+(1-2)) 
    - 2ζϑσ_1(1-)  ) 
    +β(σ_2+σ_3)(1-) )ρ]
       	+ ∂_^2  [(  ϑζ(+(1-2)) -4ζϑ(1-))ρ] 
    	- ∂_∂_[( 
    4ϑζ(1-) )ρ] 
    
    L^(2)ρ =    ∂_[( 
    ζ (+(1-2))/(2ϑ^2) - σ_1ζ(1-) )ρ] 
    +  ∂_^2 [( ζ(1-) )ρ]
       	+ ∂_[(⋯)ρ]
    + ∂_[(⋯)ρ]

As before, we introduce τ = ε^2t=t/N, expand ρ w.r.t. , 


    ρ(,,,t) = ρ^(0)(,,,ε^2t) + ερ^(1)(,,,ε^2t) + ε^2ρ^(2)(,,,ε^2t)+𝒪(ε^3)


and obtain


    L^(0)ρ^(0)=0,   L^(0)ρ^(1)=-L^(1)ρ^(0),   L^(0)ρ^(2)=_̣τρ^(0)-L^(1)ρ^(1)-L^(2)ρ^(0).

The reduced Fokker-Planck equation is given by

    _̣τf = ∫ρ^(1)d(,) + ∫ L^(2)d̊(,).

In the following computations we use that 

    L^(0)+[+(1-2)^2))]
       =   
    -(2μ+ζ)(
    +(1-2)
    )
    
    
    
    L^(0)+(+(1-2)))
       =    -(μ+ζ) (+(1-2)))
    
    
    L^(0)+[(+(1-2))^2]    =    
    -2(μ+ζ)(+(1-2))^2+8ζϑ^2(1-)

Therewith we find

    ∫ L^(2)ρ^(0) d(,zt)
       =   ζ/2ϑ^2 ∂_ ∫ (+(1-2)) ρ^(0)  d(,)
    - σ_1ζ∂_[((1-) )f(,τ)] 
          	+  ∂_^2 [( ζ(1-) )f(,τ)]
       =   
    - σ_1ζ∂_[((1-) )f(,τ)] 
    +  ∂_^2 [( ζ(1-) )f(,τ)].

Furthermore, 

    ∫ L^(1)ρ^(1)  d(,)
       =   -ζ/2ϑ∂_∫  (+(1-2)) ρ^(1)  d(,)
       =   ζ/2ϑ(μ+ζ)∂_∫  L^(0)+(+(1-2)) ρ^(1)  d(,)
       =   -ζ/2ϑ(μ+ζ)∂_∫  (+(1-2)) L^(1)ρ^(0)  d(,)
       =   -ζ/2ϑ(μ+ζ)∂_(
    T_1+T_2+T_3+T_4+T_5
    )

where T_1,…,T_5 are given by

    T_1 
       =   ∫  (+(1-2)) (- ∂_)[(
    ζ(+(1-2))/(2ϑ)
    )ρ^(0) ]  d(,)
       =    - ζ/2ϑ[∂_x ∫ (+(1-2))^2ρ^(0)  d(,) +2 ∫  (+(1-2))ρ^(0)  d(,)]
       =    - ∂_ζ/2ϑ   4ϑ^2ζ(1-)/ζ+μ   f
    =   -2ϑζ^2 /ζ+μ  ∂_((1-) f)
    
    
    
    T_2    =   ∫  (+(1-2)) (- ∂_)(
    - σ_3β(1-ß) - σ_2μϑ(1-)
    )ρ^(0)  d(,)
    
    
    
       =    (1-2)(
    - (σ_2+σ_3)β(1-ß) 
    )  f	
    
    
    
    T_3    =   ∫  (+(1-2)) (	- ∂_)(
    ζ(1/ϑ)(+(1-2)) 
    - 2ζϑσ_1(1-) +β(σ_2+σ_3)(1-)
    )ρ^(0)  d(,)
    
    
    
    
    
    
    
       =    - 2ζϑσ_1(1-)  f +β(σ_2+σ_3)(1-)f
    
    
    
    T_4    =   ∫  (+(1-2)) (∂_^2)(
    ϑζ(+(1-2)) -4ζϑ(1-))ρ^(0)  d(,)	= 0
    
    
    
    T_5 
    
    
    
       =   ∫ ∂_(
    4ϑζ(1-) 
    )ρ^(0)  d(,) 
    =   ∂_(
    4ϑζ(1-) f
    )

Adding up the corresponding terms yields the reduced Fokker-Planck equation

    f_τ   =   
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    -ζ (σ_1+σ_2+σ_3)/1+ζ/μ ∂_[(1-)f]
    + 
    ζ/(1+ζ/μ)^2 ∂_^2[(1-)f].













