Comparison between Stochastic and Deterministic Selection-Mutation Models

We present a deterministic selection-mutation model with a discrete trait variable. We show that for an irreducible selection-mutation matrix in the birth term the deterministic model has a unique interior equilibrium which is globally stable. Thus all subpopulations coexist. In the pure selection case, the outcome is known to be that of competitive exclusion, where the subpopulation with the largest growth-to-mortality ratio will survive and the remaining subpopulations will go extinct. We show that if the selection mutation matrix is reducible, then competitive exclusion or coexistence are possible outcomes. We then develop a stochastic population model based on the deterministic one. We show numerically that the mean behavior of the stochastic model in general agrees with the deterministic one. However, un like the deterministic one, if the differences in the growth-to-mortality ratios are small in the pure selection case, it cannot be determined a priori which subpopulation will have the highest probability of surviving and winning the competition.


Introduction
Deterministic selection-mutation models, also referred to as distributed rate population models with closed-open reproduction, have been studied by many researchers [1,2,3,4,8,10,11,15].These are models for the density of individuals with respect to some evolutionary discrete or continuous trait.In [4] a pure selection model (closed reproduction) with logistic type nonlinearity and a continuous 2-dimensional trait variable (growth and mortality) was studied.Therein, the authors prove that competitive exclusion occurs and the surviving subpopulation is the one with the largest growth to mortality ratio.In [3] these results were extended to a more general population model constructed on the (natural) space of measures.In this case, the limiting measure (a Dirac delta function) is an element of the state space.
In [1,2] a nonlinear size-structured population model with discrete trait (a finite number of subpopulations each having its own growth, mortality, and reproduction functions) was studied.It was shown that in the case of closed reproduction competitive exclusion occurs and the winning subpopulation is the one with the highest ratio of reproduction to mortality.In the case of open reproduction the authors show that survival of all subpopulations is possible.
In [9] a simple population model with two groups of age (juveniles and adults) for the mean age at maturity was studied.The authors show the existence and uniqueness of a globally attractive stationary solution.This simple model is then used in [10] to build a system of equations for the density of individuals with respect to age at maturity.In this model there are two birth terms one for selection (closed reproduction) with a probability (in the form of a coefficient) 1 − and another for pure mutation (open reproduction) with a probability .Therein, the authors prove the existence of L 1 stationary solutions which tend to concentrate around the evolutionary stable value of the trait when the open reproduction term coefficient → 0. In [11] a pure mutation model is considered.The authors prove the existence of stationary solutions and they investigate the behavior of these stationary solutions when the mutation is small.
It is well known that the deterministic and stochastic models may behave differently.
Hence, the goal of this paper is to present a deterministic selection-mutation model, and then formulate a stochastic differential equation model based on it and compare the dynamics of these two models.The derivation of the stochastic differential equation model is based on the method developed in [6,7,13], where the variability inherent in the system is only due to the demographic variability.This paper is organized as follows: in Section 2 we present a selection-mutation model with a 2-dimensional discrete trait variable.We recall that in the case of a pure selection birth term competitive exclusion between the different traits occurs and the winner trait is the one with highest growth to mortality ratio.In the case of a selection-mutation birth term we prove the existence of a unique stationary solution and we show that this stationary solution is globally attractive.In Section 3 we develop an Itô stochastic differential equation model which is based on the deterministic model.In Section 4 we consider two numerical examples one for a pure selection case and another for a selection-mutation case.We show that in both cases the stochastic model behaves like the deterministic model.However, the main difference is that in the pure selection case if the ratio of growth to mortality of the first subpopulation is not much larger than the other subpopulations then unlike the deterministic model the subpopulation with the larger ratio may go to extinction.Section 5 is devoted to concluding remarks.

