Analysis of the Stochastic Population Model with Random Parameters

Abstract

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.

1. 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 modeled by adding a noise term [4] or by assuming random variations for one or more parameters e.g., stochastic or random carrying capacity [5].

There are some differences between dynamics and properties of the deterministic and stochastic models. The most important one is the long-time asymptotic dynamics. Other differences include the outbreak probability and the distribution of the final population size [6]. Additionally, the noise term in the population model plays a feedback rule and forces to the system dynamics to stabilize. In the stochastic model, the noise prevents the explosion of the population and is responsible for extinction for one of the species.

The population stochastic differential equation (SDE) is used as a model for the birth–death and population growth processes. As an example, the size of tumor submitted to radiotherapy or phototherapy treatment can be modeled and interpreted mathematically with the population SDE. There are three basic factors affecting the tumor growth. The first is the natural duplication growth of the cells, the second is the self-limitation (may be due to limited available space) of the tumor and the third is due to the treatment. The noise term accounts for the uncertainty due to external influences of environmental factors and/or interaction with other species or due to internal tumor growth [7]. New interesting effects appear due to the presence of the stochastic variations. Randomness introduces more flexibility than the deterministic setting.

Although there is an analytical solution for the 1D population model, the statistical properties are not easily obtained. The problem becomes more complicated for higher dimensions and hence the numerical and spectral techniques should be used. Many numerical techniques based on realizations (sampling) are used to analyze the stochastic population model in 1D and higher dimensions. The numerical techniques are simple to implement but have a slow convergence rate. The Euler–Maruyama (EM) technique has only 0.5 convergence order, while the Milstein (linearized to Runge–Kutta) schemes has convergence of order 1. Some variance-reduction enhancements are done to increase the convergence rate of the numerical techniques. The quasi-Monte Carlo (QMC), as an example, is combined with the tau-leaping technique to enhance the convergence of the stochastic biological models [8], but still the rate is not that efficient. This causes the numerical techniques to be not practical for real applications.

On the other hand, the decomposition Fourier-like techniques such as the Wiener–Hermite expansion (WHE) [9] and the Wiener chaos expansion (WCE) [10,11] have higher orders of convergence. The convergence order of the spectral techniques can be exponential for linear and near-Gaussian processes [12].

For models with random, the polynomial chaos expansion (PCE) motivated by the work of Ghanem and Spanos in 1991 and its generalization (gPC) technique [13] are used efficiently to analyze models with random parameters. They are also efficient techniques in estimating the system sensitivity indices [14,15].

Practically, models that account for both random parameters and stochastic variations are needed for accurate analysis. In the current work, population models with both random parameters and noise variations are considered. The population SDE with random parameters is analyzed with a combination of WHE and gPC decompositions. An equivalent deterministic system is derived and analyzed numerically. The solution statistics, such as mean and variance, are obtained. The effects of the different parameters are quantified, and hence the sensitivity indices are estimated. The results obtained help to understand and control the growth of the population size.

This paper is organized as follows: Section 2 introduces the spectral/decomposition techniques commonly used in the case of noise and random model parameters. Both deterministic and stochastic population models are introduced in Section 3 and Section 4. In Section 5, the stochastic population model with noise imposed is analyzed using WHE. Section 6 extends the stochastic model to have random parameters. A numerical example is analyzed in Section 7 with computations of the sensitivity indices.

2. Stochastic Spectral Techniques

Assume that the model contains a set of random parameters $U\left(\omega \right)$ and depend on a set of $n$ real-valued random variables $\zeta =\left\{{\zeta }_1,\cdots ,{\zeta }_n\right\}$ with predefined probability distributions i.e., $U\left(\omega \right)=U\left(\zeta \left(\omega \right)\right)$. Consider the orthonormal basis functionals $\left\{{\psi }_k\left(\zeta \right);\hspace{0.17em}k\in \mathbb{N}\right\}$ for the space L² of second-order functionals in the random variables $\zeta \left(\omega \right)$. Many choices for the basis functionals ${\psi }_k$ including multiwavelets or the multivariate polynomials [12]. Any second-order random process $X\left(t,\hspace{0.17em}\omega \right)$ that depends on $U\left(\omega \right)$ can then be decomposed as follows [16]:

$$X\left(t,\hspace{0.17em}\omega \right)={\displaystyle \sum _{k\in \mathbb{N}}{x}_k\left(t\right)\hspace{0.17em}\hspace{0.17em}}{\psi }_k\left(\zeta (\omega )\right).$$

(1)

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.,

$$\underset{P\to \infty }{\lim }\hspace{0.17em}{\left(X-{\displaystyle {\displaystyle \sum _{k=0}^P}{x}_k\hspace{0.17em}{\psi }_k}\right)}^2=0.$$

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]:

