Analysis of the Stochastic Population Model with Random Parameters

The population models allow for a better understanding of the dynamical interactions with the environment and hence can provide a way for understanding the population changes. They are helpful in studying the biological invasions, environmental conservation and many other applications. These models become more complicated when accounting for the stochastic and/or random variations due to different sources. In the current work, a spectral technique is suggested to analyze the stochastic population model with random parameters. The model contains mixed sources of uncertainties, noise and uncertain parameters. The suggested algorithm uses the spectral decompositions for both types of randomness. The spectral techniques have the advantages of high rates of convergence. A deterministic system is derived using the statistical properties of the random bases. The classical analytical and/or numerical techniques can be used to analyze the deterministic system and obtain the solution statistics. The technique presented in the current work is applicable to many complex systems with both stochastic and random parameters. It has the advantage of separating the contributions due to different sources of uncertainty. Hence, the sensitivity index of any uncertain parameter can be evaluated. This is a clear advantage compared with other techniques used in the literature.


Introduction
The Verhulst logistic or population equation is a model for many practical applications. It has many applications that imply nonlinear dynamics such as sociology, biology and economics. For example, in epidemiology where the nonlinear term is scaled by a small parameter to denote the contagious rate and in physics to model diffusion with nonlinear perturbations due to unexpected material changes caused by the environment [1].
Deterministic population models have been extensively studied, including models with time variation of the carrying capacity [2,3]. The deterministic models are practical only for sufficiently large populations. Additionally, they neglect many factors that can be significant to the model such as the stochastic and random variations due to different sources.
Population modeling with stochastic variations accounts for the environmental and/or external perturbations in many ecosystems. The environmental fluctuations in the population model can be

Stochastic Spectral Techniques
Assume that the model contains a set of random parameters U(ω) and depend on a set of n real-valued random variables ζ = {ζ 1 , · · · , ζ n } with predefined probability distributions i.e., U(ω) = U(ζ(ω)). Consider the orthonormal basis functionals ψ k (ζ); k ∈ N for the space L 2 of second-order functionals in the random variables ζ(ω). Many choices for the basis functionals ψ k including multiwavelets or the multivariate polynomials [12]. Any second-order random process X(t, ω) that depends on U(ω) can then be decomposed as follows [16]: For practical computations, the Expansion (1) is truncated to its first P + 1 terms. The gPC expansion converges in the mean square sense and [17], i.e., In the case of the model with Brownian motion randomness, the stochastic process X can be decomposed using (M + 1) truncated WHE as follows [9]: The functionals x ( j) (t 1 , t 2 , · · · , t j ) are usually assumed to be symmetric in the variables t 1 , t 2 , . . . , t j . But dW(t 1 )dW(t 2 ) · · · dW(t j ) = H ( j) (t 1 , t 2 , · · · , t j ) dt 1 dt 2 . . . dt j , where H ( j) is the Hermite functional of order j. For brevity we shall write x ( j) without its independent variables. Now we can write where dτ j = dt 1 dt 2 . . . dt j and the integral R j is a j-dimensional iterated integral over the variables t 1 , t 2 , . . . , t j . The WHE details are explained in many references cited in the current work, especially the contribution in [9] and the references therein. Additionally, below a simple linear stochastic model is analyzed in more details. Consider the stochastic linear transport population model: using WHE Expansion (2) to get d dt multiply by H (0) = 1 and H (1) , respectively, and then apply the expectation to get and so on in cases where higher-order kernels and moments are required. In this linear case, analytical solutions can be obtained, for example, and hence

The Deterministic Population Model
Consider the transport population model with nonlinear losses: where the growth rate a = λ − µ, which is the difference between the birth rate λ and death rate µ.
The model parameters, a and ε, are assumed deterministic in this case. The model Equation (3) is of Bernoulli type of order 2 and has the exact solution: where ε/|a| = 1/L is the "crowding term", the reciprocal of the carrying capacity L (saturation level). For a > 0, the system will reach an infinite limit equals the carrying capacity L and the well-known sigmoid/logistic curve is obtained. For a < 0, the system decays to zero, i.e., the population will become extinct, as shown in Figure 1. In the current work we shall assume the case of an infinite limit (population persistence) but the techniques presented are also applicable for the population extinction scenario.
and so on in cases where higher-order kernels and moments are required. In this linear case, analytical solutions can be obtained, for example, x t x e = and ( )

