1. Introduction
In December 2019, after the increment in pneumonia cases in Wuhan (China), the authorities ordered respiratory tests to be carried out to find out its origin. They discovered that the new rise in these cases was caused by a novel virus, which was coincident with severe acute respiratory syndrome coronavirus (SARS-CoV) and other mammal coronaviruses (especially bat and pangolin); their similarity was more than 70% and 95%, respectively [
1,
2,
3]. Therefore, this new virus was named SARS-CoV-2 (more commonly known as COVID-19) and, in spite of being spotlighted due to the pneumonia cases, it had recently been found via survey that the most common symptoms were cough, fever, and weakness [
4,
5]. Even though many countries implemented prevention measurements to avoid the spread of COVID-19, by 9 March 2020, more than 118,000 cases and 4291 deaths were reported in 114 countries, and the World Health Organization (WHO) characterized this virus as a pandemic [
6]. This situation highlighted the importance of epidemic models, and some derived terms such as the effective reproduction number or herd immunity, which allow the epidemic situation and its future behaviour to be determined.
The basic idea behind these mathematical models is to split the population into different subpopulations depending on their medical condition [
7,
8,
9,
10,
11,
12,
13,
14,
15], so the SIR model (which divides the population into susceptible, infectious, and recovered subpopulations), SEIR model (the exposed subpopulation is included in the SIR model), etc. are built up, and the conditions to eradicate the disease are usually obtained and analyzed. For example, in [
14], the vaccination coverage level needed to eliminate Ebola from a population is given and, in [
15], different prevention policies are compared. Usually, diseases do not affect different age groups homogeneously (i.e., the elderly population are often more sensitive to diseases), so many studies consider age-group epidemic models [
16,
17,
18,
19]. While the complexity is increased, the obtained results are more precise. For instance, in [
16], data concerning the COVID-19 epidemic evolution in Shijiazhuang City, Hebei Province in China and three different prevention strategies, with respect to different age groups, are characterized by a set of epidemic features, which gives an overview of each strategy’s impact. In [
17], to reach the WHO’s target measles incidence rate in India, they concluded that it is necessary to increase the vaccination coverage rate among children of age 0–4 years. The literature regarding control techniques for epidemic models is also exhaustive; in some studies, they apply constant and pulse vaccination [
20,
21,
22,
23], and they show it to be simple and effective. Other studies use control theory techniques, such as feedback control [
24,
25,
26], and they show better convergence time and steady-state errors. Machine learning has also been applied to make more accurate models and therefore tackle more efficiently the control problem [
27,
28,
29,
30,
31].
In addition, due to the global impact caused by COVID-19, many studies have put into practice this kind of mathematical model for this specific case; the works in [
32,
33] consider a quarantine subpopulation, and the condition that makes the disease disappear from the whole population is obtained. Moreover, in both cases, the model is validated with real data regarding the epidemic evolution in Saudi Arabia [
32] and Italy [
33], so its future evolution can be predicted. During the first stage of the COVID-19 epidemic evolution, as no vaccine was available, reducing social contact was a widely adopted measurement all over the world, and the impact of this prevention method is researched in [
34,
35]. Other works included vaccination strategies such as newborn vaccination [
36], vaccinating a proportion of the susceptible subpopulation [
37,
38], and implementing dynamic vaccination and treatment with an SIR-like model [
39]. All of these analyze the local and global stability of their respective systems’ equilibrium points, among other things.
In this paper, an SIHR model has been built up where the total population has been divided into susceptible (S), non-seriously infected (I), hospitalized (H), and recovered (R) subpopulations, and a vaccination control strategy has been included. Thus, taking each subpopulation as a state variable and vaccination as the input variable, this model lays out a set of first-order non-linear differential equations. It has been assumed, as in a real-case situation, that S, I, and R are unknown, so only vaccination and hospitalized variables have been fed back to the vaccination dynamics through the gains and , respectively. Moreover, to make a flexible vaccination, a free-design time-dependent function has been included in the vaccination control dynamics. From the formulated system, the disease-free equilibrium point () and the endemic equilibrium point () have been calculated, and it has been proved that is reachable under certain conditions. In addition, considering the existence of a vaccination and making the right choice of the free-design parameters (, , and ), it is possible to turn into . The conditions for the global and local asymptotic stability of have been obtained. In the case of the endemic equilibrium point in absence of vaccination (), the Routh–Hurwitz criterion has been used to conclude that it is locally asymptotically stable whenever it is reachable and, from Rouche’s theorem, the conditions for the local asymptotic stability of have been inferred.
Finally, a value has been given to each parameter based on the background literature [
40,
41,
42,
43,
44,
45], and several simulations have been carried out to reinforce and display the achieved theoretical results. Overall, the main novelties that this paper presents are as follows:
- (i)
The implementation of a vaccination (dynamic control input), where hospitalization and the vaccination itself are fed back;
- (ii)
Including a free-design time-dependent function into vaccination;
- (iii)
Considering the vaccine stock in the vaccination strategy;
- (iv)
A strategy to choose the vaccination free-design time-dependent function and parameters such that and are locally asymptotically stable and unstable, respectively.
The first novelty, in comparison with the background literature, gives a more realistic feedback vaccination control; it includes only hospitalization and vaccination state variables as feedback of the vaccination dynamics, which are related to data that are usually known or accessible. This vaccination can easily be adapted to the desired specifications since it includes a free-design time-dependent function (second novelty). Meanwhile, in practice, the vaccines stockpile is reduced as they are provided; many epidemic models usually ignore this constrain. Therefore, the introduced vaccination has been constrained by the vaccines stockpile; the number of individuals who receive a dose cannot exceed the vaccines stock (third novelty). This characteristic, with the fact of free-design parameters (, , and ), allows a realistic vaccination strategy to be designed. The controlled reproduction number (), which depends on and , has been found to be a fundamental part of the local asymptotic stability of and its uniqueness. So, with the objective of reducing as much as possible the epidemic impact, one can consider and select a pair of and to ensure the local asymptotic stability of (fourth novelty). However, if the condition for the local asymptotic stability of can not be accomplished (i.e., there is a low vaccine reserves), one can consider a suitable and the conditions that make it locally asymptotically stable.
The paper is organized as follows:
Section 2 describes the SIHR model within the vaccination policy: what characterized the subpopulations, which parameters are taken into account, how vaccination affects each subpopulation, and how each subpopulation transforms into another. Thus, the resulting first-order non-linear differential system is formulated. Afterwards, considering non-negative initial conditions, the non-negativity of each state variable, with and without vaccination, is studied and proved. Finally, from the dynamical system,
and
are obtained, and the conditions for their existence are discussed.
Section 3 is divided into three subsections;
Section 3.1 and
Section 3.2 are focused on the stability analysis of
, whereas
Section 3.3 focuses on
. In
Section 3.1, the condition for the local asymptotic stability of
is proved with the next generation matrices and, in
Section 3.2, the solutions of the differential equations are used to derive a sufficient condition for the global stability. In
Section 3.3, the system is linearized about
, and the Jacobian matrix is inferred. Then, the eigenvalue problem of the Jacobian matrix is formulated, and the conditions that ensure a negative real part for all eigenvalues (so
is locally asymptotically stable) are obtained from the Routh–Hurwitz criterion and Rouche’s theorem. In
Section 4, several simulations have been carried out which show an accordance between the theoretical results and the numerical ones. Moreover, in a simulation where a vaccine stock function has been introduced, the free-design parameters and function have been chosen based on a desired vaccination strategy. In addition to illustrating the desired behavior, it shows that the disease is eradicated. Finally, the paper ends with
Section 5, where the results are discussed.
2. Model Description
When diseases spread through big populations, in terms of prediction, deterministic models have been proven to give good results, so a deterministic model has been built up. This type of models is usually represented by flow charts, which are composed by two main parts: blocks and arrows. Each block stands for a subpopulation characterized by its medical condition (e.g., healthy, sick, recovered, etc.), while the arrows indicate how individuals from one block transform into another.
Figure 1 shows the deterministic system which adds a new block (
H) to the SIR model. Those compartments are used to represents the following subpopulations:
Susceptible (S): Group of individuals that can catch the disease. They are not yet infected nor have immunity against it.
Non-seriously infected (I): Group of individuals that caught the disease and present symptoms. They are contagious and therefore responsible for the disease spreading. It is assumed that they do not have grave symptoms. Note that in order to simplify the technical nomenclature, the non-seriously infected population is simply referred to as the "infected subpopulation" in the sequel.
Hospitalized (H): This group of people is characterized by suffering from serious symptoms and therefore being hospitalized. Since they are hospitalized, it is assumed that prevention measurements are strict, and consequently, they are not contagious.
Recovered (R): Group of individuals with immunity on them, that is; infected or hospitalized people that already have been recovered or those susceptibles that were vaccinated.
Vaccinated (v): It defines the number of susceptible people vaccinated per unit time which are introduced in the recovered subpopulation. It is used for control purposes.
The parameters appearing in
Figure 1 are:
: newborns per unit time.
: transmission rate.
: natural mortality rate.
p: probability of being hospitalized once you catch the disease.
: recovery rate of non-seriously infected individuals.
: recovery rate of hospitalized individuals.
: death rate of hospitalized individuals.
: immunity loss rate.
These parameters are conditioned by
N is the total population size, and it is assumed to be sufficiently large. Newborns are incorporated to the class S with a rate . Considering that all inhabitants die naturally indistinctly to the subpopulaiton they belong to, a portion proportional to will be removed from all compartments. Hospitalized individuals not only will die naturally but also due to the disease; a portion will also be removed.
Assuming that the population is homogeneously mixed, the rate
indicates the probability of infecting susceptible individuals when they are in contact with infected individuals. Thus, the mathematical expressions is
where
c is the average number of close contacts per day of a member from
I, and
is the probability of infecting susceptible individuals when there is a close contact. Hence, how the susceptible individuals get infected can be defined; there is a probability
that an individual from
S contacts an individual from
I. Therefore, an individual from
S has a probability
of being infected, and consequently,
susceptible individuals will be infected. Models with this type of transmission are known as mass action models.
The probability of those who caught the disease being hospitalized is p, so individuals will be hospitalized. The rest () will go to I. Infected individuals from I and H subpopulations will recover at rates and , respectively; that is, and will be removed from I and H, respectively, and they will be introduced into R. The vaccinated susceptible individuals per day (v) are extracted from S and brought to R. However, immunity disappears after days, and they become susceptible again.
Each state variation with respect to time is equal to the added individuals minus the removed ones. Thus, the flow chart shown in
Figure 1 leads to the following time differential system,
and initial conditions
,
,
, and
. Regarding the total population,
The controller
is defined by the following equation:
where
are tuning parameters adjusted to yield a desired behavior, and
represents the free-design time-dependent function. Note that Equation (
5) is a first-order differential equation, where only
and
are taken into account, and its solution will be a feedback function of the state variables and time.
2.1. Non-Negativity of the Solution
To demonstrate the consistency of the equations given in (
3), (
4), and (
5) the sizes of the subpopulations
,
,
, and
, and the total population
must be non-negative as well as the vaccination per unit time
. It can be shown that all the solutions will be kept non-negative for a given finite non-negative initial condition under certain reasonable constrains.
Theorem 1. Assuming that there is no vaccination ( for all time), the following properties hold:
- (i)
The solution of the epidemic model defined by the set of Equation (3) is non-negative for all and for any finite initial condition such that , , , and .
- (ii)
The setwhere , is positively invariant with respect to the dynamical system (3).
Proof of Theorem 1. Since the differential equations shown in (
3) and (
4) are first-order differential equations, each solution is obtained with the superposition of the respective homogeneous part and the forcing part, so it leads to
From (
7), it is straightforward that
is non-negative for all time if the initial condition is non-negative. By introducing the expressions for
and
into
, one obtains:
From (
7) and (
8), it follows that:
After removing some terms from (
9), the next inequality can be obtained:
From the expression above, it follows that
. Therefore, we proceed by contradiction to prove that
. Assume that
is the first time instant at which
. Taking into account that
for all
since
, see the second equation in (
3), the condition
can only be fulfilled if the total population was negative in a previous time interval. Therefore, it will be assumed that there exists a time interval at which the total population becomes negative; that is,
for any
, and
for any
.
From the expression for
in (
7), it is concluded that
for
as a result of
non-negativeness, since
,
, and
are non-negative during the given time interval. As all the integrands of
are non-negative, see (
7),
will also be non-negative.
Now, considering the time derivative of
given in (
4), its solution is obtained by direct calculation:
and, by substituting the expression for
in (
7) into (
11), it follows that:
Let
be defined as
and it yields
Taking into account that
, it leads to
Previously, it was found that all the subpopulations are non-negative for any
. Therefore, it follows that
for any
and
, where
. Continuing with the contradiction, for any
, Equation (
12) can be rewritten as follows:
which contradicts the existence of a time interval where
. Therefore,
for all
, and consequently, a
at which
does not exist. Thus, from (
10), one concludes that
, and
too, for all
. It was previously proved that
,
, and
non-negativity implies
and
non-negativity, so property (
i) is proved.
Moreover, considering that nobody is dying because of the disease (
), it is possible to prove that the entire population will be bounded. The solution of the differential Equation (
4) is given by
and the limit
So, let
be a set defined as
where
. Hence, given any initial condition belonging to
, the solution of system (
3) will remain in
. Therefore, the set
is positively invariant and Property (ii) is proved. □
Theorem 2. If a vaccination exists, the following properties hold:
- (i)
The solution of the epidemic model, defined by the set of equations in (3) and (5), is non-negative for any initial condition such that , , , , and , if the following condition is accomplished:where , , and are defined as follows: and are the maximum number of susceptible individuals that can be vaccinated and the vaccination stock, respectively, where , and - (ii)
The set is positively invariant with respect to the dynamical system (3).
Proof of Theorem 2. First, the solution of
is calculated:
Considering that
and substituting into (
22), one obtains:
and, by applying the integration by parts, it follows that:
so the expression for
is
Depending on the non-negativity of
, two possibilities are found; if
, then
is necessary so that
. In contrast, if
, then
is sufficient so that
. Therefore, if the left part inequality in (
20) is accomplished,
is proved to be true. If the right part inequality in (
20) is fulfilled, then
. Then,
provided that (
20) is fulfilled.
To prove each subpopulation’s non-negativity through time, the solutions of the equations shown in (
3) have been calculated:
Taking into account the solution of
given in (
26), and the constrain
, it leads to
By contradiction, it will be demonstrated that
is positive for all
and
. Let us suppose that
becomes negative during a time interval; that is,
for any
, and
for any
. Then,
Thus, a time interval at which does not exist, so for all .
Now, as in the proof of Theorem 1, we proceed by contradiction to conclude that
for all
; let us assume that
for any
, and
for any
, then
for
since
is non-negative during
. From the solution of
given in (
26), it follows that
for any
as all the integrands are non-negative. Regarding the solution of
, the expression given in (
12) is still valid for the system (
3) when vaccination is being applied, and
for all
since all the subpopulations are non-negative at
. Thus,
so a time interval
where
does not exist, and consequently,
for all
. Property (i) is proved.
Considering that nobody dies due to the disease, the solution of the differential Equation (
4) is given by the expression (
17). Consequently, the set
, see (
19), is positively invariant with respect to the dynamical system (
3) with vaccination. Property (ii) is proved. □
2.2. Equilibrium Points
At an equilibrium point, the system variables do not change with respect to time; therefore, the time derivatives in (
3) and (
5) are set to zero (
), and its solution will give the equilibrium points. Two different equilibrium points are found:
and
.
Definition 1. is an equilibrium point of the system (3), and it is characterized by a null infected subpopulation:whereand the total population .
This result shows that the susceptible and the recovered subpopulations hold the whole population, and a population exchange exists between them due to vaccination. Moreover, vaccination effort must be limited to not contradict Theorem 2 and consequently ensure the solution’s non-negativeness. In a real situation, this implies that it is not possible to vaccinate more individuals than those in the susceptible suppopulation.
Note that has been considered as time-invariant () in order to calculate the equilibrium points.
Proposition 1. is reachable (each subpopulation has a non-negative value) if and only if the vaccination control law tuning parameters and , and f, are chosen such that:where .
Taking into account the expressions for each subpopulation given in (
31),
is defined as follows:
where
and, in this particular case, it can be observed that the
holds the whole population:
.
Definition 2. is an equilibrium point of the system (3), and it is characterized by a non-zero infected subpopulations ( and ):whereand The total population .
Remark 1. Considering Theorem 2 and all the subpopulations in (36), the following hold: - (i)
From the expression for in (36), it can be seen that a portion of the susceptible subpopulation, which is proportional to (deaths due to the disease), is removed from . In the case of a disease with a high death rate, the number of deaths could surpass the newborns, and it could cause a negative value of , which contradicts Theorem 2.
- (ii)
Other variables’ non-negativeness only depends on ’s non-negativeness; that is, the variables , and are linear with respect to , and their respective independent parameters are positive, so the condition is sufficient to ensure , , and .
- (iii)
Let , then ; when , then , , and . Considering that (equivalent to ), it follows that .
- (iv)
From the point above, it is possible to transform into with a suitable vaccination; that is, B can be modified to , or equivalently , so .
Proposition 2. The endemic equilibrium point is not reachable, in the sense that it has some negative component, if Proof of Proposition 2. Considering the expressions for
and
in (
36), the normalized variable
is calculated:
If the condition is fulfilled (i.e., ), then . This result contradicts Theorem 2 since it demonstrates that all the subpopulations are positive, and therefore all of them will be bounded by the whole population N; that is, the normalized subpopulation can not exceed 1.
If , then , which implies that the susceptible subpopulation holds the whole population (i.e., ). This particular case does not correspond with the definition of , see Definition 2, but with the definition of , see Definition 1. □
Proposition 3. Assuming that:the endemic equilibrium point is reachable if the following conditions are fulfilled: - (i)
- (ii)
Proof of Proposition 3. Considering the expression for the denominator of
in (
36):
and taking into account
and
, by direct calculations, it is obtained that:
Assuming that
, it follows that:
which leads to
If the condition (i) from Proposition 3 is fulfilled, the nominator of will be negative, and consequently, will be positive since its denominator is negative. □
Remark 2. As it was mentioned in Remark 1, part (i), could reach a negative value due to a very aggressive disease. Taking into account the expression for in (36), it is concluded that this situation cannot occur if . By doing this, one obtains the following:and considering expression (41), it follows that the inequality above is still true if the next inequality is fulfilled:which leads to The condition above is always fulfilled, so a value of α that neglects the existence of does not exist.
The expressions for each subpopulation given in (
36) include the particular case in which no vaccination is being applied. Thus, considering the expressions in (
36), one can define
as follows:
where
and the total population
.
4. Simulations
Some simulations have been carried out to validate the results obtained in the previous sections. The parameter values have been gathered from different sources; regarding the natural birth and death rate, the Spanish public source Instituto Nacional de Estadística (National Statistics Institute) [
40] has been used. In a survey where 28,503 people participated [
41], their contact rate with respect to their incomes, locations, age, etc. was analyzed, and they showed an overall value of 14.5 contacts per day. The transmission probability has been proved to vary in different settings [
42]; that is, in cases in which the contact is more prolonged, the transmission probability value can reach the
, and in working places, it decreases to
. It has been observed that the hospitalization risk among the COVID-19 patients in England changes notoriously between groups, being the greatest (
) in the elderly population (80 years or older), and in the case of a medium-age population (40–49 years), the probability decreases to
[
43]. In addition, all age groups achieve a peak value in winter months, which is a trend that may be caused by a weakener health system during this period. Regarding hospitalized people’s recovery and death rates, confirmed COVID-19 cases from Belgium are processed in [
44], and it is concluded that people around 20–60 years old stay 8.2 days in the hospital until they recover. In addition, on average, they stay 12.2 days until they die. With respect to the duration of immunity, the immunity wanes after 3–24 weeks after vaccination, and is not until the 24th week that there is an important decline in immunity [
45]. Based on available data,
Table 1 shows a summary of the parameter values used for simulation proposes within their respective references. Note that depending on the analysis, the contact rate (
c) will be modified.
Considering the parameters given in
Table 1, the condition
is fulfilled, and by substituting their values into the denominator of
, see (
36), one obtains:
which is negative as far as
, as it was seen in the
proof of Proposition 3.
One can rewrite the nominator of
in terms of
:
and in agreement with the statement exposed in Remark 3, from (
82) and (
83), it is concluded that the endemic equilibrium point is reachable (unrechable) when
(
).
The simulations have been carried out with the version R2022a of Matlab, and the differential equations shown in (
3) have been solved with the solver
ode45(), where the selected time step has been 1 day; that is, all the results regarding the epidemic evolution (any subpopulation or vaccine evolution) are given in discrete form, which are equally spaced data points in time (i.e., 1 day). For the graphical representation, the function
plot(), which generates a line plot from the data points, has been used. Note that in the interest of clarity in data visualization, the function
plot() has been applied instead of
stem(), which is used to plot discrete sequence data.
Since in
Section 3.2, it was proved that the disease-free equilibrium point is globally asymptotically stable under certain conditions, in the following simulations, this characteristic will be verified, and subsequently other simulations will be executed for the local asymptotic stability analysis. With regard to the endemic equilibrium point, with and without vaccination, the condition for its stability will be evaluated for different
values, and it will be verified that the requirements given in Theorem 5 are accomplished. In
Section 4.3, the state feedback vaccination method is compared with other common vaccination methods based on the hospitalized subpopulation and vaccination evolution. Finally, a vaccination specification is defined, and the tuning parameters are set up to match the desired specification.
4.1. Stability of the Disease-Free Equilibrium Point
4.1.1. Global Stability
Firstly, the condition (
32) shown in Proposition 1, which ensures the existence of
, has been considered (i.e.,
). Consequently, by establishing
, the condition is fulfilled and
is reachable. Then, the conditions for the global stability have been considered, see Proposition 5; taking into account that
must be fulfilled (first condition), with
, the requirement is satisfied, and it leads to the contact rate
(for the simulation carried out in this subsection, instead of using the value of the contact rate shown in
Table 1, the value
has been used). With
, it follows that
(second condition). Consequently, if one chooses
, the second condition is fulfilled and
is globally stable. With
, all parameters needed for simulation purposes have been already defined. Regarding the initial conditions, it has been assumed that there is a low quantity of non-seriously infected and hospitalized people and that nobody has been vaccinated yet; that is,
,
,
,
, and
. Then, the epidemic evolution has been simulated; see
Figure 2.
As it was expected from the theoretical results, it can be observed that all the subpopulations tend asymptotically to the disease-free equilibrium value. After 37 days the number of hospitalized people is reduced to less than one, and after 30 days, the same happens with the infected individuals. During this period, 404 people die as a consequence of the disease. This type of epidemic evolution is always desired, since the the number of infected and hospitalized individuals decay rapidly to zero regardless the outbreak magnitude. However, the social effort needed for this propose is high; considering that , the contact rate must be reduced to less than contacts per day so that , which implies self-isolation.
4.1.2. Local Stability
As it was previously seen,
is reachable with
. Now, the most common way to achieve a
lower than one, so
is locally asymptotically stable, is to modify the contact rate or/and the vaccination tuning parameters. Considering the contact rate given in
Table 1, and supposing that its value is fixed, in the following two simulations, the asymptotic stability will be analyzed with respect to the different vaccination tuning parameters. For both simulations, it has been assumed that initially,
R is big compared to
S;
,
,
,
, and
.
From the condition
, so
is unstable, it follows that
: for
,
, and
, one obtains that
and
. The epidemic evolution shown in
Figure 3 is obtained where the subpopulations do not reach
, but
:
,
,
,
and
. This is a predictable result since, for the parameters itemized in
Table 1,
is reachable when
.
From the condition
, so
is locally asymptotically stable, it follows that
. With
,
, and
,
and
are attained. The resulting epidemic evolution has been depicted in
Figure 4, where one can see that the model tends to
; after 164 days approximately, the non-seriously infected and the hospitalized subpopulations are less than one, and some time after,
S and
R reach the points
and
, respectively, as expected.
In
Table 2, some data of interest regarding previous simulations are summarized.
During the first days, when the infected subpopulation starts to increase again, even if the differences between both cases are minor, it is observed that a higher vaccination helps to avoid future deaths. At the end of 2 years, the differences between both cases are more notorious; when , the total vaccination is higher than the total vaccination in , whereas the total deaths decreases four times.
The simulations of these particular cases are significant because the initial conditions are comparable to the current COVID-19 epidemic situation, and they highlight the importance of maintaining , since it avoids a significant increment in the total number of deaths.
4.2. Stability of the Endemic Equilibrium Point
4.2.1. Local Stability without vaccination
Theorem 4 states that if
is reachable (i.e.,
), then it is locally asymptotically stable. With the objective of reinforcing this theorem, the terms from the Routh tablet (
and
) given in (
A7),
Appendix A, can be evaluated numerically for different
values and verify its stability; if the terms are positive, then
is locally asymptotically stable. Whereas, if any term is negative, then
is unstable. Considering the parameter values given in
Table 1, the terms in (
A7) have been evaluated for different
values, and the result is depicted in
Figure 5.
In
Figure 5, a red line is also shown to emphasize that only the values of
and
that correspond to
must be taken into account. As it was expected from the theoretical results,
and
are positive if
.
4.2.2. Local Stability with Vaccination
Theorem 5 states that it is sufficient if the conditions (
74) and (
75) are fulfilled to ensure the local asymptotic stability of
. In this section, the procedure to compute the condition for
, which is exposed in Remark 5, will be followed.
Let us consider the parameters given in
Table 1 and the following vaccination tuning parameters values
and
, so
(
is not locally asymptotically stable). Then, the condition (
74) has been computed:
To calculate the numerical value for the second condition, see expression (
75), the bode magnitude plot of the function
has been computed, which is shown in
Figure 6. From
Figure 6, the peak has been obtained, which corresponds to
dB. Thus,
Finally, the minimum value between the values given in (
84) and (
85) has been chosen;
.
The roots of the characteristic equation obtained from the linearized system (
3) about
, see Equation (
A12), have been computed for different
cases;
,
and
. See
Table 3 for more details.
Taking into account the results in
Table 3, variations in
barely affect
and
. In the case of
,
, and
, when
increases, their real parts decrease (the time for the response to reach the final value will be bigger), whereas the imaginary parts of
and
tend to increase (the response of the associated eigenvectors will be characterized by bigger oscillations). Regarding the stability, two of three are stable since all
, where
, have negative real parts. As it was mentioned in Remark 5, even if the condition
is not fulfilled, the system might still be stable since not fulfilling the conditions (
74) and (
75) does not imply that the equilibrium point will be unstable. However, when
,
and
have a positive real part and they are complex-conjugate poles. This instability is caused by the violation of the condition (i) given in Theorem 2 (i.e.,
).
4.3. Study of State Feedback Vaccination
To analyze the benefits of the state feedback vaccination method, its performance will be compared with that of other common vaccination control methods, such as the constant and the proportional methods. These controllers will be defined as follows.
Note that
is the simplified version of the equation given in (
5). In addition,
and
are functions of the state
H, since it has been assumed that data exist regarding the hospitalized people, which makes its implementation possible. Each vaccination function will be set up to drive the hospitalized subpopulation to the desired value
, and each performance will be evaluated based on the following features:
,
, the total number of people who receive a vaccine dose up to 365 days, and the total number of deaths up to 365 days.
Let
be the endemic equilibrium point belonging to
, where
. Taking the parameter values in
Table 1,
has been calculated, and the condition
has led to
,
, and
for
. Two values of
have been considered,
and
, thereby giving
and
, respectively. Then, with the initial conditions
,
,
,
, and
, a simulation of the epidemic evolution has been carried out, and in
Figure 7, the evolution of the hospitalized subpopulation per day and the number of vaccinated people per day has been displayed. Note that the settling time is bigger than 365 days, since none of the four cases reach the desired state
at the end of 1 year.
As it can be observed in
Figure 7, with our state feedback controller, the vaccination performance can be modified; it is possible to switch from an underdamped-like behavior to an overdamped. Meanwhile, with
and
, vaccination is totally conditioned by a constant value and
H, respectively.
Taking into account the performance measurements given in
Table 4, if the objective is to reduce as much as possible the number of deaths, then
gives the best response while
gives the worst. However, these types of responses are usually difficult to implement (e.g., there could be staff limitations or cost limits), and in such cases, big vaccination peaks (e.g., the peaks observed with
, and
when
) are not desired. By decreasing the value of
,
can be adjusted to the desired response.
4.4. Study of Vaccination Strategy Design and Implications
For the proposal of a vaccination strategy, it will be assumed that the initial conditions are known and that
vaccines are delivered every
T days:
Then, the vaccines stock
can be defined as:
where
, and
is the Heaviside function. Therefore,
is a picewise continuous with jump discontinuities.
During the period , the vaccination strategy will be characterized by two main points:
A continuous vaccination of about
: considering the expression for
in (
24), the term
is a good candidate for this propose, since
Therefore,
. It is possible to define
by imposing a condition related to the vaccination response; that is, if it is desired to reach
of the vaccination final value (i.e.,
) when
, then
and consequently,
Note that
depends on
, whose value changes every period. Therefore, for every time period,
must be updated before simulating the solution of system (
3).
As fast as there is a new vaccine delivery, an increment in vaccination is desired, and this rise must disappear with respect to time. This condition can be accomplished if one chooses
such that
. Thus, the term containing
in Equation (
22) turns into
Regarding , different values will be used to see its implications.
Supposing that
days and
, it leads to
. With the initial conditions
,
,
,
, and
, a simulation has been carried out for
, and
has been updated for each period; see
Table 5 for more information. After repeating this procedure for different
values, the results have been depicted (blue lines) in
Figure 8. To show the advantage of vaccination strategy,
Figure 8 also shows how infected and hospitalized individuals as well as the cumulative deaths evolve with respect to time when no vaccines are implemented (red lines). For each
,
, and
combination, it has been checked that the number of vaccinated individuals never exceeds the vaccines stock nor the number of susceptible individuals, so Theorem 2 is not contradicted.
Regarding
I,
H, and the total number of deaths, during the first few days, the vaccination strategies do not make a big difference. At the end, the vaccinations drive the mentioned subpopulations to zero, and consequently, the number of cumulative deaths is stabilized, since
; see
Table 5. Without vaccination
, hence the increase of the total number of deaths observed during the last period.
With the chosen and values, follows the desired behavior; when there is a vaccine delivery, there is an increment in vaccinations and thereafter it tends to . An increase of can modify this tendency when there is a rise in the number of hospitalized people; during the time period 0–30 days, when , the number of vaccinations continuous growing instead of drifting to a constant value. With respect to , at the beginning of each time period, the stock is greater than in the previous one (i.e., ), since not all the vaccines are being implemented.
5. Discussion
An epidemic model consisting of susceptibles, non-seriously infected, hospitalized, and recovered subpopulations (SIHR) has been built up. In addition, a vaccination policy has also been included where the vaccination and hospitalization state variables are fed back and controlled with the tuning parameters and , respectively. For more design freedom, another parameter () and a time-varying function () have been added. This vaccination design is advantageous, since data regarding hospitalized people and vaccinated individuals are usually known.
The positivity of the model without vaccination has been proved. Considering the model with vaccination, a very impulsive one could provoke the negativeness of either the susceptible subpopulation or the vaccination stock. Therefore, in Theorem 2, the mentioned characteristic is taken into account, and the positivity of this model is proved under certain conditions. Then, and as well as the conditions for their existence have been obtained. From the analysis of the local asymptotic stability of , the analytical expression for has been calculated, which gives broad information; when , in addition to being locally asymptotically stable, is not reachable. Contrary, when , is not locally asymptotically stable and is reachable. In addition, considering that the vaccination tuning parameters are constant with respect to time, the expression for , so the local asymptotic stability of is ensured, has been calculated. has been proved to be globally stable when is below . The local asymptotic stability of has also been proved; when a vaccination is being applied, the stability is conditioned by some constrains, whereas the model without vaccines is not.
To study how the COVID-19 epidemic might evolve under different initial conditions, vaccination strategies, or average contact rate, data concerning population such as the natural death rate, birth rates, etc. have been gathered from official sources, and some features regarding COVID-19 disease (i.e., transmission probability) have also been collected. Firstly, an outbreak has been simulated when the conditions for the global stability are met, and even if a good result regarding the total number of deaths is attained, the number of contact rates needed for that purpose (approximately an average of contacts per day) is very low and not affordable from a real point of view. Then, two simulations have been carried out to study the implications of the local asymptotic stability/instability; with a slight difference in vaccination and the same contact rates, during the first days, the differences between both is not prominent. However, after some time, one reaches while the other reaches , which implies that the deaths stop and continue constantly, respectively. Therefore, it is concluded that from the public health point of view, governments should always try to maintain , since by increasing the total vaccines administration 1.3% (from to in two years), the total number of deaths during that time decreases approximately four times (from to ).
Regarding the stability of the endemic equilibrium point, two cases have been considered: and . In the first particular case, the stability is conditioned by the positivity of two terms (the equilibrium point is locally asymptotically stable when both terms are positive) derived from the Routh–Hurwitz criterion, which has been proved to be positive if is reachable. When vaccines are being implemented, from Rouche’s theorem, it has been concluded that it is sufficient if fulfills two conditions to ensure the local asymptotic stability of . In both cases, numerical simulations have been carried out with a regular vaccination, and it has been found that both cases match with the theoretical results.
The response of the hospitalized subpopulation with respect to different vaccination functions (constant, proportional, and state feedback) has been simulated, and the results have been compared; with the state feedback vaccination, its behavior can be easily adapted to a desired performance (e.g., a vaccination without peaks).
The last simulation could stand for a real-case example; a vaccines delivery function () and vaccines stock function () have been included, and , , and have been chosen to match a vaccine specification; that is, has been chosen as a periodic pulse-type function, so an impulsive vaccination is obtained, and the selected and drive vaccination to a desired steady value after a while. Several simulations have been carried out for different values, and the results have been compared with the case in which no vaccination is applied. When vaccines are implemented, the infected subpopulation tends to zero, and consequently the number of deaths do, too. With greater values of , better results are shown since the number of vaccinated individuals at the beginning increases. For instance, the last simulation shows that with a periodic vaccines delivery, it is possible to eradicate the disease without administrating all the doses (i.e., the vaccines stock increases with respect to time), which is very important in case future outbreaks are anticipated (i.e., fall and winter are characterized by the prevalence of outbreaks).
Overall, the epidemic model with vaccination has been proved to be coherent in the sense that its solution and the equilibrium points are non-negative under certain conditions. Moreover, from the global and local asymptotic stability analysis of , the conditions that force the disappearance of the disease have been obtained. With the Routh–Hurwitz criterion, the conditions for the local asymptotic stability of have been also derived. Then, based on real data, numerical values have been assigned to the parameters, and different simulations have been carried out; it has been found that the conditions for global stability are difficult to accomplish, whereas for the local asymptotic stability, they are quite affordable. In the last simulation, it has been demonstrated that the implemented vaccination, as a result of the way in which it is built up and its free-design parameters and time-dependent function, can be easily adapted to desired requirements.
Taking into account that usually the death rate differs notoriously between age groups, our future work will consider this model as the baseline for an age-group epidemic model. Then, the optimal control problem can be formulated, where the performance measure could be defined to follow the minimum control effort (i.e., least waste of resources), to transfer a system from an arbitrary initial state to the target state (the ) in minimum time and to maximize the deviation of the hospitalized subpopulation state from the hospital’s bed limit.