Mean-square stability of the zero equilibrium of the nonlinear delay differential equation: Nicholson’s blowflies application

We are concerned about the stochastic nonlinear delay differential equation. The stochasticity arises from the white Gaussian noise, which is the time derivative of the standard Brownian motion. The main objective of this paper is to introduce a new technique using the Lyapunov functional for the study of stability of the zero solution of the stochastic delay differential system. Constructing a new appropriate deterministic system in the neighborhood of the origin is an effective way to investigate the necessary and sufficient conditions of stability in the sense of the mean square. Nicholson’s blowflies equation is one of the major problems in ecology; necessary conditions for the possible extinction of the Nicholson’s blowflies population are investigated. We support our theoretical results by providing areas of stability and some numerical simulations of the solution of the system using the Euler–Maruyama scheme, which is mean square stable Maruyama (Rendiconti del Circolo Matematico di Palermo 4(1):48, 1955), Cao et al. (Appl Math Comput 159(1):127–135, 2004).


Introduction
Functional differential equations (FDEs) have motivated many mathematical and applied statistical researches. For instance, one of the major problems is the dynamics of animal population that has been proposed by [3][4][5][6]. Time delay occurs so often; many processes involve time delays in physics, medicine, finance, ecology, chemistry, etc. Realistic models must include some of the past history of the state of the system, and this, in turn, leads to the delay differential equations (DDEs). For more details regarding processes with time delay and their qualitative behavior, see [3][4][5][7][8][9][10][11][12].
In these epidemiological and ecological mathematical models, uncertainty in the interactions between individuals of population and/or higher-frequency environmental noises develops stochasticity. Considering the stochasticity in these models reveals the mechanisms that influence the transmission and control the disease in epidemiology, dynamics of the population in ecology, etc. [13]. Stochastic models are proposed to capture the uncertainty and variations in the mathematical models by perturbing the deterministic system with a white Gaussian noise, which is ill-defined, and then change it to a Brownian motion, which is well defined by ζ(t)dt = dB(t).
Without solving these systems, we study the stability of the equilibria, which is a very effective way to have a good insight of the solutions and their properties. Constructing a suitable Lyapunov functional is a good approach for investigating the necessary and sufficient conditions of stability of the equilibrium points [10,[14][15][16].
Dedicated to the study of stability of stochastic delay differential equations (SDDEs), there are many works, including, but not limited to, use this class of equations to control the stabilization problem of the controlled inverted pendulum, model the infectious diseases and model the population dynamics of the Australian sheep-blowfly known as Nicholson's blowfly [17][18][19]. By stochastic neutral differential equations (SNDEs), the stochastic neural networks and stochastic cellular neural networks have been used to model many of the human activities in science [20][21][22].
Our work is involved in the study of stability of delayed dynamical systems. In comparison with the known methods of investigating stability in the mean square of these systems [10,18,19,23], our approach is more effective as it lays in constructing an appropriate deterministic system using appropriate Lyapunov functionals. It can be implemented in many application problems, see [4,20,22]. For instance, regarding the Nicholson's blowflies model, many works in the literature have studied the extinction of these blowflies in the deterministic and stochastic sense [18,[24][25][26]. Our concept provides more stability delay-dependent and delay-independent conditions with better stability regions compared to the work done on this equation. Moreover, by this approach, this work can be extended by studying the persistence of these species of blowflies.
We shall determine the mean-square stability conditions of the stochastic nonlinear delay differential equation and its general form using a new way of constructing a delayed-deterministic system by Lyapunov functional in the presence of the white noise term. Consider the stochastic delay differential equation The solution of this equation is where the last term is known as an Itô integral. Existence and uniqueness of this solution are considered in [27,28]. The delay in the stochastic process is τ > 0, assume C := C ([−τ, 0], R) is a Banach space of continuous functions defined on [−τ, 0] with the norm φ = sup −τ <s≤0 |φ(s)|. In population dynamics, let the initial function φ ∈ C + , where C + := C ([−τ, 0], R + ). The continuous functionals μ(t, X t ) and σ (t, X t ) are defined on [0, ∞) × C[−τ, 0] and satisfy Lipschitz condition in the second argument, i.e., for ∂μ(t, X t ) ∂ X < L and also for σ (t, X t ). The functional μ is assumed to be a slow deterministic process, and σ is assumed to be a stochastic part of a fast phenomenon. We can perform the transformation used by ( [28], p. 200), which motivates the study of the trivial equilibrium of (1), i.e., μ(t, 0) = σ (t, 0) = 0. B(t) is a one-dimensional standard Wiener process defined on the complete filtered probability space ( , F, F B t , P), where F B t is the filtration generated by it up to time t. The initial function φ is a stochastic process which is independent of the minimal σ -algebra generated by random variables B(t) − B(s) for 0 ≤ s < t < ∞. Moreover, we also consider the general stochastic delay differential equation where the functionals μ(t, X t ) and σ (t, X t ) are mdimensional continuously differentiable. B i (t), i = 1, · · · , m are m-dimensional standard Wiener processes.

Nicholson's blowflies model
One of the major problems in ecology is the population dynamics of the Australian sheep blowfly, Lucilia cuprina, which is known as Nicholson's blowfly. Author in [29] introduced the differential equation that models the population dynamics of this blowfly with fixed delay in the following form X (t) represents the population of the mature adults at time t. All parameters p, a, δ ∈ [0, ∞), where p represents the maximum daily production rate of eggs per capita, a −1 is the size at which the population reproduces at the maximum rate, δ is the adult death rate of per capita daily and τ > 0 is the delay in the production process.
The basic reproduction ratio is denoted by This model admits only the trivial solution, X 0 = 0, if the basic reproduction ratio of the system R 0 < 1, i.e., the maximum reproduction rate per capita, is less than the death rate, and for R 0 > 1, a unique positive equi- , exists and at this point, this kind of blowflies does not become extinct. The species become extinct due to the global stability of the zero solution of (4), see [30]. Many authors have studied the dynamics of this equation in the deterministic sense [23,26,[30][31][32], and stochastically [18,19,24,25,33].
This system is exposed to stochastic perturbation in the form of white noise, which is assumed to be proportional to the deviation of the current state of the system from the zero solution; therefore, the stochastic version of (4) is in the form The parameter λ represents the intensity of the noise. Necessary and sufficient conditions of mean-square stability of the trivial equilibrium of (5) are studied using a new approach, which provides better areas of stability with some numerical simulations that strengthen the theoretical findings. The plan of the paper is as follows: In Sect. 2, some important preliminaries and notations are introduced. In Sect. 3, we introduce our new concept by proving the mean-square asymptotic stability of the zero solution of (1), (3) accordingly (5). A detailed numerical example is introduced in Sect. 4, with stability regions and some numerical simulations and interpretation. Conclusions and future works are presented to close the paper in Sect. 5.
Definition 2 [38] The standard one-dimensional Brownian motion B(t) is a real-valued continuous {F t }adapted process defined on the complete probability space ( , F, F B t , P) with the properties 1.
, it is a simple or step process if there exists the partition a = t 0 < t 1 < · · · < t k = b Definition 3 [39,40] The zero solution of (1) is 2. Asymptotically mean square stable if it is mean square stable and lim t→∞ E|X (t, t 0 , φ)| 2 = 0.

Main result
Our main finding is to investigate the stability of the zero solution of the general Eqs. (1), (3) in the sense of the mean square using a suitable quadratic form of Lyapunov functional. Next theorems provide the meansquare stability under a constructed deterministic function in the neighborhood of the origin. This way is more effective and provides different stability conditions. Based on this result, the extinction of the Nicholson's blowflies species will be studied.

Theorem 1
The trivial solution of (1) is mean square asymptotically stable if there exists a negative definite functional θ(t, X t ) in the neighborhood of the origin such that Proof Define the Lyapunov functional in the quadratic form where Q is a positive-definite symmetric matrix, T is the transposition and define Q d as a diagonal matrix that has the same element of the diagonal of Q.
hold, then taking the expectation implies , and it is a negative-definite deterministic functional defined on the neighborhood of the origin; therefore, for M > 0 Therefore, the zero solution of (1) is asymptotically mean square stable.

Theorem 2 The trivial solution of (3) is mean-square asymptotically stable if there exists a negative-definite functional θ(t, X t ) in the neighborhood of the origin such that
Proof Consider the linear part of (3) as a special case by letting Then, the linear system becomes where A, C are constant m × m matrices. Now, assume the Lyapunov functional (7) dV(t, X t , − X T (t) Q X(t).

E [dV(t, X t , X (t))] = X T (t) Q A X(t) dt
Following the proof of previous theorem, assume Then, Therefore, the zero solution of (3) is asymptotically mean square stable.

Proposition 1 The trivial solution of (5) is mean square asymptotically stable if
The proof follows directly from Theorem 1, by considering the same quadratic Lyapunov functional (7) and using properties (8), and leads to the functional which should be negative definite in the neighborhood of the origin for asymptotic mean-square stability. Now, centering system (5) on the zero equilibrium solution X 0 = 0, linearizing using the notation e −a X (t−τ ) = 1−a X (t −τ )+O(X (t −τ )), and neglect

Proposition 2
The trivial solution of (11) is meansquare asymptotically stable if The proof follows from Theorem 2, by considering the same-quadratic Lyapunov functional (7) and using properties (8), and leads to the functional which should be negative definite in the neighborhood of the origin for asymptotic mean-square stability for the linear system (11) and consequently for (5).

Example
Consider the linear system (11). By introducing the quadratic Lyapunov functional (7) for Q = 1 and Then, from Proposition 1, and system (11) For the negative definiteness, choose Then According to Theorem 1, for Q = 1 which is negative definite if 2 p + λ 2 < 2δ. (12) Condition (12) is a necessary condition for the meansquare stability of the zero solution of (5); it is a delayindependent condition. Accordingly, Fig. 1 shows the regions of mean-square stability. By introducing different Lyapunov functionals, we get different stability conditions. Assume another Lyapunov functional in the form Then, Therefore, (a) (b) (c) Fig. 1 Regions of mean-square stability of (5) according to condition (12) for different values of the parameter of the intensity λ (a) (b) (c) Fig. 2 Regions of mean-square stability of (5) according to condition (14) for different values of τ This delay-dependent condition is necessary for the mean-square stability. Via this new condition, we can show the impact of delay on the stability regions with λ. Figure 2 shows the mean-square stability regions for different values of τ with fixed λ = 1.8. For fixed delay τ = 0.5, mean-square stability regions are shown in Fig. 3 for different λ.
Regarding the numerical simulation, we use the Euler-Maruyama algorithm which has been shown by [1]. According to (2), Then, the scheme of EM has the form where x(t n ) is the data required in the scheme and x(t n+1 ) is the resulting process at t n+1 . The EM scheme is strongly convergent with order 0.5. We choose the stepsize t appropriately to avoid time and errors of computation. For the initial history function φ = 2.5 cos(3s), s ∈ [−3.5, 0] and λ = 1.8, τ = 3.5, Fig. 4a shows stable trajectories of the solution of (5), and conditions (12), (14) are satisfied. If these conditions are not met, then we get unstable trajectories of the zero solution as shown in Fig. 4b. The impact of the intensity of the noise λ and the delay in the production process τ on the behavior of the solution should be noted. Figure 5 shows three different behaviors of the solution according to condition (14). The impact of τ is shown in Fig. 6. For small values of λ, τ in the light of condition (14), we get asymptotic stable solutions.

Interpretation
By introducing a suitable Lyapunov functional in the previous example, we have obtained a delayindependent condition (12) that gives areas of meansquare stability of the zero solution in ( p, δ) space of parameters for different values of λ. By choosing small values of the intensity of the noise, we get better areas of stability as shown in Fig. 1. By introducing a different Lyapunov functional, we have obtained a delaydependent condition (14), which is better as it considers the variation in the delay of the production process of eggs. We investigate the impact of λ, and we arrive at the same result of better regions by decreasing λ; this is shown in Fig. 3. Moreover, choosing small values of delay τ helps in the extinction of the Nicholson's blowflies, and better regions of stability are obtained for small τ . Some computer simulations are carried out to support the result. Choosing the point (0.1, 3.0) which is in the stability area, performing the numerical simulation of the solution with X = 0, we get 25 blue stable trajectories in Fig. 4a. Figure 4b shows 25 red unstable trajectories of the solution by choosing the coordinate (3.0, 1.0) which is out of stability region, and accordingly condition (14) is not satisfied. Figure 5 shows the impact of the environmental noise on the behavior

Fig. 6
Numerical simulation of the solution of (5) for p = 0.9, δ = 4.5, λ = 0.1, asymptotic stable solution (blue), bounded, stable but not asymptotically stable (black) and unstable solution (red) of the solution, by increasing λ, the solution becomes unstable (X (t) → ∞) even if R 0 < 1 (red trajectory), and besides, small scale of the environmental noise can lead to the possible extinction as shown by the blue trajectory, which is asymptotically stable. In the same figure, we have a stable but not asymptotically stable solution (periodic) for a specific value of λ (black trajectory). In Fig. 6, we have got the same result for the change in τ , and it is advisable to choose small enough values of the delay τ for stable solutions.

Conclusion and further directions
Our work has led us to conclude the necessary and sufficient conditions for stability of the zero solution in mean square under the influence of white noise term. This new approach can provide different stability (delay-dependent) conditions. This work can be extended to many applications in many disciplines; for instance, neural networks, infectious diseases and the generalized system of animal population such as the general equation of the Nicholson's blowflies model, which is known as the neoclassical growth model Of course, our concept of stability has some limitations, 1. It is not yet known whether this approach can be applied to systems involving the general fractional Brownian motion B H (t), 0 < H < 1, like the applications studied by [41,42], for instance. 2. A challenging point to us is to study the meansquare stability using our approach in case of systems involving distributed delays.