Dealing with variability in ecological modelling: An analysis of a random non-autonomous logistic population model

This paper presents a methodology to deal with the randomness associated to ecological modelling. Data variability makes it necessary to analyze the impact of random perturbations on the ﬁtted model parameters. We conduct such analysis for the logistic growth model with a certain sigmoid functional form of the carrying capacity, which was proposed in the literature for the study of parasite growth during infection. We show how the probability distributions of the parameters are set via the maximum entropy principle. Then the random variable transformation method allows for computing the density function of the population.


INTRODUCTION
The logistic equation models the growth of human, plants, animal and bacterial populations 1 .Developed by Pierre-François Verhulst in 1838 2 , it generalizes the exponential growth model proposed by Thomas Robert Malthus in 1798 3 (which is considered as the first law of population dynamics in the field of population ecology 4 ), by taking into account the lack of resources as the population grows.There is thus a limited capacity in the amount of population.
Classically, the carrying capacity has been considered constant.However, some works started to consider it as a function of time, under the principle that a changing environment may result in a significant change in the carrying capacity.From the viewpoint of dynamical systems theory, the logistic equation with time-dependent limiting capacity was studied in 5 .Some applications of this equation, with different functional forms for the carrying capacity, can be found in a variety of works: in 6 , for the growth of Earth's human population; in 7 , for parasite growth during infection; in 8 , for the population histories of England and Japan; and in 9 , for the total microbial biomass under occlusion of healthy human skin.In 10 , the authors conducted a mathematical analysis of the model from 9 by using an algebraic method.
These types of models, an in general any mathematical model, have three sources of error: model, data and numerical methods errors.The numerical methods error depends on the algorithms and can always be decreased, at least theoretically; moreover, when there is an exact solution, this error does not arise.On the contrary, the model error accounting for the incomplete knowledge on the phenomenon under study, and the data error due to incorrect measurements, lack or scarcity of information, etc. cannot be avoided.Thereby, the process of modelling has associated uncertainty.Accounting for data variability, the parameters of the model should be regarded as random variables with probability distributions, and the solution becomes a random variable that evolves with time, that is, a stochastic process 11 .
The maximum entropy principle infers consistent probability laws for the parameters, by maximizing the ignorance on their density functions while not violating physical principles.The ignorance on the density function is usually expressed via the Shannon entropy functional, subject to certain constraints on the statistical moments and the support 12,13,14 .
Once the model coefficients follow specific probability distributions, the probability law of the stochastic solution must be found.When a closed-form solution is available, the random variable transformation method gives the exact probability density function under mild conditions.This technique has been employed in different settings: Bertalanffy growth model 15 , SIS-type epidemiological model 16 , autonomous logistic equation 17,18 , and radiative transfer equation 19 , for instance.When there is no closed form of the solution, there are hybrid strategies that combine the random variable transformation technique and approximation methods, such as: Monte Carlo simulation, recently applied for the advection linear partial differential equation 20 ; finite difference numerical schemes, applied to some linear random differential equations 21 and to the heat partial differential equation 22 ; and spectral expansions, such as those of Karhunen-Loève type for the damped pendulum differential equation 23 and the logistic differential equation 24 , and those of polynomial chaos type 25 .These approaches capture discontinuity and non-differentiability points of the target density function correctly and are more efficient than non-parametric kernel density estimations.On the other hand, let us mention that Liouville's equation may be used to validate the densities obtained via the random variable transformation method.This equation takes the form of a partial differential equation satisfied by the density function, which, in certain cases, can be explicitly solved by employing the associated Lagrange system 26,27 .Finally, other strategies that only focus on statistics estimations, which rely on simulations and expansions, are available, see 28,29,30 for instance.
In this paper, we deal with the stochastic version of the logistic model proposed in 7 , where the carrying capacity varies with time.The structure of the paper is the following.In Section 2, we study the deterministic logistic equation (no randomness) with time-dependent carrying capacity: we solve the model by using the theory on Bernoulli ordinary differential equations and we compare the solution with the cases discussed in 7,9,10 .In Section 3, we randomize the logistic equation, we prove that the randomized solution satisfies the equation in the path-wise and the mean-square senses, and we show how to infer consistent probability distributions for the model parameters via the maximum entropy principle.In Section 4, we derive the probability density function of the solution proposed in 7 by means of a comprehensive application of the random variable transformation method.In Section 5, we carry out numerical test examples.Finally, in Section 6, we discuss the main results and present the conclusions.

