Improving Stability Conditions for Equilibria of SIR Epidemic Model with Delay under Stochastic Perturbations

: So called SIR epidemic model with distributed delay and stochastic perturbations is considered. It is shown, that the known sufﬁcient conditions of stability in probability of the equilibria of this model, formulated immediately in the terms of the system parameters, can be improved by virtue of the method of Lyapunov functionals construction and the method of Linear Matrix Inequalities (LMIs). It is also shown, that stability can be investigated immediately via numerical simulation of a solution of the considered model.


Introduction
Investigation of different versions of the mathematical model of the spread of infections diseases, so called, SIR epidemic models, has a long history, and until now these models are very actual and are very popular in research (see, for instance, ). A big bibliography of SIR epidemic models investigation one can find also in ( [31], Ch-11). During its history of development, SIR epidemic models were considered both with constant and distributed delay, both in a deterministic and in a stochastic version. We will consider the SIR epidemic model in the form of the following system of differential equations with distributed delay [2,31] S(t) = b − βS(t) ∞ 0 I(t − s)dK(s) − µ 1 S(t), Here, S(t) is the number of members of a population susceptible to the disease at time t, I(t) is the number of infective members at time t and R(t) is the number of members who have been removed from the possibility of infection at time t, through full immunity, b is the recruitment rate of the population, µ 1 , µ 2 and µ 3 are the natural death rates of the susceptible, infective and recovered individuals, respectively, β is the transmission rate, and λ is the natural recovery rate of the infective individuals. The parameters b, β, λ, µ 1 , µ 2 , µ 3 of the system (1) are positive constants, K(s) is a nondecreasing function, such that ∞ 0 dK(s) = 1, the integrals are understanding in the Stieltjes sense. The equilibria of the system (1) are defined by the conditionsṠ(t) ≡ 0,İ(t) ≡ 0,Ṙ(t) ≡ 0, i.e., by the system of algebraic equations b = βSI + µ 1 S, βSI = (µ 2 + λ)I, which has two solutions: E 0 = (bµ −1 1 , 0, 0) and E + = (S * , I * , R * ), where From (4) it follows that E + is a positive equilibrium of the system (1) by the condition Let us suppose that the system (1) is influenced by stochastic perturbations of the white noise type that are directly proportional to the deviation of the system state from the equilibrium (S * , I * , R * ). By that the SIR epidemic model (1) under stochastic perturbations is described by the system of Ito's stochastic differential equations [31,32] where w i (t), i = 1, 2, 3, are mutually independent standard Wiener processes [31,32], Note that the equilibrium (S * , I * , R * ) of the initial deterministic system (1) is the equilibrium of the stochastic system (6) too.
In [31] the following simple sufficient conditions for stability in probability of the equilibria E 0 and E + are obtained that are formulated immediately in the terms of the system (6) parameters.

Lemma 1.
If then the equilibrium E 0 of the system (6) is stable in probability.

Remark 1.
Note that the second condition (8) contradicts with (5). It means that by the conditions (8) the system (6) does not have the positive equilibrium E + .

Lemma 2.
Let be δ 3 < µ 3 and Then the equilibrium E + of the system (6) is stable in probability.
Below new sufficient conditions for stability in probability of the equilibria E 0 and E + are obtained by virtue of the method of Lyapunov functionals construction [31] in the terms of Linear Matrix Inequalities (LMIs) [33] that are essentially less conservative than the conditions (8)-(10).

Some Auxiliary Definitions and Statements
Using the new variables It is clear that stability of the equilibrium (S * , I * , R * ) of the system (6) is equivalent to stability of the zero solution of the system (11). It is known [31] that for getting sufficient conditions for stability in probability of the zero solution of the nonlinear system (11) it is enough to get sufficient conditions for asymptotic mean square stability of the zero solution of the linear part of this system

Definition 2.
The zero solution of the system (12) with the initial condition y(s) = φ(s), s ≤ 0, is called: -mean square stable if for each ε > 0 there exists a δ > 0 such that E|y(t, φ)| 2 < ε, t ≥ 0, provided that φ 2 = sup s≤0 E|φ(s)| 2 < δ; -asymptotically mean square stable if it is mean square stable and for each initial function φ the solution y(t, φ) of the system (12) satisfies the condition lim t→∞ E|y(t, φ)| 2 = 0; -exponentially mean square stable if it is mean square stable and there exists λ > 0 such that for each initial function φ there exists C > 0 (which may depend on φ) such that E|y(t, φ)| 2 ≤ Ce −λt for t ≥ 0.

Remark 2.
Understanding, that in the neighbourhood of the zero a nonlinear function with the order of nonlinearity higher that one converges to the zero quicker than a linear function, explains why asymptotic mean square stability of the zero solution of the linear part of a nonlinear system guarantees stability in probability of the zero solution of the initial nonlinear system.

Represent the system (12) in the matrix form
where the matrix C i has all zero elements besides of c ii = σ i , i = 1, 2, 3, Let y(t) be a value of the solution of the Equation (13) in the time moment t, y t = y(t + s), s < 0, be a trajectory of the solution of the Equation (13) until the time moment t. Consider a functional Denote by D the set of the functionals, for which the function V ϕ (t, y) defined in (15) has a continuous derivative with respect to t and two continuous derivatives with respect to y. Let be the sign of transpose, ∇V ϕ (t, y) and ∇ 2 V ϕ (t, y) be, respectively, the vector of the first and the matrix of the second derivatives of the function V ϕ (t, y) with respect to y, i.e., For the functionals from D the generator L of the Equation (13) Theorem 1. [31] Let there exist a functional V(t, ϕ) ∈ D, positive constants c 1 , c 2 , c 3 , such that the following conditions hold: Then the zero solution of the Equation (13) is asymptotically mean square stable.

