\documentclass[reqno]{amsart}
\usepackage{hyperref}
\usepackage{graphicx}

\AtBeginDocument{{\noindent\small
\emph{Electronic Journal of Differential Equations},
Vol. 2017 (2017), No. 121, pp. 1--18.\newline
ISSN: 1072-6691. URL: http://ejde.math.txstate.edu or http://ejde.math.unt.edu}
\thanks{\copyright 2017 Texas State University.}
\vspace{8mm}}

\begin{document}
\title[\hfilneg EJDE-2017/121\hfil Dynamics of a SIRC epidemiological model]
{Dynamics of a SIRC epidemiological model}

\author[H. Li, S. Guo \hfil EJDE-2017/121\hfilneg]
{Haijiao Li, Shangjiang Guo}

\address{Haijiao Li  \newline
School of Business Administration, and 
College of Mathematics and Econometrics,
Hunan University,
Changsha, Hunan 410082, China}
\email{601865814@qq.com}

\address{Shangjiang Guo \newline
College of Mathematics and Econometrics,
Hunan University,
Changsha, Hunan 410082, China}
\email{shangjguo@hnu.edu.cn}

\dedicatory{Communicated by Zhaosheng Feng}

\thanks{Submitted  May 18, 2016. Published May 4, 2017.}
\subjclass[2010]{34D20, 92D30}
\keywords{SIRC model; cross-immunity; global stability; 
\hfill\break\indent basic reproduction number}

\begin{abstract}
 This article concerns the SIRC epidemiological model for influenza A,
 which efficiently describes the mechanism of disease spreading,
 including the susceptible (S), the infected (I) and the recovered (R),
 along with a cross-immune class (C) that recovers after being inflected
 by different strains of the same viral subtype.
 The dynamics of the model is completely determined by the basic reproduction
 number $R_0$. If $R_0\leq 1$, the disease-free equilibrium of the SIRC
 model is globally asymptotically stable, which means influenza A will die out.
 Otherwise, the SIRC model may have exactly one endemic equilibrium which
 is globally asymptotically stable under certain parametric conditions.
 Also, numerical simulations are given to support our analytical results.
\end{abstract}

\maketitle
\numberwithin{equation}{section}
\newtheorem{theorem}{Theorem}[section]
\newtheorem{lemma}[theorem]{Lemma}
\newtheorem{remark}[theorem]{Remark}
\newtheorem{corollary}[theorem]{Corollary}
\allowdisplaybreaks


\section{Introduction}

Influenza \cite{Earn}, a member of the family Orthomyxoviridae, is a RNA 
virus which can give rise to epidemic disease between mankind and animals.
In general, influenza is primarily divided into A, B and C types and every 
type contains a wide variety of subtypes according to hemagglutinin (HA) 
and neuraminidase (NA) differences. From an epidemiologic standpoint,
influenza A is the most common and the most terrible virus among three types, 
and can result in the highest pathogenicity because of the easiest way to 
generate variation. What's worse, the virus has brought about more than century 
pandemic influenza in the past years. The pathogenicity of influenza B virus
is the same as type A, but performed studies are shown that type B virus does not 
contribute to century pandemic. Last, influenza C virus just leads to 
unconspicuous or feeble respiratory infection and almost doesn't incur a pandemic 
either.

In this article, we pay attention to influenza A virus whose surface often 
has a tiny variation, which is referred to as drift. That is, the virus
camouflages itself by a subtle change, and it thereby can elude identificaton
by the human immune system. As a result of the drift,
all kinds of flu strains appear every year. 
Therefore, human annually need to take a vaccine against influenza for prevention. 
On the other hand, shift means influenza A takes place genetic mutations causing 
a new ``subtype''. In this case, it can lead to outbreaking a global pandemic 
influenza. For example, Spanish outbroke H1N1 influenza virus in 1918-1919
about the death toll up to 20-40 million \cite{Taubenberger}. In 1957-1958, 
Asia influenza erupted due to the result of  H2N2 virus. Furthermore, influenza 
A virus which threats the lives of even the healthiest individuals brings 
serious economic loss for many countries. Some research data indicate the direct 
economic loss up to 10-30 billion dollars in America, and the potential economic 
losses at 10-15 billion dollars \cite{Thomas}. Consequently, devoting to 
studying the spread of influenza A is extremely indispensable.

Mathematically modelling, as a bridge, plays a crucial role in understanding 
of the spread of the epidemic. Researchers transform their focus
from considering threshold value of pathophoresis to prevent epidemic propagation 
by taking some strategies for control \cite{Anderson}.
The model of influenza is parallel to the traditional standard SIR model in which 
the crowd are segmented into three compartments:
susceptible (S), infectious (I), recovered (R). A number of compartmental models 
have been  established based on this idea \cite{Zhong}. Nonetheless, some scholars 
indicate that the traditional SIR model is not adequate
to describe actual situation of influenza spreading because every type can evolve 
different subtypes. Hence,  Andreasen et al. \cite{Andreasen2} initially considered 
the dynamics of system with multiple virus strains consisting of partial 
cross-immunity \cite{Smith}, concluding that the system exists stable equilibrium 
when the number of strains are less than or equal to three. 
However, in the work of Lin et al.\ \cite{Lin}, they demonstrated that 
oscillations can be sustained under a linear chain of three cocirculating 
influenza A strains. These researchers only consider a special case where 
cross-protection is symmetric to analyze complex
system because the analysis and computation of multiple viral strains are much 
intractable. Recently, Minayev and Ferguson \cite{Minayev} proposed a
deterministic model of multi-strain pathogens with symmetric equilibrium, 
self-organized strain structures, regular periodic and chaotic regimes, which
is determined by cross-immune response function.
In addition to these properties, Kooi et al.\ \cite{Kooi} revealed bifurcation 
analysis, Lyapunov exponent calculation as well as quantitative and
qualitative results by numerical simulations in three multi-strain compartment 
models.
Nu\~{n}o et al.\ \cite{Nuno} developed a general multiple strains pathogens model 
with all kinds of cross-immunity structures. They illustrated the weaker 
cross-immunity structures are more likely to appear instability in the strain 
coexistence mode. For seasonal influenza, there is another case where strong
and weak cross-immunity can result in coexistence with the following pandemic. 
On the contrary, the intermediate level one may take the place of the seasonal
subtype.
Chung and Liu \cite{Chung} extended the result of Nu\~{n}o et al. \cite{Nuno1} 
and elucidated the local asymptotic stability of a two-strain influenza model.
Also, the stability may be lost when two strains are far apart.
 All the above published works didn't  analyze global asymptotical stability
 of models with cross-immune class. In this paper, we investigate a
SIRC model with cross-immune class proposed by  Casagrandi et al.\ \cite{Renato}, 
which is represented by the following system of four ordinary differential equations,
\begin{equation}\label{eq}
\begin{gathered}
\frac{dS(t)}{dt}=\gamma(1-S)-\beta SI+\eta C,\\
\frac{dI(t)}{dt}=\beta SI+\mu \beta CI-(\gamma + \alpha)I,\\
\frac{dR(t)}{dt}=(1-\mu) \beta CI+\alpha I-(\gamma + \delta)R,\\
\frac{dC(t)}{dt}=\delta R - \beta CI-(\gamma + \eta)C,
\end{gathered}
\end{equation}
with initial conditions:
$$
(S(0),I(0),R(0),C(0))^T\in \mathbb{R}^4_+
:=\big\{(S,I,R,C)^T\in \mathbb{R}^4: S\geq 0,I\geq 0,R\geq 0,C\geq 0\big\}
$$
and positive real parameters $\gamma$, $\beta$, $\eta$, $\mu$, $\alpha$,  
and $\delta$. The parameters $\alpha$, $\delta$ and $\eta$ are the inverses 
of the average time spent by the individuals in each of the three compartments 
$I$, $R$, and $C$, respectively. The parameter $\gamma$ denotes the mortality 
rate in every compartment and is assumed to equal to the rate of newborn 
in the population. The parameter $\mu$ is interpreted as the average reinfection 
probability of a cross-immune subject, whereas the parameter $\beta$ is the 
contact rate. The state space of  model \eqref{eq} is $\mathbb{R}^4_+$.
The novelty of this model is that it takes into account the presence of 
cross-immune (C) subjects, i.e., subjects that are temporarily immune. 
It efficiently describes the mechanism for influenza A virus spreading.
In this paper, we shall not only consider existence, number, and local 
asymptotical stability of equilibrium for the model, but also  verify 
the globally asymptotical stability. By analyzing the SIRC model, 
we can gain the law of development of influenza A. In some sense, 
this paper puts forward an important theoretical approach and decision-making 
basis in order to prevent influenza A spreading.

Let
$$
R_0=\frac{\beta}{\gamma+\alpha}
$$
which is called the basic reproduction number \cite{Anderson1,Kermack} 
or the contact number that represents the average number of secondary
infections from a single infections host. We can prove that the dynamics 
of model \eqref{eq} is completely determined by the threshold value $R_0$. 
If $R_0\leq 1$ then the disease-free equilibrium $E_0$ is globally asymptotically
 stable and hence the disease will die out (see Theorem \ref{thm42}). 
If $R_0> 1$ then the unique endemic equilibrium $E^*$ is globally asymptotically 
stable so that  the disease always persists at the unique endemic level  
(see Theorem \ref{thm63}). Thus, we can conclude that the spread of the disease 
should be controlled by way of suitable protection measures of the society to
 reduce the value of $\beta$ (transmission rate of disease) when susceptible 
individuals contact with infected individuals. This can be done by adopting 
some strategies to detect early cases among the
passengers coming from the infected countries and
individuals should follow simple steps like cough etiquettes, stay away from persons
coughing or sneezing, avoid gathering and so on. Disease spreading can also be 
kept under control by increasing $\alpha$ (recovery rate of the infected population). 
It is recommended that if one feels any respiratory distress, one should report to a
nearby hospital immediately. The global stability of $E_0$ when $R_0\leq 1$ can be 
routinely proved by using a well-known Lyapunov function, but the global stability 
of $E^*$ when $R_0>1$ has been an open problem in the literature due to the high 
dimensionality for the model. Our proof is based on a theoretical approach 
developed in \cite{Li,Li2}.

From a mathematical viewpoint, it is interesting to notice that, in the absence 
of cross-immunity ($\mu=1$), the two classes $S$ and $C$ are immunologically 
indistinguishable,
since
$$
\frac{d(S+C)}{dt}=\gamma[1-(S+C)]-\beta(S+C)I+\delta R.
$$
Consequently, in the limit of $\mu\to 1$, the SIRC model reduces to the following
 classical SIRS model
\begin{equation}\label{SIRS}
\begin{gathered}
\frac{dS(t)}{dt}=\gamma(1-S)-\beta SI+\delta R,\\
\frac{dI(t)}{dt}=\beta SI-(\gamma + \alpha)I,\\
\frac{dR(t)}{dt}=\alpha I-(\gamma + \delta)R.
\end{gathered}
\end{equation}
Thus, our results generalize some earlier results on the SIRS model. 
In this paper, we shall prove the global stability of a unique endemic 
equilibrium of \eqref{SIRS} when the basic reproduction number $R_0$ 
is greater than 1.

This article is organized as follows. The positively invariant set 
is shown in section 2. Section 3 is devoted to the existence of equilibrium 
including a disease-free equilibrium and a unique endemic equilibrium. 
The stability of the disease-free equilibrium and the endemic equilibrium 
is presented in sections 4 and 6, respectively. Section 5 is devoted to 
the persistence of the model and the work of numerical simulation reveals 
in section 7. Finally, this paper ends up with a brief conclusion in section 8.

\section{Positively invariant set}

This section is devoted to proving the positivity and boundedness of solutions 
of model \eqref{eq}  with  initial conditions 
$(S(0),I(0),R(0),C(0))^T\in \mathbb{R}^4_+$.
We first introduce the following lemma.

\begin{lemma}[\cite{ref21}]\label{lem21}
Suppose $\Omega\subset \mathbb{R} \times \mathbb{C}^{n}$ is open,
 $f_{i}\in C(\Omega,\mathbb{R})$, $i=1,2,3\ldots n$.
If $f_{i}|_{x_i(t)=0,X_{t}\in \mathbb{C}_{+0}^n}\geq0$, 
$X_{t}=(x_{1t},x_{2t},\ldots, x_{nt})^T$, $i=1,2,3\ldots n$, then  
$\mathbb{C}_{+0}^n$
is the invariant domain of the following equations
$$
\frac{d x_{i}(t) }{dt}=f_{i}(t,X_{t}), \quad  t\geq\sigma,\; i=1,2\ldots n.
$$
\end{lemma}


\begin{theorem}\label{thm22}
Each solution $(S(t),I(t),R(t),C(t))$ of model \eqref{eq} with the 
non-negative initial conditions is non-negative for all $t>0$.
\end{theorem}

\begin{proof}
Let $X=(S,I,R,C)^{\mathrm{T}}$ and
 $f(X)=(f_1(X),f_2(X),f_3(X),f_4(X))^{\mathrm{T}}$ then we can
rewrite model \eqref{eq} as
$$
\dot{X}=f(X)
$$
where
$$
f(X)=\begin{pmatrix}
f_1(X)\\
f_2(X) \\
f_3(X) \\
f_4(X)
\end{pmatrix}
= \begin{pmatrix}
\gamma(1-S)-\beta SI+\eta C\\
\beta SI+\mu \beta CI-(\gamma + \alpha)I\\
(1-\mu) \beta CI+\alpha I-(\gamma + \delta)R\\
\delta R - \beta CI-(\gamma + \eta)C
\end{pmatrix}.
$$
Note that
\begin{gather*}
\frac{dS(t)}{dt}|_{S=0}=r+\eta C>0, \quad
\frac{dI(t)}{dt}|_{I=0}=0,\\
\frac{dR(t)}{dt}|_{R=0}=(1-\mu)\beta CI+\alpha I \geq 0,\quad
\frac{dC(t)}{dt}|_{C=0}=\delta R\geq 0.
\end{gather*}
Then it follows from Lemma \ref{lem21} that $\mathbb{R}^4_+$ is invariant set.
\end{proof}

\begin{theorem}\label{thm23}
$\mathbb{D}=\{(S,I,R,C)^T\in \mathbb{R}_{+}^{4}:0\leq S+I+R+C\leq 1\}$
 is a positively invariant set and also a globally attractive set of  model
 \eqref{eq}.
\end{theorem}

\begin{proof}
Consider the total population $ N(t)=S(t)+I(t)+R(t)+C(t)$.
 Direct calculation leads to
\begin{equation}\label{eq21}
\frac{dN}{dt}=\gamma-\gamma N.
\end{equation}
Solving this equation, we obtain
\begin{align*}
N(t)=1-(1-N(0))e^{-\gamma t}.
\end{align*}
It is straight forward to show that $N(t)\leq 1$ if $N(0)\leq 1$. 
This means that the set $\mathbb{D}$ is a positively invariant set for 
model \eqref{eq}.
If $N(0)> 1$ then it turns out that $\lim_{t\to\infty}N(t)=1$.
Thereby the set $\mathbb{D}$ is the globally attractive set for model \eqref{eq}.
\end{proof}

\begin{remark}\label{rm24}\rm
In view of the above two theorems, we see that the  positively invariant 
set $\mathbb{D}$ can attract every solution with initial conditions starting 
in its state space $\mathbb{R}^4_+$.
Namely, every trajectory of model \eqref{eq} with the initial conditions in 
$\mathbb{R}^4_+$ eventually stays in $\mathbb{D}$.
\end{remark}

\section{Existence of equilibria}

In this section, we consider the existence, type, and number of the equilibria.
 Obviously, model \eqref{eq} has two equilibria: one is disease-free equilibrium 
$E_0=(1,0,0,0)$ which exists for all parameter values; and the other is the endemic 
equilibrium $E^*=(S^*,I^*,R^*,C^*)$, which is a positive solution of the 
following equation
\begin{equation}\label{eqeq}
\begin{gathered}
\gamma(1-S^*)-\beta S^*I^*+\eta C^*=0,\\
\beta S^*I^*+\mu \beta C^*I^*-(\gamma + \alpha)I^*=0,\\
(1-\mu) \beta C^*I^*+\alpha I^*-(\gamma + \delta)R^*=0,\\
\delta R^* - \beta C^*I^*-(\gamma + \eta)C^*=0.
\end{gathered}
\end{equation}
It follows from the third and forth equations of \eqref{eqeq} that
\begin{equation}\label{eqeq1}
C^*=\frac{\delta\alpha I^*}{(\delta\mu+\gamma)\beta I^*
+(\gamma+\delta)(\gamma+\eta)}.
\end{equation}
Next, substituting \eqref{eqeq1} in the second equation of \eqref{eqeq}, we obtain
\begin{equation}\label{eqeq2}
S^*=\frac{\gamma+\alpha}{\beta}-\frac{\mu\delta\alpha}{(\delta\mu+\gamma)
\beta I^*+(\gamma+\delta)(\gamma+\eta)}\cdot I^*.
\end{equation}
Combining \eqref{eqeq1} with \eqref{eqeq2} and then substituting them in the 
first equation of \eqref{eqeq}, we have
\begin{equation}\label{eqeq3}
a I^{*2}+b I^*+d=0,
\end{equation}
where
\begin{gather*}
a=\beta^2(\gamma+\alpha)(\delta\mu+\gamma)-\beta^2\mu\delta\alpha
=\beta^2(\gamma\delta\mu+\gamma^2+\alpha\gamma)>0,\\
\begin{aligned}
b&=\beta\gamma(\gamma+\alpha)(\delta\mu+\gamma)+\beta(\gamma+\alpha)
 (\gamma+\delta)(\gamma+\eta) \\
&\quad -\gamma\beta\mu\delta\alpha-\gamma\beta^2(\delta\mu+\delta)
 -\beta\eta\delta\alpha,
\end{aligned}\\
d=\gamma(\gamma+\alpha)(\gamma+\delta)(\gamma+\eta)
-\gamma\beta(\gamma+\delta)(\gamma+\eta).
\end{gather*}
Therefore, we have the following result.

\begin{theorem}\label{thm31}
If $R_0\leq1$, model \eqref{eq} always has a disease-free equilibrium
$E_0=(1,0,0,0)$. If $R_0>1$, the model has exactly one endemic equilibrium
$E^*=(S^*,I^*,R^*,C^*)$.
\end{theorem}

\begin{proof} 
It is not difficult to observe that model \eqref{eq} has a disease-free 
equilibrium $E_0(1,0,0,0)$. Now, we consider the existence of the endemic
equilibrium $E^*(S^*,I^*,R^*,C^*)$ when $R_0>1$.
It follows from $R_0>1$ that $d<0$ and hence that \eqref{eqeq3} has exactly
 one positive solution $I^*$, which, together with
\eqref{eqeq1}, implies that $C^*$ is positive.
Similarly, from the first and forth equations of \eqref{eqeq}, we see that 
both $S^* $ and $ R^*$ are positive.
Therefore, model \eqref{eq} has exactly one endemic equilibrium $E^*$
when $R_0>1$.
\end{proof}


\section{Stability of disease-free equilibrium}

This section is devoted to the local and global stability of the disease-free 
equilibrium.

\begin{theorem}\label{thm41} 
The disease-free equilibrium $E_0(1,0,0,0)$ is locally asymptotically stable
if $R_0<1$.
If $R_0>1$ then $E_0$ is unstable and all solutions starting from
sufficiently close to $E_0$ in $\mathbb{D}$ ultimately leave away from $E_0$,
except those starting on $S$-axis close to $E_0$ along this.
\end{theorem}

\begin{proof}
The Jacobian matrix of  \eqref{eq} at $E_0$ is
$$
\begin{pmatrix}
-\gamma & -\beta & 0 & \eta \\
0 & \beta-(\gamma+\alpha) & 0 & 0 \\
0 & \alpha & -(\gamma+\delta) & 0\\
0 & 0 & \delta &-(\gamma+\eta)
\end{pmatrix}.
$$
Its characteristic equation is
$$
(\lambda+\gamma)[\lambda-(\beta-(\gamma+\alpha))]
(\lambda+\gamma+\delta)(\lambda+\gamma+\eta)=0.
$$
If $R_0<1$ then all the characteristic roots are less than 0.
That is, $E_0$ is locally asymptotically stable.
If $R_0>1$ then there exists a characteristic value $\beta-(\gamma+\alpha)>0$.
Therefore, $E_0$ is unstable.
Moreover, the trajectory starting from sufficiently close to $E_0$ will be
away from a neighborhood of $E_0$
except that those are on the $S$-axis, where model \eqref{eq} can translate 
into $\frac{dS(t)}{dt}=\gamma(1-S(t))$ and hence $\lim_{t\to\infty}S(t)=1$. 
This completes the proof.
\end{proof}

Next  we will prove that $E_0$ is globally asymptotically stable
by means of Lyapunov function.

\begin{theorem}\label{thm42}
The disease-free equilibrium $E_0$ for \eqref{eq} is globally asymptotically
stable when $R_0\leq 1$.
\end{theorem}

\begin{proof}
Consider a Lyapunov function $V=V(S,I,R,C)$ defined by
$$
V=\frac{1}{2}(S-1+I+R+C)^2.
$$
The time derivative of $V$ along a solution of  \eqref{eq} is
$$
\frac{dV}{dt}=-\gamma(S-1+I+R+C)^2\leq 0.
$$
We see that $\frac{dV}{dt}$=0 if and only if $S+I+R+C=1$. 
Consequently, the maximal
invariant set in $\{(S,I,R,C)^T\in \mathbb{R}^4_+:\frac{dV}{dt}=0\}$ is in 
$\{(S,I,R,C)^T\in \mathbb{R}^4_+:S+I+R+C=1\}$. Therefore,  to prove that
the disease-free equilibrium $E_0$ of \eqref{eq} is globally asymptotically
stable, it suffices to show that the equilibrium $(0,0,0)$ of the following 
model is globally asymptotically stable:
\begin{equation}\label{eq41}
\begin{gathered}
\frac{dI(t)}{dt}=\beta I (1-I-R-C)+\mu \beta CI-(\gamma + \alpha)I,\\
\frac{dR(t)}{dt}=(1-\mu) \beta CI+\alpha I-(\gamma + \delta)R,\\
\frac{dC(t)}{dt}=\delta R - \beta CI-(\gamma + \eta)C.
\end{gathered}
\end{equation}
Again define a Lyapunov function by
$$
V_1=\frac{1}{2}(R+C)^2+\frac{k_1}{2}R^2+k_2I+\frac{k_3}{2}C^2,
$$
where
$$
k_1=\frac{\mu}{1-\mu},\quad
k_2=\frac{\alpha}{(1-\mu)\beta},\quad
k_3=\frac{2\gamma+\eta}{\delta}.
$$
Using a similar argument as above, we have
\begin{align*}
\frac{dV_1}{dt}
&=(R+C)(\dot{R}+\dot{C})+k_1R\dot{R}+k_2\dot{I}+k_3C\dot{C} \\
&=-k_2\beta I^2-(\gamma+k_1(\gamma+\delta))R^2
 -(\gamma+\eta+k_3(\gamma+\eta))C^2+k_2(\beta-(\gamma+\alpha))I \\
&\quad+(-\mu+k_1(1-\mu))\beta CIR+(\alpha+k_1\alpha-k_2\beta)IR
 +(-(2\gamma+\eta)+k_3\delta)RC\\
&\quad-(\mu+k_3)\beta C^2I+(\alpha+k_2(\mu-1)\beta)CI.
\end{align*}
Therefore,
\begin{align*}
\frac{dV_1}{dt}
&=-\frac{\alpha}{1-\mu} I^2-(\gamma+\frac{\mu(\gamma+\delta)}{1-\mu}) R^2
 -(\gamma+\eta)(1+\frac{2\gamma+\eta}{\delta}) C^2\\
&\quad+\frac{\alpha}{(1-\mu)\beta}(\beta-(\gamma+\alpha)) I
 -(\mu+\frac{2\gamma+\eta}{\delta})\beta C^2I.
\end{align*}
Since $R_0\leq 1$, we obtain
$$
\frac{dV_1}{dt}\leq 0.
$$
Obviously, \eqref{eq41} implies that the largest invariant set in the set of 
$\frac{dV_1}{dt}=0$ is $(0,0,0)$. By LaSalle's invariant
principle \cite{LaSalle,Lyapunov}, we conclude that model \eqref{eq41} 
is global asymptotically stable at $(0,0,0)$.
\end{proof}

\begin{remark}\label{rm43} \rm
From the perspective of epidemiological significance, Theorem \ref{thm42} 
demonstrates that the disease ultimately dies out regardless of
the initial values of  model \eqref{eq} when $R_0\leq1$.
\end{remark}

Take $\mu=1$, then  \eqref{eq} can reduce to the classical SIRS model of the 
form \eqref{SIRS}. Hence, we can deduce that the disease-free equilibrium 
$G_0=(1,0,0)$ of \eqref{SIRS} is globally asymptotically  stable
when $R_0\leq 1$. Namely, we have the following result.


\begin{corollary}\label{cor41}
The disease-free equilibrium for  model \eqref{SIRS} is globally asymptotically  
stable when $R_0\leq$ 1.
\end{corollary}

In fact, it is straightforward to show that the disease-free equilibrium 
$G_0=(1,0,0)$ is locally asymptotically  stable.
We can also illustrate that $G_0$ is  globally asymptotically  stable by
considering the following Lyapunov function
$$
V_2=\frac{1}{2}(S-1+I+R)^2 + \frac{2\gamma}{\beta} I+\frac{\gamma}{\alpha} R^2.
$$

\section{Persistence}

In the section, we will establish the persistence theorem of disease when $R_0>1$.
Persistence implies that the infected individuals will persist in the future. 
In this paper, we first use the definition given by Butler and Waltman \cite{Butler}. 
That is,  model \eqref{eq} is said to be uniformly persistent if there 
exists a positive number $b$ such that
\begin{equation}\label{eq51}
\min\Big\{\liminf_{t\to\infty}(S(t)),\liminf_{t\to\infty}(I(t)),
\liminf_{t\to\infty}(R(t)),\liminf_{t\to\infty}(C(t))\Big\}=b
\end{equation}
for every trajectory with positive initial conditions.

For a region $E$, denote by $\partial E$ and $\mathring{E}$ the boundary and 
the interior of $E$, respectively.
Denote by $\partial \mathcal{F}$ the restriction of the flow $\mathcal{F}$ to 
$\partial E$ and note that $\partial E$ is, in general, not positively invariant.
Let $N$ be the maximal invariant set of $\partial E$. Suppose $N$ is a closed 
invariant set  and there exists a cover $\{N_{\alpha}\}_{\alpha\in A} $
of $N$, where $A$ is a nonempty index set. $N_{\alpha}\subset\partial E$, 
$N\subset \cup_{\alpha\in A}N_{\alpha}$ and $N_{\alpha} (\alpha\in A)$ 
are pairwise disjoint closed invariant sets. Furthermore, we propose the 
following hypotheses:
\begin{itemize}
  \item[(H1)] All $N_{\alpha}$ are isolated invariant sets of the flow 
$\mathcal{F}$;

  \item[(H2)] $N_{\alpha}$ $(\alpha\in A)$ is acyclic, that is, any finite subset 
of $N_{\alpha}$ $(\alpha\in A)$ does not form a cycle;

  \item[(H3)] Any compact subset of $\partial E$ contains, at most,
 finitely many sets of $N_{\alpha}$ $(\alpha\in A)$.
\end{itemize}
The following lemma plays an important role in analyzing the uniformly persistence.

\begin{lemma}[\cite{ref29}]\label{lem51} 
Let $E$ be a closed positively invariant subset of $X$ on which a continuous 
flow $\mathcal{F}$ is defined. Suppose there is a constant
$\alpha>0$ such that $\mathcal{F}$ is point dissipative on 
$S[\partial E,\alpha]\cap \mathring{E}$ and the assumptions {\rm (H1)--(H3)} hold.
 Then the flow $\mathcal{F}$
is uniformly persistent if and only if 
$W^{+}(N_{\alpha})\cap S[\partial E,\alpha]\cap \mathring{E}=\emptyset$ for all 
$\alpha\in A $,
where $W^{+}(N_{\alpha})=\{y\in X|\,\Lambda^{+}(y)\subset N_{\alpha}\}$.
\end{lemma}

\begin{theorem}\label{thm52}
Model \eqref{eq} is uniformly persistent in $\mathring{\mathbb{D}}$ if and only 
if $R_0>1$.
\end{theorem}

\begin{proof}
It is easy to prove the necessity by means of Theorems \ref{thm41} and 
\ref{thm42} because the asymptotical stability of
$E_0$ excludes all kinds of persistence. Now, we prove the sufficiency
of this theorem by using Lemma \ref{lem51}. 
Choose $X=\mathbb{R}^{4}$ and $E=\mathbb{D}$. 
We only need to prove that model \eqref{eq}
satisfies all the conditions of Lemma \ref{lem51}. Note that the maximal 
invariant set on the boundary $\partial \mathbb{D}$ only
contains a point $E_0$ which is isolated. Then the assumptions (H1)--(H3)
are satisfied. By Lemma \ref{lem51}, we observe that the
uniform persistence of model \eqref{eq} is equivalent to the instability of 
the disease-free equilibrium. Thereby, the proof is complete.
\end{proof}

\begin{remark}\label{rm54} \rm
It follows from Theorem \ref{thm52} that the uniform persistence of 
model \eqref{eq} in the bounded set $\mathring{\mathbb{D}}$ is equivalent to the
existence of a compact attractor $K\subset \mathring{\mathbb{D}}$.
\end{remark}

\section{Stability of the endemic equilibrium}

This section concerns the stability of the
endemic equilibrium $E^*(S^*,I^*,R^*,C^*)$ when $R_0> 1$.

\begin{theorem}\label{thm61}
$E^*(S^*,I^*,R^*,C^*)$ is locally asymptotically stable when $R_0> 1$.
\end{theorem}

\begin{proof}
Let $N(t)=S(t)+I(t)+R(t)+C(t)$. It is easy to see that 
$\dot{N}(t)=\gamma-\gamma N(t)$. By a change of variables, we see that
model \eqref{eq} is equivalent to the  model
\begin{equation}\label{eq61}
\begin{gathered}
\dot{N}=\gamma-\gamma N ,\\
\dot{S}=\gamma(1-S)-\beta SI+\eta C,\\
\dot{I}=\beta SI+\mu \beta CI-(\gamma + \alpha)I,\\
\dot{C}=\delta (N-S-I-C) - \beta CI-(\gamma + \eta)C\,.
\end{gathered}
\end{equation}
By Theorem \ref{thm31} model \eqref{eq61} has a unique endemic 
equilibrium  $(N^*,S^*,I^*,C^*)$,
where $N^*=S^*+I^*+R^*+C^*$. Now, it suffices to verify that
$(N^*,S^*,I^*,C^*)$ is locally asymptotically stable  to show
that $E^*(S^*,I^*,R^*,C^*)$  of
model \eqref{eq} is locally asymptotically stable.

The Jacobian matrix of model \eqref{eq61} at $(N^*,S^*,I^*,C^*)$ is
$$
\begin{pmatrix}
-\gamma & 0 & 0 & 0 \\
0 & -\gamma-\beta I^* & -\beta S^* & \eta \\
0 & \beta I^* & \beta S^*+\mu \beta C^*-(\gamma+\alpha) & \mu \beta I^* \\
\delta & -\delta & -\delta-\beta C^* & -\delta-\beta I^*-(\gamma+\eta)
\end{pmatrix}.
$$
Since $\beta S^*+\mu \beta C^*-(\gamma+\alpha)=0$, its characteristic
equation is 
$$
(\lambda+\gamma)(d_0\lambda^{3}+d_1\lambda^2+d_2\lambda+d_3)=0,$$
where
\begin{gather*}
d_0=1,\quad
d_1=\delta+\beta I^*+\gamma+\eta+\gamma+\beta I^*,\\
d_2=\mu\beta I^*(\delta+\beta C^*)+(\gamma+\beta I^*)
 (\delta+\beta I^* +\gamma+\eta)+\beta I^*\beta S^*+\delta\eta,\\
\begin{aligned}
d_3&=(\gamma+\beta I^*)\mu\beta  I^*(\delta+\beta C^*)
 +\beta I^*\beta S^*(\beta I^* +\gamma+\eta )
+\beta I^*\eta(\delta+\beta C^* ) \\
&\quad +\delta(1-\mu)\beta I^* \beta S^*.
\end{aligned}
\end{gather*}
Therefore, the conclusion of this theorem is verified if the real parts of 
all solutions of $d_0\lambda^{3}+d_1\lambda^2+d_2\lambda+d_3=0$
are negative.
From above, we have $d_0>0$, $d_1>0$, $d_2>0$, $d_3>0$, and
\begin{align*}
d_1d_2-d_0d_3
=&(\gamma+\beta I^*)(\delta+\eta+\gamma+\beta I^*)
 (\delta+\eta+\gamma+\beta I^*+\gamma+\beta I^*)\\
&+\mu\beta I^*(\delta+\beta C^*)(\delta+\eta+\gamma+\beta I^*)
 +\delta\eta(\delta+\gamma+\eta+\gamma+\beta I^*)\\
&+\delta\mu\beta I^*\beta S^*+\beta I^* \beta S^*(\gamma+\beta I^*)
 -\beta I^* \beta C^* \eta.
\end{align*}
It follows from  the second of model \eqref{eq61} that
$$ 
\eta\beta C^*=-\gamma\beta+\gamma\beta S^*+\beta S^* \beta I^*.
$$
We get
$$
\beta I^* \beta S^*(\gamma+\beta I^*)-\beta I^* \beta C^* \eta
=\gamma\beta\beta I^*,
$$
namely, $d_1d_2-d_0d_3>0. $ By applying the Routh-Hurwitz criterion,
we can verify that $(N^*,S^*,I^*,C^*)$ is locally asymptotically stable.
\end{proof}

Next, we consider the global asymptotical stability of the endemic 
equilibrium $E^*$. The routine technique of the global asymptotical stability 
of endemic equilibrium is based on Lyapunov function and the 
Poincar\'e-Bendixson trichotomy. Here we will utilize another method, which 
is developed by Li and Muldowney \cite{Li,Li2}.

Let $x\mapsto f(x)\in \mathbb{R}^{n}$ be a $C^{1}$ function for $x$ in an 
open set $D\subset \mathbb{R}^{n}$. Consider the differential equation
\begin{equation}\label{eq62}
\dot{x}=f(x).
\end{equation}
Denote by $x(t,x_0)$ the solution to \eqref{eq62} such that $x(t,x_0)=x_0$,
and introduce the following two assumptions:
\begin{itemize}
  \item[(H3)] There exists a compact absorbing set $K\subset D$;
  \item[(H4)] Equation \eqref{eq62} has a unique equilibrium $\overline{x}$ in $D$.
\end{itemize}

Let $A$ ba an $n\times n$ matrix, $A^{[2]}$ is called the second additive 
compound matrix of $A$, which is an $(_2^{n})\times(_2^{n})$ matrix.
For instance, when $n=3$,
$$
f(x)=\begin{pmatrix}
a_{11}+a_{22} & a_{23} & -a_{13} \\
a_{32} & a_{11}+a_{33} & a_{12}  \\
-a_{31} & a_{21} & a_{22}+a_{33}
\end{pmatrix}.
$$
For the detailed discussions of compound matrix and their properties we refer 
the reader to \cite{Fiedler,Muldowney}.
Let $ x\mapsto P(x)\ $ be a $(_2^{n})\times(_2^{n})$ matrix-value function
that is $C^{1}$ for $x\in D$. Assume that
$P^{-1}(x)$ exists and is continuous in $ x\in K$, where $K$ is the compact 
absorbing set. A quantity $\overline{q}_2$ is defined as
$$
\overline{q}_2={\limsup_{t\to\infty}\sup_{x_0\in K}
\frac{1}{t}\int_0^{t}\mu(B(x(s,x_0)))ds},
$$
where
$$ 
B=P_{f}P^{-1}+P\frac{\partial f}{\partial x}^{[2]}P^{-1},
$$
and the matrix $P_{f}$ is obtained by replacing each entry $p_{ij}$ of $P$ 
by its derivative in the direction of $f$, $p_{ijf}$.
The quantity $\mu(B)$ is Lozinski\u{i} measure of $B$ with respect to a 
vector norm $|\cdot|$ in $\mathbb{R}^{N}$, $N=(_2^{n})\times(_2^{n})$,
defined by
$$
\mu(B)=\lim_{h\to 0+}\frac{|I+hB|-1}{h},
$$
see \cite{Martin}. The following global stability result is 
\cite[Theorem 3.5]{Li}.

\begin{lemma}\label{lem62}
Assume that $\mathbb{D}$ is simple connected and that assumptions {\rm (H3)}
 and {\rm (H4)} hold.  Then the unique equilibrium $\bar{x}$ of \eqref{eq62}
is global stable in $\mathbb{D}$ if  $\overline{q}_2<0$.
\end{lemma}

\begin{theorem}\label{thm63}
If $R_0>1$ then the endemic equilibrium $E^*$ of  \eqref{eq} is
globally asymptotically stable when
$\eta<\gamma$ and $\delta-\eta+\alpha+\beta <\gamma$.
\end{theorem}

\begin{proof}
From the discussions of Theorem \ref{thm31} and Remark \ref{rm54}, 
we conclude that model \eqref{eq} satisfies the assumptions (H3) and (H4) 
in $\mathbb{D}$.
From the proof of Theorem \ref{thm23}, we have $\lim_{t\to\infty}N(t)$=1.
Thus, we just need to consider the following limiting equation of model \eqref{eq}:
\begin{equation}\label{eq63}
\begin{gathered}
\dot{S}=\gamma(1-S)-\beta SI+\eta C,\\
\dot{I}=\beta SI+\mu \beta CI-(\gamma + \alpha)I,\\
\dot{C}=\delta (1-S-I-C)-\beta CI-(\gamma + \eta)C.
\end{gathered}
\end{equation}
Let $f=(f_1,f_2,f_3)^{T}$, where $f_1, f_2$ and $f_3$ represent
the right-hand sides of model \eqref{eq63}, respectively.
Furthermore, let $x=(S,I,C)^{T}$, then the Jacobian matrix associated 
with a general solution $x(t)$ of model \eqref{eq63} is
$$
\frac{\partial f}{\partial x}
=\begin{pmatrix}
-\gamma-\beta I & -\beta S &\eta\\
\beta I & \beta S+\mu\beta C-(\gamma+\alpha) & \mu\beta I\\
-\delta & -\delta-\beta C & -\delta-\beta I-(\gamma+\eta)
\end{pmatrix}.
$$
The second additive compound matrix of $\frac{\partial f}{\partial x}$ is
$$
\frac{\partial f}{\partial x}^{[2]}=
\begin{pmatrix}
g_{11} & \mu\beta I & -\eta \\
-\delta-\beta C  & g_{22} & -\beta S \\
\delta & \beta I & g_{33}
\end{pmatrix},
$$
where
\begin{gather*}
g_{11}=-\gamma-\beta I + \beta S+\mu\beta C-(\gamma+\alpha),\\
g_{22}=-\gamma-\beta I -\delta - \beta I -(\gamma+\eta),\\
g_{33}=\beta S+\mu\beta C-(\gamma+\alpha)-\delta - \beta I -(\gamma+\eta).
\end{gather*}
Set the function $P(x)=P(S,I,C)
=\operatorname{diag}(\frac{S}{I},\frac{S}{I},\frac{S}{I})$, then we have
$$
P_{f}P^{-1}=\operatorname{diag} (\frac{\dot{S}}{S}-\frac{\dot{I}}{I},
\frac{\dot{S}}{S}-\frac{\dot{I}}{I},\frac{\dot{S}}{S}-\frac{\dot{I}}{I}).
$$
Therefore,
$$
B=P_{f}P^{-1}+P\frac{\partial f}{\partial x}^{[2]}P^{-1}
=
\begin{pmatrix}
B_{11} & B_{12}\\
B_{21} & B_{22}
\end{pmatrix},
$$
where
\begin{gather*}
B_{11}=g_{11}+\frac{\dot{S}}{S}-\frac{\dot{I}}{I},\quad
B_{12}=(\mu\beta I,-\eta), \\ 
 B_{21}=(-\delta-\beta C,\delta)^{T},\quad 
B_{22}=\begin{pmatrix}
g_{22}+\frac{\dot{S}}{S}-\frac{\dot{I}}{I} & -\beta S\\
\beta I & g_{33}+\frac{\dot{S}}{S}-\frac{\dot{I}}{I}
\end{pmatrix}.
\end{gather*}
Let $(u,v,w)$ denote the vector in $\mathbb{R}^{3}$, we consider the 
following norm in $\mathbb{R}^{3}$,
$$
|(u,v,w)|=\max\{|u|,|v|+|w|\},
$$
and let $\mu_1$ denote the Lozinski$\breve{\mathrm{i}}$ measure with respect to this norm. Using the method of estimating $\mu_1$
in \cite{Martin}, we have
$$
\mu_1(B)\leq \sup\{g_1,g_2\},
$$
where
$g_1=\mu_1(B_{11})+|B_{12}|$, $g_2=\mu_1(B_{22})+|B_{21}|$,
and $|B_{12}|,|B_{21}|$ are matrix norms with respect to the $L^{1}$ 
vector norm and $\mu_1$ denote
the Lozinski\v{i} measure with respect to the $L^{1}$ norm.
From  \eqref{eq61}, it implies that
$$
\frac{\dot{I}}{I}=\beta S+\mu \beta C-(\gamma + \alpha).
$$
Therefore,
\begin{gather*}
\mu_1(B_{11})=\frac{\dot{S}}{S}-\gamma-\beta I,\quad
|B_{12}|=\max\{\mu\beta I,\eta\},\quad  |B_{21}|=2 \delta+\beta C, \\
\mu_1(B_{22})
=\frac{\dot{S}}{S}-\gamma+\max\{-\delta-\beta I- \eta
 - \beta S-\mu\beta C+\alpha,\beta S- \delta-\beta I-\eta\}.
\end{gather*}
Thus, we have
\begin{gather*}
g_1=\frac{\dot{S}}{S}-\gamma+\max\{(\mu-1)\beta I,-\beta I+\eta\}
<\frac{\dot{S}}{S}-\gamma+\eta,\\
\begin{aligned}
g_2&=\frac{\dot{S}}{S}-\gamma+\max\{\delta-\beta I-\eta-\beta S
+(1- \mu)\beta C+\alpha,\beta S+\delta-\beta I-\eta+\beta C\}\\
&<\frac{\dot{S}}{S}-\gamma +\delta-\eta +\max\{\alpha+\beta C, \beta S+\beta C\}\\
&<\frac{\dot{S}}{S}-\gamma+\delta-\eta+\alpha+\beta .
\end{aligned}
\end{gather*}
This leads to
\[
\mu_1(B)\leq \frac{\dot{S}}{S}-\gamma +\max\{\eta, \delta-\eta+\alpha+\beta \},
\]
where $\omega=\max\{\eta, \delta-\eta+\alpha+\beta\}<\gamma$.
Consequently,
$$
\frac{1}{t}\int_0^{t}\mu_1(B)ds
\leq\frac{1}{t}\log\frac{S(t)}{S(0)}-(\gamma-\omega),
$$
which yields $\overline{q_2}<0$. The proof is complete.
\end{proof}

\begin{remark}\label{rm64} \rm
Theorems \ref{thm52} and \ref{thm63} describe that the disease always persists 
and becomes endemic at an endemic level, no matter how small size the 
initial value of infections has. To eradicate the disease, what we need 
to do is to reduce the key of threshold value $R_0$ to below 1.
\end{remark}

In view of Theorem \ref{thm63}, we can obtain the global asymptotical 
stability of endemic equilibrium $E^*$ of model \eqref{eq} when $R_0>1$, 
$\eta<\gamma$ and $\delta-\eta+\alpha+\beta <\gamma$. 
In the following theorem, we shall see that the constraint that $\eta<\gamma$ 
and $\delta-\eta+\alpha+\beta <\gamma$ is not necessary if we consider the 
global asymptotical stability of endemic equilibrium of model \eqref{SIRS}, 
which can be regarded as a special case of model \eqref{eq}.

\begin{theorem}\label{thm65}
If $R_0>1$ then the unique endemic equilibrium of \eqref{SIRS} is globally 
asymptotically  stable.
\end{theorem}

\begin{proof}
It is easy to see that if $R_0>1$ then  model \eqref{SIRS} has a unique
endemic equilibrium $G^*=(S_1^*,I_1^*,R_1^*)$, where
$$
S_1^*=\frac{\gamma+\alpha}{\beta},\quad
I_1^*=\frac{(\gamma+\delta)(\beta-\gamma-\alpha)}{(\gamma+\delta+\alpha)\beta},
\quad
R_1^*=\frac{\alpha (\beta-\gamma-\alpha)}{(\gamma+\delta+\alpha)\beta}.
$$
Meanwhile, we further see that $G^*$ is locally asymptotically  stable.
Next, we show that $G^*$ is globally asymptotically stable by
considering the following Lyapunov function
$$
V_3=\frac{1}{2}(S-S_1^*+I-I_1^*+R-R_1^*)^2
+k_4(I-I_1^*+I_1^*
\log\frac{I}{I_1^*})+\frac{k_5}{2}(S-S_1^*+I-I_1^*)^2,
$$
where $k_4$ and $k_5$ are positive constants to be determined later. The time derivative of $V_3$
along the solutions of  \eqref{SIRS} is
$$\frac{dV_3}{dt}=(S-S_1^*+I-I_1^*+R-R_1^*)(\dot{S}+\dot{I}+\dot{R})+k_4\frac{I-I_1^*}{I}\dot{I}+k_5(S-S_1^*+I-I_1^*)
(\dot{S}+\dot{I}).$$
Note that
\begin{gather*}
\gamma(1-S_1^*)-\beta S_1^*I_1^*+\delta R_1^*=0,\\
\beta S_1^*I_1^*-(\gamma + \alpha)I_1^*=0,\\
\alpha I_1^*-(\gamma + \delta)R_1^*=0.
\end{gather*}
Then, we have
\begin{align*}
\frac{dV_3}{dt}
&=(S-S_1^*+I-I_1^*+R-R_1^*)(-\gamma)(S-S_1^*+I-I_1^*+R-R_1^*)\\
&\quad +k_4(I-I_1^*)\beta(S-S_1^*)+k_5(S-S_1^*+I-I_1^*)(-\gamma(S-S_1^*)\\
&\quad +\delta(R-R_1^*)-(\gamma+\alpha)(I-I_1^*))\\
&=-\gamma(1+k_5)(S-S_1^*)^2-(\gamma+k_5(\gamma+\alpha))(I-I_1^*)^2-\gamma(R-R_1^*)^2\\
&\quad +(-2\gamma+k_4\beta-k_5(2\gamma+\alpha))(S-S_1^*)(I-I_1^*)\\
&\quad +(-2\gamma+k_5\delta)(S-S_1^*)(R-R_1^*)
  +(-2\gamma+k_5\delta)(I-I_1^*)(R-R_1^*).
\end{align*}
Take $k_4=2\gamma(\delta+2\gamma+\alpha)/(\delta\beta)$ and 
$k_5=2\gamma/\delta$, then we have
\[
\frac{dV_3}{dt}=-\gamma(1+\frac{2\gamma}{\delta})(S-S_1^*)^2
-(\gamma+\frac{2\gamma}{\delta}(\gamma+\alpha))(I-I_1^*)^2-\gamma(R-R_1^*)^2\leq0.
\]
Therefore, the LaSalle's invariant principle \cite{LaSalle,Lyapunov} 
implies that $G^*$ is globally asymptotically  stable.
\end{proof}

\section{Numerical simulation}

In this section, we aim to provide a numerical simulation to substantiate the 
theoretical results established in
the previous sections by using the Runge-Kutta fourth order iterative method.
Consider  model \eqref{eq} with the parameters given as follows:
 $\gamma=0.3$, $\eta=0.2$, $\mu=0.05$, $\alpha=0.5$, $\delta=0.5$, $\beta=0.4$.

In Figure~\ref{SIRC1}, we set up two sets of initial values. One case is that 
$S(0)=0.3, I(0)=0.5, R(0)=0, C(0)=0$. Another case is that 
$S(0)=20, I(0)=20, R(0)=0, C(0)=0$.
It is observed in Figure~\ref{SIRC1} that all trajectories of model \eqref{eq} 
eventually stay in the positively invariant set $\mathbb{D}$ regardless of 
whether or not the initial values are in $\mathbb{D}$ and that we can obtain 
the pivotal threshold value $R_0=0.5$ for the choice of parameters.
 In this case, it follows from Theorem \ref{thm31} that  \eqref{eq} has a 
unique equilibrium $E_0=(1, 0, 0, 0)$. Theorem \ref{thm42} means that this 
disease-free equilibrium is globally asymptotically stable.

\begin{figure}[ht]
  \begin{center}
\includegraphics[width=\textwidth]{fig1} % SIRC1
\end{center}
  \caption{Solution of \eqref{eq} with 
$\gamma=0.3$, $\eta= 0.2 $, $\mu= 0.05$, $\alpha=0.5$, $\delta=0.5$, 
$\beta= 0.4$, 
where the initial value is (a): 
$S(0)=0.3$, $I(0)=0.5$, $R(0)=0$, $C(0)=0$, and 
(b): $S(0)=20$, $I(0)=20$, $R(0)=0$.}
  \label{SIRC1}
\end{figure}

Figure~\ref{SIRC1} shows that the infected individuals are eventually eradicated 
from the crowd, while the susceptible individuals will ultimately approach the 
maximum value.
The epidemiological implication of Figure~\ref{SIRC2}(a) is that the infected
 population vanish over time. In other words, the disease will die out in the 
long time.
If we change the value of $\beta$ into 0.98, then $R_0=1.225$ and hence $E_0$ 
is unstable (see Figure~\ref{SIRC2}(b)).

\begin{figure}[ht]
  \begin{center}
\includegraphics[width=\textwidth]{fig2} % SIRC2
\end{center}
  \caption{Solution of \eqref{eq} with $S(0)=0.3$, $I(0)=0.5$, $R(0)=0$, $C(0)=0$
and 
$\gamma=0.3$, $\eta=0.2$, $\mu=0.05$, $\alpha=0.5$, $\delta=0.5$, where 
(a) $\beta=0.4$ and  (b) $\beta=0.98$.}
  \label{SIRC2}
\end{figure}

If $(\gamma,\eta,\mu,\alpha,\delta,\beta)=(0.3, 0.2, 0.05, 0.5, 0.5,0.9)$, 
then there exists a unique endemic equilibrium $E^*$ which is locally asymptotically 
stable (see Theorem \ref{thm61}). Figure~\ref{SIRC3} implies that the number 
of infected individuals persist and gradually tend to a positive constant when 
$R_0=1.125>1$. If we take the parameters $\gamma=0.4$, $\eta=0.35$, $\mu=0.05$, 
$\alpha=0.05$, $\delta=0.05$, $\beta=0.5$, then all the conditions of
Theorem \ref{thm63} are satisfied, and hence that $E^*$ is globally asymptotically
 stable. The epidemiological implication of Figure~\ref{SIRC4} is that the 
infected individuals always exist if its initial value is non-negative.

\begin{figure}[ht]
  \begin{center}
\includegraphics[width=\textwidth]{fig3} % SIRC3
\end{center}
  \caption{Solution of \eqref{eq} with $S(0)=0.3$, $I(0)=0.5$, $R(0)=0$,
$C(0)=0$ and 
$\gamma=0.3$, $\eta=0.2$, $\mu=0.05$, $\alpha=0.5$,
$\delta=0.5$, $\beta=0.9$.}
  \label{SIRC3}
\end{figure}


\begin{figure}[ht]
  \begin{center}
\includegraphics[width=\textwidth]{fig4} % SIRC4
\end{center}
  \caption{Solution of \eqref{eq} with 
$S(0)=0.3$, $I(0)=0.5$, $R(0)=0.$, $C(0)=0$ and 
$\gamma=0.4$, $\eta=0.35$, $\mu=0.05$, 
$\alpha=0.05$, $\delta=0.05$, $\beta=0.5$.}
  \label{SIRC4}
\end{figure}


In what follows, we carry out sensitivity analysises for  \eqref{eq} by the change 
of the recovery rate $\alpha$ as well as the contact rate $\beta$.
In Figure~\ref{SIRC5}, we study the effect of parameter $\alpha$ on  model
 \eqref{eq}. We see that  parameter $\alpha$ is directly proportional with the
 number of susceptible, recovered, cross-immune individuals. However, it is 
inversely proportional with the number of infectious individuals, and infectious 
individuals finally eradicate. It follows from Figure~\ref{SIRC5} that
 the threshold value $R_0$ reduces to be less than 1 by increasing $\alpha$ so 
that the endemic equilibrium vanishes.
Therefore, model \eqref{eq} just has a disease-free equilibrium.  
This  means that the larger the parameter $\alpha$
is, the smaller the basic reproduction number $R_0$ is, and hence the faster 
the disease die out. From a biological perspective, we should reduce the value of
$R_0$ as possible as we can in order that the disease dies out quickly.

Finally, we examine the influence of the contact rate $\beta$. 
In Figure~\ref{SIRC6}, we observe that  the number of infectious, recovered, 
cross-immune individuals are directly proportional with the parameter $\beta$, 
but the number of susceptible individuals are inversely proportional with the 
parameter $\beta$, and finally approach to 1. Why these phenomena happened is 
that the disease-free equilibrium becomes unstable and an endemic equilibrium 
appears in \eqref{eq} when the threshold value $R_0$ increases and passes 
through 1. This implies that the larger the parameter $\beta$
is, the larger the value of $R_0$ is, and then the higher the endemic level 
will be. As a result, more and more population contacted with infected individuals
 will make the disease persist at an endemic level.

\begin{figure}[ht]
  \begin{center}
\includegraphics[width=\textwidth]{fig5} % SIRC5
\end{center}
  \caption{Sensitivity of  model \eqref{eq} for different values of $\alpha$.}
  \label{SIRC5}
\end{figure}

\begin{figure}[ht]
  \begin{center}
\includegraphics[width=\textwidth]{fig6} % SIRC6
\end{center}
  \caption{Sensitivity of  model \eqref{eq} for different values of $\beta$.}
  \label{SIRC6}
\end{figure}


\subsection*{Conclusions}
This paper presents a mathematical study on the dynamics of an SIRC epidemiological 
model established by Casagrandi et al.\ \cite{Renato}.
The basic reproduction number $R_0$ plays a vital role in determining the global 
dynamics of \eqref{eq}.
It is noted that the model always has a
disease-free equilibrium, which is globally asymptotically stable when 
$R_0\leq1$. When $R_0>1$,
we apply the Routh-Hurwitz criterion to prove that the model has a unique endemic 
equilibrium, which is locally asymptotically stable. In this case, the
disease-free equilibrium become unstable. Based on Li-Muldowney's 
global-stability criterion \cite{Li}, we show that the unique endemic equilibrium 
can be globally asymptotically stable in a feasible region, i.e., 
influenza A becomes endemic. Although we have established the global stability 
of the unique endemic equilibrium $E^*$ when $R_0>1$, our results are obtained 
under the assumptions that $\eta<\gamma$ and $\delta-\eta+\alpha+\beta <\gamma$. 
From Theorem \ref{thm65}, we conjecture that the condition that $\eta<\gamma$ 
and $\delta-\eta+\alpha+\beta <\gamma$ is not necessary. Therefore,
the perspective of our work is to show the assumption that $R_0>1$ is a 
sufficient and necessary condition ensuring the globally asymptotical 
stability of the unique endemic equilibrium $E^*$.


\subsection*{Acknowledgments}
This work has been supported by the Natural Science Foundation of 
China (Grant No. 11671123).

\begin{thebibliography}{99}

\bibitem{Andreasen2} V. Andreasen, J. Lin, S. A. Levin; 
\emph{The dynamics of cocirculating influenza strains conferring partial 
cross-immunity}, Journal of Mathematical Biology, 35 (1997), 825--842.

\bibitem{Anderson} R. M. Anderson, R. M. Mary;
\emph{Infectious Diseases of Humans: Dynamics and Control},
 Oxford University Press, 1991.

\bibitem{Anderson1} R. M. Anderson, R. M. May;
\emph{Population biology of infectious diseases: Part I},
 Nature, 280 (1979), 361--367.

\bibitem{Butler} G. J. Butler, P. Waltman;
\emph{Persisitence in dynamical systems}, 
Journal of Differential Equations 63 (1986), 255--263.

\bibitem{Renato} R. Casagrandi, L. Bolzoni, S. A. Levin, V. Andreasen;
\emph{The SIRC model and influenza A}, Mathematical Biosciences 200 (2006), 152--169.

\bibitem{Chung} K. W. Chung, R. Lui;
\emph{Dynamics of two-strain influenza model with cross-immunity and no 
quarantine class}, Journal of Mathematical Biology, 73 (2016), 1467--1489.

\bibitem{Earn} D. J. Earn, J. Dushoff, S. A. Levin;
\emph{Ecology and evolution of the flu}, 
Trends in Ecology and Evolution, 17 (2002), 334--340.

\bibitem{Fiedler} M. Fiedler;
\emph{Additive compound matrices and an inequality for eigenvalues of symmetric 
stochastic matrices}, Czechoslovak Mathematical Journal, 99 (1974), 392--402.

\bibitem{ref29} H. I. Freedman, S. Ruan, M. Tang;
\emph{Uniform persistence and flows near a closed positively invariant set}, 
 Journal of Dynamics and Differential Equations, 6 (1994), 583--600.

\bibitem{Kermack} W. O. Kermack, A. G. McKendrick;
\emph{Contributions to the mathematical theory of epidemics: II. 
The problem of endemicity}, Proceedings of Royal Social London. Series A, 
Containing Paper of a Mathemeatical and Physiacal Character, 138 (1932), 55--83.

\bibitem{Kooi} B. W. Kooi, M. Aguiar, N. Stollenwerk;
\emph{Bifurcation analysis of a family of multi-strain epidemiology models},
Journal of Computational and Applied Mathematics, 252 (2013), 148--158.

\bibitem{LaSalle} J. P. LaSalle;
\emph{The stability of dynamical systems}, 
Society for Industrial and Applied Mathematics, 27 (1976), 1121--1130.

\bibitem{Lin} J. Lin, V. Andreasen, S. A. Levin;
\emph{Dynamics of influenza A drift: the linear three-strain model}, 
Mathematical Biosciences, 162 (1999), 33--51.

\bibitem{Li} M. Y. Li, J. S. Muldowney;
\emph{A geometric approach to global stability problems}, 
SIAM Journal on Mathematical Analysis, 27 (1996), 1070--1083.

\bibitem{Li2} M. Y. Li, J. S. Muldowney;
\emph{On R. A. Smith's autonomous convergence theorem}, 
The Rocky Mountain Journal of Mathematics,  25 (1995), 365--378.

\bibitem{Lyapunov} A. M. Lyapunov;
\emph{The general problem of the stability of motion}, 
International Journal of Control, 31 (1992), 531--534.

\bibitem{Martin} R. H. Martin, Jr.;
\emph{Logarithmic norms and projections applied to linear differential systems}, 
Journal of Mathematical Analysis and Applications, 45 (1974), 432--454.

\bibitem{Minayev} P. Minayev, N. Ferguson;
\emph{Improving the realism of deterministic multi-strain models: 
implications for modelling influenza A},
Journal of The Royal Society Interface, 6 (2009), 509--518.

\bibitem{Muldowney} J. S. Muldowney;
\emph{Compound matrices and ordinary differential equations}, 
Rocky Mountain Journal of Mathematics, 20 (1990),  857--872.

\bibitem{Nuno1} M. Nu\~{n}o, Z. Feng, M. Martcheva, C. Castillo-Chavez;
\emph{Dynamics of two-strain influenza with isolation  and partial cross-immunity}. 
SIAM Journal on Applied Mathematics, 65 (2005), 964--982.

\bibitem{Nuno}M. Nu\~{n}o, M. Martcheve, C. Castillo-Chavez;
\emph{Immune level approach for multiple strain pathogens},
             Journal of Biological Systems, 17 (2009), 713--737.

\bibitem{Smith} D. J. Smith, S. Forrest, D. H. Ackley, A. S. Perelson;
\emph{Variable efficacy of repeated annual influenza vaccination. Proceedings 
of the National Academy of Science}, 96 (1999) ,14001--14006.

\bibitem{Thomas} T. Szucs;
\emph{The socio-economic burden of influenza}, 
Journal of Antimicrobial Chemotherapy, 44 (1999), 11--15.

\bibitem{Taubenberger} J. K.Taubenberger, A. H. Reid, A. E. Krafft, 
K. E. Bijwaard, T. G. Fanning;
\emph{Initial genetic characterization of the 1918 Spanish influenza virus}, 
Science 275 (1997), 1793--1796.

\bibitem{ref21} X. Yang, L. Chen, J. Chen;
\emph{Permanence and positive periodic solution for the single-species 
nonautonomous delay diffusive models}, Computers and Mathematics with Applications, 
32 (1996) 109--116.

\bibitem{Zhong} X. Zhong, S. Guo, M. Peng;
\emph{Stability of stochastic SIRS epidemic models with saturated incidence 
rates and delay}, Stochastic Analysis and Applications, 35 (2017), 1--26.

\end{thebibliography}

\end{document}