NON-AUTONOMOUS LOGISTIC EQUATION
The autonomous logistic equation has the form where 0 > 0 is the initial condition, > 0 is the growth rate parameter and > 0 is the carrying capacity, with 0 < .The solution to (1) is well-known: Notice that, as expected, = lim →∞ ( ).As explained in the previous section, we aim at analyzing the non-autonomous logistic equation driven by a time-varying carrying capacity ( ): Different forms of ( ) have been conceived in the existing literature.For instance, 7 employed the sigmoid growth while 9 considered the logistic growth where 0 = (0) is the initial limiting capacity, = lim →∞ ( ) is the saturation (or equilibrium) level, and is the saturation constant.It is assumed 0 < 0 < and > 0. Other functional forms of ( ), which vary sinusoidally, exponentially or linearly, can be consulted in 8,10 .
Equation ( 3) is a Bernoulli ordinary differential equation.It can be solved via the change of variables = 1∕ .By differentiating ( ), a linear ordinary differential equation is derived for : When ( ) is given by ( 5) and is substituted into (6), we have the solution derived in 10 : To evaluate the integral from the denominator, the authors employed an algebraic trick, by expanding part of the integrand as a series: In practice, the series is truncated to a finite-term sum.Accurate approximations to the exact solution ( ) are obtained for small orders of truncation.Notice that this method is applicable whenever ( ) can be written as ̃ (e − ), where ̃ is a decreasing function expressed as a power series On the other hand, when ( ) is defined by (4), we have a closed-form solution (without integrals) after some simple manipulations in (6): In the development from Section 4, we will focus on the closed-form solution (7).