The Deterministic Population Model
Consider the transport population model with nonlinear losses: where the growth rate a λ μ = − , which is the difference between the birth rate λ and death rate μ . The model parameters, a and ε , are assumed deterministic in this case. The model Equation (3) is of Bernoulli type of order 2 and has the exact solution: is the "crowding term", the reciprocal of the carrying capacity L (saturation level). For 0 a > , the system will reach an infinite limit equals the carrying capacity L and the well-known sigmoid/logistic curve is obtained. For 0 a < , the system decays to zero, i.e., the population will become extinct, as shown in Figure 1. In the current work we shall assume the case of an infinite limit (population persistence) but the techniques presented are also applicable for the population extinction scenario.   A more general form of (3) is called Rechards model [19] and it takes the form A more general form of (3) is called Rechards model [19] and it takes the form Here, the exponent β > 0 is called an allometric parameter and it is a measure of curvature. Hence, it can be useful for modeling purposes but will not affect the steady-state behavior. The exact solution is then In the current work, we shall consider the original Verhulst model which is Richards model with β = 2.

The Stochastic Model
Consider the stochastic transport population model with nonlinear losses: with λ is a diffusion coefficient. In the case of deterministic parameters, a and ε, the noise is the only source of randomness in the model (4). It is more practical to consider the growth rate as a random parameter due to the environmental fluctuations or any other excitations. The nonlinear coefficient ε can be also random for different reasons. In this case, the process X will be under multiple sources of randomness.
In the case of multiple sources of randomness, the process X should locally and globally satisfy the finite variance criteria [18], i.e., where U is the vector of all random parameters in Model (4).
The SDE model (4), with x 0 > 0, has a unique positive solution X: which reduces to the deterministic solution in the case where the diffusion parameter λ = 0. The stochastic process X(t) given in (4) and for a > λ 2 /2 is recurrent and converges to the stationary probability gamma distribution γ 2a/λ 2 − 1 , λ 2 /2ε . For a < λ 2 /2, the process decays/converges almost sure to zero [7]. In the case in which a = λ 2 /2, a noise-induced transition solution is obtained. Solution (5) is not helpful in obtaining the exact mean and variance in analytical form. Only approximate analytical formulae can be obtained, see [19] for example. The problem becomes more complicated in the case where one or more of the model parameters is random. In the current work, a numerical technique is introduced that can be used to analyze the stochastic population model even in the case of random variations of the model parameters.
Approximate mean and variance analytical formulae are obtained in [19], but the variance formula is not accurate enough and does not reflect the system dynamics.

