Computation of the Domain of Attraction for Suboptimal Immunity Epidemic Models Using the Maximal Lyapunov Function Method

and Applied Analysis 3 whereQ is a fixed positive definite matrix and k ≥ 3. Then, one has the following Lyapunov functions V n (x) = R2 (x) + R3 (x) + ⋅ ⋅ ⋅ + Rn (x) 1 + Q 1 (x) + ⋅ ⋅ ⋅ + Q n−2 (x) . (3)


Introduction
The computing of domain of attraction (DOA), that is, the region where the dynamical system is asymptotically stable, is an interesting research topic in the stability analysis of nonlinear systems such as the systems for compartmental ODE epidemic models.In other words, the mathematical analysis of epidemic models often involves computing the asymptotic stability region for both the disease-free equilibrium and endemic equilibrium.The set of initial states whose corresponding trajectories converge to an asymptotically stable equilibrium point as time increases is known as the stability region or domain of attraction (DOA) of the equilibrium under study.If the initial state lies within the DOA, the disease will evolve towards an endemic state.On the contrary, if the initial state is outside the DOA, the system will converge to a disease-free state.Therefore, it is important to study the DOA of endemic equilibrium.
Lyapunov's second method (the direct method) is generally used to analyze the stability of epidemic models.In this method, the asymptotic stability of the origin can be examined if a positive definite function whose derivative along the solutions of the system is negative definite.However, it is not only difficult to construct the Lyapunov Function, but also hard to guarantee the asymptotical stability of the equilibrium.Apart from that, it is known even if the Lyapunov function exists in an autonomous ODE, it may not be unique.A maximal Lyapunov function is a special Lyapunov function on  (where  denotes the DOA) which can be used to determine the DOA for a given locally asymptotical stable equilibrium point.
Considerable work on DOA estimation and optimized DOA for epidemic dynamical models has been done.In [1], the authors had computed the DOA in epidemiological models with constant removal rates of infected individuals.An optimization approach for finding the DOA of a class of SEIR models, based on the sum of square optimization, is presented in [2].Recently, the authors in [3] had adopted a recurrence formula established by Kaslik et al. [4] by using an -analytical function and the sequence of its Taylor polynomial to construct the Lyapunov function, and solved the linear matrix inequality (LMI) relaxations of a global optimization problem to obtain the DOA.However, all the epidemic models in the papers mentioned above are limited to relatively simple epidemic models, without taking into account nonlinear incidence rates or saturated recovery rates.In this paper, we study the DOA for the suboptimal immunity models with nonlinear incidence rates and saturated recovery rates, by utilizing the maximal Lyapunov function in [5].Throughout the paper, we focus on DOA for the suboptimal immunity models.This kind of suboptimal immunity models is more appropriate for the study of microparasite infections which usually occurs during childhood.After a primary infection, one may get temporary immunity (namrly, immune protection that will wane over time) or partial immunity (namrly, immunity that is not fully protective).Examples of this kind of diseases include Pertussis (temporary immunity) and Influenza (partial immunity).
The rest of the paper is organized as follows.In Section 2, we will establish a theorem based on [5] and an iterative procedure for the construction of the Maximal Lyapunov Function.In Section 3, we will briefly explain the suboptimal immunity model.In Section 4, two examples of suboptimal immunity models are given to demonstrate the validity of the procedure.A conclusion is then given in Section 5.

Maximal Lyapunov Function
Consider the following system ẋ =  (x) , where  : R  → R  is an analytical function with the following properties: (a) (0) = 0, that is, x = 0 is an equilibrium point of system (1); (b) all the eigenvalues of the Jacobian matrix at x = 0, that is, (/x)(0), have negative real parts, namely, x = 0 is an asymptotically stable equilibrium point.It is well known that in the Lyapunov sense, if there exists a Lyapunov function for the equilibrium point x = 0 of the system (1), then x = 0 is asymptotically stable.
Definition 1 (Lyapunov function).Let (x) be a continuously differentiable real-valued function defined on a domain (0) ⊆ R  containing the equilibrium point x = 0.The function (x) is a Lyapunov function of the equilibrium x = 0 of the system (1) if the following conditions hold: (a) (x) is positive definite on (0); (b) the time derivative of ( ẋ ) is negative definite on (0).
If (x) is a Lyapunov function which fulfills the conditions in the Definition 1, the estimation of DOA is given by the following definition.Definition 2. Given an autonomous system (1), where x ∈ R  and (0) = 0, the domain of attraction (DOA) of x = 0 is   = {x 0 ∈ R  : lim  → ∞ x(, x 0 ) = 0}, where x(⋅, x 0 ) denotes the solution of the autonomous system corresponding to the initial condition x(0) = x 0 .The Lyapunov function is not unique.A maximal Lyapunov function, (x), is a special Lyapunov function on  (where  denotes the DOA) which can be used to determine the DOA for a given locally asymptotically stable equilibrium point.
We have the following definition for the DOA of an asymptotically stable equilibrium point which derives from the maximal Lyapunov function.From the above definition, based on the work in [5], we derive the following theorem.