Randomization
Model (3) is randomized.Mathematically, there is an underlying complete probability space (Ω,  , ℙ), where Ω is the sample space,  ⊆ 2 Ω is the -algebra of events, and ℙ ∶  → [0, 1] is the probability measure.The parameters and 0 are random variables and the carrying capacity ( ) is a stochastic process.In this manner, ( ), 0 ( ) and ( , ) depend on the outcome ∈ Ω of an experiment.Notice that ( ) is a stochastic process by taking , 0 and as random variables, which evaluate at the outcomes of experiments as ( ), 0 ( ) and ( ).The general solution ( ) to (3), given by ( 6), is a stochastic process ( , ).It is assumed that, given a population datum at a certain time , it is just a possible realization of the model solution over all possible realizations { ( , ) ∶ ∈ Ω}.The outcome will be implicitly assumed without being systematically written.
Let L (Ω), 1 ≤ < ∞, be the Lebesgue space of random variables where is the expectation operator.When = ∞, L ∞ (Ω) is defined as the Lebesgue space of almost surely bounded random variables, and ‖ ⋅ ‖ ∞ stands for the least upper-bound.The set L (Ω), 1 ≤ ≤ ∞, is a Banach space, and for = 2, it is a Hilbert space with the inner product ( , Convergence in L 2 (Ω) is referred to as mean-square convergence.Convergence in L (Ω) preserves the convergence of the statistical moments up to order ; in particular, mean-square convergence implies convergence of the expectation and the variance.Mean-square calculus considers the limits from the definitions of continuity, differentiability, Riemann integrability, etc. in the mean-square sense.Under this alternative calculus, one may solve differential equations with random parameters 11,Ch. 4 , 31,32 .
We point out that, as a consequence of 35, pp. 440-441,properties 1-2 , the integral from the denominator of ( 6), ∫ 0 e ( ) d , can be considered in the mean-square sense (the Riemann sums have mean-square convergence), apart from path-wise.

Inverse parameter estimation
In practice, one must set consistent probability distributions for the parameters.Let us see how the Shannon entropy measure is useful here 12,13 .Let be any random parameter.Its Shannon entropy is expressed via the functional where is the probability density function of , [ , ] denotes its support (of course, it could be = −∞ and/or = ∞), and log is the natural logarithm (in base e).It is understood here that 0 log 0 = 0. Prior information on is expressed via the following restrictions: Here denotes the -th statistical moment of .For instance, 0 = 1 by definition of density function, 1 is the mean of , 2 is the variance of plus its squared mean, etc.The maximum entropy principle maximizes the objective functional (8) subject to (9).This is usually done in the literature via calculus of variations 36 .The Lagrange multipliers method says that are the functional derivatives with respect to at , determined by means of the Euler-Lagrange equation, and ⟨⋅⟩ denotes the linear span.This gives the solution for certain Lagrange constants 0 , 1 , … , ∈ ℝ.These constants are determined by solving the nonlinear system of equations that appears when substituting (10) into (9).
It is worth pointing out that a derivation of (10) using just elementary calculus was proposed in 37 .The idea is simple.Let be given by ( 10) and let ̃ be any other non-negative function satisfying the corresponding restrictions (9).From the inequality log ≤ − 1 for > 0, one derives and we are done.
Suppose that no information on is known, only a bounded support [ , ].Then = 0 and (10) becomes the density function − .On the other hand, suppose that > 0 and its mean ̂ is known.In this case, = 1, and follows the Exponential law, where the rate parameter is 1∕ ̂ : well-known case is when is defined on the real line, with a certain mean value and standard deviation ( = 2), which leads to the Gaussian distribution.Other cases were investigated in 12 .In general, one obtains 0 , … , numerically.In practice, determining the prior information of a parameter is not an easy task.Suppose that we have a single pointwise deterministic estimator of , say 1 .This value could be used as the mean of (Exponential distribution for ), or to define a certain support for by adding and subtracting a certain variability around 1 (Uniform distribution for ).But, in any case, the information is vague.A better scenario arises when we have several pointwise deterministic estimators of , say 1 , … , , ≥ 2. By employing this sample, the mean and the variance of can be estimated.Moreover, the support of can also be estimated via Markov's inequality (which states that ℙ where [ ] is the standard deviation of ).In this case, = 2, and the Lagrange constants 0 , 1 and 2 satisfying (9) are calculated numerically.

PROBABILISTIC SOLUTION: DENSITY FUNCTION
In this section, we assume that the random parameters of the model (3) with carrying capacity (4) have specific probability distributions, and we aim at computing the probability density function of the solution ( ) given by (7).We employ the random variable transformation method, see 38, Th. 2.1.5 .

Random variable transformation method
The random variable transformation technique gives the density function of a response under the relation output-input = ( ), where is a deterministic map called the transformation mapping and is a random quantity.It is assumed that the dimensions of and are equal.When is smooth on an open set containing ( ) = { ( ) ∶ ∈ Ω}, one-to-one and with non-vanishing Jacobian ( ) = det( ( )), the density function of is given by where ℎ = −1 is the inverse of on its domain.Indeed, for any Borel set contained in the image of , we have The requirements may be overcome when the domain of the transformation mapping is divided into sub-domains where the conditions hold 38, Th. 2.1.8 .

Application
Notice that the solution ( 7) is a closed-form transformation of the random input parameters 0 , , , 0 and .Then, using (11), the density function of ( ), ( ) ( ), can be expressed in terms of the joint density function of ( 0 , , , 0 , ), ( 0 , , , 0 , ) ( 0 , , , 0 , ).The transformation mapping is This manner of proceeding is usual: one of the components of the transformation mapping is ( ), while the rest of components are fixed; the parameter corresponding to the component ( ) (in this case 0 ) must be easily isolated in the equation ( ) = .
The inverse mapping is given by where the positivity comes from the fact that > 0 and ( − )(e ( − ) − 1) > 0 for all values of and (except when = , but this happens with zero probability when and are independent and absolutely continuous random variables).By (11) and marginalizing, where ( 0 , , , ) = {( 0 ( ), ( ), ( ), ( )) ∶ ∈ Ω} is the image of ( 0 , , , ).To apply the random variable transformation method, we must have the denominator of ℎ( 0 ( ), ( ), ( ), ( , ), ( )) distinct from 0 almost surely; this always holds because the solution ( ) is absolutely continuous.Notice that, given an equation of the form ( ) = , the parameters and 0 are also easily isolated.Only and cannot be explicitly isolated, as they appear within and outside an exponential function.Hence (11) is applicable by isolating and 0 , instead of 0 .First, for , the transformation mapping is Therefore, by (11) and marginalizing, Now we justify that the term e − 1 − − (e ( − ) − 1) coming from the numerator of the Jacobian, does not need to be written with absolute value as indicated in the general formula (11).For this goal, let > 0 be fixed, and and be positive random variables.Observe that the above-mentioned term can be expressed as Moreover, the function ( ) = (e − 1)∕ , defined for all ∈ ℝ − {0}, satisfies that lim →0 ( ) = 1, ( ) > 0 and ′ ( ) > 0 (i.e. is increasing).Now we distinguish two cases with respect to the two positive random variables and : > > 0 and > > 0. Notice that other situations have null probability to occur, by absolute continuity.In both cases, we shall show that expression ( 14) is positive.Assume that > > 0. Then > − and > ( − ) for > 0. Since is increasing, ( ) > (( − ) ), i.e. e − 1 − e ( − ) − 1 ( − ) > 0.
Multiplying this last expression by > 0, from ( 14) one gets the stated result.Now, let us assume that 0 < < .Then, > − and > ( − ) for > 0. Now the reasoning follows exactly as before, and we are done.For 0 , the transformation mapping is defined as By (11) and marginalizing, Two general considerations must be made here.Firstly, when some input random parameter is independent of the rest, then the joint density function can be factored as a product.In fact, when applying the maximum entropy principle, the input random parameters that are distinct are assumed to be independent.Secondly, when some input parameter is a constant 0 , then its probability density function is a Dirac delta function centred at the constant, = 0 = ∞1 { 0 } .This generalized function is defined via the heuristic property ∫ ∞ −∞ ( ) 0 ( ) d = ( 0 ), for any real function .In practice, one chooses between (12), ( 13) and ( 15) depending on whether 0 , and 0 have proper density function, respectively.
Given the density function ( ) , any statistic of ( ) can be determined by using the relation , where is any deterministic function.Also, given a probability ∈ (0, 1), a confidence interval [ , ] is constructed by using ≈ ∫ ( ) ( ) d .An analogous procedure serves to calculate the density function of the carrying capacity ( ) given by (4).Consider the transformation mapping , .By (11) and marginalizing, we obtain

Its inverse mapping, computed by isolating
We can also consider the transformation mapping This selection corresponds to isolating when computing the inverse of from the equation ( ) = .The inverse mapping is then 0 − e − 2 > 0. By (11) and marginalizing, we arrive at Finally, in contrast to the equation ( ) = , we can isolate in the equation ( ) = ; consider Then Let us summarize the formulae derived in this section.Depending on which parameter we can isolate in the equations ( ) = and ( ) = and on whether the density function exists, one chooses the appropriate expression of ( ) ( ) between ( 12) (when = 0 ), ( 13) (when = ) and ( 15) (when = 0 ), and the appropriate expression of ( ) ( ) between ( 16) (when = 0 ), ( 17) (when = ) and ( 18) (when = ).Notice that 0 does not play any role in formula (18), since it refers to the carrying capacity.On the other hand, for ( ) ( ), the situations = and = are not contemplated because they cannot be isolated in the equation ( ) = (they appear within and outside an exponential).
All of the integrals that have appeared in this section can be computed via tensorized Gauss quadratures.If the integrand has significant jump discontinuities, a parametric Monte Carlo method may be employed to estimate the integral, instead 39 .

NUMERICAL EXAMPLES
In this section we undertake numerical test examples for specific input parameters.In Example 1, we will consider a single pointwise deterministic estimator for each parameter.While in Example 2, we will assume several pointwise deterministic estimators for each parameter.
Example 1.Let us consider the parameters values = 0.15, = 0.1, 0 = 0.3, = 1 and 0 = 0.2.The coefficients and are measured in time −1 , while 0 , 0 and represent proportions of counts, being the maximum.The parameters values may be fitted experimentally or via least-squares minimization procedures.By using the exact deterministic solution (7), we can forecast the population ( ) for ≥ 0. The results are reported in Figure 1, where the fat line refers to ( ) defined by ( 4), the solid thin line represents ( ) defined by (7), and the dots correspond to the numerical solution employing the classical Runge Kutta scheme for validation.We observe that, as expected, expression (7) and the numerical solution agree.Also, the solution ( ) tends to the carrying capacity ( ) as increases, at the time when the growth of both functions is decelerated.The slope of the sigmoid trajectory of ( ) depends on the saturation constant .On the other hand, the slope of the sigmoid trajectory of ( ) depends on both the saturation constant and the growth rate .4), the solid thin line represents ( ) defined by (7), and the dots correspond to the numerical solution employing the classical Runge Kutta scheme.This figure corresponds to Example 1.
Next we study the effect of randomness on the model output.Let us suppose that, by some experimental evidence, 0 lies within [0.19, 0.21], is bounded on [0.13, 0.17], > 0 with mean equal to 0.1, and 0 is bounded on [0.26, 0.34].According to the maximum entropy principle, 0 has a Uniform distribution on [0.19, 0.21], follows a Uniform distribution on [0.13, 0.17], has an Exponential distribution with rate parameter 1∕0.1 = 10, and 0 is Uniform on [0.26, 0.34].These distributions maximize the ignorance on the random behaviour of 0 , , and 0 , while not violating the restrictions on their supports and statistics, if any.The random variables are assumed to be independent.We keep with its constant value 1, as it represents the maximum proportion; its density function may be considered in terms of the Dirac delta function as 1 .The saturation level being 1 may be seen as a general fact, as it can usually be scaled out of the problem for simplicity.The joint density function ( 0 , , , 0 , ) is factorized in terms of the marginal density functions:

As
does not have a proper density function, we must consider (12) or (15) for the density ( ) ( ).We take (12).On the other hand, for ( ) ( ), we must consider (16) or (18).We pick (16).Thus, for both density functions, the inverse of the transformation mapping is defined by isolating 0 .The density function ( 12) is a triple integral in terms of d 0 d d , while ( 16) is a onedimensional integral with respect to d .In Figure 2, we plot the density functions ( 12) and ( 16) for different times = 10 (solid line), 20 (dashed line), 30 (dotted dashed line) and 40 (dotted line).For (16), the plot has been truncated to [0, 8] in the vertical axis for = 20, = 30 and = 40, otherwise the density functions take too large values near = 1.
The support of ( ) is defined from the value of (7) when = 0.13, = 0, 0 = 0.26, 0 = 0.19 (these values correspond to the left endpoints of the domains of their corresponding distributions) and = 1, which is 0.236 for = 10, 0.253 for = 20, 0.258 for = 30, and 0.259 for = 40.Regarding the support of ( ), its left endpoint is 0.26 for any ≥ 0. The upper bound for the supports is 1.Notice that small dispersions of the parameters give wide supports for ( ) and ( ): this means that, for small perturbations in the experimental values of the parameters, the response may be sensitive and present large variation.Hence the need of appropriate tools to understand the random behaviour of the response.For standard probability distributions for the
The evolution of the density function of ( ) can be better understood by means of Figure 3.It shows the contour plot of the densities for locations in [0, 1] and non-negative discretized times.Light colours correspond to large values, being the "outliers" represented in white colour.From 0 at = 0, the densities spread as increases and, for large , they tend to 1 .In Table 1, we tabulate the deterministic value of ( ) together with its main statistical information (mean value, standard deviation and 95% confidence range), for times = 10, 20, 30 and 40.The parameter is set to 1, as it stands for the maximum proportion in any case.Compared to Example 1, now there is much more information about the parameters, and the conclusions will be robust.The mean and the variance of each parameter are estimated via the corresponding sample of length ten.The support of each parameter is defined by adding and subtracting three times the standard deviation to the mean (this is motivated by the 3 rule and the Vysochanskii-Petunin inequality, see 40 ).When applying the maximum entropy principle, we have = 2, and the Lagrange constants 0 , 1 and 2 satisfying the complex nonlinear system (9) are determined numerically by employing Newton's method, with starting value (−1, 1, 1), for each parameter.In Table 2 2 For each parameter, we tabulate the mean and the standard deviation estimated from the corresponding sample of length ten.The support of each parameter is defined by adding and subtracting three times the standard deviation to the mean.The Lagrange constants 0 , 1 and 2 are determined numerically.This table corresponds to Example 2.
Figure 4 plots the resulting density functions (10) for 0 (solid line), (dashed line), (dotted dashed line) and 0 (dotted line).As explained in 12,41 based on theory and numerical evidence, whenever the mean, the variance and a compact support are available, the entropy distribution is approximately truncated Gaussian.The density function of is the Dirac delta function centred at 1, 1 .All of these random variables are assumed to be independent, so their joint density function factorizes.As does not have a proper density function, we must consider (12) or (15) for the density ( ) ( ), and ( 16) or ( 18) for ( ) ( ).We take ( 12) and ( 16), which come from isolating 0 when defining the inverse of the transformation mapping.In Figure 5, we show the density functions ( 12) and ( 16) for different times = 10 (solid line), 20 (dashed line), 30 (dotted dashed line) and 40 (dotted line).On the other hand, in Figure 6 we can see the contour plot of the densities ( ) ( ) for positions ∈ [0, 1] and non-negative discretized times .Light tones are associated to large values, and the "outliers" are painted in white colour.
In Table 3, we document the deterministic value of ( ) together with its main statistical information (mean value, standard deviation and 95% confidence range), for times = 10, 20, 30 and 40.