$$X={\displaystyle {\displaystyle \sum }_{j=0}^M}{\displaystyle {\displaystyle \underset{-\infty }{\overset{\infty }{\int }}}{\displaystyle {\displaystyle \underset{-\infty }{\overset{t_j}{\int }}}{\displaystyle {\displaystyle \underset{-\infty }{\overset{t_{j-1}}{\int }}}\cdots {\displaystyle {\displaystyle \underset{-\infty }{\overset{t_2}{\int }}}{x}_{}^{(j)}(t_1,t_2,\cdots ,t_j)\hspace{0.17em}\hspace{0.17em}dW(t_1)dW(t_2)\cdots dW(t_j)}}}}\hspace{0.17em}.$$

The functionals $x_{}^{(j)}(t_1,t_2,\cdots ,t_j)$ are usually assumed to be symmetric in the variables $t_1,\hspace{0.17em}{t}_2,\hspace{0.17em}\dots ,\hspace{0.17em}{t}_j$. But $dW(t_1)dW(t_2)\cdots dW(t_j)=H^{(j)}(t_1,t_2,\cdots ,t_j)\hspace{0.17em}d{t}_{1}d{t}_2\dots \hspace{0.17em}d{t}_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

$$X={\displaystyle {\displaystyle \sum }_{j=0}^M}\underset{R^j}{{\displaystyle \int }}\hspace{0.17em}{x}_{}^{(j)}\hspace{0.17em}{H}^{(j)}\hspace{0.17em}d{\tau }_j,$$

(2)

where $d{\tau }_j=d{t}_{1}d{t}_2\dots \hspace{0.17em}d{t}_j$ and the integral ${\displaystyle {\int }_{R^j}}$ is a j-dimensional iterated integral over the variables $t_1,\hspace{0.17em}{t}_2,\hspace{0.17em}\dots ,\hspace{0.17em}{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:

$$\frac{dX}{dt}=a\hspace{0.17em}X\hspace{0.17em}dt+\lambda \hspace{0.17em}X\frac{dW}{dt}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em};\hspace{0.17em}\hspace{0.17em}x(0)=x_0\hspace{0.17em},$$

using WHE Expansion (2) to get

$$\frac{d}{dt}\left({\displaystyle {\displaystyle \sum }_{j=0}^M}\hspace{0.17em}\underset{R^j}{{\displaystyle \int }}\hspace{0.17em}{x}_{}^{\left(j\right)}\hspace{0.17em}{H}^{\left(j\right)}\hspace{0.17em}d{\tau }_j\right)=a\hspace{0.17em}{\displaystyle {\displaystyle \sum }_{j=0}^M}\hspace{0.17em}\underset{R^j}{{\displaystyle \int }}\hspace{0.17em}{x}_{}^{\left(j\right)}\hspace{0.17em}{H}^{\left(j\right)}\hspace{0.17em}d{\tau }_j\hspace{0.17em}+\lambda \hspace{0.17em}{H}^{(1)}{\displaystyle {\displaystyle \sum }_{j=0}^M}\hspace{0.17em}\underset{R^j}{{\displaystyle \int }}\hspace{0.17em}{x}_{}^{\left(j\right)}\hspace{0.17em}{H}^{\left(j\right)}\hspace{0.17em}d{\tau }_j$$

multiply by $H^{\left(0\right)}=1$ and $H^{\left(1\right)}$, respectively, and then apply the expectation to get

$$\begin{gathered} \frac{d{x}_{}^{\left(0\right)}}{dt}=a\hspace{0.17em}{x}_{}^{\left(0\right)}, \\ \frac{d{x}_{}^{\left(1\right)}}{dt}=a\hspace{0.17em}{x}_{}^{\left(1\right)}\hspace{0.17em}+\lambda \hspace{0.17em}{x}_{}^{\left(0\right)}\hspace{0.17em}\delta \left(t-t_1\right), \end{gathered}$$

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_{}^{\left(0\right)}(t)=x_0\hspace{0.17em}{e}^{a\hspace{0.17em}t}\mathrm{and}{x}_{}^{\left(1\right)}(t,t_1)=\lambda x_0\hspace{0.17em}{e}^{-a\hspace{0.17em}t};\hspace{0.17em}\hspace{0.17em}t\ge t_1$$

and hence

$$E[X]=x_{}^{\left(0\right)}(t)=x_0\hspace{0.17em}{e}^{a\hspace{0.17em}t}\mathrm{and}Var\left[X\right]={\displaystyle {\displaystyle \underset{t_1=0}{\overset{t_1=t}{\int }}}{\left(x_{}^{\left(1\right)}\right)}^{2}d{t}_1}=t{\left(\lambda x_0\hspace{0.17em}{e}^{-a\hspace{0.17em}t}\right)}^2.$$

3. The Deterministic Population Model

Consider the transport population model with nonlinear losses:

$$\frac{dX}{dt}=a\hspace{0.17em}X\hspace{0.17em}-\hspace{0.17em}\epsilon \hspace{0.17em}{X}^2\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em};\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}x(0)=x_0\hspace{0.17em}\hspace{0.17em}\hspace{0.17em},\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}t>0,$$

(3)

where the growth rate $a=\lambda -\mu$, which is the difference between the birth rate $\lambda$ and death rate $\mu$. The model parameters, $a$ and $\epsilon$, are assumed deterministic in this case. The model Equation (3) is of Bernoulli type of order 2 and has the exact solution:

$$X(t)=\frac{a\hspace{0.17em}{x}_0}{\epsilon x_0+\left(a-\epsilon x_0\right)\hspace{0.17em}{e}^{-at}}=\frac{L\hspace{0.17em}}{1+\left(L/x_0-1\right)\hspace{0.17em}{e}^{-at}},$$

where $\epsilon /\left|a\right|=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.

A more general form of (3) is called Rechards model [19] and it takes the form

$$\frac{dX}{dt}=a\hspace{0.17em}X\hspace{0.17em}-\hspace{0.17em}\epsilon \hspace{0.17em}{X}^{\beta }\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em};\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}x(0)=x_0\hspace{0.17em}\hspace{0.17em}\hspace{0.17em},\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}t>0.$$

Here, the exponent $\beta >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

$$X(t)=\frac{L\hspace{0.17em}}{{\left(1+\left(L/x_0-1\right)\hspace{0.17em}{e}^{-at}\right)}^{\frac{1}{\beta -1}}}.$$

In the current work, we shall consider the original Verhulst model which is Richards model with $\beta =2$.

4. The Stochastic Model

Consider the stochastic transport population model with nonlinear losses:

$$dX=\left(a\hspace{0.17em}X\hspace{0.17em}-\hspace{0.17em}\epsilon \hspace{0.17em}{X}^2\right)dt+\lambda \hspace{0.17em}dW\hspace{0.17em}X\hspace{0.17em}\hspace{0.17em};\hspace{0.17em}\hspace{0.17em}x(0)=x_0\hspace{0.17em}\hspace{0.17em}\hspace{0.17em},\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}t>0,$$

(4)

with $\lambda$ is a diffusion coefficient. In the case of deterministic parameters, $a$ and $\epsilon$, 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 $\epsilon$ 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.,

$$E\left[{\left(X|U\right)}^2\right]<\infty ,E\left[{\left(X|W\right)}^2\right]<\infty \mathrm{and}E\left[X^2\right]<\infty ,$$

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$:

$$X(t)=\frac{e^{\left(a\hspace{0.17em}-\hspace{0.17em}\frac{1}{2}{\lambda }^2\right)\hspace{0.17em}t+\lambda \hspace{0.17em}B(t)}}{\frac{1}{x_0}+\epsilon {\displaystyle {\displaystyle \underset{0}{\overset{t}{\int }}}{e}^{\left(a\hspace{0.17em}-\hspace{0.17em}\frac{1}{2}{\lambda }^2\right)\hspace{0.17em}t+\lambda \hspace{0.17em}B(t)}\hspace{0.17em}dt}},$$

(5)

which reduces to the deterministic solution in the case where the diffusion parameter $\lambda =0$. The stochastic process $X(t)$ given in (4) and for $a>{\lambda }^2/2$ is recurrent and converges to the stationary probability gamma distribution $\gamma \left(2a/{\lambda }^2-1\hspace{0.17em},\hspace{0.17em}\hspace{0.17em}{\lambda }^2/2\epsilon \right)$. For $a<{\lambda }^2/2$, the process decays/converges almost sure to zero [7]. In the case in which $a={\lambda }^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.

5. Analysis Using WHE

Use WHE Expansion (2) in SDE (4) to get

$$\begin{array}{ll}\frac{d}{dt}\left({\displaystyle {\displaystyle \sum }_{j=0}^M}\hspace{0.17em}\underset{R^j}{{\displaystyle \int }}\hspace{0.17em}{x}_{}^{\left(j\right)}\hspace{0.17em}{H}^{\left(j\right)}\hspace{0.17em}d{\tau }_j\right)=a\hspace{0.17em}{\displaystyle {\displaystyle \sum }_{j=0}^M}\hspace{0.17em}\underset{R^j}{{\displaystyle \int }}\hspace{0.17em}{x}_{}^{\left(j\right)}& H^{\left(j\right)}\hspace{0.17em}d{\tau }_j\hspace{0.17em}-\epsilon \hspace{0.17em}{\left({\displaystyle {\displaystyle \sum }_{j=0}^M}\hspace{0.17em}\underset{R^j}{{\displaystyle \int }}\hspace{0.17em}{x}_{}^{\left(j\right)}\hspace{0.17em}{H}^{\left(j\right)}\hspace{0.17em}d{\tau }_j\right)}^2\\ & +\lambda \hspace{0.17em}{H}^{(1)}{\displaystyle {\displaystyle \sum }_{j=0}^M}\hspace{0.17em}\underset{R^j}{{\displaystyle \int }}\hspace{0.17em}{x}_{}^{\left(j\right)}\hspace{0.17em}{H}^{\left(j\right)}\hspace{0.17em}d{\tau }_j.\end{array}$$

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_{}^{\left(2\right)}$, apply WHE algorithm to get

