Coronavirus disease 2019 (COVID-19) is a serious global health threat. On March 11, 2020, the World Health Organization (WHO) declared COVID-19 a pandemic, as it spreads worldwide rapidly following the logistic growth pattern. Reference Kamrujjaman, Mahmud and Islam1 More than 103 million cases were reported worldwide by January 31, 2021, where the United States of America (USA) has over 26.5 million cases alone. The virus that causes the COVID-19 disease is known as severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2). SARS-CoV-2 is spreading very quickly in the human population. The main route of spreading this virus is close person-to-person contact, and the spread is sustainable, as it goes from person-to-person without stopping. The ongoing COVID-19 pandemic is spreading more efficiently than influenza, but not as efficiently as highly contagious measles. 2 Another complicating factor is that asymptomatic patients are still able to spread the virus. According to the Centers for Disease Control and Prevention (CDC), common symptoms include but are not limited to fever or chills, cough, shortness of breath or difficulty breathing, sore throat, and many more.
As infected populations and the death tolls continue to rise rapidly, governments worldwide are trying to control the pandemic by reducing people’s close contact, such as shutting down public places, schools, colleges, universities, restaurants, playgrounds, and the list continues. Due to the lack of proper viral medicine or vaccine, travel bans from highly infected areas, social distancing, lock-down policy, isolation of an infected person, self-quarantine of exposed individual, mandatory use of face masks or face coverings, and strictly following all social-conscious and prevention strategies Reference Mahmud, Kamrujjaman and Jubyrea3 have been widely used strategies to contain the virus.
Infectious disease modeling is one of the most critical parts of understanding the current pandemic’s ongoing scenario. Mathematical modeling techniques, particularly the compartmental modeling technique, can help us estimate crucial parameters, such as disease transmission dynamics, number of infected people, number of hospitalized people, and number of recovered/dead individuals. Consequently, it helps to forecast outbreak timelines and the overall dynamics of the disease. The knowledge gained from mathematical modeling helps determine strategies to mitigate the outbreak and determine whether the interventions taken are useful or not.
The modeling community took the challenge to face the COVID-19 pandemic together. It started modeling to understand the dynamics of COVID-19 and help determine the intervention to slow the spread of the virus. Most of the models related to COVID-19 use population-based, Susceptible-Infected-Removal (SIR) models with differential equations or stochastic differential equations Reference Calafiore, Novara and Possieri4–Reference Kamrujjaman, Jubyrea and Islam7 or susceptible-Exposed-Infectious-Removal (SEIR) and extended SEIR compartmental models. Reference Chowdhury, Heng and Shawon9–Reference Kamrujjaman, Mahmud and Ahmed18 Some agent-based models developed related to COVID-19 use structured networks to connect individuals and model infection exchanges stochastically. Reference Biswas, Khaleque and Sen19–Reference Kamrujjaman, Mahmud and Islam23
Mathematical models often made redundant use of parameters and equations. Our primary focus is to keep the model as easy and straightforward as possible, so it will be understandable to the nonmath major population for a broader perspective. In this study, we develop a deterministic
$S{I_a}{I_s}R$
(modified SEIR) model to deepen our understanding of COVID-19’s dynamic. Theoretically, we analyze the model to show the existence and positivity invariance of the system’s solutions and determine the Disease-Free Equilibrium (DFE) and Endemic Equilibrium (EE). We calculate the basic reproduction number,
${{\cal R}_0}$
by using the next-generation matrix approach. We also analyze the local stability of the DFE and EE fixed points. Finally, we estimate the probability of an outbreak using the Continuous-time Markov Chain (CTMC) model for better prediction to control the disease outbreak.
We parameterize our model for Texas, United States, by using the data from March 6, 2020, to January 31, 2021, available at https://www.dshs.texas.gov/coronavirus/.
24
Using numerical simulations and data analyses, we predict the daily confirmed case and the cumulative case of Texas and compare the model prediction with the existing data.
Reference Murray25
We provide a sensitivity and elasticity index analysis of
${{\cal R}_0}$
to understand better the most fluctuating (sensitive) parameter of the model. Furthermore, in this study, we developed the CTMC model. We estimate the probability of a disease outbreak depending on the most critical parameter of the model to determine better control measures of the diseases. The main objectives of the study are: (a) We will work with real-life available discrete data of Texas to understand the cases and project the control of the infection; and (b) we expect that the prediction of controlling measure of COVID-19 will be able to validate the dynamics of the
$S{I_a}{I_s}R$
and CTMC models to obtain more accurate results.
We organized this article as follows: Mathematical Model discussed elaborately with positivity and boundedness of solution in the Mathematical Model and Existence of Solutions section. The fixed points, auxiliary results are described in the Determination of Fixed Points section. The local stability analysis, parameter estimation, and sensitivity analysis described in the Stability Analysis, Parameter Estimation, and Sensitivity section. The data analysis compared with the model solution with further prediction to control the epidemic, as a case study in Texas accomplished in the Numerical Simulation and Results section. The probability of disease outbreak with CTMC analysis presented in the Probability of a Disease Outbreak section. Finally, the Concluding Remarks section outlines the summary and concluding remarks of the results.
Mathematical Model and Existence of Solutions
The classical SIR model predicts the dynamics of infectious disease. We started with the following compartmental SIR model proposed in Murray Reference Murray25 :
$\left\{ {\matrix{ {S'(t) = \Lambda - \beta {{IS} \over N},} \cr {I'(t) = \beta {{IS} \over N} - \gamma I - \mu I,} \cr {R'(t) = \gamma I,} \cr } } \right.$
for
$t \in (0,\infty )$
with initial conditions
$S(0) = {S_0},\;\;I(0) = {I_0}\;\;{\rm{and}}\;R(0) = {R_0}.$
and for total population,
$N(t) = S(t) + I(t) + R(t).$
Here
$S(t),\;I(t),\;R(t)$
are the number of individuals in the susceptible, infected and removed compartments, respectively at time
$t$
with a day unit. The parameter
$\beta $
denotes the infection rate/disease transmission rate, and
$\gamma $
and
$\mu $
are the removal and disease induced mortality rate, respectively. The solution and detailed analysis of Equation (2.1) are available in Murray.
Reference Murray25
Compared with the SIR epidemic model, the next updated and advanced model is SEIR, which is biologically more feasible in many pandemics and infectious diseases. In this study, we consider the following 4 compartments
$S{I_a}{I_s}R$
mathematical model; a modified version of the typical SEIR model:
$\left\{ {\matrix{ {{{dS} \over {dt}} = \Lambda - ({\beta _1}{I_a} + {\beta _2}{I_s})S - {\mu _1}S,\quad t \ge 0,} \cr {{{d{I_a}} \over {dt}} = ({\beta _1}{I_a} + {\beta _2}{I_s})S - (\sigma + {\gamma _1} + {\mu _1}){I_a},\quad t \ge 0,} \cr {{{d{I_s}} \over {dt}} = \sigma {I_a} - ({\gamma _2} + {\mu _1} + {\mu _2}){I_s},\quad t \ge 0,} \cr {{{dR} \over {dt}} = {\gamma _1}{I_a} + {\gamma _2}{I_s} - {\mu _1}R,\quad t \ge 0} \cr } } \right.$
with initial conditions
$S(0) = {S_0},\;\;{I_a}(0) = {I_a}_0,\;\;{I_s}(0) = {I_s}_0\;\;{\rm{and}}\;\;R(0) = {R_0}\;,$
and
$N(t) \equiv S(t) + {I_a}(t) + {I_s}(t) + R(t)\;,$
Here,
$S(t),\;{I_a}(t),\;{I_s}(t)$
and
$R(t)$
are the number of individuals in the susceptible, asymptomatically infected, symptomatically infected (for simplicity, we will call the symptomatically infected population as an infected population), and removed compartments, respectively, at time
$t$
per day unit.
$\Lambda $
is the recruitment number in the susceptible compartment. Natural and disease induced deaths are denoted by
${\mu _1}$
and
${\mu _2}$
, respectively.
${\beta _1}$
and
${\beta _2}$
are the diseases transmission rates of susceptible individuals with asymptomatically infected and infected ones, respectively, which may cause the transmission of the infection,
$\sigma $
is the transition rate from asymptomatically infected to infected compartment and
${\gamma _i},\;i = 1,2$
are the recovery rates from asymptomatically infected and infected to removal compartment, respectively. After transmission, a susceptible individual initially becomes asymptomatically infected. As the disease progresses, asymptomatically infected individuals may develop symptoms, and transition from
${I_a}$
to
${I_s}$
at rate
$\sigma $
or they may never develop symptoms and recover at rate
${\gamma _1}$
. Because the model monitors dynamics of population, it follows that all its dependent variables and parameters, for example,
$\Lambda ,{\mu _i},{\beta _i},{\gamma _i}$
and
$\sigma $
must be non-negative along with
${\mu _1} \gt 0$
as in the model Equations (2.3)-(2.5). The definition of all parameters are elaborated in Table 1 and the flow diagram of the main model Equation (2.3) is shown in Figure 1.
Model parameters and their descriptions
Compartmental diagram for model.
Remark 1 In this study, we consider the modified SEIR model for SARS-CoV-2 dynamics to cover both infectious and exposed-asymptomatic individuals compared to the classical SIR and SEIR compartment models. The model has a combined compartment for exposed and asymptomatic classes, which reduces the number of parameters and eventually reduces the model’s complexity.
The following result ensures the existence and positivity of solutions of Equation (2.3).
Theorem 1
The closed region
${\rm \Omega} = \{ (s,{I_a},{I_s},R)\varepsilon \,{R}_ + ^4:0 \lt N \le {\Lambda \over {{\mu _1}}}\} $
is positively invariant set for the system in Equation (2.3).
The proof of theorem 1 can be found in the Appendix.
Determination of Fixed Points
To find the equilibrium points
$(\widetilde S,\widetilde {{I_a}},\widetilde {{I_s}},\widetilde R)$
of the system Equation (2.3), we set the derivatives equal to zero. So, at equilibrium states, we get
$\left\{ {\matrix{ {\Lambda - ({\beta _1}\widetilde{I_a} + {\beta _2}\widetilde{I_s}) - {\mu _1} = 0,} \cr {({\beta _1}\widetilde{I_a} + {\beta _2}\widetilde{I_s})\widetilde S - (\sigma + {\gamma _1} + {\mu _1})\widetilde{I_a} = 0,} \cr {\sigma \widetilde{I_a} - ({\gamma _2} + {\mu _1} + {\mu _2})\widetilde{I_s} = 0,} \cr {{\gamma _1}\widetilde{I_a} + {\gamma _2}\widetilde{I_s} - {\mu _1} = 0.} \cr } } \right.$
DFE Point
For the DFE, we replace the variables as
$(\widetilde S,\widetilde {{I_a}},\widetilde {{I_s}},\widetilde R) \equiv ({S_0},{I_{a0}},{I_{s0}},{R_0}).$
This gives,
$\left\{ {\matrix{ {\Lambda - ({\beta _1}{I_{s0}} + {\beta _2}{I_{s0}}){S_0} - {\mu _1}{S_0} = 0,} \cr {({\beta _1}{I_{a0}} + {\beta _2}{I_{s0}}){S_0} - (\sigma + {\gamma _1} + {\mu _1}){I_{a0}} = 0,} \cr {\sigma {I_{a0}} - ({\gamma _2} + {\mu _1} + {\mu _2}){I_{s0}} = 0,} \cr {{\gamma _1}{I_{a0}} + {\gamma _2}{I_{s0}} - {\mu _1}{R_0} = 0.} \cr } } \right.$
Therefore, the DFE point can easily be found as
$$\left( {{S_0},{I_{a0}},{I_{s0}},{R_0}} \right) \equiv {\rm{ }}\left( {{\Lambda \over {{\mu _1}}},0,0,0} \right).$$
EE Point
For the EE, we replace the variables as
$(\widetilde S,\widetilde {{I_a}},\widetilde {{I_s}},\widetilde R) \equiv ({S^*},I_a^*,I_s^*,{R^*}),$
where,
$I_s^* \gt 0$
. And we have the following system
$\left\{ {\matrix{ {\Lambda - ({\beta _1}I_a^* + {\beta _2}I_s^*){S^*} - {\mu _1}{S^*} = 0,} \cr {({\beta _1}I_a^* + {\beta _2}I_s^*){S^*} - (\sigma + {\gamma _1} + {\mu _1})I_a^* = 0,} \cr {\sigma I_a^* - ({\gamma _2} + {\mu _1} + {\mu _2})I_s^* = 0,} \cr {{\gamma _1}I_a^* + {\gamma _2}I_s^* - {\mu _1}{R^*} = 0.} \cr } } \right.$
Then the third equation of the system Equation (3.3) gives,
$$I_s^* = {\rm{ }}{{\sigma I_a^*} \over {{\gamma _2} + {\mu _1} + {\mu _{2.}}}}$$
Similarly, the fourth equation yields
$${R^*} = {{[{\gamma _1}({\gamma _2} + {\mu _1} + {\mu _2}) + {\gamma _2}\sigma ]} \over {{\mu _1}({\gamma _2} + {\mu _1} + {\mu _2})}}I_a^*\;.$$
Next, the first equation of Equation (3.3) yields
$${S^*} = {\rm{ }}{{\Lambda ({\mu _1} + {\mu _2} + {\gamma _2})} \over {\sigma {\beta _2}I_a^* - ({\mu _1} - {\beta _1}I_a^*)({\mu _1} + {\mu _2} + {\gamma _2})}}\;.$$
Finally, the second equation of Equation (3.3) gives
$I_a^* = {{a_1}\over{a_2}}$
where,
${a_2} = (\sigma + {\mu _1} + {\gamma _1})(\sigma {\beta _2} + {\mu _1}{\beta _1} + {\mu _2}{\beta _1} + {\gamma _2}{\beta _1})$
$\eqalign{{a_1} & = \sigma \mu _1^2 + \mu _1^2{\mu _2} + \mu _1^2{\gamma _1} + \mu _1^2{\gamma _2} + \mu _1^3 + \Lambda \sigma {\beta _2} + \Lambda {\mu _1}{\beta _1} \cr & + \Lambda {\mu _2}{\beta _1} + \Lambda {\gamma _2}{\beta _1} + \sigma {\mu _1}{\mu _2} + \sigma {\mu _1}{\gamma _2} + {\mu _1}{\mu _2}{\gamma _1} + {\mu _1}{\gamma _1}{\gamma _2}}$
Hence the endemic steady state is completely depending on
$I_a^*$
.
Basic Reproduction Number Using Next-Generation Matrix
In this section, we calculated the basic reproduction number, which is a crucial threshold in analyzing infectious disease modeling. It regulates whether the disease will die out or persist in the population.
26,Reference Diekmann, Heesterbeek and Roberts27
The basic reproduction number, denoted
${{\cal R}_0}$
, “the expected number of secondary cases produced, in a completely susceptible population, by a typical infective individual”.
Reference Martcheva28
If
${{\cal R}_0} \gt 1$
, the DFE is unstable, which means 1 primary infection can produce more than 1 secondary infection and epidemic breaks out. If
${{\cal R}_0} \lt 1$
, the DFE is locally asymptotically stable, the disease cannot persist in the population, and the situation is sustainable.
In this manuscript, we have used the next generation matrix method Reference Van den Driessche and Watmough8 to find basic reproduction number of the system Equation (2.3). We obtain 2 following matrix from the system Equation (2.3), which are F and V, they are given below
$F = \left( {\matrix{ {{\beta _1}{S_0}} & {{\beta _2}{S_0}} \cr 0 & 0 \cr } } \right)$
and
$V = \left( {\matrix{ {\sigma + {\gamma _1} + {\mu _1}} & 0 \cr { - \sigma } & {{\gamma _2} + {\mu _1} + {\mu _2}} \cr } } \right).$
Therefore, the
${V^{ - 1}}$
matrix is
${V^{ - 1}} = \left( {\matrix{ {{1 \over {\sigma + {\gamma _1} + {\mu _1}}}} & 0 \cr {{\sigma \over {(\sigma + {\gamma _1} + {\mu _1})({\mu _1} + {\mu _2} + {\gamma _2})}}} & {{1 \over {{\gamma _2} + {\mu _1} + {\mu _2}}}} \cr } } \right)$
Thus, the next-generation matrix
$F{V^{ - 1}}$
is
$F{V^{ - 1}} = \left( {\matrix{ {{{{\beta _1}({\gamma _2} + {\mu _1} + {\mu _2}){S_0} + {\beta _2}\sigma {S_0}} \over {(\sigma + {\gamma _1} + {\mu _1})({\gamma _2} + {\mu _1} + {\mu _2})}}} & {{{{\beta _2}{S_0}} \over {{\gamma _2} + {\mu _1} + {\mu _2}}}} \cr 0 & 0 \cr } } \right).$
Hence, the basic reproduction number
${{\cal R}_0}$
is
${{\cal R}_0} = \rho (F{V^{ - 1}}) = {{{{\beta _1}({\gamma _2} + {\mu _1} + {\mu _2}){S_0} + {\beta _2}\sigma {S_0}}}\over{{(\sigma + {\gamma _1} + {\mu _1})({\gamma _2} + {\mu _1} + {\mu _2})}}},$
$$\;\;\;\;\;\;\;\; = {{{\beta _1}\Lambda } \over {{\mu _1}(\sigma + {\gamma _1} + {\mu _1})}} + {{{\beta _2}\sigma \Lambda } \over {{\mu _1}(\sigma + {\gamma _1} + {\mu _1})({\gamma _2} + {\mu _1} + {\mu _2})}}.$$
We also have determined the Jacobian matrix of the system Equation (2.3) at any equilibrium point
$(\widetilde S,\widetilde {{I_a}},\widetilde {{I_s}},\widetilde R)$
which will be used for further analysis. which will be used for further analysis. The Jacobian matrix of the system Equation (2.3) is given by
${J = \left( {\matrix{ { - ({\beta _1}\widetilde{I_a} + {\beta _2}\widetilde{I_s}) - {\mu _1}} & { - {\beta _1}} & { - {\beta _2}} & 0 \cr {{\beta _1}\widetilde{I_a} + {\beta _2}\widetilde{I_s}} & {{\beta _1} - (\sigma + {\gamma _1} + {\mu _1})} & {{\beta _2}} & 0 \cr 0 & \sigma & { - ({\gamma _2} + {\mu _1} + {\mu _2})} & 0 \cr 0 & {{\gamma _1}} & {{\gamma _2}} & { - {\mu _1}}} } \right)}.$
Stability Analysis, Parameter Estimation, and Sensitivity
Initiallly, we have studied the stability analysis at the DFE point and the EE point and the statement of the results are as follows.
Theorem 2
The DFE
$({S_0},{I_{a0}},{I_{s0}},{R_0})$
of Equation (2.3) is locally stable if
${{\cal R}_0} \lt 1$
and unstable if
${{\cal R}_0} \gt 1$
.
The proof of theorem 2 can be found in the Appendix.
Theorem 3
The EE
$({S^*},I_a^*,I_s^*,{R^*})$
of the system Equation (2.3) is locally stable if
${{\cal R}_0} \gt 1$
.
The proof of theorem 3 is available in the Appendix.
Parameter Estimation
We used the Texas data from March 6, 2020, to January 31, 2021 (available at https://www.dshs.texas.gov/coronavirus/)
24
to estimate the parametric values. The transmission rates
${\beta _1}$
and
${\beta _2}$
are estimated as piece-wise function values according to the mutation behavior of COVID-19 virus using the source data in https://www.dshs.texas.gov/coronavirus/
24
and https://www.usapopulation.org/texas-population,
Reference Diekmann, Heesterbeek and Metz29
given as a possible interval in Table 2, initially making an assumption on the asymptomatically infected class population. The recruitment rate in S class, daily natural deaths, and the total population of Texas are collected from https://www.usapopulation.org/texas-population.
Reference Diekmann, Heesterbeek and Metz29
We used the average incubation period, - to -d interval, to estimate the disease transition from
${I_a}$
to
${I_s}$
class,
$\sigma $
. The recovery rates
${\gamma _1}$
and
${\gamma _2}$
from
${I_a}$
and
${I_s}$
compartments are calculated using the
$S,\;{I_a},\;$
and
${I_s}$
class source data (available in https://www.dshs.texas.gov/coronavirus/
24
and https://www.usapopulation.org/texas-population
Reference Diekmann, Heesterbeek and Metz29
), with total recovery and total cases. The disease induced death rate is estimated following the formula of WHO.
30
Model parameters values and sensitivity index
Sensitivity and Elasticity Index of
$\bf\it{{\cal R}_0}$
Because
$\bf\it{{\cal R}_0}$
provides qualitative information of an infectious disease modeling, the sensitivity and elasticity of
${{\cal R}_0}$
can play an important role in determining a disease’s control strategy. The sensitivity index of
${{\cal R}_0}$
with respect to any parameter, u is defined by
$$\partial {R_0}\partial u$$
. The elasticity index, also known as the normalized sensitivity index of
${{\cal R}_0}$
, measures the relative change of
${{\cal R}_0}$
with respect to a parameter. The elasticity index of
${{\cal R}_0}$
with respect to any parameter, u is defined by
$\Upsilon _u^{{{\cal R}_0}} = {{{\partial {{\cal R}_0}}}\over{{\partial u}}} \times {{u}\over{{{{\cal R}_0}}}}.$
${{\cal R}_0}$
increases with the positive sign of the elasticity index of the parameter and decreases with the negative sign. The magnitude of the elasticity index tells us the importance of the parameter. These measures used to determine the control of the parameters of an epidemic model. More examples can be found in Van den Driessche.
31
The numerical elasticity index of
${{\cal R}_0}$
for the baseline parameter values is provided in Table 2.
We observe that the models’ most critical parameters are the
${\beta _1}$
and
${\gamma _1}$
with the average value of
$0.8985$
and
$ - 0.771$
(see sensitivity index in Table 2), respectively. An effective way to control the outbreak could be controlling the parameters
${\beta _1}$
and
${\gamma _1}$
. As
${\beta _1}$
represents the transmission rate from S to
${I_a}$
class with the contact of
${I_a}$
, the reduction of contact between the asymptomatic individual to the general population can play an important role. One of the proven measures is to use a face mask uniformly for all populations,
Reference Van den Driessche32
significantly reducing the contact between the classes and reducing disease transmission. If doctors can find some effective viral medicine to improve the recovery rate, it could be an effective way to control the disease.
Numerical Simulation and Results
To reach the first 100,000 cases, Texas took 105 d, the second hundred thousand cases confirmed withing the next 18 d, but third, fourth, and fifth hundred thousand cases took only 11, 12, and 13 d, respectively. It is noticeable that in the month of January 2021, Texas reported over a half million confirmed cases of COVID-19 and approximately 4000 deaths. 24
The exact data also show that the second wave came almost after day 90 of the first infection wave, and during the second wave, the peak of daily case confirmation reached up to 14,916 on day 134 of pandemic. Again, the third wave had hit the locals approximately at day 250, with a peak of 28,020 on day 299. We assumed 2 major and highly effective mutations of the malign virus all-around Texas, and so imposed 3 impulsed values for the transmission rates
${\beta _1}$
and
${\beta _2}$
.
Current Epidemic Situation of COVID-19 at Texas State
In this subsection, we provide the numerical results of our proposed model and compare the results with the data of Texas. 24 Figure 2 presents the current epidemic situation of COVID-19 in Texas. The highest number of daily new confirmed cases of COVID-19 reported during December 29, 2020, and January 26, 2021 (Figure 2a). Our model predicted the same number of daily new cases compared with the daily new cases reported by https://www.dshs.texas.gov/coronavirus/. 24 As of January 31, 2021, Texas has a total of 2,059,143 cases, and our model predicted 2,081,753 total cases for the same date (Figure 2b).
Comparative solutions between data and model prediction of Equation (2.3) for (a) daily cases vs model, and (b) cumulative vs model.
Projecting the Epidemic Situation of COVID-19
The proposed model predicts another wave (third wave) to happen at the beginning of May 2021, as depicted in Figure 3a, which will have its peak at the end of June 2021. This prediction may come true if the current situation persists, ie, no clinically effective vaccine is available. But the hope is, vaccination has been started in Texas at the beginning of this Spring 2021. Reference Van den Driessche32 Moreover, if this vaccination program fails, there will be 3,439,804 cases of infection after July 2021 (Figure 3b).
Comparative solutions between data and model prediction of Equation (2.3) for (a) daily cases vs model solution, and (b) cumulative data vs. model solution.
When
$\sigma = 0.2083$
, we observe the predicted peak of daily new cases is
$19,197$
on January 03, 2021, day 303 of the pandemic. If we increase
$\sigma $
by
$10\% $
to
$0.2292,$
then the peak for the daily new cases moved to December 05, 2020, as 13,757 which also makes the third wave happen earlier. When we decrease
$\sigma $
by
$10\% $
to
$0.1875,$
then the peak for the daily new cases moved to January 14, 2021, as 27,149. This
$\sigma $
reduction also delays the next wave as depicted in Figure 4a. We see similar dynamics for the cumulative cases (Figure 4b). So
$\sigma $
plays a vital role in the dynamics of the epidemics.
Forecasting due to the effect of using model solution of Equation (2.3) for (a) daily and (b) cumulative cases.
It forecasts that on the baseline value
$\sigma = 0.2083$
by the end of July 2021, Texas may have 3,439,804 million infected cases of COVID-19. If we increase
$\sigma $
by
$10\% $
to
$0.2292,$
then the total infected cases increases from
$3.44$
to approximately
$3.51$
million (3,509,384). If we decrease the value of
$\sigma $
by
$10\% ,$
then it reduces to
$3.30$
million (3,304,451) infections (see Figure 4b).
By simulating the transmission rate from the susceptible individuals to the asymptomatic individuals,
${\beta _1}$
, we notice a little bit of change of
${\beta _1}$
plays a vital role in the forecast of the entire epidemic. For example, an increase in
${\beta _1}$
by only
$1\% $
delays the next third wave for almost 20 d, but costs more lives than the base value of
${\beta _1}$
. The base line shows 122,492 more cases than the base line solution outcomes. On the other hand, deduction of
$1\% $
in
${\beta _1}$
from the base value predicts predicts 135,103 fewer cases in total than when
${\beta _1} = 2.8361 \times {10^{ - 8}}$
(see Figure 5b). Again, it is clear that these changes in
$\beta $
cannot resist the third wave at any cost; the comparison is proved in Figure 5a. Hence, because at the time of this writing, the COVID-19 vaccine has yet to distributed to a significant number of people, the most effective way to prevent the virus is to maintain social distancing and to use a face mask (recommended by the Eikenberry et al.
Reference Eikenberry, Mancuso and Iboi33
).
Forecasting due to the only 1% changing effect of using model solution of Equation (2.3) for (a) daily cases, and(b) cumulative cases.
Similarly, decreasing
${\beta _2}$
by
$10\% $
will shorten the interval between the concurrent waves. However, the total number of cases does not change significantly because of the difference between the wave heads and bottom-lands. It shows fewer cases (159,450) than the base model fitting simulation (see Figure 6b). Again, a
$10\% $
raise at
${\beta _2} = 1.8395 \times {10^{ - 8}}$
makes the waves a little bit out-lying than the other values in Figure 6a and makes the peaks sharper. It predicts 127,946 more cases in total until July 2021 (Figure 6b).
Forecasting due to the 10% changing effect of using model solution of Equation (2.3) for (a) daily cases, and (b) cumulative cases.
From above, it is now clear that, to get a similar type of epidemic change, we needed to change
${\beta _2}$
almost 10 times compared to
${\beta _1}$
. The contact between S to
${I_a}$
class is random instead of S to
${I_s}$
class.
Figure 7 presents the daily and cumulative death rate in Texas with the best data fitting results. On the peak day, the model solution predicted the maximum daily death count was people on the last week of December 2020 (see Figure 7a). The total number of dead was approximately 41,000 by February 2021 (Figure 7b), and predicted to be by July 2021 if the spread is not under control already.
Numerical solutions and data fitting for (a) daily deaths, and (b) cumulative deaths.
Probability of a Disease Outbreak
CTMC Model
To derive a continuous-time Markov chain model assumed 4 random variables for 4 states
$(S,{I_a},{I_s},R)$
of the deterministic
$S{I_a}{I_s}R$
mathematical epidemic model. The variables are discrete-valued and time is continuous,
$t \in [0,\infty )$
,
$S(t),{I_a}(t),{I_s}(t),R(t) \in \{ 0,1,2,3, \ldots ,N\} .$
For the simplification of the notations, we use the same notations as we used in the mathematical model to define the infinitesimal transition probabilities, given in Table 3. Let
$\Delta t \gt 0$
be sufficiently small so that at most one event occurs during the
$\Delta t$
time interval. Let
$\vec X(t) = (S(t),{I_a}(t),{I_s}(t),R(t))$
and
$\Delta \vec X(t) = (\Delta S(t),\Delta {I_a}(t),\Delta {I_s}(t),\Delta R(t))$
, where
$\Delta {I_s}(t) = {I_s}(t + \Delta t) - {I_s}(t)$
, etc.
Infinitesimal transition probabilities for the
$S{I_a}{I_s}R$
mathematical model
For example, event 2 is a new infection and the probability of a new infection in time
$\Delta t$
is,
${\rm{P\{ }}\Delta X(t) = ( - 1,1,0,0)|X(t){\rm{\} = (}}{\beta _1}{{\rm{I}}_a}{\rm{ + }}{\beta _2}{{\rm{I}}_s}{\rm{)}}S\Delta t + o(\Delta t).$
To estimate the probability of a disease outbreak from the CTMC model, the effect of transmission rates from S to
${I_a}$
and from S to
${I_s}$
, which is
${\beta _1}$
and
${\beta _2}$
are considered. To approximate the probability of an outbreak, 5000 sample paths are simulated for 1 initial infected individuals, and the simulation will stop if either
${I_a} + {I_s} = 0$
or
${I_a} + {I_s} = 100$
is reached. If the total asymptomatically infected and infected population reaches 100, it is counted as an outbreak. When it reaches zero, it is assumed that there is a probability of extinction. The calculated probability of extinction is a proportion of 5000 sample paths as
$P/5000$
. Then the probability of an outbreak will be
$1 - (P/5000)$
. Parameter values are shown in Table 2. The value of
${{\cal R}_0}$
is computed from the expression of Equation (3.8).
Table 4 records the basic reproduction numbers and the probability of an outbreak for a set of values of
${\beta _1}$
and
${\beta _2}$
. The outbreak value increases when the transmission rate from the susceptible to asymptomatically infected class and from the susceptible to infected class increases. It is obvious that, when the transmission from susceptible to infected populations increase, the risk of outbreak will increase. But for the transmission from the susceptible class to the asymptomatically infected class, there is a threshold value for
${\beta _1}$
, not for
${\beta _2}$
. Here,
${\beta _1} = 6.7357 \times {10^{ - 8}}$
serves as a threshold parameter value for the disease outbreak of the ordinary differential equation system Equation (2.3) when all other parameter are fixed. Also, the probability of an outbreak depends on the initial number of infected populations.
Reference Nipa and Allen34
For our model parameter values, the basic reproduction number is
${{\cal R}_0} = 2.0087$
and the probability of an outbreak is
${{\rm{P}}_{outbreak}} = 0.4094$
.
Basic reproduction number and the probability of an outbreak are computed from the CTMC model for different value of
${\beta _1}$
and
${\beta _2}$
. Initial number of infected and asymptomatically infected populations are
${I_s}(0) = 1$
and
${I_a}(0) = 0$
Note: Parameter values are given in Table 2.
Discussion
Now we are in a world where interdisciplinary research is the most critical mechanism to live a healthy and standard life for all of humankind. Over the decades, mathematical prediction modeling used to be one of the essential tools to study and understand any epidemics’ behavior. It helps the policy-maker to make a crucial decision and get prepared for the coming days. This study’s aim was to develop a mathematical model with real-time data 24 and predict the dynamics of the COVID-19 in the state of Texas.
We used the
$S{I_a}{I_s}R$
(modified SEIR) modeling framework to design our model. We estimated the model parameters using the data from Texas
24
from March 6, 2020, to January 31, 2021. The average basic reproduction number,
${{\cal R}_0}$
, calculated as
$2.65$
for the system Equation (2.3). A study by Liu et al. 202035 calculates the median of
${{\cal R}_0}$
is
$2.79$
with a range from
$1.4$
to
$6.49$
.
Our model predicted that the net peak of daily new cases occurred at the first week of January, 2021, with new infections on the peak day (see Figure 2b). The epidemic will have another wave after April of 2021 infecting people in total during this period (Figure 3b).
According to the CDC, the incubation period of the COVID-19 is 2-14 d.
2
Lauer et al.
Reference Lauer, Grantz and Bi36
estimated the COVID-19 mean incubation period as
$5.1$
days. We considered the baseline incubation period of
$4.8$
d in our model. The incubation period may be different for the nature of the population or the individual patient. A
$10\% $
larger incubation period forces the peak epidemic size bigger and earlier with an
${{\cal R}_0}$
of
$2.6835$
and a
$10\% $
smaller incubation period delayed peak epidemic in a smaller size with
${{\cal R}_0}$
of
$2.6178$
(see Figure 4).
Scientists worldwide are trying to find the proper control measures to slow the spread of the COVID-19 virus; social distancing, quarantine of the exposed individuals, isolation of infected persons, and wearing face masks are widely used strategies. The main target for containing the virus is to reduce the contact between susceptible and exposed/infected persons, which eventually will reduce the transmission rate from the susceptible class (S) to the infected classes (
${I_a}{\kern 1pt} {I_s}$
). Texas closed schools, colleges, nonessential businesses during March 2020 and finally put the stay-at-home order at the beginning of April until the end of April 2020. This restriction delayed the earlier predicted peak on June 15 by Cooper et al.
Reference Cooper, Mondal and Antonopoulos5
We concluded that it is more important to control the asymptomatically exposed individual’s contact with the population than the symptomatic infected person, as the symptomatic infected person already develops some symptoms. Hence, they know that they are contagious and restricted them from mixing with others. The asymptomatically exposed individual who does not have any symptoms may spread the virus even without knowing that they are contagious.
Our model suggests that the sensitivity of the transmission rate due to contact between susceptible and asymptomatically infected persons is 5 times greater than the transmission rate caused by the contact between susceptible to infected individuals. A
$1\% $
increase of the transmission rate from the susceptible population to the asymptomatically exposed population increases the
${{\cal R}_0}$
from
$2.65$
to
$2.6736$
, and
$1\% $
decrease the
${{\cal R}_0}$
to
$2.6262$
(see Figure 5 for epidemic size and the total number of infection). To get a similar change on
${{\cal R}_0}$
, we had to increase/decrease the transmission rate from susceptible person to infected person,
${\beta _2}$
by
$10\% $
(see also Figure 6 for change in epidemic size and total infection). Due to the unavailability of the data, this article did not track the dynamics of the asymptomatic carrier. Tracking asymptomatic infection could provide additional insight into the dynamics of the disease.
A total of 250 of 254 counties reported SARS-CoV-2 infections, with nearly 2 million cases in all of Texas by January 31, 2021. Fifty-eight counties reported more than 1000 cases, 22 counties reported more than 5000 cases, and 11 counties reported more than 10,000 cases. In contrast, 89 counties reported less than 100 confirmed cases of coronavirus. It is noticeable that the top 6 infected counties (Harris, Dallas, Bexar, Tarrant, Travis) reported more than 50% of the total number of infections in Texas. 24
County-wise daily case analysis for the top 15 (70% of Texas) infected counties (see Supplementary Tables A1, A2, A3, A4) shows that some counties (Harris, Dallas, Tarrant, Bexar, El Paso Travis) may pass their second epidemic peak. However, other counties (such as Montgomery, Williamson) still have an uptrend in daily new infections. Proper implementation of the health expert’s suggestions, such as using face coverings, maintaining social distance, and frequent hand washing, could play an important role in epidemic dynamics.
Additionally, the under-reporting of numerous daily deaths due cuased by coronavirus in different counties within the state again goes in line with the results of our
$S{I_a}{I_s}R$
model. According to this simple yet effective exploratory prediction model, the early relaxation of lockdown (in April) within the state (overall in the United States) has potentially contributed to the exponential rise in COVID-19 cases. If the uptrend continues, the government and policy-makers have to rigorously prepare at a highly rapid pace for a large nationwide surge in patients presenting to intensive care units, also anticipating the daily rate at which front-line workers (both medical and law enforcement professionals) may succumb to death while serving the nation. The Texas governor may need to introduce draconian measures such as stay-at-home order and massive fines to limit population mobility and gatherings immediately.
The mathematical model has some limitations; our model has some as well. Texas has a different kind of epidemic scenario for different cities and counties; for instance, at a particular time, some cities and counties have stay-at-home order, other does not have the same restriction. The diversity of the population changes the contact behavior between population and changes the dynamics. Texas has has a mandatory mask policy since July 2, 2020. 37 We did not consider those variations in our modeling approach. The next iteration of the modeling approach could implement the discussed idea to reveal additional insights into disease dynamics. Furthermore, SARS-CoV-2 is a novel coronavirus, and scientists all over the world discovered new information every day, which could lead to a new understanding of the disease dynamics.
Conclusions
In conclusion, our
$S{I_a}{I_s}R$
prediction model has successfully exhibited findings in line with the current practical situation regarding the ongoing and forthcoming COVID-19 epidemic transmission in Texas, USA. The state is already in an extremely critical situation, and the peak time near the door or passing. Learning from the experiences of other high-index countries with robust health systems, the government and policy-makers of Texas need to be extra vigilant and look for alternate, effective and aggressive measures for limiting the catastrophic impacts of the contagion and protecting the front-line workers. They risk their lives to provide essential services.
Supplementary material
To view supplementary material for this article, please visit https://doi.org/10.1017/dmp.2021.151
Acknowledgments
The authors are grateful to Dr. Angela Peace for her constructive suggestions on the primary version of the manuscript.
Funding statement
M. Kamrujjaman’s research was partially supported by Bose Centre for Advanced Study and Research in Natural Sciences, 2021-2022
Competing interests
The authors declared no competing interests exist.
Disclaimer
All data are provisional and subject to change. Probable cases are not included in the total case numbers Murray. Reference Murray25
Appendix
All the supportive and supplementary results are presented in this section.
Proof of Theorem 1
Proof. Equation (2.5) gives
$\eqalign {{{{{\rm{d}}N}}\over{{{\rm{d}}t}}} & = {{{{\rm{d}}S}}\over{{{\rm{d}}}t}} +{{{{\rm{d}}{I_a}}}\over{{{\rm{d}}t}}} + {{{{\rm{d}}{I_s}}}\over{{{\rm{d}}t}}} + {{{{\rm{d}}R}}\over{{{\rm{d}}t}}} \cr & = \Lambda - {\mu _1}(S + {I_a} + {I_s} + R) - {\mu _2}{I_s} \cr & = \Lambda - {\mu _1}N - {\mu _2}{I_s}.}$
This implies
$${{{\rm{d}}N} \over {{\rm{d}}t}} \le \Lambda - {\mu _1}N.$$
Hence,
${{{{\rm{d}}N}}\over{{{\rm{d}}t}}} \lt 0$
whenever
$${\Lambda \over {{\mu _1}}} \equiv \Gamma {\rm{ < }}N(t).$$
Hence, this inequality claims that
${{{\rm{d}}N}}\over{{{\rm{d}}t}}$
is bounded by
$\rm \Gamma $
.
Now, integrating the inequality in (A.1) and using the initial condition, we get
$N(t) \le N(0){e^{ - {\mu _1}t}} + (1 - {e^{ - {\mu _1}t}})\rm \Gamma .$
Now, when
$t \to \infty $
, we obtain
$N(t) \le \Gamma $
asymptotically.
Thus it has been established that, all components of the solution of Equation (2.3) are positive and bounded in the closed region.
Proof of Theorem 2
Proof. The Jacobian matrix of the system Equation (2.3) at the DFE point is
${J_0} = \left( {\matrix{ { - {\mu _1}} & { - {\beta _1}{S_0}} & { - {\beta _2}{S_0}} & 0 \cr 0 & {{\beta _1}{S_0} - (\sigma + {\gamma _1} + {\mu _1})} & {{\beta _2}{S_0}} & 0 \cr 0 & \sigma & { - ({\gamma _2} + {\mu _1} + {\mu _2})} & 0 \cr 0 & {{\gamma _1}} & {{\gamma _2}} & { - {\mu _1}} } } \right).$
Then the characteristic equation gives the eigenvalues as
${\lambda _1} = - {\mu _1}\;,\;\;\;{\lambda _2} = - {\mu _1}\;,$
and the other two eigenvalues are the roots of the following quadratic equation
$\Psi ({\lambda _3}) \equiv \lambda _3^2 + A{\lambda _3} + B = 0,$
where,
$A = \;(\sigma + {\gamma _1} + {\mu _1}) + ({\gamma _2} + {\mu _1} + {\mu _2}) - {\beta _1}{S_0}$
$B = \;(\sigma + {\gamma _1} + {\mu _1})({\gamma _2} + {\mu _1} + {\mu _2}) - {\beta _1}({\gamma _2} + {\mu _1} + {\mu _2}){S_0} - {\beta _2}\sigma {S_0}$
$\eqalign{ & = \;(\sigma + {\gamma _1} + {\mu _1})({\gamma _2} + {\mu _1} + {\mu _2}) \left( {1 - {{{{\beta _1}({\gamma _2} + {\mu _1} + {\mu _2}){S_0} + {\beta _2}\sigma {S_0}}}\over{{(\sigma + {\gamma _1} + {\mu _1})({\gamma _2} + {\mu _1} + {\mu _2})}}}} \right) \cr & = \;(\sigma + {\gamma _1} + {\mu _1})({\gamma _2} + {\mu _1} + {\mu _2})\left( {1 - {{\cal R}_0}} \right)}$
Here, two possible scenarios can happen:
1. If
${{\cal R}_0} \lt 1$
, then B obviously positive. Again
${{\cal R}_0} \lt 1$
indicates that
$\eqalign {{{\cal R}_0} \lt 1& \cr & \Rightarrow \;{{{{\beta _1}({\gamma _2} + {\mu _1} + {\mu _2}){S_0} + {\beta _2}\sigma {S_0}}}\over{{(\sigma + {\gamma _1} + {\mu _1})({\gamma _2} + {\mu _1} + {\mu _2})}}} \lt 1 \cr & \Rightarrow \;{{{{\beta _1}{S_0}}}\over{{\sigma + {\gamma _1} + {\mu _1}}}} + {{{{\beta _2}\sigma {S_0}}}\over{{(\sigma + {\gamma _1} + {\mu _1})({\gamma _2} + {\mu _1} + {\mu _2})}}} \lt 1 \cr & \Rightarrow \;{\beta _1}{S_0} + {{{{\beta _2}\sigma {S_0}}}\over{{{\gamma _2} + {\mu _1} + {\mu _2}}}} \lt \sigma + {\gamma _1} + {\mu _1} \cr & \Rightarrow \;{{{{\beta _2}\sigma {S_0}}}\over{{{\gamma _2} + {\mu _1} + {\mu _2}}}} \lt (\sigma + {\gamma _1} + {\mu _1}) - {\beta _1}{S_0} \cr & \Rightarrow \;{{{{\beta _2}\sigma {S_0}}}\over{{{\gamma _2} + {\mu _1} + {\mu _2}}}} + ({\gamma _2} + {\mu _1} + {\mu _2}) \lt (\sigma + {\gamma _1} + {\mu _1}) \cr & \ \ \ \ + ({\gamma _2} + {\mu _1} + {\mu _2}) - {\beta _1}{S_0} \equiv A \cr & \Rightarrow \;0 \lt {{{{\beta _2}\sigma {S_0}}}\over{{{\gamma _2} + {\mu _1} + {\mu _2}}}} + ({\gamma _2} + {\mu _1} + {\mu _2}) \lt A. }$
Thus,
$A,B \gt 0$
. Then Routh-Hurwitz criterion for polynomials implies that DFE is stable.
2. If
${{\cal R}_0} \gt 1$
, then B obviously negative. Thus,
$\Psi (0) = B \lt 0$
. Again,
$\Psi ({\lambda _3}) \to \infty $
as
${\lambda _3} \to \infty $
. Since
$\Psi ({\lambda _3})$
is a continuous function of
${\lambda _3}$
, hence by Bolzano’s theorem on continuous function we have
$\Psi ({\lambda _i}) = 0$
for some
${\lambda _i} \gt 0$
. Therefore, at least one eigenvalue of the Jacobian matrix is positive. Hence, DFE point is unstable equilibrium point.
This completes the proof.
Proof of Theorem 3
Proof. The characteristic equation of the Jacobian matrix (3.9) at the EE is
$|{J^*} - \lambda {\cal I}| = 0$
where,
${\cal I}$
is an
$4 \times 4$
identity matrix; which gives,
$( - {\mu _1} - \lambda )|{\lambda ^3} + A{\lambda ^2} + B\lambda + C| = 0$
where,
$A = \;\sigma + {\gamma _1} + {\gamma _2} + 3{\mu _1} + {\mu _2} + {\beta _1}I_a^* + {\beta _2}I_s^* - {\beta _1}{S^*}$
$\eqalign{ & B = \sigma {\gamma _2} + {\gamma _1}{\gamma _2} + 2\sigma {\mu _1} + \sigma {\mu _2} + 2{\gamma _1}{\mu _1} + {\gamma _1}{\mu _2} + 2{\gamma _2}{\mu _1} \cr & \ \ \ \ \ + 2{\mu _1}{\mu _2} + 3\mu _1^2 + \sigma {\beta _1}I_a^* + \sigma {\beta _2}I_s^* \cr & \ \ \ \ \ + {\gamma _1}{\beta _1}I_a^* + {\gamma _2}{\beta _1}I_a^* + {\gamma _1}{\beta _2}I_s^* + {\gamma _2}{\beta _2}I_s^* \cr & \ \ \ \ \ + 2{\mu _1}{\beta _1}I_a^* + {\mu _2}{\beta _1}I_a^* - \sigma {\beta _2}{S^*} + 2{\mu _1}{\beta _2}I_s^* \cr & \ \ \ \ \ + {\mu _2}{\beta _2}I_s^* - {\gamma _2}{\beta _1}{S^*} - 2{\mu _1}{\beta _1}{S^*} - {\mu _2}{\beta _1}{S^*} \cr & \ \ = \;2{\mu _1}A + \sigma {\gamma _2} + {\gamma _1}{\gamma _2} + \sigma {\mu _2} + {\gamma _1}{\mu _2} + \sigma {\beta _1}I_a^* \cr & \ \ \ \ \ + \sigma {\beta _2}I_s^* + {\gamma _1}{\beta _1}I_a^* + {\gamma _2}{\beta _1}I_a^* + {\gamma _1}{\beta _2}I_s^* \cr & \ \ \ \ \ + {\gamma _2}{\beta _2}I_s^* + {\mu _2}{\beta _1}I_a^* - \sigma {\beta _2}{S^*} + {\mu _2}{\beta _2}I_s^* \cr & \ \ \ \ \ - 3\mu _1^2 - {\gamma _2}{\beta _1}{S^*} - {\mu _2}{\beta _1}{S^*} \cr} $
$\eqalign{ & C = \;\sigma \mu _1^2 + {\gamma _1}\mu _1^2 + {\gamma _2}\mu _1^2 + \mu _1^2{\mu _2} + \mu _1^3 + \sigma {\gamma _2}{\mu _1} \cr & \ \ \ \ \ + {\gamma _1}{\gamma _2}{\mu _1} + \sigma {\mu _1}{\mu _2} + {\gamma _1}{\mu _1}{\mu _2} + \mu _1^2{\beta _1}I_a^* \cr & \ \ \ \ \ + \mu _1^2{\beta _2}I_s^* - \mu _1^2{\beta _1}{S^*} - {\mu _1}{\mu _2}{\beta _1}{S^*} + \sigma {\gamma _2}{\beta _1}I_a^* \cr & \ \ \ \ \ + \sigma {\gamma _2}{\beta _2}I_s^* + {\gamma _1}{\gamma _2}{\beta _1}I_a^* + \sigma {\mu _1}{\beta _1}I_a^* \cr & \ \ \ \ \ + \sigma {\mu _2}{\beta _1}I_a^* + {\gamma _1}{\gamma _2}{\beta _2}I_s^* + \sigma {\mu _1}{\beta _2}I_s^* + \sigma {\mu _2}{\beta _2}I_s^* \cr & \ \ \ \ \ + {\gamma _1}{\mu _1}{\beta _1}I_a^* + {\gamma _1}{\mu _2}{\beta _1}I_a^* + {\gamma _2}{\mu _1}{\beta _1}I_a^* \cr & \ \ \ \ \ + {\gamma _1}{\mu _1}{\beta _2}{I^*} + {\gamma _1}{\mu _2}{\beta _2}I_s^* + {\gamma _2}{\mu _1}{\beta _2}I_s^* + {\mu _1}{\mu _2}{\beta _1}I_a^* \cr & \ \ \ \ \ - \sigma {\mu _1}{\beta _2}{S^*} + {\mu _1}{\mu _2}{\beta _2}I_s^* - {\gamma _2}{\mu _1}{\beta _1}{S^*} \cr & = \;\mu _1^2A + \sigma {\gamma _2}{\mu _1} + {\gamma _1}{\gamma _2}{\mu _1} + \sigma {\mu _1}{\mu _2} + {\gamma _1}{\mu _1}{\mu _2} \cr & \ \ \ \ \ + \sigma {\gamma _2}{\beta _1}I_a^* + \sigma {\gamma _2}{\beta _2}I_s^* + {\gamma _1}{\gamma _2}{\beta _1}I_a^* \cr & \ \ \ \ \ + \sigma {\mu _1}{\beta _1}I_a^* + \sigma {\mu _2}{\beta _1}I_a^* + {\gamma _1}{\gamma _2}{\beta _2}I_s^* + \sigma {\mu _1}{\beta _2}I_s^* \cr & \ \ \ \ \ + \sigma {\mu _2}{\beta _2}I_s^* + {\gamma _1}{\mu _1}{\beta _1}I_a^* + {\gamma _1}{\mu _2}{\beta _1}I_a^* \cr & \ \ \ \ \ + {\gamma _2}{\mu _1}{\beta _1}I_a^* + {\gamma _1}{\mu _1}{\beta _2}I_s^* + {\gamma _1}{\mu _2}{\beta _2}I_s^* + {\gamma _2}{\mu _1}{\beta _2}I_s^* \cr & \ \ \ \ \ + {\mu _1}{\mu _2}{\beta _1}I_a^* + {\mu _1}{\mu _2}{\beta _2}I_s^* - 2\mu _1^3 \cr & \ \ \ \ \ - {\mu _1}{\mu _2}{\beta _1}{S^*} - \sigma {\mu _1}{\beta _2}{S^*} - {\gamma _2}{\mu _1}{\beta _1}{S^*}\;. \cr} $
Here, the first eigenvalue is negative, then the system will be locally asymptotically stable if other 3 eigenvalues are all negative or their real parts are negative. The other 3 eigenvalues will be negative or will have negative real parts if the Routh-Hurwitz criterion is satisfied.
From Routh-Hurwitz criterion, we can say that the EE point is stable if
$A \gt 0,\;C \gt 0$
and
$AB \gt C$
. Therefore the EE of Equation (3), which exists if
${{\cal R}_0} \gt 1$
is locally asymptotically stable.
Open access