Lemma 3.
[34] Let R ∈ R n×n be a positive definite matrix, z = Q y(s)µ(ds), where z, y(s) ∈ R n , µ(ds) is some measure on Q such that µ(Q) < ∞ and the integral is defined in the Lebesgue sense. Then

Stability of Equilibria
In this section we obtain sufficient conditions for asymptotic mean square stability of the zero solution of the Equation (13). Via Remark 2 these conditions are also sufficient conditions for stability in probability of the appropriate equilibrium of the system (6).
Below the sign * in a matrix means the symmetric element of a matrix, the inequality A > 0 (A < 0) for a symmetric matrix A means that it is a positive (negative) definite matrix.

The First Stability Condition
Theorem 2. Let there exist positive definite 3 × 3-matrices P and R satisfying the LMI Then the equilibrium (S * , I * , R * ) of the system (6) is stable in probability.
Proof. Let L be the generator of the Equation (13) and J(y t ) as it defined in (7). Using the general method of Lyapunov functionals construction [31] we will construct the Lyapunov functional for the Equation (13) in the form Let us choose the additional functional V 2 in the form Using the inequality (18) and (2), (7), we obtain So, for the functional V 2 we have From (20), (21) for the functional V = V 1 + V 2 it follows that where η(t) = (y (t), J (y t )) and the matrix Ψ 0 < 0 is defined in (19). So, the constructed above Lyapunov functional V satisfies the conditions (17) of Theorem 1. Therefore, the zero solution of the Equation (13) is asymptotically mean square stable. The proof is completed.

The Second Stability Condition
To get the second stability condition note that, using the equality d dt (13) can be represented in the form of a neutral type equation Similarly to [36,37], it can be shown that if there exists a positive definite matrix Q ∈ R n×n such that the LMI holds then the integral Equation (24) is exponentially stable.
A simpler, but, generally speaking, more rough than (25) sufficient condition for exponential stability of the integral Equation (24) is the inequality k 1 B < 1, where B is the matrix norm of a matrix B [31]. Note, however, that in the scalar case (n = 1) both these conditions coincide.
Theorem 3. Let for some positive definite matrices P, Q, R ∈ R n×n the LMIs (25) and hold. Then the equilibrium (S * , I * , R * ) of the system (6) is stable in probability.
Proof. Let L be the generator of the Equation (23). Via (16) and (23) for the functional V 1 (z t ) = Z (t)PZ(t), P > 0, we have From the definition of G(t) (23) and the inequality (18) it follows that So, for the additional functional Via (27), (28) for the functional V = V 1 + V 2 we have where η(t) = (z (t), G (t)) and the matrix Ψ 1 < 0 is defined in (26). So, the constructed above Lyapunov functional V satisfies the conditions (17) of Theorem 1. Therefore, the zero solution of the Equation (13) is asymptotically mean square stable. The proof is completed.

The Third Stability Condition
Let us complement the Equation (23) by the following way Note that for the LMI (19) the matrix A must be the Hurwitz matrix, for the LMI (26) the matrix A + B must be the Hurwitz matrix. Here we also assume that A + B is the Hurwitz matrix. The standard approach to stability analysis of the Equation (30) includes construction of a Lyapunov functional V(x t ) with the conditions that hold for some positive constants c 1 and c 2 , provided the integral equation Z(t) = 0 is asymptotically stable [36,37]. In the approach, proposed in this section, similarly to [33,38,39] so called augmented Lyapunov functional V(z t ) = V(z t , G(t)) is used, that satisfies to the conditions (17) of the Theorem 1, which is a classical theorem of the type of Lyapunov-Krasovskii. In this case there is no need to verify the stability of the integral equation Z(t) = 0.

Theorem 4.
Let there exist 3 × 3-dimensional matrices P 1 , P 2 , P 3 , R > 0 and Q > 0 that satisfy the following LMIs: and Then the equilibrium (S * , I * , R * ) of the system (6) is stable in probability.
Proof. Let L be the generator of the Equation (30). Via (30) for the functional we have Using the additional functional with Q, R > 0, we have Via the inequality (18) From (37), (38) it follows that As a result for the functional V( (33).
Note also that via (34) and via (36), (38) From this and (32) it follows with c 1 > 0 since Φ > 0. Therefore, the constructed functional V satisfies the conditions (17) of Theorem 1 and the zero solution of the system (30) is asymptotically mean square stable. The proof is completed.

Remark 5.
For numerical simulation of the solution of the SIR epidemic model (6) the algorithm of numerical simulation of trajectories of the Wiener process was used described in [31] and the Euler-Maruyama discretization [31,40,41].

Conclusions
On the example of the SIR epidemic model that is very popular in research it is shown how a sequence of stability conditions for equilibria of this model can be obtained in the presence of stochastic perturbations. The method of Lyapunov functionals construction and the method of Linear Matrix Inequalities are used for getting of these stability conditions. Besides it is shown that stability of equilibria can be investigated immediately via numerical simulation of a solution of the considered system. The proposed research can be applied to many other nonlinear mathematical models in different applications, for instance, for social epidemic models: model of alcohol consumption and model of obesity epidemic [31]. From the other hand nobody can say that investigation of the SIR epidemic model can be considered as a complete one. In connection with the current world pandemic it could be interesting and important to establish improving stability conditions for a more realistic SIR epidemic model, which takes into account the death rate and the contact rates where the contact rates depend on the total number of infections; and which under general assumptions, the recovery and death rates then become functions of the total number of infections.