$$\frac{d{x}_{}^{\left(0\right)}}{dt}=a\hspace{0.17em}{x}_{}^{\left(0\right)}\hspace{0.17em}-\hspace{0.17em}\epsilon \hspace{0.17em}\hspace{0.17em}\left[{\left(x_{}^{\left(0\right)}\right)}^2+{\displaystyle \underset{R}{\int }{\left(x_{}^{\left(1\right)}\right)}^2}d{t}_1+2{\displaystyle \underset{R^2}{\int }{\left(x_{}^{\left(2\right)}\right)}^2}d{t}_{1}d{t}_2\right],$$

(6)

$$\frac{d{x}_{}^{\left(1\right)}}{dt}=a\hspace{0.17em}{x}_{}^{\left(1\right)}\hspace{0.17em}-\hspace{0.17em}2\epsilon \hspace{0.17em}{x}_{}^{\left(0\right)}{x}_{}^{\left(1\right)}-\hspace{0.17em}4\epsilon {\displaystyle \underset{R}{\int }{x}_{}^{\left(1\right)}\left(t_1\right)x_{}^{\left(2\right)}\left(t_1,t_2\right)}\hspace{0.17em}d{t}_2+\lambda \hspace{0.17em}\delta (t-t_1)\hspace{0.17em}{x}_{}^{\left(0\right)},$$

(7)

$$\frac{d{x}_{}^{\left(2\right)}}{dt}=a\hspace{0.17em}{x}_{}^{\left(2\right)}\hspace{0.17em}-\hspace{0.17em}2\epsilon \hspace{0.17em}{x}_{}^{\left(0\right)}{x}_{}^{\left(2\right)}-\epsilon \hspace{0.17em}{x}^{(1)}\left(t_1\right)x^{(1)}\left(t_2\right)+\lambda \hspace{0.17em}\delta (t-t_2)\hspace{0.17em}{x}^{(1)}\left(t_1\right),$$

(8)

where $\delta (.)$ is the Dirac delta function. The initial conditions, after applying WHE, will be:

$$x_{}^{\left(0\right)}(0)=x_0,x_{}^{\left(1\right)}(0,\hspace{0.17em}{t}_1)=zero\mathrm{and}{x}_{}^{\left(2\right)}(0,\hspace{0.17em}{t}_1,\hspace{0.17em}{t}_2)=zero.$$

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 $\epsilon$ and $\lambda$ 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_{}^{\left(0\right)}$ Equation (6) to get

$$x_{{}_{i+1}}^{\left(0\right)}=\left(1+a\hspace{0.17em}\Delta t\right)x_{{}_i}^{\left(0\right)}-\epsilon \hspace{0.17em}\Delta t\left[{\left(x_i^{\left(0\right)}\right)}^2+I_1+I_2\right],$$

(9)

where $I_1={\displaystyle \underset{R}{\int }{\left(x_{}^{\left(1\right)}\right)}^2}d{t}_1$ and $I_2=2{\displaystyle \underset{R^2}{\int }{\left(x_{}^{\left(2\right)}\right)}^2}d{t}_{1}d{t}_2$ which are usually small compared with the square of the mean kernel ${\left(x_i^{\left(0\right)}\right)}^2$. In this section, the subscripts $i\hspace{0.17em},\hspace{0.17em}j$ and $k$ are used as indices for the time variables $t\hspace{0.17em},\hspace{0.17em}\hspace{0.17em}{t}_1$ and $t_2$, respectively.

For the first kernel $x_{}^{\left(1\right)}$ Equation (7), we get

$$x_{i+1,j}^{\left(1\right)}=\left(1+a\Delta t\right)x_{i,j}^{\left(1\right)}-2\epsilon \Delta t{x}_i^{\left(0\right)}{x}_{i,j}^{\left(1\right)}-4\epsilon \Delta t{I}_{12}+\lambda {\delta }_{ij}{x}_i^{\left(0\right)},$$

(10)

where $I_{12}={\displaystyle \underset{R}{\int }{x}_{}^{\left(1\right)}\left(t_1\right)x_{}^{\left(2\right)}\left(t_1,t_2\right)}\hspace{0.17em}d{t}_2$. In the current work, the integrals $I_1\hspace{0.17em},\hspace{0.17em}\hspace{0.17em}{I}_2$ and $I_{12}$ are computed using the midpoint integration rule.

For the second kernel $x_{}^{\left(2\right)}$ Equation (8), we get

$$x_{i+1,j,k}^{\left(2\right)}=x_{i,j,k}^{\left(2\right)}+a\Delta t{x}_{i,j,k}^{\left(2\right)}-2\epsilon \Delta t{x}_i^{\left(0\right)}{x}_{i,j,k}^{\left(2\right)}-\epsilon \Delta t{x}_{i,j}^{(1)}{x}_{i,k}^{(1)}+\lambda {\delta }_{ik}{x}_{i,j}^{(1)}.$$

(11)

The time step $\Delta 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_{}^{\left(0\right)}$ after neglecting I₁ and I₂ compared with ${\left(x_i^{\left(0\right)}\right)}^2$ to get