Theorem 5. Consider the nonlinear system of equations
, where   (⋅) is a homogeneous function of degree .Suppose that the linearized system ẋ =  1 (x) = x is asymptotically stable at x = 0. Let   ,   be homogeneous functions of degree , and the functions   and   satisfy the following recursive equations: where  is a fixed positive definite matrix and  ≥ 3.Then, one has the following Lyapunov functions Proof.Rewrite which satisfy condition (a) in the Definition 3. Differentiating (4) with respect to x, we have One choice to ensure that V  (x) is negative definite is V  (x) = −x  x.Then from (5), we have Equating the coefficients of the same degree  of the two sides of (6), we get the following recursive relations: where  is a fixed positive definite matrix, and  ≥ 3.
Based on Theorem 5, the procedure for obtaining the maximal Lyapunov function and calculating the DOA is established as follows.
Step 1. From the linearized system,  1 = x, find  > 0 such that    +  = −, then set where  2 =  1  2 +  2  +  3  2 .In this case,  is a fixed positive definite matrix.Hence, one of the good choices for  is the identity matrix.
Step 2. For  = 3, we have where Equating the coefficients of same degree in (10), we will obtain a system of linear equations in terms of  1 ,  2 ,  3 ,  4 ,  1 and  2 .The solution of these linear equations will be used as constraints in the minimization problem to get   () in the later steps.
Step 3.For  = 4, we have where Then, we solve the system of linear equations as in Step 2.
Step 4 (optional).For  = 5, we have where Then, the system of linear equations is solved.We should address that equations similar to (12) can also be obtained for  > 5.
For each of the Steps 2 and 4, one will lead to a number of choices for the value of the coefficients for   and  −2 .Consider where   () is the squared 2-norm of the coefficients of the terms with degree greater than or equal to  + 1 in the expression of V  .This ensures that V  (x) is negative definite over a neighbourhood of the origin.To make it as similar as possible to V 2 (x) = −x  x, we take   () as small as possible.Hence, it creates a new condition that can be formulated as a minimization problem where the constraints are obtained from the recursive relations in each of the steps above.
Step 5. Once we get   () sufficiently small, say at Step 3, we can have the maximal Lyapunov function as To obtain the DOA, one needs to find the largest possible value  * when  4 (x) =  * such that the interior of the resulting ellipsoid is entirely bounded within the region given by Ω = {x : V  (x) ≤ 0}.In this case, one can determine  * by solving an optimization problem: subject to the constraints V 4 (x) = 0. (15) Then, the set   = {x :  4 (x) <  * } is contained in the DOA .Appropriate  * also can be determined manually as suggested in [6].In this case, one can choose the largest positive value  * such that the sublevel set   = {x :  4 (x) <  * } is contained in the region given by {x : V 4 (x) < 0}.Hence, we obtain the DOA in the form of   .

Suboptimal Immunity Epidemic Models
In this paper, we estimate the domain of attraction (DOA) for suboptimal immunity epidemic models with saturated treatment/recovery rate and nonlinear incidence rate.Apart from using the saturated treatment/recovery rate, an additional parameter  is used to form the suboptimal immunity model as in [7].The new model lies in between the SIS and SIR models.
The suboptimal immunity model with nonlinear incidence rate and saturated treatment/recovery rate is as follows: where all the parameters are positive.We assume that the population is fixed, namely, / = () + () + () and where  denotes susceptible population,  represents infective population, and  is the recovered population. is the recruitment rate of susceptible population,  is the disease transmission rate,  is the natural death rate, and () is the recovery rate.We take () =  + /(1 + ) as the recovery rate function in which /(1 + ) and  are, respectively, the recovery rate of the infected population with and with no treatment.The function (, ) denotes the incidence rate.
In comparison with previous models, our model presented here has various new features and contributions.Firstly, it is more general and includes some previous models as special cases.For example, if we take (, ) as  2 , we will have the nonlinear incidence rate.If we take (, ) as a bilinear function, then it reduces to the suboptimal immunity model [7], while it reduces to the nonlinear SIR model if () is taken to be zero.It should also be addressed that  = 1 corresponds to the SIS model in which immunity is assumed not to protect against reinfection, while  = 0 corresponds to the SIR model in which immunity is assumed to be fully protective and prevents any reinfection.The suboptimal immunity model where  ∈ [0, 1] is more appropriate for the study of microparasite infections which usually occur during childhood.After a primary infection, one may get temporary immunity (namely, immune protection that wanes over time) or partial immunity (namely, immunity that is not fully protective).Examples of this kind of diseases include Pertussis (temporary immunity) and Influenza A (partial immunity) [8].Secondly, due to the combination of the nonlinearities in both incidence rate and recovery rate, it is hard to obtain exact solution.Hence, estimating its domain of attraction using the numerical procedure is important in order to know whether the disease in a particular state will evolve towards an endemic state or converge to a disease free state.
Example 2. We consider the following reduced system for the suboptimal immunity model with nonlinear incidence rate and saturated treatment rate: For this example, we consider the nonlinear incidence rate  2 .According to [9], one of the reasons to consider the nonlinear incidence rate,     , is to represent heterogeneous mixing.Take (, , , , , , ) = (6, 0.5, 1.27, 2, 4, 1, 0.5).
By applying the same procedure, we calculate the DOA when  = 12, and we obtain the DOA as in Figure 5.It is clear that with (, , , , , , ) = (6, 0.5, 1.27, 2, 12, 1, 0.5), increasing the value of  will increase the DOA for the suboptimal model given by (20).

Concluding Remarks
In this paper, we deal with the problem of estimating the domain of attraction (DOA) for the suboptimal epidemic model.We have successfully established a procedure to determine the maximal Lyapunov function in the form of rational functions and compute the domain of attraction for epidemic models.Determination of the DOA is extremely important in order to understand the dynamic behaviour of the transmission of disease as a function of the initial population distribution.In our first example, we show that, for certain values of the parameter, a larger  value (i.e., the model is more toward the SIR model) leads to a smaller DOA.In our second example, we show that within certain values of the parameter, decreasing the  value will yield a smaller DOA.

Definition 4 .
Suppose we can find a set  ⊆ R  containing the origin in its interior and a continuous function (x) :  → R + such that (a) (x) is positive definite on ; (b) V(x) is negative definite on ; (c) (x) → ∞ as x →  and/or as ‖x‖ → ∞.Then  = , where  denotes the DOA.