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
and depend on a set of
real-valued random variables
with predefined probability distributions i.e.,
. Consider the orthonormal basis functionals
for the space
L2 of second-order functionals in the random variables
. Many choices for the basis functionals
including multiwavelets or the multivariate polynomials [
12]. Any second-order random process
that depends on
can then be decomposed as follows [
16]:
For practical computations, the Expansion (1) is truncated to its first
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
can be decomposed using (
) truncated WHE as follows [
9]:
The functionals
are usually assumed to be symmetric in the variables
. But
, where
is the Hermite functional of order
j. For brevity we shall write
without its independent variables. Now we can write
where
and the integral
is a j-dimensional iterated integral over the variables
.
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
multiply by
and
, 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
3. The Deterministic Population Model
Consider the transport population model with nonlinear losses:
where the growth rate
, which is the difference between the birth rate
and death rate
. The model parameters,
and
, are assumed deterministic in this case. The model Equation (3) is of Bernoulli type of order 2 and has the exact solution:
where
is the “crowding term”, the reciprocal of the carrying capacity
(saturation level). For
, the system will reach an infinite limit equals the carrying capacity
and the well-known sigmoid/logistic curve is obtained. For
, 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
Here, the exponent
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 .
4. The Stochastic Model
Consider the stochastic transport population model with nonlinear losses:
with
is a diffusion coefficient. In the case of deterministic parameters,
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
will be under multiple sources of randomness.
In the case of multiple sources of randomness, the process
should locally and globally satisfy the finite variance criteria [
18], i.e.,
where
is the vector of all random parameters in Model (4).
The SDE model (4), with
> 0, has a unique positive solution
:
which reduces to the deterministic solution in the case where the diffusion parameter
. The stochastic process
given in (4) and for
is recurrent and converges to the stationary probability gamma distribution
. For
, the process decays/converges almost sure to zero [
7]. In the case in which
, 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
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
, 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
Equation (6) to get
where
and
which are usually small compared with the square of the mean kernel
. In this section, the subscripts
and
are used as indices for the time variables
and
, respectively.
For the first kernel
Equation (7), we get
where
. In the current work, the integrals
and
are computed using the midpoint integration rule.
For the second kernel
Equation (8), we get
The time step
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
after neglecting
I1 and
I2 compared with
to get
The condition 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 and .
The mean and variance will be computed as
where
The integral represents the Gaussian part of the model variance and is the nonGaussian part of the total variance.
In the current work, the numerical simulations will consider the case of
which should converge, as descried above in
Section 4, to the stationary gamma distribution. Other values of the parameter
are analyzed similarly as will be shown below.
Figure 2 shows the time variation of the mean kernel
and the variance components: total variance, Var_1 =
and Var_2 =
. The convergence of the kernels to steady-states is clear in the figure for the given parameters (
= 0.05,
= 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 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
while the steady-state nonGaussian variance component is proportional to
.
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
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
are decomposed as
with
as the mean and
as the variance of
.
Substitute in
Equation (6) to get
Multiply Equation (12) by
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
and
. From the orthogonality property of functionals
, we use
, where
is the Kronecker delta function. Details about computing expectations
can be found in [
20]. Similar expressions for
and
are straight-forward.
In the case where the parameter
depends only one random variable
, we can write
. This yields the following for
:
Similarly, for
:
and for
:
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
can be analyzed/decomposed into three components: variance
due to random parameters, variance
due to noise and variance
due to the mix between noise and the random parameters, i.e.,
where
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].
7. Numerical Example with Combined Randomness
Assume that the parameter
a fluctuates uniformly around the mean with deviation 2% of the mean value
, i.e.,
with
. This is equivalent to a uniform distribution
. Using
= 0.02,
= 0.01 and zero initial conditions for all kernels except for
= 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
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
where the population extinction is expected.
Figure 10 shows the mean and total variance in the case of
,
= 0.02,
= 0.01 and zero initial conditions for all kernels except for
= 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.