$$\begin{array}{l}1+a\hspace{0.17em}\Delta t-2\epsilon \hspace{0.17em}\Delta t\hspace{0.17em}{x}_i^{\left(0\right)}<0\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\to \hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\left|a\hspace{0.17em}\Delta t-2\epsilon \hspace{0.17em}\Delta t\hspace{0.17em}{x}_i^{\left(0\right)}\right|<1\\ \therefore \hspace{0.17em}\hspace{0.17em}\Delta t<{\left|a\hspace{0.17em}-2\epsilon \hspace{0.17em}{x}_i^{\left(0\right)}\right|}^{-1}\to \Delta t<\max {\left|a\hspace{0.17em}-2\epsilon \hspace{0.17em}{x}_i^{\left(0\right)}\right|}^{-1}\to \Delta t<{\epsilon }^{-1}{\left|\frac{a}{\epsilon }\hspace{0.17em}-2\hspace{0.17em}{x}_i^{\left(0\right)}\right|}^{-1}\\ \to \Delta t<{\epsilon }^{-1}{\left|L\hspace{0.17em}-2\hspace{0.17em}{x}_i^{\left(0\right)}\right|}^{-1}\to \Delta t<{\epsilon }^{-1}{L}^{-1}\to \hspace{0.17em}\Delta t<1/a\\ \mathrm{Note}:\hspace{0.17em}\max \left[x_i^{\left(0\right)}\right]=a/\epsilon =L:\mathrm{Carrying}\hspace{0.17em}\mathrm{Capacity}\end{array}$$

The condition $\Delta 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_{}^{\left(1\right)}$ and $x_{}^{\left(2\right)}$.

The mean and variance will be computed as

$$E\left[X\right]=x^{(0)}\mathrm{and}Var\left[X\right]={\displaystyle {\displaystyle \sum _{k=1}^M}\left(k!\right){\displaystyle \underset{R^k}{\int }{\left(x_{}^{\left(k\right)}\right)}^{2}d{\tau }_k}}=I_1+2I_2,$$

where

$$I_1={\displaystyle \underset{R}{\int }{\left(x_{}^{\left(1\right)}\right)}^2}d{t}_1\mathrm{and}\hspace{0.17em}{I}_2=2{\displaystyle \underset{R^2}{\int }{\left(x_{}^{\left(2\right)}\right)}^2}d{t}_{1}d{t}_2.$$

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>{\lambda }^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_{}^{\left(0\right)}$ 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 ($\Delta t$ = 0.05, $a$ = 0.5, $\epsilon$ = 0.01, $\lambda$ = 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⁻⁶, 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 $\lambda$ on the total variance is shown in Figure 4. We can note the direct proportional relation between $\lambda$ and the total variance.

To quantify the effect of the parameter $\lambda$, Figure 5 shows the variance, Gaussian only and Gaussian with nonGaussian, for $\lambda$ = 0.01 and $\lambda$ = 0.02. We can note that as $\lambda$ 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 $\lambda$. From Figure 5 we can estimate the nonGaussian effect compared with the Gaussian response of the model. For $\lambda$ = 0.01, the nonGaussian variance contribution is only 1.4% of the total variance, while for $\lambda$ = 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 $\lambda$. 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 $\lambda$ are shown in Figure 7 and

Table 1. Steady-state variance components with $\lambda$.

$\mathit{\lambda }$	Total Var.	NonGaussian Var.
0.0000	0.000	0
0.0100	0.257	4.59 × 10−6
0.0150	0.593	2.50 × 10−5
0.0175	0.835	5.10 × 10−5
0.0200	1.220	1.21 × 10−4

. We can note, by interpolation, that the steady-state total variance is proportional to ${\lambda }^2$ while the steady-state nonGaussian variance component is proportional to ${\lambda }^4$.

6. 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\le 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

$$a\left(\omega \right)={\displaystyle {\displaystyle \sum _{k=0}^P}{a}_k\hspace{0.17em}{\psi }_k},$$

where the subscript k is the index of gPC term. Then, the kernels $x^{(j)};j\le M$ are decomposed as

$$x^{(j)}={\displaystyle {\displaystyle \sum _{k=0}^P}{x}_k^{(j)}\hspace{0.17em}{\psi }_k}\hspace{0.17em}\hspace{0.17em};\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}j\le M$$

with $x_0^{(j)}$ as the mean and ${\displaystyle {\displaystyle \sum _{k=1}^P}{\left[x_k^{(j)}\right]}^2\hspace{0.17em}}$ as the variance of $x^{(j)};j\le M$.

Substitute in $x_{}^{(0)}$ Equation (6) to get

$$\begin{array}{ll}\frac{d}{dt}{\displaystyle {\displaystyle \sum _{k=0}^P}{x}_k^{(0)}{\psi }_k}=a\hspace{0.17em}{\displaystyle {\displaystyle \sum _{k=0}^P}{x}_k^{(0)}{\psi }_k}\hspace{0.17em}-\hspace{0.17em}\epsilon \hspace{0.17em}{\left({\displaystyle {\displaystyle \sum _{k=0}^P}{x}_k^{(0)}{\psi }_k}\right)}^2& \hspace{0.17em}-\hspace{0.17em}\epsilon {\displaystyle \underset{R}{\int }{\left({\displaystyle {\displaystyle \sum _{k=0}^P}{x}_k^{(1)}{\psi }_k}\right)}^2}d{t}_1\\ & -\hspace{0.17em}2\epsilon \hspace{0.17em}{\displaystyle \underset{R^2}{\int }{\left({\displaystyle {\displaystyle \sum _{k=0}^P}{x}_k^{(2)}{\psi }_k}\right)}^2}d{t}_{1}d{t}_2\end{array}$$