Deterministic Model
The selection-mutation deterministic model that we consider is given as follows: (2.1) Here the total population X = n j=1 x j , and the i-th subpopulation is identified by the 2- , where a i is a scaled per capita growth rate and b i a scaled per capita mortality rate.The parameter p ij denotes the probability of an individual in the jth subpopulation reproducing an individual in the ith population, 0 ≤ p ij ≤ 1 and For convenience, we will denote by P = (p ij ) n,n the matrix whose (i, j) entry is given by p ij .Clearly P is a column stochastic matrix since each of the column sums of P is equal to one.
Remark: For the special case where P is identity matrix, i.e., p ii = 1 and p ij = 0 for i = j (2.1) reduces to the following pure selection model: (2.2) Here the quotient a i b i can be thought of as a scaled reproductive ratio, i.e. it is a measure of the average amount of offspring an individual of trait (a i , b i ) produces during its lifetime.
The actual reproductive ratio at population density X is given by T ∈ R n , then the system of (2.1) can be written in the following vector form: where x = [x 1 , x 2 , . . ., x n ] T and x 0 = [x 0 1 , x 0 2 , . . ., x 0 n ] T .We now establish equilibrium and stability results for the model (2.1).In the next theorem, we show that under some conditions there exists a unique positive equilibrium for system (2.3).
Theorem 2.1 Suppose that P is an irreducible matrix, then there exists a unique positive equilibrium for system (2.3).Note that since BR = A and R −1 A = B, we see that where r = Re.In order to show that there exists a unique positive equilibrium for (2.3), we only need to show that it is true for (2.4).Note that since a i /b i > 0 for i = 1, 2, . . ., n and P is a column stochastic irreducible matrix, P R is a nonnegative irreducible matrix.By Perron-Frobenius theorem, we know that there exists a unique positive normalized eigenvector ȳ corresponding to the spectral radius ρ y such that Clearly e T ȳ = 1, e T P = e T and e T R = r T .Therefore, we find that By substituting ρ y = r T ȳ into (2.5),we have Hence, ȳ is an equilibrium of system (2.4).Suppose that there exists another positive equilibrium ū for (2.4), then ū is a positive eigenvector of P R corresponding to the positive eigenvalue r T ū.By Perron-Frobenius theorem, we have ρ y = r T ū.Since P is column stochastic, we get This implies that ū is a positive normalized eigenvector corresponding to ρ y .Hence, ū = ȳ.
Thus, there exists a unique positive equilibrium ȳ for the system (2.4), which implies that Rȳ is the unique positive equilibrium for the system (2.3).
In the next theorem we prove the global asymptotic stability of the unique interior equilibrium under some conditions on the model parameters.
Theorem 2.2 Assume that P is an irreducible matrix, and b i = β for i = 1, 2, . . ., n, where β is a positive constant.Then the unique positive equilibrium is globally asymptotically stable.
Proof.Since P is an irreducible matrix, AP is an irreducible matrix.By Perron-Frobenius theorem, there exists a unique positive normalized eigenvector x correponding to the spectral radius ρ * such that AP x = ρ * x, and this eigenvalue which is equal to ρ * is a simple eigenvalue of AP .Let M = AP − ρ * I with I being the identity matrix.It is known that there exists a positive constant c such that lim t→∞ exp(M t)x 0 = cx (e.g., [5]).Substituting 3) with g being a scalar function of t, we find that Hence, the limiting equation of (2.8) is given by the following logistic model which has a unique positive equilibrium ḡ∞ = ρ * βc and it is globally asymptotically stable.
Furthermore, it is easy to see from (2.8) that g is bounded.Therefore, it follows from the theory on asymptotically autonomous systems that lim t→∞ g(t) = ḡ∞ (e.g., [16]), which implies the theorem.
Remark: By a similar process, we can show that Theorem 2.2 still holds for the more general system ẋ = AP x − βf (e T x)x, where f is a Lipschitz continuous strictly increasing function with f (0) = 0 and lim
The probability of a single birth and death in subpopulation i is given by X j ∆t, respectively, i = 1, 2, . . ., n.Then the probability that no change occurs in any subpopulation in the time interval ∆t is X j )∆t.Therefore, the infinitesimal mean is approximated as follows: Notice that E(∆X)E(∆X) T is of order ∆t 2 , and X j ∆t, i = 1, 2, . . ., n, and the off-diagonal elements of V are 0.

Now we have
where i = 1, 2, . . ., n, and ξ 1 , ξ 2 , . . ., ξ n are independent random variables following a normal distribution with mean 0 and variance ∆t.Clearly the above equation is just an Euler approximation of X i (t) with time step ∆t.Notice that from [12,14], it follows that as ∆t → 0, X(t) converges in the mean square sense to the solution of the following stochastic system where i = 1, 2, . . ., n, and

Numerical Results
In this section, two cases will be studied.In Section 4.1 we study a selection birth term where for the deterministic model the subpopulation with the highest ratio max i=1,2,...,n will survive and all the other subpopulations will die out.In Section 4.2 we study a selectionmutation case with P being an irreducible matrix (where all subpopulations survive in the deterministic model).
We use Euler's method to numerically approximate the sample paths for the stochastic model (3.3).In all the simulations for the stochastic model, we choose the time mesh size ∆t = 0.001.The mesh points are given by: t k = k t, k = 0, 1, . . ., m. Denote by X k i the numerical solution of X i (t k ).Then we have the following numerical scheme: where ξ k 1 , ξ k 2 , . . ., ξ k n are independent random variables following a normal distribution with mean 0 and variance ∆t.In the simulation, if X k+1 i ≤ 0, then we set In all of the examples below, we simulated N sample paths X i,l (t) with l = 1, 2, . . ., N and i = 1, 2, . . ., n for the corresponding stochastic model.The mean conditioned on nonextinction for the i-th subpopulation µ i (t) is calculated by averaging all the sample paths conditioned on nonextinction of the total population The sample standard deviation for the i-th subpopulation is represented by σ i (t).We use π i (0, t) = Prob{X i (t) = 0}, i = 1, 2, . . ., n to denote the probability of extinction for the ith subpopulation at time t, and π(0, 0, . . ., 0, t) = Prob{X 1 (t) = 0, X 2 (t) = 0, . . ., X n (t) = 0} to denote the probability that all the subpopulations are driven to extinction at time t.We denote by π c (t) = Prob{X 1 (t) > 0, X 2 (t) > 0, . . ., X n (t) > 0} the probability of coexistence at time t.

