Resident-invader dynamics of similar strategies in fluctuating environments

We study resident-invader dynamics in fluctuating environments when the invader and the resident have close but distinct strategies. First we focus on a class of continuous-time models of unstructured populations of multi-dimensional strategies, which incorporates environmental feedback and environmental stochasticity. Then we generalize our results to a class of structured population models. We classify the generic population dynamical outcomes of an invasion event when the resident population in a given environment is non-growing on the long-run and stochastically persistent. Our approach is based on the series expansion of a model with respect to the small strategy difference, and on the analysis of a stochastic fast-slow system induced by time-scale separation. Theoretical and numerical analyses show that the total size of the resident and invader population varies stochastically and dramatically in time, while the relative size of the invader population changes slowly and asymptotically in time. Thereby the classification is based on the asymptotic behavior of the relative population size, and which is shown to be fully determined by invasion criteria (i.e., without having to study the full generic dynamical system). Our results extend and generalize previous results for a stable resident equilibrium (particularly, Geritz in J Math Biol 50(1):67–82, 2005; Dercole and Geritz in J Theor Biol 394:231-254, 2016) to non-equilibrium resident population dynamics as well as resident dynamics with stochastic (or deterministic) drivers.


Introduction
Two important issues in the framework of adaptive dynamics are which mutant strategies can invade a population of given resident strategies, and what would be the population dynamical outcomes of an invasion event. The long-term growth rate of a newly arrived and initially rare mutant with strategy y in the environment generated by a population of resident strategy x (i.e., invasion fitness of mutant y in resident x) determines whether an invasion event may occur or not (Metz et al. 1992(Metz et al. , 1996Dieckmann and Law 1996;Geritz et al. 1997Geritz et al. , 1998. If y can invade x but not vice versa, does it mean that y will eventually take over x and becomes the new resident? In general, this need not happen as examples of unprotected coexistence and the "resident strikes back" phenomenon show (see e.g., Doebeli 1998;Parvinen 1999;Mylius and Diekmann 2001;. However, if y is close to, but not identical to x, Geritz (2005) and Dercole and Geritz (2016) have shown that invasion without back-invasion generically implies substitution, and mutual invasion generically leads to coexistence in a broad class of population models. They found that the generic population dynamical outcomes of an invasion event of similar strategies are determined by invasion criteria alone, i.e., without having to study the full generic dynamical system.
The focus of Geritz (2005) and Dercole and Geritz (2016) is on resident-invader dynamics of similar strategies in a constant environment, providing that the resident population in a large ecological community is at a hyperbolic attracting steady state. Recently, there has been considerable interest in questions related to evolution of phenotypic diversity in fluctuating environments (Kussell and Leibler 2005;Kisdi and Liu 2006;Geritz et al. 2007;Schreiber 2012b;Wakano and Iwasa 2013;Ripa and Dieckmann 2013;Melbinger and Vergassola 2015;Saether and Engen 2015;Ferris and Best 2019). The purpose of this paper is to extend and generalize the results of Geritz (2005) and Dercole and Geritz (2016) to environmental fluctuations due to nonequilibrium population dynamics (e.g., cycle, quasiperiodic trajectory and chaos) or environmental stochasticity. This requires (i) considering an explicit model that takes account of environmental feedback and environmental stochasticity, and (ii) specifying the population dynamics of residents in fluctuating environments.
A population affects its environment, and the environment in turn affects the population. Such an environmental feedback loop characterizes the interaction of populations with their environments and plays a central role in their ecological and evolutionary dynamics (Metz et al. 1996Meszéna and Metz 1999;Kisdi and Geritz 2016;Lion 2018b). Examples of feedback variables are: interacting species (e.g., predator, prey and food); different classes in a population who is structured by age, sex, spatial location, etc; and physical variables (e.g., temperature and humidity) which have been shown that they may provide a feedback between ecological communities and local physical patterns in some concrete examples (see the Introduction of Benaïm and Schreiber (2019) and references therein). These feedback variables need to be modelled with their own dynamics. Environmental stochasticity, however, generally refers to effects of fluctuations in external factors (which can be biotic or abiotic) on model parameters (Nisbet and Gurney 1982;Kliemann 1983;Chesson 1986;Tuljapurkar 1990). These external factors drive the dynamics of the aforementioned feedback variables but not vice versa. This is what differentiates environmental stochasticity from environmental feedback.
A population of a resident strategy in a fluctuating environment generally requires to be non-growing on the long-run and stochastically persistent. The non-growing property means that the long-term growth rate of the resident population in a given environment is zero. Stochastic persistence is an important concept to characterize population dynamics in fluctuating environments, which captures the property that if initially present, even if only in small size, then it will be present over arbitrary long periods of time in a stochastic sense. There are serval different (not necessarily equivalent) definitions of stochastic persistence (see e.g., a review in Schreiber 2012a). In this paper, we use the definition of Schreiber et al. (2011) and Benaïm (2018) which asserts that the probability of a population being near extinction is arbitrarily small.
The rest of this paper is organized as follows. In Sect. 2, we focus on the analysis of resident-invader dynamics of a class of unstructured population models. In Sect. 2.1, a class of polymorphic unstructured population models for multi-dimensional strategies and an explicit formulation of environmental feedback and environmental stochasticity is presented. In Sect. 2.2, we review some basic concepts used in this paper. In Sect. 2.3, we specify when a population of a given strategy is a resident. Section 2.4 gives the definition of invasion fitness. In Sect. 2.5, we present the basic analysis of invasion dynamics of similar strategies. Particularly, in Sects. 2.5.2 and 2.5.3, we classify the generic population dynamical outcomes of an invasion event. In Sect. 3, we illustrate how to generalize the results of the unstructured population models to a class of structured population models. In Sect. 4, we apply our results to examples of evolving bacteria in a chemostat, Lotka-Volterra competition, structured SIRS epidemic dynamics, and the evolution of timidity of the prey in a prey-predator model. The first two examples are designed to highlight different ways that our results can be used. The third example is designed to illustrate how a concrete structured population model can be reformulated into our framework. The last example is designed to demonstrate that our results are applicable to the evolution with the non-equilibrium resident dynamics. In Sect. 5, we conclude with discussing how our results relate to the existing literature and what possible generalizations are. Proofs are given in the "Appendix".
Before proceeding with the analysis, we introduce some default notations used in this paper. Let X be the strategy space, N be the space of non-negative population sizes, E be the space of feedback variables and Θ be the space of external factors. Each of these spaces is assumed to be a subset of a normed vector space. f (i) is the ith derivative of function f with respect to its first argument, and f (i, j) is the ith derivative of function f with respect to its first argument and the jth derivative with respect to its second argument. The superscript means transpose of a vector or a matrix.

The polymorphic population model
Consider a population of strategies x 1 , . . . , x k ∈ X whose sizes at time t are given by n 1,t , . . . , n k,t ∈ N , respectively. The population of x 1 , . . . , x k interact with e t ∈ E that describes the feedback environment. Generally, e t may correspond to interacting species, different classes of the population, physical factors, or any combination thereof. n 1,t , . . . , n k,t as well as e t may be influenced by a stochastic driver θ t ∈ Θ that describes external factors. The per-capita growth rate f of individuals with strategy x i depends on the current environmental condition (e t , θ t ). Under these considerations, the dynamics of the polymorphic unstructured population are given bẏ where the dot denotes differentiation with respect to time variable t. Notice that for a given environmental condition (e t , θ t ), (2.1a) is linear in the population size. All the nonlinearity of (2.1a) comes from how the environment e t depends on the population sizes.
Before introducing the formulation of the dynamics of e t , let us look at how the environment e t depends on the population sizes in the following two deterministic models.
Example 1 Lehtinen and Geritz (2019) studied the evolution of timidity in a prey species whose predator has cannibalistic tendencies, in which the population dynamics is given by where the reproduction of prey i is limited by competition (of resources, territories, breeding sites, etc.) among the foraging prey n F i,t = n i,t 1+x i p t with the per-capita birth rate g j n F j,t = a − c j n F j,t , a searching adult predator p S t = p t 1+βh j n F j,t captures the foraging prey at the rate β, the maturation of juvenile predators J t = γβ p S t j n F j,t σ +(1−λ)α p S t with the mean time T depends on the cannibalistic pressure they experience throughout their juvenile period, and the per-capita death rate of the adult predator and the predation-independent per-capita death rate of the prey are δ and μ, respectively. The prey strategies differ only in their level of timidity x i ∈ R + := X . Other parameters are all positive constant. Let (adult predators) e 1,t = p t , (foraging preys) e 2,t = j n F j,t = j 1 1+x j e 1,t n j,t , Since e 1,t is an implicit function of the prey sizes and e 2,t is an explicit function of e 1,t , the environmental variable e 2,t becomes nonlinear in the prey sizes. Thus, the per-capita population growth rate f of prey type i in the nonlinear environment (e 1,t , e 2,t ) := e t is given by Example 2  and Dercole (2003) studied the evolution of cannibalistic strategies, in which the interactions between k cannibalistic consumer sub-populations with sizes n i,t and strategies x i ∈ R + := X , i = 1, . . . , k, are described byṅ where the conversion factor γ , the competition coefficient c, and the size of a common resource n 0 are assumed to be positive constant, while the attack rate a (and a 0 ) and the handling time h (and h 0 ) depend upon the strategies. In the bracket of the equation, the first term is reproduction through the harvesting of the common resource and cannibalism toward the sub-population j, the second term is mortality due to cannibalism, and the last term is mortality due to competition. For a given strategy x i , let Here the environmental variables e 21,t , e 22,t and e 24,t are all linear in the population sizes, but the environmental variable e 23,t is a nonlinear and explicit function of e 22,t which leads to the nonlinearity of e 23,t in the population sizes. Thus, the per-capita population growth rate f of sub-population i with a given strategy x i in the nonlinear environment (e 21,t , e 22,t , e 23,t , e 24,t ) := e t is given by 24,t (N.B., e 21,t , e 22,t and e 23,t are functions of strategy x i and the strategy of subpopulation type i has infinitely many choices, which causes an infinite-dimensional environment e t to appear. We will provide more interpretations for the case of infinitedimensional environment below).
Following Diekmann et al. (2001) and Geritz (2005), the environment e t is assumed to be a linear function of the population sizes where I (x j ) describes the environmental impact of a single individual with strategy x j . However, this assumption excludes many interesting models with nonlinear environments such as the previous two examples. In Example 1, the predator size necessarily is an environmental variable for the preys, and the fraction of time individual prey spend foraging depends on the adult predators. In Example 2, the harvested food, the total handling times, and the compete-caused and cannibalize-caused death are environmental variables for the sub-population i, in which the cannibalize-caused death depends on the total handling times. We note that different environmental variables may have different forms even in the same model. For instance, in Example 1, e 1,t is govern by a differential equation but e 2,t is a function of e 1,t . For these reasons, we describe the environment by implicit (differential) equations.
In this paper, we consider an explicit formulation of the dynamics of the feedback environment e t = (e 1,t , e 2,t ) as illustrated in Fig. 1, which is given bẏ (2.1b) The dynamics of e 1,t is given by a differential equation, while e 2,t is an implicit function of the current state of the system. The functions G 1 and G 2 describe the intrinsic dynamics of the virgin environment (the environment unaffected by the population of x 1 , . . . , x k ). Given the current environment, the functions H 1 and H 2 describe the environmental impact of a single individual with a given strategy. Consequently, the second terms of (2.1b) gives the total environmental impact of the entire population of x 1 , . . . , x k . These functions all involve the current environmental condition (e t , θ t ).
Here we would like to highlight that e 2,t is assumed to be unique defined by the implicit equation. As Example 2 and "Appendix A" show, H 2 can be a function of e 2,t , but e 2,t is still unique defined. Alternatively, by allowing e t to be infinite-dimensional, we can also incorporate it into (2.1b) (Geritz 2005). For instance, in Example 2 and the model in Sect. 4.2, where the feedback variables at any time are a function of strategy and where f (x i , e t , θ t ) involves the evaluation of e t in x i (ref to "Appendix A" and Sect. 4.2, respectively).
To consider external factors, we assume that {θ t } t≥0 is a continuous time Markov process which is defining on an underlying probability space and taking value in Θ. From the definition of Markov process, the future state is independent of the past and depends only the present, which is a natural consideration in reality and contains a large number of ecologically interpretable factors. There are many ways to construct a Markov process, e.g., a pure jump process and a process generated by stochastic differential equations. In this paper, we assume that the process {θ t } t≥0 is governed by the following stochastic differential equation: (2.1c) whereẆ t is a white noise, namely W t is a Brownian motion defined on a complete probability space (Ω, F , P) with a filtration {F t } t≥0 satisfying the usual conditions (i.e., it is right continuous and increasing while F 0 contains all P-null sets). (2.1c) can be calculated using one of two methods: Itô or Stratonovich, but in fact our results do not depend on the calculation methods used. Model (2.1) is an extension of the deterministic population models used in Geritz (2005), Dercole and Rinaldi (2008) and Dercole and Geritz (2016), which is intuitive and interpretable in biology and covers a large class of ecological models with an explicit formulation of environmental feedback and environmental stochasticity. Meanwhile, model (2.1) can be viewed as a revision of the stochastic population model used in Kliemann (1983), but now takes account of strategies and feedback variables. From the application point of view, any complicated considerations of biological mechanisms can be formulated into functions f , G 1 , H 1 , G 2 , H 2 , A and B, and a certain level of smoothness is assumed to make them mathematically meaningful.
Regarding (2.1), we make the following assumptions:

A1
The Markov process {θ t } t≥0 is ergodic with invariant probability measure ν on Θ.
A2 The partial derivatives Assumption A1 indicates the process {θ t } t≥0 is completely realized before we start to observe the dynamics of (n 1,t , . . . , n k,t , e t ). Both meanings of ergodicity and invariance are given in Sect. 2.2. Process {n 1,t , . . . , n k,t , e t } t≥0 is usually not Markovian. By Assumption A1, however, process {n 1,t , . . . , n k,t , e t , θ t } t≥0 is Markovian under weak assumptions (ref to Arnold and Kliemann (1983, Lemma 2.1)). Thereby it allows us to utilize the excellent works in the theory of Markov processes for the analysis of the process {n 1,t , . . . , n k,t , e t , θ t } t≥0 . Assumption A2 is a technical requirement to guarantee the existence of chain rule and justify the series expansion.

Some basic concepts
In order to state our results, we review some basic concepts in Markov process and ergodic theory.
Throughout Sect. 2, let which lies in space For any Borel set C ⊂ Z, let be the associated semigroup of Markov process {Z t } t≥0 on Z defined by conditional expectations: for bounded and measurable function h : Z → R. Assumption A2 implies that the Markov semigroup {T t } t≥0 is C b (Z)-Feller (Benaïm 2018), i.e., operator T t takes bounded continuous functions h on Z to bounded continuous functions for all bounded and measurable functions h : Z → R and all t ≥ 0. An invariant probability measure μ is ergodic if it can not be written as a non-trivial convex combination of other distinct invariant probability measures. If Z t initially follows the distribution of an ergodic probability measure μ, then Birkhoff's ergodic theorem tells us that, in simple terms, the time average of a function of the process along the trajectories exists with probability one and equals the space average. More precisely, for all measurable functions h : with probability one.

Resident and resident system
From (2.1a) with initial value Z 0 = z ∈ Z, a straightforward calculation shows that the species x i tends to increase in its population size if f (x i , e s , θ s )ds < 0, provided that the limit exists. Let μ be an invariant probability measure for process with probability one and for μ-almost every z, where λ x i (μ) is the per-capita growth rate of species x i with respect to μ: We now introduce two notations for further analysis. We define the extinction set of a population of strategies x 1 , . . . , x k Z 0 = z ∈ Z : min i n i = 0 is the set which at least one species is absent, and the η-neighborhood of the extinction set is the set which at least one species has a size less than η.
Next, we specify when a population of strategies x 1 , . . . , x k is a resident.
In Definition 1, the supports of required invariant probability measures should exclude the extinction set Z 0 . The first reason is that it is the technical requirement to guarantee the properties of non-growing and stochastically persistence. The second reason is from a probabilistic perspective that if a population initially follows an invariant probability measure μ with μ(Z 0 ) = 0, then the extinction of the population is a non-zero probabilistic event on the long-term. From Definition 1, a resident is nongrowing on the long-run and stochastically persistent, where the concept of persistence asserts that the probability of population size being near the extinction state is very small. In this paper, a resident (resp. a resident system) is characterized by an invariant probability measure defined on the state space. Actually, there might be more than one invariant probability measures, but we focus our attention on one of them. Each adequate invariant probability measure corresponds to a resident (resp. a resident system).

Invasion fitness
Assume that (x 1 , . . . , x k , n 1,t , . . . , n k,t , e t , θ t ) is a resident system with an invariant probability measure μ satisfying Definition 1, which further determines the ecological environment for a mutant type y ∈ X with initially infinitesimal size m t ∈ N . Thus, the population dynamics of a newly arrived and initially rare mutant y is described by the following linear differential equatioṅ (2.2) Notice that (2.2) holds as long as the mutant population remains rare. The tuple (y, m t ) is referred to an invader. Figure 1 shows the relationship of the invader (y, m t ) and the resident system (x 1 , . . . , x k , n 1,t , . . . , n k,t , e t , θ t ) where arrows with special indexes illustrate different impacts or contributions between them.
Formally, the invasion fitness of mutant y in resident (x 1 , . . . , x k ) is given by with probability one and for μ-almost every z and for every positive size m 0 (Metz et al. 1992;Ferriere and Gatto 1995). The invader y dies out if S x 1 ,...,x k (y) < 0 and can spread if S x 1 ,...,x k (y) > 0. From Definition 1, if y = x i for any i ∈ {1, . . . , k}, then we have which is so-called the principle of selective neutrality of residents in the framework of adaptive dynamics.

Population dynamics of similar strategies
Consider now a population of two similar strategies (i.e., two close but distinct strategies) x 1 , x 2 ∈ X with corresponding sizes n 1,t , n 2,t ∈ N at time t. Without loss of generality, let strategies where ξ 1 = ξ 2 with ξ 2 − ξ 1 < +∞ and where small > 0 (N.B., here is different from the ε used in Definition 1) and where x is a reference strategy. Biologically, can be viewed as the intensity of mutation, and ξ i is the direction of mutation with respect to the reference strategy x. Model (2.1) applied to the case of k = 2 we further have a dimorphic population model. Obviously, n 1,t , n 2,t and e t are functions of x 1 and/or x 2 . Thus they are functions of , which are denoted by n 1,t , n 2,t and e t , respectively. To study the dimorphic population dynamics of x 1 and x 2 in fluctuating environments, following Meszéna et al. (2005) and Dercole and Geritz (2016), we introduce the total population size N t of x 1 and x 2 at time t and the relative population size P t of x 2 at time t, i.e., Under these assumptions, we arrive at the following equations: (2.5e) Next, we unfold the dimorphic population dynamics of x 1 and x 2 based on the series expansion of model (2.5) in powers of , and on time-scale separation. In the following analysis, we will see that (N t , e t , θ t ) are fast variables with dynamics acted on ttimescale and P t is a slow variable with dynamics expounded on i t-timescales for i = 1, 2, . . . .

Fast dynamics
(2.6) A comparison of (2.6) and the dynamics of the reference monomorphic system (x, n t , e t , θ t ) shows that (N 0 t , e 0 t , θ t ) has the same dynamics as (n t , e t , θ t ). Assume that (x, n t , e t , θ t ) is a resident system with associated invariant probability measures μ. Then, in the limit = 0, (N t , e t , θ t ) enter the dynamics of the observables through μ, and P t becomes irrelevant. Thus, (2.5) is a stochastic fast-slow system with fast variables (N t , e t , θ t ) and slow variable P t .
For positive but small , the trajectories of (N t , P t , e t , θ t ) can be viewed as a small variation of that of (N 0 t , P 0 t , e 0 t , θ t ) on t-timescale. Under some mild assumptions, the dynamics of (2.6) are arbitrarily good approximation of the dynamics of (2.5) in terms of conditional expectations when → 0 (ref to "Appendix B").

Slow dynamics on t-timescale
To study the dynamics of the relative population size on slow timescales, we have to take account of the high-order terms in the series expansion.
Consider the first-order term in the series expansion of the right side of (2.5b) and let t 1 = t, we havė (2.7) Here the term e 0 t 1 / in the series expansion has been replaced by e t 1 / because they have the same dynamics on t-timescale (see the interpretation in Sect. 2.5.1). Notice that f (1,0) (x, e, θ) becomes a row vector if x is a multi-dimensional strategy. Since (x, n t , e t , θ t ) is assumed to be a resident system with associated invariant probability measures μ satisfying μ(Z \ Z 0 ) = 1, we can apply Birkhoff's ergodic theorem to an observation f (1,0) of (x, n t , e t , θ t ), i.e., we have with probability one and for μ-almost Metz et al. 1996;Geritz et al. 1997Geritz et al. , 1998Geritz et al. , 1999.
Throughout this section, assuming that Consider the following deterministic systeṁ From the averaging principle for fast-slow systems with stochastic and stationary fast dynamics, solutions of (2.7) are well approximated by solutions of (2.9) on every finite time interval when → 0, where (2.9) is so-called the averaged system with respect to (2.7) (ref to Freidlin and Wentzell (2012 Figure 2a gives all generic cases. Under the stability of the averaged system (2.9), we show the property of solutions of (2.7) in arbitrarily long time intervals.
In words, for > 0 sufficiently small and at some later time, the trajectories of P t 1 will remain δ-close to the stable equilibrium of the averaged system (2.9) for arbitrarily long time with probability one.

Remark 1
In the hypotheses of Lemma 1, the boundedness ensures that the invasion fitness S x (y) for a mutant strategy y close to the reference resident strategy x and the selection gradient ∂ y S x (y)| y=x are well defined, and it further guarantees the approximation of solutions of (2.7) and (2.9) in arbitrarily long time intervales.
Expansions of invasion fitnesses S x 1 (x 2 ) and S x 2 (x 1 ) up to the first order terms in give Using (2.11), a stochastic analogy of the invasion implies fixation theorem of Geritz (2005, Proposition 1) immediately follows from Lemma 1, which is called as invasion implies substitution theorem in the present paper.
Theorem 1 ("Invasion implies substitution theorem") Assume that the hypotheses of Lemma 1 are satisfied. Let the strategy pair (x 1 , x 2 ) and the reference strategy x are For X ⊂ R, provided that ∂ y S x (y)| y=x = 0 (i.e., x isn't an evolutionarily singular strategy), the dimorphic population dynamics of x 1 and x 2 is completely determined by the invasion fitnesses. Since there are only two equilibria that are at the boundary for the averaged system (2.9), invasion implies substitution in terms that the dynamics of the relative population size satisfies (2.10). In other words, similar strategies cannot coexist when both of them are away from an evolutionarily singular strategy.
For X ⊂ R d with d ≥ 2, geometrically, the deviation ξ 2 − ξ 1 (which is equivalent to x 2 − x 1 ) shall be away from the null-set of the linear form α → ∂ y S x (y)| y=x α for all α = 0 and ∂ y S x (y)| y=x non-vanishes, so that (2.9) is non-degenerate (i.e., the right side isn't equal to zero). Consequently, in a multi-dimensional strategy space, invasion implies substitution provided that not only the two strategies are away from an evolutionarily singular strategy, but also their deviation isn't orthogonal to the selection gradient at the reference strategy.
The averaged system (2.9) gives a geometric interpretation of the population dynamics of similar strategies in terms of their deviation and the selection gradient at a reference strategy, which indicates the sufficient condition of substitution of similar strategies. Thus, on t 1 -timescale, any invader that can spread in a given resident system will ultimately oust the resident.

Slow dynamics on 2 t-timescale
Once ∂ y S x (y)| y=x (ξ 2 − ξ 1 ) = 0, the right side of the averaged system (2.9) becomes zero. Then the high-order term O( ) in (2.7) plays a role on slower timescales.
Consider the second-order term in the series expansion of the right side of (2.5b) and let t 2 = 2 t, we have the associated averaged system: with matrixes For the reader's convenience, we give the detailed derivation of (2.12) in "Appendix D".
Throughout this section, assuming that don't vanish at the same time. If (ξ 2 − ξ 1 ) (C 22 + C 11 )(ξ 2 − ξ 1 ) = 0, then the averaged system (2.12) has an interior equilibrium. In case of (ξ 2 −ξ 1 ) (C 22 +C 11 )(ξ 2 −ξ 1 ) = 0, only boundary equilibria exist. Figure 2b and c conclude all possible equilibria and their stability. From the series expansions of invasion fitnesses S x 1 (x 2 ) and S x 2 (x 1 ), the secondorder terms in give (2.14) where the derivation of the last equality can be found in "Appendix D". Using (2.14), the averaged system (2.12) have at most three equilibria in [0, 1]: with stabilities dominated by the combination of signs of S x 1 (x 2 ) and S x 2 (x 1 ). Under the stability of the averaged system (2.12), we have a similar consequence as Lemma 1.
Lemma 2 Let P 0 =P 0 = p and Z 0 = z. Assume that for a small σ > 0, Then for every δ ∈ (0, 1), there exists a constant T δ > 0 such that (i) if the averaged system (2.12) has a unique stable equilibrium, then the solution of (2.12) satisfies ii) if the averaged system (2.12) has two stable equilibria (i.e., 0 and 1), then we have In words, for > 0 sufficiently small and at some later time, the trajectories of P t 2 will remain δ-close to the attractor (either a unique stable equilibrium or the bistable equilibria) of the averaged system (2.12) for arbitrarily long time with probability one.

Remark 2
In the hypotheses of Lemma 2, the boundedness ensures that (2.13) are well defined, and it further guarantees the approximation of solutions of P t 2 andP t 2 in arbitrarily long time intervals. The boundedness is a technical requirement for Lemmas 1 and 2, which meets by many models.
Under the stability of the averaged system (2.12) with equivalent form (2.14), a stochastic analogy of the classification theorem of Geritz (2005, Proposition 2) immediately follows from Lemma 2.
Theorem 3 ("Special invasion-substitution theorem") Assume that the hypotheses of Lemma 2 are satisfied. Let the strategy pair (x 1 , x 2 ) and the reference strategy x are For X ⊂ R, Theorem 2 corresponds to the invasion dynamics of similar strategies close to a generic evolutionarily singular strategy that satisfies C 22 + C 11 = 0, while Theorem 3 corresponds to the invasion dynamics of similar strategies with ξ 1 + ξ 2 = 0 close to an evolutionarily singular strategy satisfying C 22 + C 11 = 0 and C 22 + C 21 = 0. In the neighborhood of a generic evolutionarily singular strategy, the invasion dynamics is essentially "Lotka-Volterra"-dominance of one strategy, coexistence of two strategies, and mutual exclusion of two strategies-in sense that the relative Fig. 3 The geometric interpretation of the coexistence condition in a one-dimensional (or a one-dimensional parameterization of a multi-dimensional) strategy space. The "+" in the Pairwise Invasibility Plot (PIP) indicates the area where the mutant can invade the resident, and the "−" indicates the area where the mutant cannot invade population size remains δ-close to the attractor of the averaged system for arbitrarily long time with probability one. In the neighborhood of an evolutionarily singular strategy satisfying C 22 + C 11 = 0 and C 22 + C 21 = 0, only substitution happens for the invasion dynamics of any two similar but non-opposite strategies (i.e., ξ 1 + ξ 2 = 0). In Fig. 3, the dashed line indicates the scenario C 22 + C 11 = 0 associated with the local configuration of the Pairwise Invasibility Plot (Matsuda 1985;van Tienderen and de Jong 1986). Invasion implies substitution if (C 11 , C 22 ) is on the dash line. Coexistence of similar strategies may occur if (C 11 , C 22 ) is away from the dash line. For The correspondingly geometric interpretation is complicate. However, if let ξ 1 = 0 (i.e., x 1 = x), the two hypotheses are equivalent to ξ 2 (C 22 + C 11 )ξ 2 = 0 and ξ 2 (C 22 + C 21 )ξ 2 = 0. From a conservation law C 11 + C 22 + C 21 + C 12 = 0 (ref to (d.9) in "Appendix D"), we have is symmetric and C 21 − C 12 is skew-symmetric. Then by the property of skew-symmetric matrices, , the higher order terms in of the series expansion are imperative to reveal the invasion dynamics of similar strategies. The degeneracy of a higher degree opens the possibility of unprotected coexistence for two strategies (ref to Priklopil (2012) and Dercole and Geritz (2016) for population models with point equilibria).
Consider a structured population of strategies x 1 , . . . , x k ∈ X with corresponding sizes n 1,t , . . . , n k,t ∈ N at time t where n i,t = (n 1 i,t , . . . , n i,t ) for ≥ 2 and for i = 1, . . . , k. Under the same considerations in Sect. 2.1, our polymorphic structured population model isṅ where F(x i , e t , θ t ) is a × matrix whose entries are transition rates between different states of individuals with strategy x i and depending on the current environmental condition (e t , θ t ). Model (3.1) is an extension of the deterministic models used in Durinx et al. (2008) and Priklopil and Lehmann (2019), but now with an explicit formulation of environmental feedback and environmental stochasticity.
To place (3.1) into our framework, let n i,t 1 be the total population size of x i , and let v i,t be the proportion of each component in the population of x i , i.e., for i = 1, . . . , k, From (3.1), we obtain that for i = 1, . . . , k, where 1 = (1, . . . 1) and I is the identity matrix. Although (3.2) cannot be written in the form of (2.1), the results in Sect. 2 still can be generalized to model (3.2) where population structures v 1,t , . . . , v k,t are viewed as auxiliary variables that lie in the space Δ = v ∈ [0, 1] : v 1 = 1 . Throughout this section, denote Z t = n 1,t 1 , . . . , n k,t 1 , v 1,t , . . . , v k,t , e t , θ t . The state space of the dynamics of Z t is Z = R k + × Δ k × E × Θ, the extinction set is Z 0 = z ∈ Z : min i n i 1 = 0 , and the η-neighborhood of the extinction set is Z η = z ∈ Z : min i n i 1 ≤ η . Like Assumption A2, the partial derivative F (i, j) requires to exit and to be locally Lipschitz continuous and measurable in (e, θ) for all strategies and for all nonnegative integers i, j ≤ 2.
Next, we focus on the analysis of dimorphic population dynamics of similar strategies. Like the discussion in Sect. 2.5, let x 1 = x + ξ 1 , x 2 = x + ξ 2 ∈ X with corresponding population sizes n 1,t , n 2,t ∈ N and population structures v 1,t , v 2,t ∈ Δ at time t. Let x, n t 1 , v t , e t , θ t is the reference monomorphic resident system with invariant probability measure μ. We now are interested in the dynamics of the total population size N t of x 1 and x 2 and the relative population size P t of x 2 , where N t = n 1,t 1 + n 2,t 1 , P t = n 2,t 1 N t .
Expanding (3.2) around = 0 (i.e., x 1 = x 2 = x) and then taking = 0, we havė Recalling the analysis in Sect. 2.5.1, the difference here is that P 0 t is not longer a constant on t-timescale. Instead, P 0 t evolves with N 0 t on t-timescale. The following we derive that (i) the dynamics of the monomorphic system n t 1 , v t , e t , θ t still represents the dynamics of (N t , v i,t , e t , θ t ) on t-timescale when = 0 for i = 1, 2, and (ii) P t evolves on t-timescale but eventually returns to its initial state when = 0. Denotev which is the averaged proportion of population structures v 0 1,t and v 0 2,t at time t. From (3.3a) and the dynamics of (v 0 1,t , v 0 2,t , e 0 t ), we havė (3.4) It shows that (N 0 t ,v 0 t , e 0 t , θ t ) satisfies the same equations as n t 1 , v t , e t , θ t . In words, the dynamics of n t 1 , v t , e t , θ t represents that of (N t ,v t , e t , θ t ) on ttimescale when = 0. Turning to the dynamics of v 0 1,t and v 0 2,t , with e 0 v t ). Thus, we claim that (N t , v i,t , e t , θ t ) has the same dynamics as n t 1 , v t , e t , θ t on t-timescale when = 0 for i = 1, 2. Now we return to (3.3b), a straightforward calculation shows that provided that the initial state is (n, v 1 , v 2 , e, θ, p). Using Birkhoff's ergodic theorem and the principle of selective neutrality of residents (i.e., zero growth rate of the reference resident x), we obtain lim t→+∞ P 0 t = p with probability one and for all p ∈ [0, 1] and for μ-almost every (n, v 1 , v 2 , e, θ) (N.B., initial states v 1 and v 2 do not have to be identical). It means that P 0 t will sufficiently close to its initial state some time later with probability one. Thus, we need consider the dynamics of P t on slow timescales to revel the dimorphic population dynamics of x 1 and x 2 .
The remaining generalization of our theorems to model (3.2) can be done by using of arguments similar to Sects. 2.5.2 and 2.5.3.

Applications
To illustrate the utility of our theorems, we apply them to four concrete examples. The first example we consider the evolving bacteria in a chemostat model. We will see that substitution is the unique outcome after an invasion event. The second example we consider a Lotka-Volterra competition model. Theorems 1, 2 and 3 are applied to predict the population dynamical outcomes of an invasion event when the strategy pairs satisfy the corresponding hypotheses of these theorems. The third example we consider a structured SIRS epidemic model, which is designed to illustrate how a concrete structured population model can be reformulated into our framework. In the last example, we consider the evolution of timidity of the prey in a predator-prey model, which is designed to demonstrate that our results are applicable to the evolution with the non-equilibrium resident dynamics.
In these four examples, we start from a monomorphic population with a certain strategy. To verify the existence of an invariant probability measure for a resident system, one might consider the tightness of one of the following two families of measures: (a) the statistics associated with a single realization of process {Z t } t≥0 , i.e. measure with a Dirac measure δ Z s at Z s (ref to the argument presented in Schreiber et al. (2011) or Benaïm (2018)); (b) the statistics associated with transition probabilities {P t } t≥0 of process {Z t } t≥0 , i.e. measure which is a common tool in the literature on Markov process (see e.g. Arnold and Kliemann 1983, Proposition 3.15).
Different ways to study the existence of invariant probability measures give different interpretations for the asymptotic behavior of process {Z t } t≥0 . Assume that the limit exists as time goes to infinity, the (a) gives the long-term frequency of a typical realization of process {Z t } t≥0 visiting a particular configuration. Since a transition probability corresponds to the frequency of observing a particular event across many realizations of process {Z t } t≥0 , the (b) gives the long-term frequency of probabilities for process visiting a particular configuration. Following the strategies of proofs of Benaïm (2018, Theorem 4.4), Schreiber et al. (2011) and Arnold and Kliemann (1983, Proposition 3.15), one can show that all weak limit points of Π z t and π z t are candidate μ of Definition 1.
In each example, we present the simulated trajectories of (i) the resident and invader population in the corresponding phase plane, (ii) the total population size N t , and (iii) the relative population size P t of the invader. The numerical analysis shows that N t varies stochastically and dramatically in time, while P t changes slowly and asymptotically in time. Since we don't perform the complete time-scale separation in numerical simulations, one can find that there are some small fluctuations remain in the simulated trajectories of P t in Figs. 4, 6, 8, and 10 for different examples, respectively. With complete time-scale separation (i.e., → 0), the trajectories of P t will be smooth.

Chemostat model
Consider the evolving bacteria with strategy x = (β, γ, δ) ∈ R 3 + := X and concentration n t at time t in a chemostat model, in which the dynamics is given bẏ where R t is the concentration of nutrient at time t. Let R in t be the time-varying concentration of nutrient at the input. Positive constants D, β, γ and δ are dilution rate, nutrient uptake rate, conversion efficiency and death rate of bacteria, respectively. In this model, we assume R in t = ρ 1 − ρ 2 2θ t 1 + θ 2 t with constants ρ 1 > ρ 2 > 0 and {θ t } t≥0 is the stationary Ornstein-Uhlenbeck process generated by a linear stochastic differential equatioṅ θ t = − aθ t + bẆ t , θ 0 = 0, (4.2) with positive constants a and b. Then R in t takes value in the interval [ρ 1 − ρ 2 , ρ 1 + ρ 2 ] for all time t and tends to peaks around ρ 1 ± ρ 2 , which is a random switching scenario. Different ways that an unbounded noise induces bounded fluctuations on model parameters refer to d'Onofrio (2013) and Caraballo and Han (2016). Denote Then (4.1) can be written as the general form (2.1). From the boundedness of R in t , one can show that R t ≤ ρ 1 + ρ 2 and n t ≤ γβ(ρ 1 +ρ 2 ) δ+D for t sufficiently large. Thus, the state space of the dynamics of (n t , e t , θ t ) is The partial derivatives f (i, j) , G (i, j) 1 and H (i, j) 1 exist for all strategies x ∈ X and for all nonnegative integers i, j ≤ 2, which satisfy Assumption A2.
Once bacteria x becomes a resident, the invasion fitness of an initially rare mutant y = (β , γ , δ ) in resident x is For similar strategies (x, y) satisfying ∂ y S x (y)| y=x (y − x) = 0, Theorem 1 can be employed to predict the invasion dynamics. Instead we apply Theorems 2 or 3 if ∂ y S x (y)| y=x (y − x) = 0. However, we will see that substitution is the unique outcome after an invasion event whatever strategy pairs are. In fact, from (2.4), we get that E[R] = δ+D γβ . Hence, invasion happens if and only if δ +D γ β < δ+D γβ . It further implies that for all x, y ∈ X , S y (x) < 0 if S x (y) > 0. Figure 4 shows how a specific population trajectory develops from n res -axis to n invaxis in the phase plane provided that x = (2, 0.8, 0.095) and y = (1.99, 0.85, 0.09) with S x (y) > 0 and S y (x) < 0. The initially rare bacteria y will actually take over the bacteria x and becomes the new resident. The total concentration N t of the resident x and the invader y varies stochastically and dramatically in time. However, the relative concentration P t of the invader y changes slowly in time and asymptotically closes to 1. The dynamics of P t illustrates that the resident x is substituted by the invader y eventually.

Lotka-Volterra competition model
Consider a stochastic Lotka-Volterra competition model where n i,t is the population size of strategy x i ∈ R := X . Assume that the intrinsic growth rate and the carrying capacity of the i-th species are influenced by an external factor θ t such that they fluctuate around strategy-related values. For simplicity, let where d 1 , d 2 , ρ 1 and ρ 2 are scaling parameters and where the {θ t } t≥0 is a stationary Ornstein-Uhlenbeck process generated by (4.2). Then the expected growth rate and the expected carrying capacity of the i-th species are r (x i ) exp 2a , respectively. Let the competitive coefficient between strategy x i and x j be of the form where α(x i , x i ) = 1 for all i. Let E be the set of all functions e : X → R 2 of the form e t (·) = (e 1,t , e 2,t )(·) = 0, j α(·, x j )n j,t for n j,t ≥ 0 and x j ∈ X . Since the competition is modelled as a direct interaction, we can define the per-capita environmental impact H 2 : which is the competitive impact of a single individual with strategy x j on its competitor with any strategy. Since the competitor's strategy has infinitely many choices, e 2,t has infinite dimensions. In addition, define the map G 2 : E × R → R by G 2 (e t , θ) = 0.
Further, define the map f : Thus, (4.4) can be written as the general form (2.1). Furthermore, the state space of the dynamics of (n 1,t , . . . , and H (i, j) 2 exist for all strategies and for all nonnegative integers i, j ≤ 2, which satisfy Assumption A2.
We now start from a monomorphic population of strategy x with size n t at time t. Similar to the previous example, one can show that species x successfully establishes with population dynamics satisfying Definition 1. The only difference that should be pointed out is that the population is unbounded in this example. Instead the population is ultimately bounded in mean (Miyahara 1972(Miyahara , 1973, i.e., lim sup t→+∞ E[n t ] ≤ E lim sup t→+∞ n t = K (x)E lim sup t→+∞ exp(ρ 2 θ t ) = K (x) < +∞ for all x ∈ X . Further, this stochastic boundedness combined with the Feller property guarantees the existence of invariant probability measures for process {n t , θ t } t≥0 . Now, from (4.4), the invasion fitness of an initially rare mutant y in resident x is given by (4.5) Here we have used the principle of selective neutrality of residents (i.e., S x (x) = 0 for all x ∈ X ) to derive (4.5). Figure 5a shows the PIP and the MIP corresponding to (4.5) with parameters (c 0 , c 1 , c 2 , c 3 , c 4 , c 5 ) = (1, −11, 11, −4, 1, 0.5). The features of these two plots are similar to that of Kisdi et al. (2001) which also studies the evolutionary of strategy x but for a deterministic Lotka-Volterra model withr t (x) = K t (x) = 1 for all x ∈ X and for all t ≥ 0. Away from evolutionarily singular strategies x * er (evolutionary repeller) and x * bp (evolutionary branching point), it follows from Theorem 1 that every successful invasion of mutant y in a sufficiently small neighborhood of resident x will takeover the population. Further, directional evolution proceeding by successive invasions and substitutions with small phenotypic effect leads to an evolutionary branching point x * bp provided that the initial strategy value is above x * er , otherwise is towards lowvalued strategies if the initial strategy value is below x * er . Both x * er and x * bp satisfy C 22 + C 11 = 0, so that we can apply Theorem 2 to all strategies in a neighborhood of x * er or x * bp . Having approached x * bp , the population undergoes evolutionary branching that gives rise to two district strategies coexistence. Figure 6a gives a simulated population trajectory in the phase plane for the coexistence of strategies x = 0.48 and y = 0.51 with corresponding invasion fitnesses S x (y) > 0 and S y (x) > 0. The dynamics of P t shows that the relative population size of the invader y asymptotically tends to an interior value of (0, 1).
To find the case of mutual exclusion, we choose two strategies lie in the white area of the MIP in Fig. 5a and close to the x * er . Figure 6b gives two simulated population trajectories in the phase plane for strategies x = − 0.53 and y = − 0.47 with corresponding invasion fitnesses S x (y) < 0 and S y (x) < 0, provided that initial population states are identical. In the phase plane, we see that one simulated population trajectory eventually reaches n res -axis, while the other trajectory eventually reaches n inv -axis. The corresponding dynamics of P t shows a simulated trajectory asymptotically tends to 0 but the other one asymptotically tends to 1. Figure 5b shows the PIP and the MIP corresponding to (4.5) with parameters (c 0 , c 1 , c 2 , c 3 , c 4 , c 5 ) = (1, −2, 0, −2, 1, 0.5). Now the lower singularity is still an evolutionary repeller, while the upper singularity becomes an evolutionarily stable strategy. What's more, the associated C 22 + C 11 = 0 but C 22 = 0 and C 11 = 0 for both two evolutionarily singular strategies x * er and x * ess . Therefore, we apply Theorem 3 to all strategies in a neighborhood of x * er or x * ess , in which substitution is the unique outcome of an invasion event. Figure 6c gives a simulated population trajectory in the phase plane for strategies x = 0.55 and y = 0.51 with corresponding invasion finesses S x (y) > 0 and S y (x) < 0. The dynamics of P t verifies that the invader y substitutes the resident x eventually.

A structured SIRS epidemic model
Consider a structured SIRS epidemic model for the evolution of viral type. For a given moment t and i = 1, . . . , k, let S t be the numbers of susceptible individuals, I i,t be the numbers of individual infected by a virus x i ∈ R + := X , and R i,t be the number of recovered individual whose is infected by the virus x i . Let M t be the total number of susceptible, infected and recovered individuals at time t, i.e., Assume that all individuals are born free of the disease at a constant rate Λ > 0, and all individuals naturally die at a density-dependent rate δ M t with constant δ > 0. In a well-mixed population, the probability of a contact being with an infected individual of x i is given by The probability of that contact giving rise to an infection of x i is given byβ t (x i ), which is commonly called transmission rate. Infectious individuals of x i recover at a rateγ t (x i ) and disease-caused die at a rateα t (x i ). Recovered individuals of x i lose their protection against a reinfection and become susceptible again at a ratẽ ζ t (x i ). Hereβ t (x i ),γ t (x i ),α t (x i ) andζ t (x i ) are all positive, time-varying and virusdependent. For simplicity, we assume that the transmission rate fluctuates around a virus-dependent value in the way:β t (x i ) = β(x i ) exp(ρ 1 θ t ) with β(x i ) > 0 for all x i ∈ X , where the {θ t } t≥0 is a stationary Ornstein-Uhlenbeck process generated by (4.2) and ρ 1 is a scaling parameter to measure the effect of θ t on the transmission rate. Likewise,γ t (x i ),α t (x i ) andζ t (x i ) fluctuate in time in the same way asβ t (x i ) but the scaling parameters are different, i.e.,γ t ( Assume further that one disease-free individual is infected with a certain virus, it cannot be infected by any other virus. Under these assumptions, a structured SIRS epidemic model is given by the following differential equations: where the state space of v i,t is Δ = v ∈ [0, 1] 2 : v 1 = 1 . Then the dynamics of (4.6) are equivalent to which can be written as the form specified in (3.2) with e t = (e 1,t , e 2,t ) = (S t , M t ), It is easy to show that M t ≤ Λ δ for t sufficiently large. Thus, the state space of the dynamics of Z t = n 1,t 1 , . . . , n k,t 1 , v 1,t , . . . , v k,t , e t , θ t is Z = 0, Λ δ k × Δ k × 0, Λ δ 2 × R, the extinction set Z 0 = z ∈ Z : min i n i 1 = 0 , and the η-neighborhood of the extinction set Z η = z ∈ Z : min i n i 1 ≤ η . The partial and H (i, j) 2 exist and are locally Lipschitz continuous and measurable in (e, θ) for all strategies and for all nonnegative integers i, j ≤ 2.
We now start from a monomorphic population of viral type x with the total population size n t 1 and the population structure v t . To proceeding the further analysis, we need to know when a given viral type x becomes a resident, i.e., it successfully establishes in the virgin environment. The dynamics on the boundary of n t 1 , S t , M t -space globally converges to a unique stable equilibrium 0, Λ δ , Λ δ , provided that S 0 = 0. Hence, the ergodic measure on the boundary of n t 1 , S t , M tspace is the Dirac measure at 0, Λ δ , Λ δ . Substitution of this stable equilibrium into the equation Providing that the expectation the dynamics of (4.7) consist of a Dirac measure at the unstable equilibrium v = (0, 1) and invariant probability measuresv supported on Δ \ {(0, 1)}. Let μ be the product of the Dirac measure at n 1 , S, M = 0, Λ δ , Λ δ and thev. Then corresponds to the per-capita growth rate of n t 1 with respect to μ. From the dynamics of v 1 t we have that with probability one lim t→+∞ Thus, if (4.8) holds, a similar argument as the first example implies that the population of viral type x is non-growing on the long-run and stochastically persistent. Simulated population trajectories in the phase plane and associated dynamics of the total population size N t = (I res,t + R res,t ) + (I inv,t + R inv,t ) and the relative population size P t = (I inv,t + R inv,t )/N t for strategy pairs (a) (x, y) = (8.7, 8.68) with corresponding invasion fitnesses S x (y) ≈ 0.00708633 > 0 and S y (x) ≈ −0.00709391 < 0, and (b) (x, y) = (4.1158, 4.1238) with corresponding invasion fitnesses S x (y) ≈ 2.27275 × 10 −6 > 0 and S y (x) ≈ 2.27315 × 10 −6 > 0, where parameter values refer to Fig. 7 Once virus x becomes a resident, the invasion fitness of an initially rare mutant y in resident x is given by (4.9) Here we have used the principle of selective neutrality of residents (i.e., S x (x) = 0 for all x ∈ X ) to derive (4.9). Generally, there isn't an explicit expression of (4.9) in terms of strategies x and y. Thus the numerical PIP and MIP are based on a small-noise approximation of (4.9) (see e.g., Vilar and Rubi 2018). Figure 7 shows the PIP and the MIP corresponding to a small-noise approximation of (4.9) (ref to the legend). Away from the two evolutionarily singular strategies x * er and x * bp , every successful invasion of mutant y in a sufficiently small neighborhood of resident x will takeover the population. Figure 8a gives a simulated population trajectory in the phase plane for strategies x = 8.7 and y = 8.68 with corresponding invasion fitnesses S x (y) > 0 and S y (x) < 0. The dynamics of P t shows that the invader y will eventually oust the resident x and become the new resident. Close to the two evolutionarily singular strategies, Figure 8b gives a simulated population trajectory in the phase plane for strategies x = 4.1158 and y = 4.1238 with corresponding invasion fitnesses S x (y) > 0 and S y (x) > 0. The dynamics of P t shows that the relative population size asymptotically tends to an interior value of (0, 1), which implies that the two viral types can coexist.

A prey-predator model: continuation of Example 1
In order to show that the deterministic model of Example 1 belongs to the general class of models considered in this paper, let then the deterministic model of Example 1 can be written as the general form (2.1).
The partial derivatives f (i, j) , G and H (i, j) 2 exist for all strategies and for all nonnegative integers i, j ≤ 2, which satisfy Assumption A2.
The evolutionary dynamics of the deterministic model of Example 1 has been well studied by Lehtinen and Geritz (2019). We are interesting in the evolution of timidity of the prey with the non-equilibrium resident dynamics. When only a single prey type x is present, the system has a stable interior equilibrium at x = 1. Decreasing x destabilises the equilibrium through a supercritical Hopf bifurcation at x Hopf = 0.6289, after which periodic attractors are present (ref to Lehtinen and Geritz (2019, Fig. 3B)). In a periodic resident environment set by a single resident type x, the invasion fitness of an initially rare mutant y is given by where τ (x) is the period (ref to Metz et al. 1992). Figure 9 shows the PIP and the MIP corresponding to (4.10) with the same parameters as Lehtinen and Geritz (2019, Fig . 4), in which the monomorphic resident environment is an equilibrium environment if x ≥ x Hopf and becomes a periodic environment if x < x Hopf . In the periodically resident environment, away from the evolutionarily singular strategy x * ess , it follows from Theorem 1 that successful invasion of mutant y in resident x will takeover the population. Figure 10a gives a simulated population trajectory in the phase plane for strategies x = 0.27 and y = 0.3 with corresponding invasion fitnesses S x (y) > 0 and S y (x) < 0. The total population size of N t = n res,t + n inv,t varies significantly in time. However, the relative population size P t = n inv,t N t changes slowly in time and asymptotically increases from 0 to 1. The dynamics of P t implies that the invader y will oust the resident x and becomes a new resident. Figure 10c shows the periodic attractor of the monomorphic model with the new resident strategy x = 0.3. In the neighborhood of x * ess , we can apply Theorem 2 to predict the population dynamical outcomes of an invasion event. If both strategies of the resident and the mutant are in the gray area of MIP, it follows from Theorem 2 that the resident and the mutant can coexist eventually (notice that all coexistence strategies are not evolutionarily stable). Figure 10b gives a simulated population trajectory in the phase plane for the coexistence strategies (x, y) = (0.3764, 0.4314) with corresponding invasion fitnesses S x (y) > 0 and S y (x) > 0. The total population size N t still fluctuates significantly in time. but the relative population size P t slowly and asymptotically tends to an interior value of (0, 1). The dynamics of P t verifies that the invader y eventually coexists with the resident x. Figure 10d shows the periodic attractor of the dimorphic model with coexistence strategies (x 1 , x 2 ) = (0.3764, 0.4314).

Discussion
Our main result is the complete classification of generic population dynamical outcomes of resident-invader dynamics in fluctuating environments when the invader and the resident have similar strategies. The outcomes are essentially "Lotka-Volterra": (i) one strategy ousts the other if the one can invade a population of the other but not the other way around, i.e., invasion implies substitution; (ii) the two strategies coexist if they can mutual invade; (iii) they mutually exclude one another if neither can invade a population of the other. Which of the four Lotka-Volterra-type outcomes occurs depends only on the signs of the invasion fitnesses. In a one-dimensional strategy space or a one-dimensional parameterization of a multi-dimensional strategy space (see e.g., Kisdi 2015), invasion implies substitution away from an evolutionarily singular strategy, while all four Lotka-Volterra-type outcomes are possible close to an evolutionarily singular strategy. In a multi-dimensional strategy space, however, the situation is generally complicated because all four Lotka-Volterra-type outcomes may occur also away from an evolutionarily singular strategy. We extend and generalize previous results of resident-invader dynamics of similar strategies to models incorporating (i) explicit feedback environments with their own dynamics, (ii) scalar-valued as well as vector-valued strategies, (iii) unstructured as well as structured populations, (iv) monomorphic as well as polymorphic resident populations, and (v) non-equilibrium resident population dynamics as well as resident dynamics with stochastic (or deterministic) drivers. Although we show all results for models with a monomorphic resident population (i.e., single resident phenotype), the generalization of them to polymorphic resident populations (i.e., multiple resident phenotypes) is straight forward because of the way we modelled the environment feedback loop. Arbitrarily polymorphic resident populations can be accounted for by treating the corresponding population sizes of the extra resident phenotypes as environmental feedback variables. In the next paragraphs we focus on the differences between our work and previous studies.
For a class of unstructured population models with point equilibria, Geritz (2005) has shown that the resident-invader dynamics generically behaves like the Lotka-Volterra competition model by using of the Poincaré-Bendixson theorem. The results hold for vector-valued strategies but were proved only for a feedback environment that is given as an explicit linear function of the resident and invader population sizes. However, environmental feedback variables such as resources, predators, or competitors are often given implicitly (e.g., by differential equations) as opposed to explicitly. Our implicit representation is more general than the explicit and liner environmental feedback used in Geritz (2005). Dercole and Geritz (2016) has the same aims and results as Geritz (2005), but it is more general in some aspects but less general in others, and it uses a basically different mathematical approach, i.e., time-scale arguments. The results only apply to unstructured population models with point equilibria and scalar-valued strategies, but the formulation of environmental feedback is very general-probably as general as in our approach (in their formulation, they allow the resident and the invader to interact with finitely many other populations whose corresponding sizes are packed in the feedback variable e t . The dynamics of e t is govern by a differential equation, in which the growth rate is an implicit function of the resident, the invader, and the environmental variable itself. Given four structural properties of the growth rate of the e t as well as the per-capita growth rates of the resident and the invader, their formulation generalizes the law of mass action (see Dercole 2016). The generalization takes into account that pairwise interactions can depend on the concomitant activities of the encountered individuals, which leads a nonlinear density dependence in the per-capita growth rates of the resident and the invader). The paper also studies some mathematically degenerate cases which confirm the results of Priklopil (2012) on unprotected coexistence of two strategies near an evolutionarily singular strategy.
The "invasion implies substitution" outcome for unstructured population models with point equilibria also can be found in Dercole and Rinaldi (2008, Section 3.4) and Oba and Kigami (2018).
Our paper is the first to study resident-invader dynamics in fluctuating environments. We study the resident-invader dynamics as a stochastic fast-slow system where the total population size of the resident and the invader is the fast variable and the relative population size of the invader is the slow variable. We show that trajectories of the slow variable on slow timescale are well approximated by that of an associated averaged system, and the stability of the averaged system depends on invasion criteria alone. From these results, Theorems 1, 2 and 3 give the complete classification of generic population dynamical outcomes of an invasion event.
For a class of structured population models with point equilibria, recently, Priklopil and Lehmann (2019) have shown that invasion implies substitution for ecological communities with finite-class-structured populations of scalar-valued strategies. Their approach is based on the analysis of a weighted average of the relative invader population sizes in each class, where the weighting coefficient is the class-specific reproductive value (see also Lion 2018a). However, it is not clear how their method can be used to give a full classification of the generic population dynamical outcomes of an invasion event. Moreover, it is not clear how their approach can be applied to fluctuating environments. Cantrell et al. (2017) has extended the tube theorem of Geritz et al. (2002) and the invasion implies fixation theorem of Geritz (2005) to a class of reaction-diffusion models for understanding evolution of dispersal in space. Their focus is on a dimorphic system (i.e., one resident and one invader) of infinite-dimensional structured populations of scalar-valued strategies.
When populations are in a fluctuating environment, we illustrate how to generalize the results of unstructured population models to a class of structured population models, but only for finite-class-structured populations. The extension of our results to populations with a more general structure (e.g., size distributions, continuous age and spatial diffusion) we leave for future work.
As mentioned in the Introduction, there are serval different (not necessarily equivalent) definitions of stochastic persistence for a resident population. Most of them can be summed up from the "ensemble" perspective in terms of transition probabilities (see e.g., Chesson 1982;Chesson and Ellner 1989;Li and Mao 2009) or the "typical trajectory" perspective in terms of empirical measures (i.e., how typical sample trajectories of the population process are distributed in time, see e.g., Benaïm andSchreiber 2009, 2019;Schreiber et al. 2011;Roth and Schreiber 2014;Benaïm 2018) for the long-term population dynamics. Different definitions often give different interpretations of the population dynamics. Extending our results to models with different kinds of resident dynamics would be useful for the applicability of the theory of adaptive dynamics. Our definition of resident persistence follows Schreiber et al. (2011) and Benaïm (2018) is a rather weak one. This means that our results also apply to more restrictive definitions.
Lemma 3 Let Z 0 = Z 0 0 = z. Then, for all bounded and measurable function h : Z → R, we have for all t ≥ 0, and is uniformly in z ∈ Z.
Proof For all t ≥ 0 Let Using the Cauchy-Schwartz inequality we have By Assumption A3, there exist some positive constants K 0 , K 1 and K 2 such that Then, by the Gronwall's lemma, it follows that for all t ∈ [0, T ] where K 0 = 3T K 0 and K 1 = 3T 2 (K 1 ξ 1 + K 2 ξ 2 ). Let h : Z → R is a bounded and measurable function. By uniform continuity, there exists δ > 0 such that for all z 1 , z 2 ∈ Z with z 1 − z 2 ≤ δ, we have |h(z 1 )−h(z 2 )| ≤ . Using the Markov's inequality, we obtain that for t ∈ [0, T ] This inequality holds for any T > 0, we thus get the desired result for all t ≥ 0.
The proof is complete.

Appendix C Proof of Lemma 1
From now on, we fix initial values P 0 =P 0 = p and Z 0 = z. The proof of Lemma 1 will be separated into the following three steps. Firstly, we show that the trajectory ofP t 1 is sufficiently close to that of P t 1 for > 0 small in every finite time interval.

Proof
Denote where z = (n, e, θ). Leth be the expectation of h(z) with respect to μ, i.e., For all t 1 ≥ 0, ignoring the term O( ) in (2.7), Fixed T > 0, using the Cauchy-Schwartz inequality we have By the Gronwall's lemma, it follows that We now focus on the limit of α and β as → 0, respectively. By Birkhoff's ergodic theorem, we obtain that for all p ∈ [0, 1] and for μ-almost every z ∈ Z \ Z 0 Here, we have used that Z\Z 0 | f (1,0) (x, e, θ) (ξ 2 − ξ 1 )| 2 μ(de, dθ) < +∞ due to the hypothesis of the lemma and where the positive constant K σ satisfies ξ 2 − ξ 1 = √ K σ (y − x) for y ∈ x ∈ X : x − x ≤ σ . In fact, for given ξ 1 , ξ 2 and x, we can find a σ to be such that the hypothesis of the lemma is satisfied, and then determine y and K σ .
To estimate the term β , we divide [0, s] into some intervals depending of a given size and define a functionP s byP where s is the integer part of s . Then, using the Cauchy-Schwartz inequality we have Following the proof of Liu and Krstic (2012, (4.110) and ( Secondly, following the proof of Liu and Krstic (2012, the assertion (i) of Theorem 4.1, see pp. 65-66), the finite-time approximation (c.1) can be extended to arbitrarily long time intervals as the following lemma shows. For the sake of self-containedness, we give the proof.
The proof is complete.

Appendix E Proof of Lemma 2
The first assertion can be proven in the same manner as in Lemma 1, in which the associated averaged system becomes (2.12). The details are left to the reader.
Next, we focus on the proof of the second assertion. Let P 0 =P 0 = p and Z 0 = z. Denote , p * 1 = 1.