DISCUSSION AND CONCLUSION
This paper addresses two important issues in ecological modelling: the mathematical analysis of the logistic growth model for time-varying carrying capacity and the incorporation of randomness into the model formulation.The former is motivated by the principle that a changing environment may result in a significant change in the limiting capacity, while the latter is justified accounting for data errors.The study of randomness on the non-autonomous logistic growth model does not seem to have been explored, at least to our knowledge.
With time-dependent carrying capacity, the deterministic logistic equation belongs to the class of Bernoulli ordinary differential equations, which can be solved by a change of variables.The constant parameters (growth rate, saturation rate, saturation level, etc.) can be fitted experimentally or via least-squares minimization procedures.The resulting model can then be used for forecast.In this paper we focus on the model from 7 , which described parasite growth during infection.
To analyze the effect of randomness on the model output, first one must set appropriate probability distributions for the input parameters.This can be done via the maximum entropy principle.Given a parameter, its density function is chosen by maximizing the Shannon entropy measure, restricted to certain statistical moments and support usually devised from deterministic fittings and experimental data.For example, if only information on the support of the parameter is available, then the Uniform distribution is selected; if the parameter is positive and only its mean is known, then the Exponential distribution is chosen.In general, the selection of the density function is made numerically.The foundation of these selections is the calculus of variations and Lagrange multipliers.Depending on the amount of data at each time instant of interest, the prior information on the parameters is either more vague or more complete.When there are replicated data at each time instant, the information on the parameters is more significant and the entropy distribution is more reliable; this "multi-data" case is the most desirable scenario in practice.The density function of the exact model solution can be determined via the random variable transformation method.This technique is based on the application of a formula when the transformation mapping between the random quantities is injective (one-to-one) with non-vanishing Jacobian.These requirements may be overcome when the domain of the transformation mapping is divided into sub-domains where these conditions hold.
In this manner, we have a probabilistic solution to the model.Apart from pointwise estimations, other relevant information such as confidence intervals, dispersion and statistical moments can be derived.This supplies a more faithful and complete description of the ecological process.
In the near future, we aim at conducting a similar analysis for the logistic growth model proposed in 10 , where the carrying capacity is itself described via a logistic curve.In that paper, the authors modelled the total bacterial biomass during occlusion of healthy human skin.The solution to that model is expressed algebraically by truncating an infinite series.Thus, when the input parameters are random variables, the application of the random variable transformation method shall be analyzed.

FIGURE 1
FIGURE 1Given the parameters values = 0.15, = 0.1, 0 = 0.3, = 1 and 0 = 0.2, the fat line refers to ( ) defined by (4), the solid thin line represents ( ) defined by(7), and the dots correspond to the numerical solution employing the classical Runge Kutta scheme.This figure corresponds to Example 1.

Example 2 .
Let us suppose several listed values of the parameters 0 (in proportion of counts), (in time −1 ), (in time −1 ) and 0 (in proportion of counts), obtained experimentally or via least-squares minimization procedures from laboratory data that deterministic mean standard deviation 95% confidence interval = 10 0.

FIGURE 4
FIGURE 4 Density functions of 0 (solid line), (dashed line), (dotted dashed line) and 0 (dotted line) determined via the maximum entropy principle.This figure corresponds to Example 2.
we report the results, where the values of 0 , 1 and 2 have been stopped at the sixth significant digit.
This figure corresponds to Example 2.