(12)

Multiply Equation (12) by ${\psi }_k\hspace{0.17em};\hspace{0.17em}\hspace{0.17em}k\ge 0$ and take the average with respect to gPC basis to get

$$\begin{array}{ll}\frac{d}{dt}{\displaystyle \sum _{k=0}^P}{x}_k^{(0)}& \langle {\psi }_k{\psi }_j\rangle =a\hspace{0.17em}{\displaystyle {\displaystyle \sum _{k=0}^P}{x}_k^{(0)}\langle {\psi }_k{\psi }_j\rangle }\hspace{0.17em}-\hspace{0.17em}\epsilon {\displaystyle {\displaystyle \sum _{k=0}^P}{\displaystyle {\displaystyle \sum _{i=0}^P}{x}_i^{(0)}{x}_k^{(0)}\langle {\psi }_i{\psi }_k{\psi }_j\rangle }}\hspace{0.17em}\\ & -\hspace{0.17em}\epsilon {\displaystyle \underset{R}{\int }{\displaystyle {\displaystyle \sum _{k=0}^P}{\displaystyle {\displaystyle \sum _{i=0}^P}{x}_i^{(1)}{x}_k^{(1)}\langle {\psi }_i{\psi }_k{\psi }_j\rangle }}}\hspace{0.17em}d{t}_1-\hspace{0.17em}2\epsilon \hspace{0.17em}{\displaystyle \underset{R^2}{\int }{\displaystyle {\displaystyle \sum _{k=0}^P}{\displaystyle {\displaystyle \sum _{i=0}^P}{x}_i^{(2)}{x}_k^{(2)}\langle {\psi }_i{\psi }_k{\psi }_j\rangle }}}\hspace{0.17em}d{t}_{1}d{t}_2\end{array}$$

To derive the equivalent deterministic system, we need the expected values for the products $\langle {\psi }_k{\psi }_j\rangle$ and $c_{ijk}=\langle {\psi }_i\hspace{0.17em}{\psi }_j\hspace{0.17em}{\psi }_k\rangle$. From the orthogonality property of functionals ${\psi }_k\hspace{0.17em};\hspace{0.17em}\hspace{0.17em}k\ge 0$, we use $\langle {\psi }_k{\psi }_j\rangle ={\delta }_{ij}$, where $\delta$ 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 ${\psi }_1$, we can write $a\left(\omega \right)=a_0+a_1{\psi }_1$. This yields the following for $x^{(0)}$:

$$\begin{gathered} \frac{d{x}_0^{(0)}}{dt}=\left(a_0{x}_0^{(0)}+a_1{x}_1^{(0)}\hspace{0.17em}-\hspace{0.17em}\epsilon \hspace{0.17em}\left(\begin{array}{l}\hspace{0.17em}{\left(x_0^{\left(0\right)}\right)}^2\hspace{0.17em}+{\displaystyle \underset{R}{\int }{\left(x_0^{\left(1\right)}\right)}^2}d{t}_1+2\hspace{0.17em}{\displaystyle \underset{R^2}{\int }{\left(x_0^{\left(2\right)}\right)}^2}d{t}_{1}d{t}_2\\ +{\left(x_1^{\left(0\right)}\right)}^2\hspace{0.17em}+{\displaystyle \underset{R}{\int }{\left(x_1^{\left(1\right)}\right)}^2}d{t}_1+2\hspace{0.17em}{\displaystyle \underset{R^2}{\int }{\left(x_1^{\left(2\right)}\right)}^2}d{t}_{1}d{t}_2\end{array}\right)\right), \\ \frac{d{x}_1^{(0)}}{dt}=\left(a_0{x}_1^{(0)}+a_1{x}_0^{(0)}\hspace{0.17em}-\hspace{0.17em}2\epsilon \hspace{0.17em}\left(\hspace{0.17em}{x}_0^{\left(0\right)}{x}_1^{\left(0\right)}\hspace{0.17em}+{\displaystyle \underset{R}{\int }{x}_0^{\left(1\right)}{x}_1^{\left(1\right)}}d{t}_1+2\hspace{0.17em}{\displaystyle \underset{R^2}{\int }{x}_0^{\left(2\right)}{x}_1^{\left(2\right)}}d{t}_{1}d{t}_2\right)\right). \end{gathered}$$

Similarly, for $x^{(1)}$:

$$\begin{gathered} \frac{d{x}_0^{(1)}}{dt}=\left(a_0{x}_0^{(1)}+a_1{x}_1^{(1)}\hspace{0.17em}-\hspace{0.17em}\epsilon \hspace{0.17em}\left(\begin{array}{l}2\hspace{0.17em}{x}_0^{\left(0\right)}{x}_0^{\left(1\right)}+\hspace{0.17em}4{\displaystyle \underset{R}{\int }{x}_0^{\left(1\right)}\left(t_1\right)x_0^{\left(2\right)}\left(t_1,t_2\right)}\hspace{0.17em}d{t}_2\\ +2\hspace{0.17em}{x}_1^{\left(0\right)}{x}_0^{\left(1\right)}+\hspace{0.17em}4{\displaystyle \underset{R}{\int }{x}_1^{\left(1\right)}\left(t_1\right)x_1^{\left(2\right)}\left(t_1,t_2\right)}\hspace{0.17em}d{t}_2\end{array}\right)\right)+\lambda \hspace{0.17em}\delta (t-t_1)x_0^{(0)} \\ \frac{d{x}_1^{(1)}}{dt}=\left(a_0{x}_1^{(1)}+a_1{x}_0^{(1)}\hspace{0.17em}-\hspace{0.17em}2\epsilon \hspace{0.17em}\left(\begin{array}{l}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}{x}_0^{\left(0\right)}{x}_1^{\left(1\right)}+\hspace{0.17em}2{\displaystyle \underset{R}{\int }{x}_0^{\left(1\right)}\left(t_1\right)x_1^{\left(2\right)}\left(t_1,t_2\right)}\hspace{0.17em}d{t}_2\\ +\hspace{0.17em}{x}_1^{\left(0\right)}{x}_0^{\left(1\right)}+\hspace{0.17em}2{\displaystyle \underset{R}{\int }{x}_1^{\left(1\right)}\left(t_1\right)x_0^{\left(2\right)}\left(t_1,t_2\right)}\hspace{0.17em}d{t}_2\end{array}\right)\right)+\lambda \hspace{0.17em}\delta (t-t_1)x_1^{(0)} \end{gathered}$$

and for $x^{(2)}$:

$$\begin{gathered} \frac{d{x}_0^{(2)}}{dt}=\left(a_0{x}_0^{(2)}+a_1{x}_1^{(2)}\hspace{0.17em}-\hspace{0.17em}\epsilon \hspace{0.17em}\left(\begin{array}{l}\hspace{0.17em}2x_0^{\left(0\right)}{x}_0^{\left(2\right)}+x_0{}^{(1)}\left(t_1\right)x_0{}^{(1)}\left(t_2\right)\\ +2x_1^{\left(0\right)}{x}_1^{\left(2\right)}+x_1{}^{(1)}\left(t_1\right)x_1{}^{(1)}\left(t_2\right)\end{array}\right)\right)+\lambda \hspace{0.17em}\delta (t-t_2)x_0^{(1)} \\ \frac{d{x}_1^{(2)}}{dt}=\left(a_0\hspace{0.17em}{x}_1^{(2)}+a_1\hspace{0.17em}{x}_0^{(2)}\hspace{0.17em}-\hspace{0.17em}2\epsilon \left(x_0^{(0)}{x}_1^{(2)}+x_1^{(0)}{x}_0^{(2)}+x_0{}^{(1)}\left(t_1\right)x_1{}^{(1)}\left(t_2\right)\right)\right)+\lambda \hspace{0.17em}\delta (t-t_2)x_1^{(1)} \end{gathered}$$

When considering only the Gaussian part of the solution with only one random parameter, the system will be reduced to the following:

$$\begin{gathered} \frac{d{x}_0^{(0)}}{dt}=\left(a_0{x}_0^{(0)}+a_1{x}_1^{(0)}\hspace{0.17em}-\hspace{0.17em}\epsilon \hspace{0.17em}\left(\hspace{0.17em}{\left(x_0^{\left(0\right)}\right)}^2\hspace{0.17em}+{\left(x_1^{\left(0\right)}\right)}^2+{\displaystyle \underset{R}{\int }\left[{\left(x_0^{\left(1\right)}\right)}^2+{\left(x_1^{\left(1\right)}\right)}^2\right]}d{t}_1\right)\right), \\ \frac{d{x}_0^{(1)}}{dt}=\left(a_0{x}_0^{(1)}+a_1{x}_1^{(1)}\hspace{0.17em}-\hspace{0.17em}2\epsilon \hspace{0.17em}\left(x_0^{\left(0\right)}{x}_0^{\left(1\right)}+\hspace{0.17em}{x}_1^{\left(0\right)}{x}_0^{\left(1\right)}\right)\right)+\lambda \hspace{0.17em}\delta (t-t_1)x_0^{(0)}, \\ \frac{d{x}_1^{(0)}}{dt}=\left(a_0{x}_1^{(0)}+a_1{x}_0^{(0)}\hspace{0.17em}-\hspace{0.17em}2\epsilon \hspace{0.17em}\left(\hspace{0.17em}{x}_0^{\left(0\right)}{x}_1^{\left(0\right)}\hspace{0.17em}+{\displaystyle \underset{R}{\int }{x}_0^{\left(1\right)}{x}_1^{\left(1\right)}}d{t}_1\right)\right), \\ \frac{d{x}_1^{(1)}}{dt}=\left(a_0{x}_1^{(1)}+a_1{x}_0^{(1)}\hspace{0.17em}-\hspace{0.17em}2\epsilon \hspace{0.17em}\left(x_0^{\left(0\right)}{x}_1^{\left(1\right)}+\hspace{0.17em}{x}_1^{\left(0\right)}{x}_0^{\left(1\right)}\right)\right)+\lambda \hspace{0.17em}\delta (t-t_1)x_1^{(0)}. \end{gathered}$$