Analysis Using WHE
Use WHE Expansion (2) in SDE (4) to get As it is known in WHE [9], the first (zero-order) term is the mean and the second (first-order) term accounts for the Gaussian part of the stochastic solution. Higher-order terms are the nonGaussian part of the stochastic solution with dominance of the second-order term. So, in most applications, even nonlinear, the second-order WHE is sufficient to capture most of the important system dynamics.
Consider only up to second-order kernel x (2) , apply WHE algorithm to get where δ(.) is the Dirac delta function. The initial conditions, after applying WHE, will be: The system of Equations (6)- (8) can be solved simultaneously using the fixed-point iterative scheme. The linear terms are moved to the left-hand side and are computed at the new time level while the nonlinear terms, and the forcing terms will be in the right-hand side and assumed at the old time level. This requires conditions on the values of ε and λ to obtain a convergent solution.
Explicit numerical solutions can be also obtained. For example, using first-order (in time) finite-difference approximation (FDM) in mean kernel x (0) Equation (6) to get where x (2) 2 dt 1 dt 2 which are usually small compared with the square of the mean kernel x In this section, the subscripts i , j and k are used as indices for the time variables t , t 1 and t 2 , respectively. For the first kernel x (1) Equation (7), we get where In the current work, the integrals I 1 , I 2 and I 12 are computed using the midpoint integration rule. For the second kernel x (2) Equation (8), we get Entropy 2020, 22, 562 The time step ∆t should be used to guarantee convergence of System (9)- (11). Applying the convergence criteria of the fixed-point iteration by differentiating right-hand side of Equation (9) with respect to x (0) after neglecting I 1 and I 2 compared with x The condition ∆t < 1/a is a sufficient condition for the time step used in FDM for convergence. Due to similarity in the other two Equations (10) and (11), the same condition can be also used in the numerical FDM approximations for x (1) and x (2) .
The mean and variance will be computed as The integral I 1 represents the Gaussian part of the model variance and I 2 is the nonGaussian part of the total variance.
In the current work, the numerical simulations will consider the case of a > λ 2 /2 which should converge, as descried above in Section 4, to the stationary gamma distribution. Other values of the parameter a are analyzed similarly as will be shown below. Figure 2 shows the time variation of the mean kernel x (0) and the variance components: total variance, Var_1 = I 1 and Var_2 = I 2 . The convergence of the kernels to steady-states is clear in the figure for the given parameters (∆t = 0.05, a = 0.5, ε = 0.01, λ = 0.01). The steady-state values for mean, total variance, Var_1 and Var_2 are 49.98, 0.2566, 0.2565 and 4.6 × 10 −6 , respectively. The nonGaussian part of the variance is small for the given parameters. The decay of the kernels is a known property and one of the known advantages of WHE.
To validate the results, EM technique with 10,000 samples is considered. Figure 3 shows a comparison between EM and second-order WHE. The EM requires large number of samples to get a smooth solution and hence it will be of low efficiency. The convergence rate of EM is inversely proportional to the square root of number of samples. This issue makes EM insufficient when analyzing nonlinear and/or nonGaussian stochastic processes compared with WHE. The second-order WHE is used efficiently to simulate the dynamics in this case.
The effect of λ on the total variance is shown in Figure 4. We can note the direct proportional relation between λ and the total variance. To validate the results, EM technique with 10,000 samples is considered. Figure 3 shows a comparison between EM and second-order WHE. The EM requires large number of samples to get a smooth solution and hence it will be of low efficiency. The convergence rate of EM is inversely proportional to the square root of number of samples. This issue makes EM insufficient when analyzing nonlinear and/or nonGaussian stochastic processes compared with WHE. The secondorder WHE is used efficiently to simulate the dynamics in this case.  The effect of λ on the total variance is shown in Figure 4. We can note the direct proportional relation between λ and the total variance.  The effect of λ on the total variance is shown in Figure 4. We can note the direct proportional relation between λ and the total variance. To quantify the effect of the parameter λ , Figure 5 shows the variance, Gaussian only and Gaussian with nonGaussian, for λ = 0.01 and λ = 0.02. We can note that as λ increases, the nonGaussian variance and hence the total variance increases as well. This reflects the nonlinearity, and hence the nonGaussian nature of the population model. The model is sensitive to the value of λ . From Figure 5 we can estimate the nonGaussian effect compared with the Gaussian response of the model. For λ = 0.01, the nonGaussian variance contribution is only 1.4% of the total variance, while for λ = 0.02, the nonGaussian variance contribution is 20.4% of the total variance. Estimating the nonGaussian variance contribution is an advantage in WHE over the EM and sampling techniques. To quantify the effect of the parameter λ, Figure 5 shows the variance, Gaussian only and Gaussian with nonGaussian, for λ = 0.01 and λ = 0.02. We can note that as λ increases, the nonGaussian variance and hence the total variance increases as well. This reflects the nonlinearity, and hence the nonGaussian nature of the population model. The model is sensitive to the value of λ. From Figure 5 we can estimate the nonGaussian effect compared with the Gaussian response of the model. For λ = 0.01, the nonGaussian variance contribution is only 1.4% of the total variance, while for λ = 0.02, the nonGaussian variance contribution is 20.4% of the total variance. Estimating the nonGaussian variance contribution is an advantage in WHE over the EM and sampling techniques. The nonGaussian part of the variance is shown in Figure 6 against λ . We can note that the nonGaussian part of the variance increases to a peak value near the inflection of the mean population and then start to decay to a uniform value. This can be explained as the system dynamics is nonlinear and severe no-Gaussian around the curve inflection. The steady-state variance components with different values of λ are shown in Figure 7 and Table 1. We can note, by interpolation, that the steady-state total variance is proportional to 2 λ while the steady-state nonGaussian variance component is proportional to 4 λ . The nonGaussian part of the variance is shown in Figure 6 against λ. We can note that the nonGaussian part of the variance increases to a peak value near the inflection of the mean population and then start to decay to a uniform value. This can be explained as the system dynamics is nonlinear and severe no-Gaussian around the curve inflection. The nonGaussian part of the variance is shown in Figure 6 against λ . We can note that the nonGaussian part of the variance increases to a peak value near the inflection of the mean population and then start to decay to a uniform value. This can be explained as the system dynamics is nonlinear and severe no-Gaussian around the curve inflection. The steady-state variance components with different values of λ are shown in Figure 7 and Table 1. We can note, by interpolation, that the steady-state total variance is proportional to 2 λ while the steady-state nonGaussian variance component is proportional to 4 λ . The steady-state variance components with different values of λ are shown in Figure 7 and Table 1. We can note, by interpolation, that the steady-state total variance is proportional to λ 2 while the steady-state nonGaussian variance component is proportional to λ 4 . The steady-state variance components with different values of λ are shown in Figure 7 and Table 1. We can note, by interpolation, that the steady-state total variance is proportional to

