#### Abstract

We are concerned with the estimation of the domain of attraction (DOA) for suboptimal immunity epidemic models. We establish a procedure to determine the maximal Lyapunov function in the form of rational functions. Based on the definition of DOA and the maximal Lyapunov function, a theorem and subsequently a numerical procedure are established to determine the maximal Lyapunov function and the DOA. Determination of the domain of attraction for epidemic models is very important for understanding the dynamic behaviour of the disease transmission as a function of the state of population distribution in different categories of disease states. We focus on suboptimal immunity epidemic models with saturated treatment rate and nonlinear incidence rate. Different from classical models, suboptimal immunity models are more realistic to explain the microparasite
infection diseases such as Pertussis and Influenza A. We show that, for certain values of the parameter, larger *k* value (i.e., the model is more toward the SIR model) leads to a smaller DOA.

#### 1. 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.

#### 2. Maximal Lyapunov Function

Consider the following system where is an analytical function with the following properties: (a), that is, is an equilibrium point of system (1); (b)all the eigenvalues of the Jacobian matrix at , that is, , have negative real parts, namely, 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 of the system (1), then is asymptotically stable.

*Definition 1 (Lyapunov function). *Let be a continuously differentiable real-valued function defined on a domain containing the equilibrium point . The function is a Lyapunov function of the equilibrium of the system (1) if the following conditions hold: (a) is positive definite on ; (b)the time derivative of is negative definite on .

If 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 and , the domain of attraction (DOA) of is , where denotes the solution of the autonomous system corresponding to the initial condition .

The Lyapunov function is not 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 asymptotically stable equilibrium point.

*Definition 3 (maximal Lyapunov function [5]). *A function is called a maximal Lyapunov function for the system (1) if (a), , for all ; (b) if and only if ; (c) as and/or ; (d) is well defined and negative definite over , where denotes the DOA.

We have the following definition for the DOA of an asymptotically stable equilibrium point which derives from the maximal Lyapunov function.

*Definition 4. *Suppose we can find a set containing the origin in its interior and a continuous function such that (a) is positive definite on ; (b) is negative definite on ; (c) as and/or as . Then , where denotes the DOA.

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 is asymptotically stable at . Let , be homogeneous functions of degree , and the functions and satisfy the following recursive equations:
**
where is a fixed positive definite matrix and . Then, one has the following Lyapunov functions
*

*Proof. *Rewrite
which satisfy condition (a) in the Definition 3. Differentiating (4) with respect to , we have
One choice to ensure that is negative definite is . 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 .

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, , find such that , then set
where . In this case, is a fixed positive definite matrix. Hence, one of the good choices for is the identity matrix.

*Step 2. *For , we have
where and . Equating the coefficients of same degree in (10), we will obtain a system of linear equations in terms of , , , , and . The solution of these linear equations will be used as constraints in the minimization problem to get in the later steps.

*Step 3. *For , we have
where and . Then, we solve the system of linear equations as in Step 2.

*Step 4 (optional). *For , we have
where and . Then, the system of linear equations is solved. We should address that equations similar to (12) can also be obtained for .

For each of the Steps 2 and 4, one will lead to a number of choices for the value of the coefficients for and . Consider
where is the squared 2-norm of the coefficients of the terms with degree greater than or equal to in the expression of . This ensures that is negative definite over a neighbourhood of the origin. To make it as similar as possible to , 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 such that the interior of the resulting ellipsoid is entirely bounded within the region given by . In this case, one can determine by solving an optimization problem:
Then, the set 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 is contained in the region given by . Hence, we obtain the DOA in the form of .

#### 3. 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 as the recovery rate function in which 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 , 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 corresponds to the SIS model in which immunity is assumed not to protect against reinfection, while corresponds to the SIR model in which immunity is assumed to be fully protective and prevents any reinfection. The suboptimal immunity model where 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.

#### 4. Numerical Examples

*Example 1. *We consider the following reduced system for the suboptimal immunity model with bilinear incidence rate and saturated treatment rate:
where and . Rewrite and for and , and translate the equilibrium point to the origin by using and . By using the numerical procedure in Section 2, we have the following result:
and . Thus, is an estimate of for the system (17) when . This estimate and its phase portrait are given in Figure 1.

Consider smaller , which means the model is more toward the SIS model. Let . Rewrite and for and , and translate the equilibrium point to the origin by using and . By using the numerical procedure as given in Section 2, we have the following result: and . Thus, is an estimate of for the system (17) when . This estimate and its phase portrait are given in Figure 2.

*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 . According to [9], one of the reasons to consider the nonlinear incidence rate, , is to represent heterogeneous mixing. Take . Rewrite and for and , and translate the equilibrium point to the origin by using and . By using the numerical procedure in Section 2, we have the following result:
and . Thus, is an estimate of for the system (20) when . This estimate and its phase portrait are given in Figure 3.

**(a)**

**(b)**

Here, we study the effect of the value on the DOA. With the value of which is ten fold of previous calculation, let . Rewrite and for and , and translate the equilibrium point to the origin by using and . By using the numerical procedure as in Section 2, we have the following result: and . Thus, is an estimate of the DOA for the system (20) when . This estimate and its phase portrait are given in Figure 4.

By applying the same procedure, we calculate the DOA when , and we obtain the DOA as in Figure 5. It is clear that with , increasing the value of will increase the DOA for the suboptimal model given by (20).

#### 5. 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.

#### Acknowledgments

The authors are grateful to the anonymous reviewers for their helpful suggestions and comments. The third author wishes to acknowledge the support of the Faculty of Science, Mahidol University.