Similar expressions can be obtained for the random variation in the parameter $\epsilon$.

We can summarize the above decomposition in the case of combined noise and random parameters with the expression

$$X(t;\omega )={\displaystyle {\displaystyle \sum _{j=0}^P}{\displaystyle {\displaystyle \sum _{k=0}^M}{\displaystyle \underset{Rk}{\int }{x}_j^{\left(k\right)}(t)\hspace{0.17em}{\psi }_j{H}^{\left(k\right)}d{\tau }_k}}}.$$

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:

$$E[X]=x_0^{(0)}\mathrm{and}Var\left[X\right]={\displaystyle {\displaystyle \sum _{j=0}^P}{\displaystyle {\displaystyle \sum _{\begin{array}{l}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}k=0\\ \hspace{0.17em}\left(j,k\right)\ne \left(0,0\right)\end{array}}^M}k!{\displaystyle \underset{R^j}{\int }{\left(x_j^{\left(k\right)}\right)}^{2}d{\tau }_k}}}.$$

The total variance $Var\left[X\right]$ can be analyzed/decomposed into three components: variance $Va{r}_{par}$ due to random parameters, variance $Va{r}_{noise}$ due to noise and variance $Va{r}_{mix}$ due to the mix between noise and the random parameters, i.e.,

$$Var\left[X\right]=Va{r}_{par}+Va{r}_{noise}+Va{r}_{mix},$$

where

$$Va{r}_{par}={\displaystyle {\displaystyle \sum _{k=1}^P}{\left[x_k^{\left(0\right)}\right]}^2},Va{r}_{noise}={\displaystyle {\displaystyle \sum _{k=1}^M}k!{\displaystyle \underset{R^k}{\int }{\left[x_0^{\left(k\right)}\right]}^{2}d{\tau }_k}},Va{r}_{mix}={\displaystyle {\displaystyle \sum _{j=1}^P}{\displaystyle {\displaystyle \sum _{\hspace{0.17em}k=1}^M}k!{\displaystyle \underset{R^k}{\int }{\left[x_j^{\left(k\right)}\right]}^{2}d{\tau }_k}}}.$$

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:

$$Va{r}_{par}={\left[x_1^{\left(0\right)}\right]}^2,Va{r}_{noise}={\displaystyle \underset{R^{}}{\int }{\left[x_0^{\left(1\right)}\right]}^{2}d{t}_1},Va{r}_{mix}={\displaystyle \underset{R^{}}{\int }{\left[x_1^{\left(1\right)}\right]}^{2}d{t}_1}.$$

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].

7. 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\hspace{0.17em}\hspace{0.17em}{\psi }_1$ with ${\psi }_1\sim U[-1,\hspace{0.17em}1]$. This is equivalent to a uniform distribution $a\sim U[0.49,\hspace{0.17em}\hspace{0.17em}0.51]$. Using $\lambda$ = 0.02, $\epsilon$ = 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.

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

$$S_{noise}=Va{r}_{noise}/Var[X].$$

The sensitivity indices due to random parameter a and due to noise are shown in

Table 2. Sensitivity indices due to noise and the random parameter a.

Parameter a	Noise Sensitivity Index	Parameter a Sensitivity Index
$a\left(\omega \right)=0.5$ (deterministic)	100%	0.0%
$a\left(\omega \right)=0.5+0.005\hspace{0.17em}\hspace{0.17em}{\psi }_1$	79.8%	20.2%
$a\left(\omega \right)=0.5+0.010\hspace{0.17em}\hspace{0.17em}{\psi }_1$	49.1%	49.9%
$a\left(\omega \right)=0.5+0.015\hspace{0.17em}\hspace{0.17em}{\psi }_1$	29.7%	70.3%
$a\left(\omega \right)=0.5+0.020\hspace{0.17em}\hspace{0.17em}{\psi }_1$	18.9%	81.1%
$a\left(\omega \right)=0.5+0.025\hspace{0.17em}\hspace{0.17em}{\psi }_1$	12.8%	87.2%
$a\left(\omega \right)=0.5+0.030\hspace{0.17em}\hspace{0.17em}{\psi }_1$	9.1%	91.9%
$a\left(\omega \right)=0.5+0.035\hspace{0.17em}\hspace{0.17em}{\psi }_1$	6.7%	93.3%

. 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<{\lambda }^2/2$ where the population extinction is expected. Figure 10 shows the mean and total variance in the case of $a\left(\omega \right)=0.0001+0.00004\hspace{0.17em}\hspace{0.17em}{\psi }_1$, $\lambda$ = 0.02, $\epsilon$ = 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.

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.

8. 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.