Stochastic Model with Random Parameters
In the case of random variation of one or more of the model parameters, the variance will increase when compared with the case of only noise as the source of randomness. The kernels x ( j) ; j ≤ M can be further expanded/decomposed using gPC. For example, we assume the parameter a is a random parameter that depends on a set of standard random variables with known distribution. This means we can write where the subscript k is the index of gPC term. Then, the kernels x ( j) ; j ≤ M are decomposed as  Multiply Equation (12) by ψ k ; k ≥ 0 and take the average with respect to gPC basis to get To derive the equivalent deterministic system, we need the expected values for the products ψ k ψ j and c ijk = ψ i ψ j ψ k . From the orthogonality property of functionals ψ k ; k ≥ 0, we use ψ k ψ j = δ ij , where δ is the Kronecker delta function. Details about computing expectations c ijk can be found in [20]. Similar expressions for x (1) and x (2) are straight-forward.
In the case where the parameter a depends only one random variable ψ 1 , we can write a(ω) = a 0 + a 1 ψ 1 . This yields the following for x (0) : Similarly, for x (1) : and for x (2) : When considering only the Gaussian part of the solution with only one random parameter, the system will be reduced to the following: Similar expressions can be obtained for the random variation in the parameter ε. We can summarize the above decomposition in the case of combined noise and random parameters with the expression Which can be applied in any order either WHE-gPC or gPC-WHE. Using WHE-gPC is more efficient and results in simpler deterministic system.
In the above derivation, we assumed independency between noise and the random parameters. This happens usually in cases of uncertainties occurring from different uncorrelated sources. In this case, the statistical properties, mean and total variance, are computed as follows: The total variance Var[X] can be analyzed/decomposed into three components: variance Var par due to random parameters, variance Var noise due to noise and variance Var mix due to the mix between noise and the random parameters, i.e., In the case where only one random parameter in the model SDE and only the Gaussian solution is to be considered, the variance components will be as follows: In the case of dependency between the random parameters and/or the noise, the same above derivation can be extended after considering the covariance between the dependent sources. Alternatively, the WCE technique can be used at which the noise is approximated by a set of random variables. In this case, the system will be affected only by uncertainties due to a set of random parameters [21]. The drawback of using WCE is the low convergence due to the noise approximation by a few number of random variables. Using WHE is efficient compared with many other techniques [22].