Pure Selection Case
In this section, we consider a pure selection case, that is, P = I, which means that individuals in the ith subpopulation can only reproduce individuals in the ith subpopulation.In what follows, we choose n = 3, T = 10, and we use ode23 in Matlab to numerically solve the deterministic model (2.1).
The probability of coexistence π c (t) is plotted in Figure 3 (right), which indicates that the subpopulations no longer coexist after t = 3.5.Hence, for this example the stochastic model behavior is similar to that of the deterministic model.We point out that when the highest ratio a 1 /b 1 is significantly larger than the ratios a 2 /b 2 and a 3 /b 3 as in the above example, the stochastic (in the mean sense) and deterministic model behave similarly.In order to test if the same outcome occurs when the highest ratio is close to the other ratios, we give the following example: the initial population size is chosen to be (20,20,20).The values for all the other parameters are chosen to be a 1 = 3λ, a 2 = 6λ,   The mean conditioned on nonextinction (µ 1 (t), µ 2 (t), µ 3 (t)) and the solution to the deterministic model are plotted in Figure 5.The results suggest that the mean conditioned on nonextinction is far away from the deterministic equilibrium (270,0,0).The sample standard deviation (σ 1 (10), σ 2 (10), σ 3 (10)) is approximately equal to (124.3057, 112.9720, 134.2170).The probabilities of extinction π i (0, t) and π(0, 0, 0, t) are plotted in Figure 6 (left), from which we can see that these probabilities are approximately constant for a large range of times.Their numerical approximations at time t = 10 are given as follows: π 1 (0, 10) ≈ 0.6937, π 2 (0, 10) ≈ 0.7691, π 3 (0, 10) ≈ 0.5371, π(0, 0, 0, 10) = 0.
The probability of coexistence π c (t) is plotted in Figure 6 (right), which indicates that the populations no longer coexist after t = 1.We performed many other numerical simulations with different λ values.In these simulations, we obtain many cases where the stochastic model differs from the deterministic one in picking the winning trait or subpopulation.The results in Table 1 and Figure 7 are obtained by simulating 1000 (N = 1000) sample paths.The probability of extinction π 1 (0, 10), π 2 (0, 10) and π 3 (0, 10) for some λ values less than 10 and bigger than 200 with initial conditions (20,20,20) are tabulated in Table 1.We can see that if λ is chosen very small or very big then all the subpopulations are driven to extinction, which is expected for the stochastic model because of the variance.Table 1: Probability of extinction π 1 (0, 10), π 2 (0, 10) and π 3 (0, 10) for λ values less than 10 and bigger than 200 with initial conditions (20,20,20).
The probability of extinction π 1 (0, 10), π 2 (0, 10) and π 3 (0, 10) with λ values between 10 and 200 are plotted in Figure 7, where the left part is the plot for results with initial conditions (20,20,20) and the right one is the plot for results with initial conditions (20,40,80).We can see from Figure 7 that the results remain similar even if we choose the initial conditions to be (20,40,80), scaled according to the choice of parameters.Both plots indicate that the stochastic model disagrees with the deterministic model much more in the middle than in the tails.
In order to get a better understanding of the stochastic model dynamics, we consider three values of λ, λ = 20, 100, 180, with initial conditions (20,20,20), and plot the frequency histograms for X 1 (10), X 2 (10) and X 3 (10) out of these 1000 sample paths in Figures 8, 9

Concluding Remarks
In this paper we compared a stochastic and deterministic selection-mutation models.For the most part the behavior of the mean conditioned on nonextinction of stochastic model follows that of the deterministic model.In particular, in the case of open reproduction (selectionmutation) all subpopulation survive in both the stochastic and deterministic models.In the case of closed reproduction (pure selection) one population survives in both models.However, a main difference is that if subpopulations have growth to mortality ratios that are close to each other then in the deterministic case the one with the highest ratio will win the competition but it is not so in the stochastic model.Because of stochasticity the population with the lower ratio may out compete the one with the higher ratio.Therefore, it is not possible to determine apriori which subpopulation will have the highest probability of winning the competition.

3 Figure 2 :
Figure 2: The solution to the model (2.1) and the mean conditioned on nonextinction.

Figure 3 :
Figure 3: (left) The probability of extinction for each subpopulation and the probability that all subpopulations are driven to extinction.(right) The probability of coexistence.

Figure 4 :
Figures 4, 5 and 6 are obtained with λ = 90.We simulate 7000 (N = 7000) sample paths for the corresponding stochastic model.Three randomly chosen stochastic realizations and

Figure 5 :
Figure 5: The solution to the model (2.1) and the mean conditioned on nonextinction.

Figure 6 :
Figure 6: (left) The probability of extinction for each population and the probability that all subpopulations are driven to be extinct.(right) The probability of coexistence.

3 Figure 13 :
Figure 13: The solution to the model (2.1) and the mean conditioned on nonextinction.