Numerical Example with Combined Randomness
Assume that the parameter a fluctuates uniformly around the mean with deviation 2% of the mean value a 0 , i.e., a = 0.5 + 0.01 ψ 1 with ψ 1 ∼ U[−1, 1]. This is equivalent to a uniform distribution a ∼ U[0.49, 0.51]. Using λ = 0.02, ε = 0.01 and zero initial conditions for all kernels except for x (0) 0 (t = 0) = 0.5. The total variance and its components, due to noise and due to random parameters, are shown Figure 8. The variance due to mixed contributions are very small compared with other variance components. It is four orders of magnitude smaller than other variance components and hence will be neglected in the analysis. The system sensitivity indices [15], due to different sources of randomness, can be estimated from the variance components. The sensitivity index is the ratio of the variance component to the total variance. For example, the sensitivity index noise S due to noise is calculated as The sensitivity indices due to random parameter a and due to noise are shown in Table 2. The model is 100% sensitive to noise in the case of the parameter a, which is deterministic. Any deviations in the parameter a will affect the system sensitivity indices. As the deviations in a increase, the system response deviates further (i.e., variance increases) and becomes more sensitive toward the deviations in a. This result is also shown in Figure 9, which compares the sensitivity indices for a wide range of a deviations.   For 2% deviations, the steady-state total variance is 1.982, while it is 1.013 in the case of the deterministic parameter a,. This means that 2% deviations in a result in 95% deviation in the steady-state total variance. For 3% deviations, the steady-state variance is 3.22, while it is 1.013 in the case of the deterministic parameter a. This means that 3% deviations in a result in 218% deviation in the steady-state total variance. The model is very sensitive to the deviations in parameter a.
The system sensitivity indices [15], due to different sources of randomness, can be estimated from the variance components. The sensitivity index is the ratio of the variance component to the total variance. For example, the sensitivity index S noise due to noise is calculated as The sensitivity indices due to random parameter a and due to noise are shown in Table 2. The model is 100% sensitive to noise in the case of the parameter a, which is deterministic. Any deviations in the parameter a will affect the system sensitivity indices. As the deviations in a increase, the system response deviates further (i.e., variance increases) and becomes more sensitive toward the deviations in a. This result is also shown in Figure 9, which compares the sensitivity indices for a wide range of a deviations.
To test different scenarios, the developed system with mixed uncertainties is simulated in the case of a < λ 2 /2 where the population extinction is expected. Figure 10 shows the mean and total variance in the case of a(ω) = 0.0001 + 0.00004 ψ 1 , λ = 0.02, ε = 0.01 and zero initial conditions for all kernels except for x (0) 0 (t = 0) = 0.5. In this case, the main population decays slowly to zero and the variance will be mainly due to noise (more than 99.9% of the total variance). This means the model is not sensitive to the random parameter variations in this case.  To test different scenarios, the developed system with mixed uncertainties is simulated in the case of 2 / 2 a λ < where the population extinction is expected. Figure 10 shows the mean and total variance in the case of ( ) The above analysis can be extended to many random parameters in addition to the noise imposed to the model. In all cases, we shall need only to solve a system of few number of deterministic equations that can be analyzed with the available well-known analytical and/or numerical techniques. This gives a great advantage over the sampling techniques that require a huge number of runs to estimate reliable statistics of the model under consideration.

Conclusions
In this paper, a decomposition technique based on spectral random and stochastic-combined expansion is considered. The combined expansion is a tool that can be applied to models under noise and random parameter variations such as the population growth model. The introduced technique can be applied sequentially in any order to derive an equivalent system that is deterministic and can be analyzed using the well-known analytical and/or numerical techniques.  To test different scenarios, the developed system with mixed uncertainties is simulated in the case of 2 / 2 a λ < where the population extinction is expected. Figure 10 shows the mean and total variance in the case of ( ) x t = = 0.5. In this case, the main population decays slowly to zero and the variance will be mainly due to noise (more than 99.9% of the total variance). This means the model is not sensitive to the random parameter variations in this case. The above analysis can be extended to many random parameters in addition to the noise imposed to the model. In all cases, we shall need only to solve a system of few number of deterministic equations that can be analyzed with the available well-known analytical and/or numerical techniques. This gives a great advantage over the sampling techniques that require a huge number of runs to estimate reliable statistics of the model under consideration.

Conclusions
In this paper, a decomposition technique based on spectral random and stochastic-combined expansion is considered. The combined expansion is a tool that can be applied to models under noise and random parameter variations such as the population growth model. The introduced technique can be applied sequentially in any order to derive an equivalent system that is deterministic and can be analyzed using the well-known analytical and/or numerical techniques. The above analysis can be extended to many random parameters in addition to the noise imposed to the model. In all cases, we shall need only to solve a system of few number of deterministic equations that can be analyzed with the available well-known analytical and/or numerical techniques. This gives a great advantage over the sampling techniques that require a huge number of runs to estimate reliable statistics of the model under consideration.

Conclusions
In this paper, a decomposition technique based on spectral random and stochastic-combined expansion is considered. The combined expansion is a tool that can be applied to models under noise and random parameter variations such as the population growth model. The introduced technique can be applied sequentially in any order to derive an equivalent system that is deterministic and can be analyzed using the well-known analytical and/or numerical techniques. The stochastic population model is analyzed using WHE, and the model statistics are obtained using the statistical properties of WHE basis functions. The gPC decomposition technique is then applied to analyze the effect of the random parameters such as the growth rate. The proposed technique allows us to decompose the total variance into components due to noise and due to random parameters. This allows for estimating the model sensitivity due to different sources of randomness. The proposed technique can be extended to allow for analyzing models with many sources of randomness such as the population models.