An SIRS Epidemic Model Supervised by a Control System for Vaccination and Treatment Actions Which Involve First-Order Dynamics and Vaccination of Newborns

Abstract

This paper analyses an SIRS epidemic model with the vaccination of susceptible individuals and treatment of infectious ones. Both actions are governed by a designed control system whose inputs are the subpopulations of the epidemic model. In addition, the vaccination of a proportion of newborns is considered. The control reproduction number $R_c$ of the controlled epidemic model is calculated, and its influence in the existence and stability of equilibrium points is studied. If such a number is smaller than a threshold value ${\overline{R}}_c$, then the model has a unique equilibrium point: the so-called disease-free equilibrium point at which there are not infectious individuals. Furthermore, such an equilibrium point is locally and globally asymptotically stable. On the contrary, if $R_c>{\overline{R}}_c$, then the model has two equilibrium points: the referred disease-free one, which is unstable, and an endemic one at which there are infectious individuals. The proposed control strategy provides several free-design parameters that influence both values $R_c$ and ${\overline{R}}_c$. Then, such parameters can be appropriately adjusted for guaranteeing the non-existence of the endemic equilibrium point and, in this way, eradicating the persistence of the infectious disease.

1. Introduction

The propagation of epidemic diseases within a host population has been studied since several decades ago. Kermack and McKendrick developed one of the first works in the subject [1]. They proposed an SIR epidemic model where the host population is split in three categories depending on the status of the individuals with respect to the disease. In such a context, S, I, and R denote, respectively, the susceptible, infectious, and recovered subpopulations. Later, a lot of models have been proposed and analysed in the specialised literature. Such models consider some additional population categories and/or control actions for eradicating or, at least, diminishing the effects of the disease in the host population [2,3,4,5]. In this sense, the models can include the exposed subpopulation E composed by individuals without symptoms of the disease and without the capacity of transmitting the infection to a susceptible individual after a contagion. Such a model is referred to as the SEIR model [6]. The more usual control actions are the vaccination of susceptible individuals and the treatment of infectious individuals by antivirals, antibiotics, or other medicaments [7,8]. Additionally, other control actions such as quarantine of infectious individuals have also been proposed [9]. These control actions give place to include other subpopulations in the models as vaccinated, treated, and quarantined ones [10,11,12,13]. For instance, the authors in [13] propose a SVEIR model to analyse the impact of vaccination in the control of spread of poliomyelitis. Moreover, the control actions can be applied in continuous-time or impulsive ways during the vaccination campaigns and/or treatment procedures [14,15,16]. The social distancing is another alternative way to fight against the propagation of infectious diseases [17]. Such a measure may be crucial when there are asymptomatic infectious individuals and/or neither vaccines nor medicaments are available to maintain the disease propagation under control. The early phase of the COVID-19 pandemic is a clear example of such a fact [18]. In such situations, the susceptible subpopulation can be split in several categories according to its preventive care and responsible behaviour to avoid the disease contagion. For instance, the study in [19] considers two susceptible subpopulations according to its risk aversion. One of the categories contains individuals with self-protection awareness, and the other one is composed by people without such an awareness. Recently, new approaches from game theory have been proposed for analysing the spreading of infectious diseases [20,21]. The studies in [20] introduce the concept of a vaccination game so as to evaluate provisions other than vaccination, including protective measures as mask wearing. They also analyse the efficiency of quarantine compared with that of isolation policies or the efficiency of preventive versus late vaccination. The work in [21] proposes a double-layer game structure of vaccination and treatment. The vaccination game considers whether the vaccine is accepted or declined by the individual, while the treatment game depends on the antiviral resistance evolution and prescribing behaviour. In this way, the vaccination game deals with the presence of anti-vaccine behaviour, while the treatment game considers the antibiotic overuse.

In this paper, a controlled epidemic model with vaccination of newborns, vaccination of the susceptible, and treatment of the infectious individuals, as control actions, is proposed. The model is composed by a basic SIRS epidemic model and a first-order continuous-time control subsystem, which is based on the feedback of the state variables. Namely, the inputs of such a control subsystem are the susceptible, infectious, and recovered subpopulations so that the control subsystem acts under the knowledge of the state of the disease at each time. The control subsystem provides several free-design parameters that can be adjusted to eradicate the persistence of the disease within the host population. Namely, an appropriately adjustment of the control parameters guarantees the non-existence of the endemic equilibrium (EE) point of the controlled SIRS model. Such a fact implies the existence of a unique equilibrium point, namely, the disease-free equilibrium (DFE) point, which is a very advantageous tool to asymptotically eradicate the disease. Moreover, the local and global asymptotic stability of this DFE point is analytically proved under appropriate conditions relative to the adjustment of the control parameters. Under such conditions, the disease is eradicated from the host population after a transient period of propagation from the initial state until the DFE point is reached. Several modified SIRS-type models have been proposed and analysed in the epidemiological research. Some models study the dynamics of such models under the influence of control actions such as vaccinations of susceptible and/or the treatment of infectious individuals. Both control actions are applied either in a time continuous or impulsive way, and the intensity of each control action is usually proportional to the susceptible and/or infectious individuals. The SIRS model in our paper also proposes the combined application of vaccination to the susceptible and the treatment to the infectious subpopulation. The main novelties of the current paper related to the background literature are (i) the use of first-order dynamics in the vaccination and treatment controls and (ii) the availability of additional free-design control parameters to shape such vaccination and treatment actions. The first novelty provides some additional parameters derived from the fact that the vaccination and treatment actions are provided by a control subsystem instead of being proportional to the susceptible and/or infectious individuals. Concretely, two parameters arise from the dynamics of the control subsystem, and then, the control actions can be designed with two more free-design parameters than in the conventional research. Moreover, the control subsystem uses the information of the state of the propagation of the disease since the subpopulation variables are used as the inputs of such a subsystem. In this way, a feedback control strategy is used for obtaining the vaccination and treatment variables. In summary, the control designer can make use of the additional free-design parameters to shape the vaccination and treatment actions in a desired way, which is the second main novelty, to achieved the required objectives. One of the main objectives is to achieve the eradication of the disease or, at least, to minimise its influence within the host population. In this context, the control parameters governing the intensity of the vaccination of newborns and susceptible individuals influence on the control reproduction number $R_c$ of the controlled SIRS model. Moreover, those parameters shaping the intensity of the treatment of the infectious individuals influence on a threshold value ${\overline{R}}_c$ of interest since if $R_c<{\overline{R}}_c$, then the EE point does not exist, so that the DFE point is the unique one. Furthermore, such a DFE point is locally and globally asymptotically stable under that condition. In this way, an appropriate choice of the control parameters associated with the vaccination actions allows the control reproduction number to reduce with respect to its value in absence of vaccination. Namely, $R_c$ decreases by increasing the intensity of the vaccination. On the other hand, an appropriate choice of the control parameters associated with the treatment action allows the threshold value ${\overline{R}}_c$ to be strictly larger than the $R_c$ value of the controlled SIRS model to guarantee the inexistence of the EE point and then the existence of a unique globally asymptotically stable DFE point. Namely, ${\overline{R}}_c$ increases by increasing the intensity of the treatment action. In this context, the infectious disease can be extinguished under an appropriate choice of the design parameters of the control laws. In summary, two options are available for such a purpose. The first one is increasing the intensity of vaccination in order to reduce the control reproduction number $R_c$, and the second one is increasing the intensity of treatment in order to augment the threshold value ${\overline{R}}_c$. Obviously, an appropriate solution to guarantee the extinction of the disease can be obtained by combining both options so that $R_c<{\overline{R}}_c$.

The rest of the paper is organized as follows. Section 2 describes the basic epidemic SIRS model as well as the subsystem providing the both proposed control actions, namely, the vaccination of the susceptible and treatment of the infectious subpopulations, respectively. The positivity of the controlled model, composed by combining the basic SIRS model with the control subsystem, is analysed. In addition, the study of the equilibrium points, as the existence as the stability properties, of the controlled model is dealt with in this section. Concretely, the conditions to be satisfied by the free-design control parameters in order to guarantee the non-existence of the EE point and then, the local and global asymptotic stability of the DFE point, which is the unique equilibrium point under such conditions, are established and mathematically proved. Finally, Section 3 illustrates the theoretical results by some simulation examples. An extended study of the influence of the control parameters in the dynamics of the controlled SIRS model is done. The results obtained with the proposed controlled SIRS model are compared with those obtained by an uncontrolled SIRS model. In addition, the results of the proposed model are compared with an SIRS model with only a control action, either the vaccination of the susceptible subpopulation or the treatment of the infectious one. These comparisons are interesting from the viewpoint of the available resources relative to the existence of vaccines and/or medicaments to fight against the propagation of the disease.

2. SIRS Epidemic Model under Vaccination and Treatment Controls

An SIRS epidemic model described by the following equations:

(1)$\dot{S}\left(t\right)=b\left(1-q\right)-\beta \frac{S\left(t\right)I\left(t\right)}{N\left(t\right)}+\rho R\left(t\right)-\mu S\left(t\right)-v\left(t\right)$

(2)$\dot{I}\left(t\right)=\beta \frac{S\left(t\right)I\left(t\right)}{N\left(t\right)}-\left(\mu +\alpha +\gamma \right)I\left(t\right)-t_r\left(t\right)$

(3)$\dot{R}\left(t\right)=bq+\gamma I\left(t\right)-\left(\mu +\rho \right)R\left(t\right)+v\left(t\right)+t_r\left(t\right)$

with an initial condition given by $S\left(0\right)\ge 0$, $I\left(0\right)\ge 0$ and $R\left(0\right)\ge 0$ is proposed. The model considers the whole host population split into three categories depending on the state of the individuals with respect to the infectious disease, namely, susceptible (S), infectious (I), and recovered (R) subpopulations. Moreover, the model considers a constant recruitment rate for the population, which is represented by b [22,23]. The model includes three control actions: (i) a constant vaccination of newborn individuals defined by the control parameter $q\in \left[0, 1\right]$, (ii) a vaccination of the susceptible subpopulation into a vaccination rate $v\left(t\right)$, and (iii) a treatment of the infectious subpopulation into a treatment rate $t_r\left(t\right)$.

The variables $v\left(t\right)$ and $t_r\left(t\right)$ are the output variables of a first-order control system whose inputs are the SIRS model state variables. Then, both control actions are synthesised by the feedback of the epidemic model state. The equations of the control system providing both actions are as follows:

(4)$\dot{v}\left(t\right)=-c_{1}v\left(t\right)+c_2 S\left(t\right)+c_3 I\left(t\right)+c_4\frac{S\left(t\right)I\left(t\right)}{N\left(t\right)}$

(5)$\dot{t_r}\left(t\right)=-c_5 t_r\left(t\right)+c_6 I\left(t\right).$

Concretely, $v\left(t\right)$ and $t_r\left(t\right)$ are the number of vaccinated and treated individuals per time unit, respectively. Equations (4) and (5) regulate respectively the amount of vaccines and medicaments to be daily applied according to the time evolution of the disease propagation. In this context, if $\dot{v}\left(t\right)>0$ ($\dot{t_r}\left(t\right)>0$), then the amount of vaccines (medicaments) to be applied in the current day is larger than those applied in the previous one. Otherwise, if $\dot{v}\left(t\right)<0$ ($\dot{t_r}\left(t\right)<0$), then the amount of vaccines (medicaments) to be applied in the current day is smaller than those applied in the previous one. The time evolution of the amount of vaccines and medicaments daily applied to control the propagation of the disease depends on the current state of the epidemics, as it can be observed in (4) and (5) from the fact that the subpopulations $S\left(t\right)$ and $I\left(t\right)$ act as input in the control subsystem. Equations (1)–(5) compose the controlled epidemic model defined by two sets of parameters. The first set includes the parameters associated with the transmission of the disease and the features of the host population, namely:

• b is the natural birth rate of the host population;

• $\mu$ is the natural death rate of the host population;

• $\beta$ is the transmission rate of the disease within the host population;

• $\rho$ is the immunity loss rate within the recovered subpopulation, whose individuals become susceptible to the disease after losing the immunity;

• $\alpha$ is the death rate by causes related to the disease;

• $\gamma$ is the recovery rate of the infectious subpopulation.

The second set of parameters is associated to the control actions, namely:

• $q\in \left[0, 1\right]$ is the proportion of newborn individuals who are vaccinated;

• $c_i$, for $i\in \left\{1, 2, 3, 4\right\}$, are the parameters for designing the vaccination of the susceptible subpopulation. Such control parameters allow us to weight the vaccination rate according to the state of the disease propagation considering the number of susceptible individuals S, infectious ones I, and/or the probability of contacts between them $\frac{S\left(t\right)I\left(t\right)}{N\left(t\right)}$ at each time t. Such parameters provide the availability of giving more importance to one term of (4) against the other ones in the design of the law for the vaccination $v\left(t\right)$. The unit of the parameter $c_1$ is (time)⁻¹, usually (day)⁻¹, and that of the parameters $c_i$, for $i\in \left\{2, 3, 4\right\}$ is (time)⁻² for coherence in (4).

• $c_i$, for $i\in \left\{5, 6\right\}$, are the parameters for designing the treatment of the infectious subpopulation. Such control parameters shape the law for the treatment $t_r\left(t\right)$ according to the number of infectious individuals. However, the values allowed for such parameters are constrained for the potency of the available medicaments. In this sense, a larger $\frac{c_6}{c_5}$ is less the recovery time interval for the treated infectious individuals. The unit of the parameter $c_5$ is (time)⁻¹, usually (day)⁻¹, and that of the parameters $c_6$ is (time)⁻² for coherence in (5).

All parameters are non-negative. The dynamics of the whole host population $N\left(t\right)=S\left(t\right)+I\left(t\right)+R\left(t\right)$ can be obtained by summing up Equations (1)–(3). In this way:

(6)$\dot{N}\left(t\right)=b-\mu N\left(t\right)-\alpha I\left(t\right).$

The nature of the epidemic models requires the non-negativity of their solutions, so an analysis of the model positivity is developed in the following.

2.1. Positivity of the Controlled SIRS Epidemic Model

The result below establishes the non-negativity of all the subpopulations of the controlled model as well as the non-negativity of the vaccination and treatment control efforts under a set of sufficient conditions on the control parameters. The proof is written in Appendix A.

Theorem 1 (positivity of the model).

The solution of the model (1)–(6) is non-negative for all time and for any initial condition such that$S\left(0\right)\ge 0$,$I\left(0\right)\ge 0$,$R\left(0\right)\ge 0$,$v\left(0\right)=0$, and$t_r\left(0\right)=0$provided that the control parameters$c_i\ge 0$, for$i\in \left\{1, 2, \dots , 6\right\}$, and$q\in \left[0, 1\right]$are chosen such that:

$\begin{array}{lll}\left(i\right) c_1>\mu +\beta +2\sqrt{c_2+c_4} & ;& \left(ii\right) q\ge q_0=1-\frac{\mu +\beta }{b}S\left(0\right)\\ \left(iii\right) c_3\le \frac{b\left(1-q\right)\left(\mu +\beta \right)}{I_{max}}& ;& \left(iv\right) c_5>\mu +\alpha +\gamma +2\sqrt{c_6} ,\end{array}$

where$I_{max}=\underset{0\le t<\infty }{max}\left\{I\left(t\right)\right\}$.

Remark 1.

• (i) 

  The conditions (i)–(iv) of Theorem 1 are only sufficient conditions, since the solutions of the model can be non-negative even if some of these conditions are not satisfied by the control parameters.

• (ii) 

  The birth rate of a host population $b$ is close to its mortality rate $\mu N\left(t\right)$ for any $t\ge 0$ under normal conditions (in absence of a lethal disease). Then, $b-\mu N\left(0\right)$ is close to zero at the beginning of the propagation of an infectious disease. Then, typically $b-\left(\mu +\beta \right)N\left(0\right)\approx -\beta N\left(0\right)<0$. Moreover, in the first stage of any epidemic disease propagation, the infectious and recovered subpopulations are much smaller than the susceptible one. Then, $N\left(0\right)\approx S\left(0\right)$ and, as a consequence, the condition (ii) is satisfied for any $q\in \left[0, 1\right]$.

• (iii) 

  The condition (iii) of Theorem 1 depends on the maximum value reached by the infectious subpopulation during the propagation of the disease. Such a value cannot be known ‘a priori’, and then, one cannot appropriately choose the values of the parameters $c_3$ and $q$ to satisfy it. However, such a condition is fulfilled if $c_3=0$ for any $q\in \left[0, 1\right]$. Then, from continuity arguments, the condition is fulfilled for $c_3\in \left[0, {\overline{c}}_3\right]$ with ${\overline{c}}_3={\overline{c}}_3\left(q,I_{max}\right)=\frac{b\left(1-q\right)\left(\mu +\beta \right)}{I_{max}}$. In summary, a value for $c_3$ small enough has to be chosen in order to satisfy the condition (iii) if a vaccination provided by (4) is proposed. Note that if $q=1$, which means a vaccination of all the newborns, then $c_3=0$ has to be taken to guarantee the non-negativity of all the model variables for all the time.

• (iv) 

  In a real situation, the control actions to fight against an epidemic outbreak are taken after the disease is detected in the infectious individuals. Then, the constraint $v\left(0\right)=t_r\left(0\right)=0$ in the initial condition for guaranteeing the non-negativity of the variables of the controlled epidemic model is coherent with such a fact.

Corollary 1.

The feasible region $\Gamma$ defined as:

(7) $\Gamma =\left\{\left[S\left(t\right) I\left(t\right) R\left(t\right) v\left(t\right) t_r\left(t\right)\right]\in {ℝ}_{0+}^5|\begin{array}{c}0\le min\left\{N\left(0\right),\frac{b}{\mu +\alpha }\right\}\le N\left(t\right)\le max\left\{N\left(0\right),\frac{b}{\mu }\right\}\\ v\left(t\right)\ge 0\\ t_r\left(t\right)\ge 0\end{array}\right\}$

is positively invariant for models (1)–(5) under the conditions of Theorem 1, where $N\left(t\right)=S\left(t\right)+I\left(t\right)+R\left(t\right)$ and ${ℝ}_{0+}^5$ denotes the non-negative five-dimension hyper-plane.

Proof. 

From (6), it follows that:

(8)$b-\left(\mu +\alpha \right)N\left(t\right)\le \dot{N}\left(t\right)\le b-\mu N\left(t\right)$

from the fact that $0\le I\left(t\right)\le N\left(t\right)$ $\forall t\ge 0$ as a consequence of Theorem 1. Then, by direct calculations from (8), one obtains:

(9)$\frac{b}{\mu +\alpha }+\left(N\left(0\right)-\frac{b}{\mu +\alpha }\right)e^{-\left(\mu +\alpha \right)t}\le N\left(t\right)\le \frac{b}{\mu }+\left(N\left(0\right)-\frac{b}{\mu }\right)e^{-\mu t}.$

From (4), it follows that $\dot{v}\left(t\right)\ge -c_{1}v\left(t\right)$ by taking into account that $S\left(t\right)\ge 0$, $I\left(t\right)\ge 0$ and $N\left(t\right)\ge 0$ $\forall t\ge 0$ as a consequence of Theorem 1. Then, $v\left(t\right)\ge v\left(0\right)e^{-c_{1}t}\ge 0$. Finally, from (5), it follows that $\dot{t_r}\left(t\right)\ge -c_5{t}_r\left(t\right)$ by taking into account that $I\left(t\right)\ge 0$ $\forall t\ge 0$ as a consequence of Theorem 1. Then, $t_r\left(t\right)\ge t_r\left(0\right)e^{-c_{5}t}\ge 0$. Equation (9) and the results $v\left(t\right)\ge 0$ and $t_r\left(t\right)\ge 0$ $\forall t\ge 0$ lead to the conclusion that solutions for the model (1)–(5) with any initial condition within $\Gamma$ persist in $\Gamma$ for all time. □

From the positivity result of Theorem 1, the following assumption is established for the rest of the paper.

Assumption 1.

A choice of the control parameters satisfying the conditions of Theorem 1 is supposed in order to guarantee the positivity of the controlled SIRS epidemic model.

Remark 2.

Assumption 1 mathematically guarantees that the model subpopulations as well as the control efforts do not take negative values for all time, as the nature of an epidemic model needs to be well defined. Such a fact justifies the adoption of such an assumption.

2.2. Control Reproduction Number and Equilibrium Points of the Controlled SIRS Model

The controlled epidemic model given by (1)–(5) possesses two equilibrium points. One of them is a DFE point, and the other one is an EE point. They are obtained by setting $\dot{S}\left(t\right)=\dot{I}\left(t\right)=\dot{R}\left(t\right)=\dot{v}\left(t\right)=\dot{t_r}\left(t\right)=0$, since the equilibrium points are those at which the model variables do not change with time. In this way and by direct calculations, one obtains that the subpopulations and the values of the vaccination and treatment efforts at the DFE point are given by:

(10)$$\begin{gathered} S_{DFE}=\frac{b{c}_1\left[\mu \left(1-q\right)+\rho \right]}{\mu \left[c_1\left(\mu +\rho \right)+c_2\right]} ; I_{DFE}=0 ; R_{DFE}=\frac{b\left(c_{1}q\mu +c_2\right)}{\mu \left[c_1\left(\mu +\rho \right)+c_2\right]} \\ v_{DFE}=\frac{b{c}_2\left[\mu \left(1-q\right)+\rho \right]}{\mu \left[c_1\left(\mu +\rho \right)+c_2\right]} ; t_r{}_{DFE}=0. \end{gathered}$$

One can define the control reproduction number $R_c$ for the model as the number of new infections that an infectious individual produces in a population at the DFE point. Such a number can be calculated by using the next-generation method [24]. For such a purpose, Equations (1)–(3) related to the dynamics of the epidemiological compartments of the controlled model (1)–(5) can be, equivalently, rewritten as:

(11)$\dot{x}=ℱ\left(x\right)-\mathcal{W}\left(x\right)$

where $x={\left[\begin{array}{ccc}S& I& R\end{array}\right]}^T$ and $ℱ\left(x\right)$ and $\mathcal{W}\left(x\right)$ are given by:

(12)$ℱ={\left[\begin{array}{ccc}0& \beta \frac{SI}{N}& 0\end{array}\right]}^T ; \mathcal{W}=\left[\begin{array}{c}-b\left(1-q\right)+\beta \frac{SI}{N}-\rho R+\mu S+v\\ \left(\mu +\alpha +\gamma \right)I+t_r\\ -bq-v-\gamma I+\left(\mu +\rho \right)R-t_r\end{array}\right].$

Each component of $ℱ\left(x\right)$ is the rate of appearance of new infections in the corresponding compartment. The new infections only appear in the infectious compartment; therefore, only the second component of $ℱ\left(x\right)$ is nonzero. On the other hand, the components of $\mathcal{W}\left(x\right)$ are the difference between the output and input transition rates corresponding, respectively, to the susceptible, infectious, and recovered compartments. The control reproduction number can be calculated by analysing the dynamics of the subsystem composed by the infected compartments. The proposed controlled model only has an infected compartment, namely, the infectious subpopulation $I$. The dynamics of such an infected subsystem around the DFE point is given by the both entries (2,2) of the matrices $F={\frac{\partial ℱ}{\partial x}|}_{x=x_{DFE}}\in {ℝ}^{3x3}$ and $W={\frac{\partial \mathcal{W}}{\partial x}|}_{x=x_{DFE}}\in {ℝ}^{3x3}$, where $x_{DFE}={\left[\begin{array}{ccc}{S}_{DFE}& I_{DFE}& R_{DFE}\end{array}\right]}^T$. By direct calculations, one obtains that such entries are, respectively, $f_{22}=\beta \frac{S_{DFE}}{N_{DFE}}$ and $w_{22}=\mu +\alpha +\gamma$. Then, the next generation matrix, which results a scalar for the current models, is $f_{2,2}·{\left(w_{22}\right)}^{-1}$. Finally, the control reproduction number under control efforts is obtained as the spectral radius of such a matrix, namely:

(13)$R_c=\sigma \left(f_{2,2}·{\left(w_{22}\right)}^{-1}\right)=\frac{f_{22}}{w_{22}}=\frac{\beta S_{DFE}}{\left(\mu +\alpha +\gamma \right)N_{DFE}}=\frac{\beta c_1\left[\mu \left(1-q\right)+\rho \right]}{\left(\mu +\alpha +\gamma \right)\left[c_1\left(\mu +\rho \right)+c_2\right]}$

where $\sigma \left(M\right)$ denotes the spectral radius of the matrix $M$ and the fact that $N_{DFE}=S_{DFE}+I_{DFE}+R_{DFE}=\frac{b}{\mu }$ and expressions in (10) have been used. Note that if the vaccination of newborn individuals is not applied, i.e., $q=0$, while the vaccination of susceptible individuals is given by (4) with $c_2=0$, i.e., without the forced term depending on S, then the control reproduction number of the current controlled model is equal to the basic reproduction number $R_0=\frac{\beta }{\mu +\alpha +\gamma }$ of a basic SIRS epidemic model. On the contrary, if $c_2\ne 0$ and/or $q\ne 0$, then the control reproduction number is smaller than the basic reproduction number, i.e., $R_c<R_0$, since the control parameters are non-negative by definition. Such a fact is key to achieve the global stability of the DFE point of the controlled epidemic model (1)–(5), by means of an appropriate choice of the control parameters $c_2$ and $q$, in a situation where the DFE point of an SIRS model without vaccination is not globally stable. Note also that the control reproduction number depends on neither $c_5$ nor $c_6$, i.e., the application of a treatment to the infectious subpopulation does not have influence on such a number.

On the other hand, one obtains that the subpopulations and the values of the vaccination and treatment efforts at the EE point are given by:

(14)$$\begin{gathered} \begin{array}{ccc}{S}_{EE}=\frac{b{k}_1\left(k_2{k}_3{R}_c+k_1{c}_4\right)}{k_2{c}_5{R}_c\left(k_2{k}_4{R}_c+k_1{k}_5\right)}& ;& I_{EE}=\frac{b{k}_6\left[k_7\left(R_c-1\right)-c_6\right]}{k_2{k}_4{R}_c+k_1{k}_5}\\ R_{EE}=\frac{b\left(k_2^2{k}_8{c}_5{R}_c^2+k_1{k}_2{k}_9{R}_c-k_1^2{c}_4\right)}{k_2{c}_5{R}_c\left(k_2{k}_4{R}_c+k_1{k}_5\right)}& ;& t_r{}_{EE}=\frac{b{k}_6{c}_6\left[k_7\left(R_c-1\right)-c_6\right]}{c_5\left(k_2{k}_4{R}_c+k_1{k}_5\right)}\end{array} \\ {v}_{EE}=\frac{b\left(k_2^2{c}_3{c}_5^2\left[\mu \left(1-q\right)+\rho \right]R_c^2+k_1{k}_2{k}_{10}{R}_c-k_1^2{c}_4\left(\mu +\rho \right)\right)}{k_2{c}_5{R}_c\left(k_2{k}_4{R}_c+k_1{k}_5\right)} \end{gathered}$$

where:

(15)$$\begin{gathered} \begin{array}{ll}{k}_1=c_1\left[\mu \left(1-q\right)+\rho \right]\left[c_5\left(\mu +\alpha +\gamma \right)+c_6\right] ; & k_2=\left(\mu +\alpha +\gamma \right)\left[c_1\left(\mu +\rho \right)+c_2\right]\\ k_3=c_1\left[c_5\left(\mu +\gamma +\rho +q\alpha \right)+c_6\right]+c_3{c}_5 ; & k_4=c_1\left(c_5\left[\mu \gamma +\left(\mu +\alpha \right)\left(\mu +\rho \right)\right]+\mu c_6\right)+\mu c_3{c}_5\\ k_5=\mu c_4-\alpha \left[c_1\left(\mu +\rho \right)+c_2\right] ; & k_6=c_1\left[\mu \left(1-q\right)+\rho \right]\left[c_1\left(\mu +\rho \right)+c_2\right]\\ k_7=c_5\left(\mu +\alpha +\gamma \right) ; & k_8=c_1\left(c_5\left[\gamma +q\left(\mu +\alpha \right)\right]+c_6\right)+c_3{c}_5\\ k_9=c_5\left[c_2-c_3+c_4-c_1\left(\alpha q+\gamma \right)\right]-c_1{c}_6& \end{array} \\ {k}_{10}=c_2\left[c_5\left(\mu +\gamma +\rho +q\alpha \right)+c_6\right]-c_3{c}_5\left(\mu +\rho \right)+c_4{c}_5\left[\mu \left(1-q\right)+\rho \right]. \end{gathered}$$

The following very important result from the disease eradication viewpoint is proved.

Theorem 2 (Non-existence of the EE point).

If the control parameters are chosen such that:

(16) $R_c<{\overline{R}}_c=1+\frac{c_6}{c_5\left(\mu +\alpha +\gamma \right)} \iff R_0<{\overline{R}}_0=\frac{\left[c_5\left(\mu +\alpha +\gamma \right)+c_6\right]\left[c_1\left(\mu +\rho \right)+c_2\right]}{c_1{c}_5\left(\mu +\alpha +\gamma \right)\left[\mu \left(1-q\right)+\rho \right]}$

then the EE point does not exist.

Proof. 

The condition (16) can be written as $R_c<1+\frac{c_6}{k_7}$ from the expression of $k_7$ in (15). Then, one obtains that $N_I=b{k}_6\left[k_7\left(R_c-1\right)-c_6\right]<0$, where $N_I$ is the numerator of $I_{EE}$. Then, $I_{EE}<0$ unless that $D_I=k_2{k}_4{R}_c+k_1{k}_5<0$, where $D_I$ is the denominator of $I_{EE}$. Suppose that $D_I<0$ so that $I_{EE}>0$. Then, $D_S=k_2{c}_5{R}_c\left(k_2{k}_4{R}_c+k_1{k}_5\right)=k_2{c}_5{R}_c{D}_I<0$ from the fact that $c_5$, $k_2$, and $R_c$ are strictly positive, where $D_S$ denotes the denominator of $S_{EE}$. Then, $S_{EE}<0$, since $N_S=b{k}_1\left(k_2{k}_3{R}_c+k_1{c}_4\right)>0$ since $b$, $c_4$, $k_1$, $k_2$, $k_3$ and $R_c$ are strictly positive, where $N_S$ denotes the numerator of $S_{EE}$. In summary, one of $S_{EE}$ or $I_{EE}$ is negative under the condition (16), which implies the non-existence of the EE point. □

Remark 3.

The proposed SIRS epidemic model defined by (1)–(3) has a globally asymptotically stable EE point whenever $R_0>1$ , as is always the case in SIR-like epidemic models with demography [25]. The proposed control actions, namely, the vaccination of newborns, the vaccination of the susceptible subpopulation, given by (4), and the treatment of infectious individuals, given by (5), eliminate the existence of such an EE point under a choice of the control parameters, satisfying the condition (16). In this way, the SIRS epidemic model coupled with the control actions only has the DFE point, which is a crucial result in this paper.

2.3. Local Stability of the Disease-Free Equilibrium Point

The study of the local stability of a nonlinear system around an equilibrium point can be done by analysing the eigenvalues of the Jacobian matrix corresponding to the linearisation of the system around such a point. For such a purpose, the controlled model (1)–(5) can be, equivalently, rewritten as ${\dot{x}}_c=f\left(x_c\right)$ where $x_c={\left[\begin{array}{cc}\begin{array}{cc}S& I\end{array}& \begin{array}{ccc}R& v& t_r\end{array}\end{array}\right]}^T$ is the extended state vector of the controlled SIRS model after including the control dynamics and:

(17)$f=\left[\begin{array}{c}\begin{array}{c}b\left(1-q\right)-\beta \frac{SI}{N}+\rho R-\mu S-v\\ \beta \frac{SI}{N}-\left(\mu +\alpha +\gamma \right)I-t_r\end{array}\\ \begin{array}{c}bq+\gamma I-\left(\mu +\rho \right)R+v+t_r\\ -c_{1}v+c_{2}S+c_{3}I+c_4\frac{SI}{N}\\ -c_5{t}_r+c_{6}I\end{array}\end{array}\right].$

The Jacobian matrix of the model around the DFE point is directly obtained as:

(18)$J_{DFE}={\frac{\partial f}{\partial x_c}|}_{x_c=x_{c,DFE}}=\left[\begin{array}{ccccc}-\mu & j_{12}& \rho & -1& 0\\ 0& j_{22}& 0& 0& -1\\ 0& \gamma & j_{33}& 1& 1\\ c_2& j_{42}& 0& -c_1& 0\\ 0& c_6& 0& 0& -c_5\end{array}\right]$

where $x_{c,DFE}={\left[\begin{array}{cc}\begin{array}{cc}{S}_{DFE}& I_{DFE}\end{array}& \begin{array}{ccc}{R}_{DFE}& v_{DFE}& t_r{}_{DFE}\end{array}\end{array}\right]}^T$ $j_{12}=-\left(\mu +\alpha +\gamma \right)R_c$, $j_{22}=\left(\mu +\alpha +\gamma \right)\left(R_c-1\right)$, $j_{33}=-\left(\mu +\rho \right)$, and $j_{42}=c_3+\frac{c_4}{\beta }\left(\mu +\alpha +\gamma \right)R_c$, and the expressions (10) and (13) for $R_c$ have been used. The following result is proved.

Theorem 3 (Local stability of the DFE point).

If the control parameters are chosen such that:

(19) $R_c<{\overline{R}}_c \iff R_0<{\overline{R}}_0$

then the DFE point is locally asymptotically stable under Assumption 1, where ${\overline{R}}_c$ and ${\overline{R}}_0$ are given in (16).

Proof. 

The characteristic equation of the linearised model around the DFE point is given by $\left|\lambda {\mathbb{I}}_5-J_{DFE}\right|=0$, where ${\mathbb{I}}_5$ denotes the 5th-order unity matrix. By direct calculations, one obtains that:

(20)$\left|\lambda {\mathbb{I}}_5-J_{DFE}\right|=\left(\lambda +\mu \right)\left({\lambda }^2+p_1\lambda +p_0\right)\left({\lambda }^2+q_1\lambda +q_0\right)$

where:

(21)$$\begin{gathered} p_1=c_5-\left(\mu +\alpha +\gamma \right)\left(R_c-1\right) ; p_0=c_6-c_5\left(\mu +\alpha +\gamma \right)\left(R_c-1\right) \\ {q}_1=\mu +\rho +c_1 ; q_0=(\mu +\rho )c_1+c_2. \end{gathered}$$

Then, the eigenvalues of the matrix $J_{DFE}$ are $-\mu$ and the roots of $P\left(\lambda \right)={\lambda }^2+p_1\lambda +p_0$ and $Q\left(\lambda \right)={\lambda }^2+q_1\lambda +q_0$. By applying the Routh–Hurwitz method to $Q\left(\lambda \right)$, one obtains that the real part of its roots is negative if $q_1>0$ and $q_2>0$. From (21), both conditions are satisfied, since the control parameters $c_1$ and $c_2$ as well as the model parameters $\mu$ and $\rho$ are positive by definition. On the other hand, the roots of $P\left(\lambda \right)$ are given by:

(22)$r_{1,2}=\frac{1}{2}\left(\left(\mu +\alpha +\gamma \right)\left(R_c-1\right)-c_5\pm \sqrt{{\left(c_5+\left(\mu +\alpha +\gamma \right)\left(R_c-1\right)\right)}^2-4c_6}\right)$

from (21), where the sign ‘+’ corresponds to the root $r_1$ and the sign “$-$” corresponds to $r_2$. The control parameters $c_5$ and $c_6$ fulfil the condition (iv) of Theorem 1 provided Assumption 1. Then, it follows that ${\left(c_5+\left(\mu +\alpha +\gamma \right)\left(R_c-1\right)\right)}^2-4c_6>0$, since $R_c>0$ from (13). In this way, both roots $r_1$ and $r_2$ are real, and $r_2<r_1$. By direct calculations from (22):

(23)$$\begin{gathered} r_1<\frac{1}{2}\left(\left(\mu +\alpha +\gamma \right)\left({\overline{R}}_c-1\right)-c_5+\sqrt{{\left(c_5+\left(\mu +\alpha +\gamma \right)\left({\overline{R}}_c-1\right)\right)}^2-4c_6}\right) \\ =\frac{1}{2}\left(\frac{c_6}{c_5}-c_5+\sqrt{{\left(c_5+\frac{c_6}{c_5}\right)}^2-4c_6}\right)=\frac{1}{2c_5}\left(c_6-c_5^2+\sqrt{{\left(c_6+c_5^2\right)}^2-4c_5^2{c}_6}\right)=0 \end{gathered}$$

by using (19), the definition of ${\overline{R}}_c$ in (16) and the fact that $c_6<c_5^2$ from the condition (iv) of Theorem 1. In summary, both roots of $P\left(\lambda \right)$ are strictly negative. Then, all the roots of the characteristic equation of the linearised model around the DFE point have a negative real part under Assumption 1, provided that the control parameters satisfy the condition (19). Then, the DFE point is locally asymptotically stable. □

Remark 4.

Note that the DFE point can be locally asymptotically stable although the control reproduction number$R_c$is larger than 1. Namely, such a property is proved if$R_c<{\overline{R}}_c$independently of the particular value of$R_c$. Indeed, the local asymptotic stability of the DFE point is guaranteed if the transition rate from the infectious subpopulation to the recovered one is larger than the transition rate from the susceptible subpopulation to the infectious one. In this context, the number$R_c$is directly proportional to the transition rate from the susceptible subpopulation to the infectious one, while${\overline{R}}_c$is directly proportional to the transition rate from the infectious subpopulation to the recovered one. In this sense, such a transition rate from the infectious to the recovered subpopulation depends on several factors: potency of medicaments, amount of material and human resources in the health-care centres to treat the infectious individuals, and so on. Then, the availability of enough resources is crucial to avoid the persistence of a disease or, at least, diminish its effects on the population by increasing the value of${\overline{R}}_c$.

2.4. Global Stability of the Disease-Free Equilibrium Point

For the purpose of analysing the global stability of the DFE point, the controlled epidemic model is rewritten as [26,27]:

(24)$\{\begin{array}{l}{\dot{X}}_n\left(t\right)=A\left(X_n\left(t\right)-X_{n,DFE}\right)+B\left(t\right)X_I\left(t\right)\hfill \\ {\dot{X}}_I\left(t\right)=C\left(t\right)X_I\left(t\right)\hfill \end{array}$

where $X_n\left(t\right)={\left[S\left(t\right)R\left(t\right)v\left(t\right)\right]}^T$, $X_I={\left[I\left(t\right)t_r\left(t\right)\right]}^T$, $X_{n,DFE}=\left[S_{DFE}{R}_{DFE}{v}_{DFE}\right]$ is the vector $X_n\left(t\right)$ at the DFE point, with its components given in (10), and:

(25)$$\begin{gathered} A=\left[\begin{array}{ccc}-\mu & \rho & -1\\ 0& -\left(\mu +\rho \right)& 1\\ c_2& 0& -c_1\end{array}\right] ; B\left(t\right)=\left[\begin{array}{cc}-\beta \frac{S\left(t\right)}{N\left(t\right)}& 0\\ \gamma & 1\\ c_3+c_4\frac{S\left(t\right)}{N\left(t\right)}& 0\end{array}\right] \\ C\left(t\right)=\left[\begin{array}{cc}\beta \frac{S\left(t\right)}{N\left(t\right)}-\left(\mu +\alpha +\gamma \right)& -1\\ c_6& -c_5\end{array}\right]. \end{gathered}$$

The following result about the global stability of the DFE point is proved.

Theorem 4 (Global stability of the disease-free equilibrium point).

The DFE point is globally asymptotically stable under Assumption 1 provided that the control parameters are chosen such that$R_c<{\overline{R}}_c$or, equivalently,$R_0<{\overline{R}}_0$and$q\in \left[0, q_c\right){\displaystyle \cap }\left[0,1\right]$so that$F\left(q\right)={\lambda }_1+\frac{\left(c_5+{\lambda }_1\right)\mu c_{1}q\beta }{\left[c_1\left(\mu +\rho \right)+c_2\right]\left({\lambda }_1-{\lambda }_2\right)}<0$, where$q_c$is such that$F\left(q_c\right)=0$while ${\lambda }_1$and${\lambda }_2$are the eigenvalues of the matrix$C_0$given by:

(26)$C_0=\left[\begin{array}{cc}\left(\mu +\alpha +\gamma \right)\left(R_c-1\right)& -1\\ c_6& -c_5\end{array}\right].$

Proof. 

The subsystem ${\dot{X}}_I\left(t\right)=C\left(t\right)X_I\left(t\right)$ can be rewritten as:

(27)${\dot{X}}_I\left(t\right)=C_0{X}_I\left(t\right)+C_1\left(t\right)X_I\left(t\right)$

where $C_0$ is given by (26), and the matrix $C_1\left(t\right)$ is:

(28)$C_1\left(t\right)=\left[\begin{array}{cc}\beta \left(\frac{S\left(t\right)}{N\left(t\right)}-\frac{S_{DFE}}{N_{DFE}}\right)& 0\\ 0& 0\end{array}\right]$

and (13) has been used. From (27), it follows that:

(29)$\left[\begin{array}{c}I\left(t\right)\\ t_r\left(t\right)\end{array}\right]=\varphi \left(t\right)\left[\begin{array}{c}I\left(0\right)\\ t_r\left(0\right)\end{array}\right]+{{\displaystyle \int }}_0^t\varphi \left(t-\tau \right)C_1\left(\tau \right)X_I\left(\tau \right)d\tau \forall t\ge 0,$

with $\varphi \left(t\right)=e^{C_{0}t}=L^{-1}\left\{{\left(s{\mathbb{I}}_2-C_0\right)}^{-1}\right\}$. By direct calculations, one obtains that:

(30)$$\begin{gathered} {\varphi }_{11}\left(t\right)=\frac{\left(c_5+{\lambda }_1\right)e^{{\lambda }_{1}t}-\left(c_5+{\lambda }_2\right)e^{{\lambda }_{2}t}}{{\lambda }_1-{\lambda }_2} ; {\varphi }_{12}\left(t\right)=-\frac{e^{{\lambda }_{1}t}-e^{{\lambda }_{2}t}}{{\lambda }_1-{\lambda }_2} ; {\varphi }_{21}\left(t\right)=\frac{c_6\left(e^{{\lambda }_{1}t}-e^{{\lambda }_{2}t}\right)}{{\lambda }_1-{\lambda }_2} \\ {\varphi }_{22}\left(t\right)=\frac{\left[{\lambda }_1-\left(\mu +\alpha +\gamma \right)\left(R_c-1\right)\right]e^{{\lambda }_{1}t}-\left[{\lambda }_2-\left(\mu +\alpha +\gamma \right)\left(R_c-1\right)\right]e^{{\lambda }_{2}t}}{{\lambda }_1-{\lambda }_2}, \end{gathered}$$

where ${\lambda }_1$ and ${\lambda }_2$ are the eigenvalues of $C_0$. Note that these eigenvalues are, respectively, the roots $r_1$ and $r_2$, given in (22), of the polynomial $P\left(\lambda \right)$ defined in the proof of Theorem 3. Under Assumption 1, the condition (iv) of Theorem 1, and $R_c<{\overline{R}}_c$, it follows that the eigenvalues ${\lambda }_1$ and ${\lambda }_2$ are real and ${\lambda }_2<{\lambda }_1<0$. From (29) and (30), it follows that:

(31)$I\left(t\right)={\varphi }_{11}\left(t\right)I\left(0\right)+{\varphi }_{12}\left(t\right)t_r\left(0\right)+{{\displaystyle \int }}_0^t{\varphi }_{11}\left(t-\tau \right)\beta \left(\frac{S\left(\tau \right)}{N\left(\tau \right)}-\frac{S_{DFE}}{N_{DFE}}\right)I\left(\tau \right)d\tau \forall t\ge 0.$

The fact that ${\lambda }_2<{\lambda }_1<0$ implies that $e^{{\lambda }_{1}t}\ge e^{{\lambda }_{2}t}\ge 0$ $\forall t\ge 0$. Moreover, $c_5+{\lambda }_2\ge 0$, since $c_6\ge 0$ from its definition. Then, ${\varphi }_{11}\left(t\right)\le \frac{\left(c_5+{\lambda }_1\right)e^{{\lambda }_{1}t}}{{\lambda }_1-{\lambda }_2}$, ${\varphi }_{12}\left(t\right)\le 0$ $\forall t\ge 0$, and one obtains from (31) that:

(32)$I\left(t\right)\le \frac{\left(c_5+{\lambda }_1\right)I\left(0\right)}{{\lambda }_1-{\lambda }_2}{e}^{{\lambda }_{1}t}+\frac{\left(c_5+{\lambda }_1\right)\beta }{{\lambda }_1-{\lambda }_2}{{\displaystyle \int }}_0^t{e}^{{\lambda }_1\left(t-\tau \right)}\left(1-\frac{S_{DFE}}{N_{DFE}}\right)I\left(\tau \right)d\tau \forall t\ge 0,$

where the fact that $\frac{S\left(t\right)}{N\left(t\right)}\le 1$ $\forall t\ge 0$ from the positivity of the model has been used. Moreover, from the definition of $S_{DFE}$ and $N_{DFE}$ given in (10), it follows:

(33)$I\left(t\right)\le \frac{c_5+{\lambda }_1}{{\lambda }_1-{\lambda }_2}{e}^{{\lambda }_{1}t}\left(I\left(0\right)+\frac{\mu c_{1}q\beta }{c_1\left(\mu +\rho \right)+c_2}{{\displaystyle \int }}_0^t{e}^{-{\lambda }_1\tau }I\left(\tau \right)d\tau \right) \forall t\ge 0.$

By applying the Bellman–Gronwall Lemma I [28] in (33), one obtains:

(34)$I\left(t\right)\le \frac{c_5+{\lambda }_1}{{\lambda }_1-{\lambda }_2}{e}^{{\lambda }_{1}t}\left(I\left(0\right)+\frac{\mu c_{1}q\beta }{c_1\left(\mu +\rho \right)+c_2}\underset{0}{\overset{t}{{\displaystyle \int }}}\frac{c_5+{\lambda }_1}{{\lambda }_1-{\lambda }_2}I\left(0\right)e^{{{\displaystyle \int }}_{\tau }^t\frac{\left(c_5+{\lambda }_1\right)\mu c_{1}q\beta }{\left[c_1\left(\mu +\rho \right)+c_2\right]\left({\lambda }_1-{\lambda }_2\right)}ds}d\tau \right)$

that leads by direct calculations to:

(35)$I\left(t\right)\le \frac{c_5+{\lambda }_1}{{\lambda }_1-{\lambda }_2}I\left(0\right)e^{\left({\lambda }_1+\frac{\left(c_5+{\lambda }_1\right)\mu c_{1}q\beta }{\left[c_1\left(\mu +\rho \right)+c_2\right]\left({\lambda }_1-{\lambda }_2\right)}\right)t}=\frac{c_5+{\lambda }_1}{{\lambda }_1-{\lambda }_2}I\left(0\right)e^{F\left(q\right)t}$

with $F\left(q\right)={\lambda }_1+\frac{\left(c_5+{\lambda }_1\right)\mu c_{1}q\beta }{\left[c_1\left(\mu +\rho \right)+c_2\right]\left({\lambda }_1-{\lambda }_2\right)}$. Note that the eigenvalues ${\lambda }_1$ and ${\lambda }_2$ depend on $R_c$, which depends on $q$ from its definition in (13). Furthermore, one obtains that $F\left(0\right)={\lambda }_1\left(0\right)<0$ since ${\lambda }_1\left(q\right)<0$ $\forall q\in \left[0,1\right]$ while:

(36)$\frac{dF}{dq}=\frac{2{\left(\mu \beta c_1\right)}^2{c}_{6}q}{{\left[c_1\left(\mu +\rho \right)+c_2\right]}^2{\left({\left[c_5+\left(\mu +\alpha +\gamma \right)\left(R_c-1\right)\right]}^2-4c_6\right)}^{\frac{3}{2}}}.$

Then, $\frac{dF}{dq}=0$ for $q=0$ and $\frac{dF}{dq}>0$ $\forall q>0$, since $c_6<\frac{{\left[c_5+\left(\mu +\alpha +\gamma \right)\left(R_c-1\right)\right]}^2}{4}$ from the condition (iv) of Theorem 1. Such a fact implies that there exists a critical value $q_c>0$ such that $F\left(q\right)<0$ $\forall q\in \left[0, q_c\right)$ by continuity of the function $F\left(q\right)$ with respect to $q$. As a consequence, one obtains that $\underset{t\to \infty }{\lim }\left\{I\left(t\right)\right\}=0$ if $q\in \left[0, q_c\right){\displaystyle \cap }\left[0,1\right]$. From (29) and (30), it follows that:

(37)$t_r\left(t\right)={\varphi }_{21}\left(t\right)I\left(0\right)+{\varphi }_{22}\left(t\right)t_r\left(0\right)+{{\displaystyle \int }}_0^t{\varphi }_{21}\left(t-\tau \right)\beta \left(\frac{S\left(\tau \right)}{N\left(\tau \right)}-\frac{S_{DFE}}{N_{DFE}}\right)I\left(\tau \right)d\tau \forall t\ge 0.$

The fact that ${\lambda }_2<{\lambda }_1<0$ implies that ${\varphi }_{21}\left(t\right)\le \frac{c_6{e}^{{\lambda }_{1}t}}{{\lambda }_1-{\lambda }_2}$ $\forall t\ge 0$. Moreover, from (22), it follows that ${\lambda }_2-\left(\mu +\alpha +\gamma \right)\left(R_c-1\right)<{\lambda }_1-\left(\mu +\alpha +\gamma \right)\left(R_c-1\right)<0$ and then:

(38)${\varphi }_{22}\left(t\right)\le -\frac{{\lambda }_2-\left(\mu +\alpha +\gamma \right)\left(R_c-1\right)}{{\lambda }_1-{\lambda }_2}{e}^{{\lambda }_{2}t} \forall t\ge 0.$

One obtains from (37) that:

(39)$t_r\left(t\right)\le \frac{c_{6}I\left(0\right)}{{\lambda }_1-{\lambda }_2}{e}^{{\lambda }_{1}t}+\frac{\left(\mu +\alpha +\gamma \right)\left(R_c-1\right)-{\lambda }_2}{{\lambda }_1-{\lambda }_2}{t}_r\left(0\right)e^{{\lambda }_{2}t}+\frac{c_6\beta }{{\lambda }_1-{\lambda }_2}\left(1-\frac{S_{DFE}}{N_{DFE}}\right)\underset{0}{\overset{t}{{\displaystyle \int }}}{e}^{{\lambda }_1\left(t-\tau \right)}I\left(\tau \right)d\tau$

where the fact that $\frac{S\left(t\right)}{N\left(t\right)}\le 1$ $\forall t\ge 0$ from the positivity of the model has been used. Moreover, from the definition of $S_{DFE}$ and $N_{DFE}$ given in (10), it follows:

(40)$t_r\left(t\right)\le \frac{c_6}{{\lambda }_1-{\lambda }_2}{e}^{{\lambda }_{1}t}\left(I\left(0\right)+\frac{\mu c_{1}q\beta }{c_1\left(\mu +\rho \right)+c_2}{{\displaystyle \int }}_0^t{e}^{-{\lambda }_1\tau }I\left(\tau \right)d\tau \right)+\frac{\left(\mu +\alpha +\gamma \right)\left(R_c-1\right)-{\lambda }_2}{{\lambda }_1-{\lambda }_2}{t}_r\left(0\right)e^{{\lambda }_{2}t}.$

By introducing (35) in (40), one obtains:

(41)$t_r\left(t\right)\le \frac{c_{6}I\left(0\right)}{{\lambda }_1-{\lambda }_2}{e}^{{\lambda }_{1}t}\left(1+\frac{\mu c_{1}q\beta \left(c_5+{\lambda }_1\right)}{\left[c_1\left(\mu +\rho \right)+c_2\right]\left({\lambda }_1-{\lambda }_2\right)}\underset{0}{\overset{t}{{\displaystyle \int }}}{e}^{\left(F\left(q\right)-{\lambda }_1\right)\tau }d\tau \right)+\frac{\left(\mu +\alpha +\gamma \right)\left(R_c-1\right)-{\lambda }_2}{{\lambda }_1-{\lambda }_2}{t}_r\left(0\right)e^{{\lambda }_{2}t}$

and finally:

(42)$t_r\left(t\right)\le \frac{c_{6}I\left(0\right)e^{{\lambda }_{1}t}+\left[\left(\mu +\alpha +\gamma \right)\left(R_c-1\right)-{\lambda }_2\right]t_r\left(0\right)e^{{\lambda }_{2}t}}{{\lambda }_1-{\lambda }_2}+\frac{\mu c_1{c}_{6}q\beta \left(c_5+{\lambda }_1\right)I\left(0\right)}{\left[c_1\left(\mu +\rho \right)+c_2\right]{\left({\lambda }_1-{\lambda }_2\right)}^2\left(F\left(q\right)-{\lambda }_1\right)}\left(e^{F\left(q\right)t}-e^{{\lambda }_{1}t}\right).$

Then, one obtains that $\underset{t\to \infty }{\lim }\left\{t_r\left(t\right)\right\}=0$ if $q\in \left[0, q_c\right){\displaystyle \cap }\left[0,1\right]$ since ${\lambda }_1<0$ and $F\left(q\right)<0$ under such a condition. On the other hand, if one applies the variable change $Y_n\left(t\right)=X_n\left(t\right)-X_{n,DFE}$ in the first equation of (24), then:

(43)${\dot{Y}}_n\left(t\right)=A{Y}_n\left(t\right)+B\left(t\right)X_I\left(t\right)$

whose solution is:

(44)$Y_n\left(t\right)={\varphi }^A\left(t\right)Y_n\left(0\right)+{{\displaystyle \int }}_0^t{\varphi }^A\left(t-\tau \right)B\left(\tau \right)X_I\left(\tau \right)d\tau \forall t\ge 0$

where ${\varphi }^A\left(t\right)=e^{At}=L^{-1}\left\{{\left(s{\mathbb{I}}_3-A\right)}^{-1}\right\}$. The eigenvalues of $A$ are ${\lambda }_1^A=-\mu$ and the roots of the Hurwitz polynomial $Q\left(\lambda \right)$, which are defined in the proof of Theorem 3, namely:

(45)${\lambda }_{2,3}^A=\frac{1}{2}\left(-\left(\mu +\rho +c_1\right)\left(R_c-1\right)\pm \sqrt{{\left(c_1-\left(\mu +\rho \right)\right)}^2-4c_2}\right)$

where the sign ‘+’ corresponds to the root ${\lambda }_2^A$ and the sign “$-$” corresponds to ${\lambda }_3^A$. Since the real part of such roots is negative, there exists some norm-dependent real constant $K_A\ge 1$ such that:

(46)$\Vert Y_n\left(t\right)\Vert \le K_A{e}^{{\sigma }_{A}t}(\Vert Y_n\left(0\right)\Vert +{{\displaystyle \int }}_0^t{e}^{{\sigma }_A\tau }\Vert B\left(\tau \right)\Vert \Vert X_I\left(\tau \right)\Vert d\tau ) \forall t\ge 0,$

where ${\sigma }_A=max\left\{-\mu ,{\lambda }_2^A\right\}<0$ is the stability abscissa of $A$. From (35) and (42), one obtains that:

(47)$\Vert X_I\left(t\right)\Vert =\sqrt{I^2\left(t\right)+t_r{}^2\left(t\right)}\le I\left(t\right)+t_r\left(t\right)\le g_1{e}^{{\lambda }_{1}t}+g_2{e}^{{\lambda }_{2}t}+g_3{e}^{F\left(q\right)t} \forall t\ge 0,$

where the positivity of the model has been used, and:

(48)$$\begin{gathered} g_1=\frac{c_{6}I\left(0\right)}{{\lambda }_1-{\lambda }_2}\left(1-\frac{\mu c_{1}q\beta \left(c_5+{\lambda }_1\right)}{\left[c_1\left(\mu +\rho \right)+c_2\right]\left({\lambda }_1-{\lambda }_2\right)\left(F\left(q\right)-{\lambda }_1\right)}\right) ; g_2=\frac{\left[\left(\mu +\alpha +\gamma \right)\left(R_c-1\right)-{\lambda }_2\right]t_r\left(0\right)}{{\lambda }_1-{\lambda }_2} \\ {g}_3=\frac{c_5+{\lambda }_1}{{\lambda }_1-{\lambda }_2}I\left(0\right)\left(1+\frac{\mu c_1{c}_{6}q\beta }{\left[c_1\left(\mu +\rho \right)+c_2\right]\left({\lambda }_1-{\lambda }_2\right)\left(F\left(q\right)-{\lambda }_1\right)}\right). \end{gathered}$$

From (25), it follows that:

(49)$\Vert B\left(t\right)\Vert =max\left\{1,c_3+\gamma +\left(\beta +c_4\right)\frac{S\left(t\right)}{N\left(t\right)}\right\}\le max\left\{1,c_3+\gamma +\beta +c_4\right\}=K_B \forall t\ge 0$

where $\frac{S\left(t\right)}{N\left(t\right)}\le 1$ has been used. By introducing (47) and (49) in (46), one obtains:

(50)$\Vert Y_n\left(t\right)\Vert \le K_A{e}^{{\sigma }_{A}t}\left(\Vert Y_n\left(0\right)\Vert +K_B{{\displaystyle \int }}_0^t{e}^{-{\sigma }_A\tau }\left(g_1{e}^{{\lambda }_1\tau }+g_2{e}^{{\lambda }_2\tau }+g_3{e}^{F\left(q\right)\tau }\right)d\tau \right).$

By direct calculation from (50), it follows that:

(51)$$\begin{gathered} \Vert Y_n\left(t\right)\Vert \le K_A\Vert Y_n\left(0\right)\Vert e^{{\sigma }_{A}t} \\ +K_B\left[\frac{g_1}{{\lambda }_1-{\sigma }_A}\left(e^{{\lambda }_{1}t}-e^{{\sigma }_{A}t}\right)+\frac{g_2}{{\lambda }_2-{\sigma }_A}\left(e^{{\lambda }_{2}t}-e^{{\sigma }_{A}t}\right)+\frac{g_3}{F\left(q\right)-{\sigma }_A}\left(e^{F\left(q\right)t}-e^{{\sigma }_{A}t}\right)\right]. \end{gathered}$$

Finally, from (51), one obtains that $\underset{t\to \infty }{\lim }\left\{\Vert Y_n\left(t\right)\Vert \right\}=0$ if $q\in \left[0, q_c\right){\displaystyle \cap }\left[0,1\right]$, since ${\lambda }_1<0$, ${\lambda }_2<0$, ${\sigma }_A<0$, and $F\left(q\right)<0$ under such a condition. Then, $\underset{t\to \infty }{\lim }\left\{\Vert X_n\left(t\right)-X_{n,DFE}\Vert \right\}=0$ and, also, $\underset{t\to \infty }{\lim }\left\{X_n\left(t\right)\right\}=X_{n,DFE}$. In summary, the DFE point is globally asymptotically stable. □

Remark 5.

Note the following features:

• (i) 

  From Theorems 1–4, the model is positive and only has an equilibrium point, namely, the DFE point, which is locally and globally asymptotically stable provided that the control parameters fulfil the conditions of Theorem 1 and that $R_c<{\overline{R}}_c$ or, equivalently, $R_0<{\overline{R}}_0$.

• (ii) 

  The control parameters $c_1$ and $c_2$ influence the components of both the DFE and EE points, see (10) and (14), (15) respectively.

• (iii) 

  The control parameters $c_3$ , $c_4$ , $c_5$ , and $c_6$ influence the components of the EE point.

• (iv) 

  The control parameters $c_1$ , $c_2$ , $c_5$ , and $c_6$ influence the stability of the DFE point according to (19) and, also, they can imply the non-existence of the EE point under an appropriate choice of their values according to (16).

• (v) 

  The threshold value ${\overline{R}}_c$ given in (16) depends on the control parameters $c_5$ and $c_6$ associated with the treatment effort, while the control reproduction number $R_c$ given in (13) depends on the control parameters $q$, $c_1$, and $c_2$ associated with the vaccination efforts so that the non-existence of the EE point can be guaranteed by a treatment strategy adapted to a designed vaccination campaign.

• (vi) 

  Neither$R_c$nor${\overline{R}}_c$depend on the control parameters$c_3$and $c_4$associated with the effort of the vaccination of susceptible individuals so that such parameters are not relevant for eradicating the disease. Such parameters affect the values of the subpopulations at the EE point if such a point is reached, which is intended to be avoided.

• (vii) 

  The expression for $R_c$ can be equivalently written as $R_c=\frac{\beta \left[\mu \left(1-q\right)+\rho \right]}{\left(\mu +\alpha +\gamma \right)\left[\mu +\rho +\frac{c_2}{c_1}\right]}=\frac{\mu \left(1-q\right)+\rho }{\mu +\rho +\frac{c_2}{c_1}}{R}_0$. Then, one can see that $R_c$ is inversely proportional to $\frac{c_2}{c_1}$ so that an increment in the value of $\frac{c_2}{c_1}$ results in a decrement of $R_c$, implying a small incidence of the infectious disease. In this context, a large value for $\frac{c_2}{c_1}$ is interesting for reducing the incidence of the disease on the host population. A large value for the relation $\frac{c_2}{c_1}$ can be obtained by considering small values for $c_1$. However, the condition (i) of Theorem 1 requires a lower bound for $c_1$, namely, $c_1>\mu +\beta +2\sqrt{c_2+c_4}$, for some prescribed values for $\mu$, $\beta$, $c_2$, and $c_4$, in order to guarantee the non-negativity of the solutions of the controlled model under any non-negative initial condition. Then, the only way of increasing the relation $\frac{c_2}{c_1}$ is by means of an increment of $c_2$, which also implies an increment of $c_1$ to guarantee the condition $c_1>\mu +\beta +2\sqrt{c_2+c_4}$ for some prescribed values for $\mu$, $\beta$, and $c_4$. In summary, the only practical way of increasing the relation $\frac{c_2}{c_1}$ is via increasing simultaneously the value of both parameters $c_1$ and $c_2$. However, a large value for $c_2$ can imply large values for the vaccination control effort, since it affects directly the forced term of Equation (4) for the dynamics of the vaccination law. In fact, the vaccination can be constrained to a number of available vaccines in a practical situation, which implies upper-bounds for the control parameters $c_2$, $c_3$, and $c_4$ of the forced terms of (4).

• (viii) 

  One can see that ${\overline{R}}_c$ is directly proportional to $\frac{c_6}{c_5}$ so that an increment in the value of $\frac{c_6}{c_5}$ implies an increase of ${\overline{R}}_c$. Moreover, note that the EE point of the controlled model does not exist if $R_c<{\overline{R}}_c$ so that an increment of ${\overline{R}}_c$ can be interesting in order to guarantee the non-existence of such an EE point for a prescribed value for $R_c$ adjusted by values for $q$, $c_1$, and $c_2$ adapted to the number of available vaccines.

• (ix) 

  The influence of the parameter $q\in \left[0,1\right]$ on the value of the control reproduction number $R_c$ is negligible if the value of the natural death rate $\mu$ is very small, as it happens in the case of humans. In such a case, the influence of $q\in \left[0,1\right]$ on the eigenvalues ${\lambda }_1$ and ${\lambda }_2$ of the matrix $C_0$ as well as on the function $F\left(q\right)$, both defined in the proof of Theorem 4, is also negligible. Such a fact implies that the DFE point is globally asymptotically stable $\forall q\in \left[0,1\right]$, since $q_c>1$, provided that the control parameters $c_i$, for $i\in \left\{1, 2, \dots , 6\right\}$, are chosen such that Assumption 1 and $R_c<{\overline{R}}_c$ are satisfied.

3. Simulation Results

Some simulation examples based on a high infectivity disease illustrate the efficacy of the proposed control strategy. The examples have been developed by using the version R2020b of MATLAB. They have been simulated by using the solver “ode3 (Bogacki-Shampine)” with a fixed step of $\frac{1}{24}$ days, i.e., one hour. First, a SIRS model without vaccination and treatment is considered. Later, the same model with the proposed control actions is analysed to show their impact in the propagation of the disease within the host population. As a consequence, a drastic mortality reduction is exhibited due to the application of the proposed feedback control actions. In addition, a set of examples are used to analyse the influence of certain control parameters on the evolution of the disease spreading. Finally, the influence of the immunity loss rate is also studied.

3.1. Example 1: SIRS Model without Vaccination and Treatment

The model (1)–(3) with the values $b=0.0384 d^{-1}$, $\mu =3.653\times {10}^{-5} d^{-1}$, $\beta =0.65 d^{-1}$, $\gamma =\frac{1}{22}=0.0455 d^{-1}$, $\rho =\frac{1}{120}=0.0083 d^{-1}$ and $\alpha =5\times {10}^{-5} d^{-1}$, where $d$ means the unit time “day”, is analysed. In addition, the parameter $q=0$ is considered, since vaccination on the newborn individuals is not applied. The main objective is to show the time evolution of the subpopulations and that of the whole population under the influence of the infectious disease. The basic reproduction number of this model is $R_0=14.2728$; that is, each primary infectious individual transmits the disease to more than 14 healthy individuals. Then, this example describes the propagation of a high infectivity disease. The DFE point of the model is unstable, while its EE point is globally asymptotically stable, as Remark 3 points out. Figure 1 shows the time evolution of the susceptible, infectious, and recovered subpopulations during the first 100 days, while Figure 2 displays that of the whole population in a very long time period. The considered initial condition is given by $S\left(0\right)=985$, $I\left(0\right)=10$ and $R\left(0\right)=5$.

In the longer simulation, one can see that the model tends to an EE point where the number of individuals is $N_{EE}=877$ with $S_{EE}=62$, $I_{EE}=127$ and $R_{EE}=688$ accordingly to (14), (15) with $c_i=0$, for $i\in \left\{1, 2, \dots , 6\right\}$. Although the mortality rate due to the disease $\alpha =5\times {10}^{-5} d^{-1}$ is not very high, the influence of the disease is very important, since the whole population notably decreases from 1000 individuals to 877. Moreover, the percentage of infectious individuals with respect to the whole population in the EE point is 14.48%. In this situation, the application of control measures is crucial in order to eliminate the infection or, at least, to diminish its effect within the host population while reducing drastically the mortality.

3.2. Example 2: SIRS Model with Vaccination and Treatment

Models (1)–(3) with the same values for the parameters $b$, $\mu$, $\beta$, $\gamma$, $\rho$ and $\alpha$ as those given in Example 1 are considered. Moreover, a vaccination for newborn individuals together with the control actions (4) and (5) are applied in order to diminish the effects of the infectious disease. Concretely, the control actions are (i) a vaccination of a proportion $q=0.1$ of the newborn individuals; i.e., 10% of newborns are vaccinated, (ii) a vaccination of the susceptible population given by (4) with the values $c_1=25$, $c_2=1$, $c_3=0.0001$ and $c_4=0.07$ for the associated parameters and (iii) a treatment of the infectious population given by (5) with the values $c_5=2$ and $c_6=0.27$ for the corresponding parameters. Such values satisfy the conditions of Theorem 1 so that the model subpopulations and the vaccination and treatment efforts are non-negative for all time. The initial condition is that of Example 1 together with $v\left(0\right)=t_r\left(0\right)=0$ for the initial control efforts. Note that dynamics of the treatment effort, given by the control law (5), is of first order with a gain that reaches the value $\frac{c_6}{c_5}=0.135$ in the stationary regime. Such a gain points out a transition rate from the infectious subpopulation to the recovered one via applied treatment to be added to the transition rate $\gamma =0.0455$ between such subpopulations from the natural response of the immunity system of the individuals against the disease. This fact points out that an infectious individual under treatment needs an average time of 5.54 days to overcome the disease. In this way, the average time of recovering is reduced in 16.46 days, since such an interval time is ${\gamma }^{-1}=22$ days if no treatment is applied. Figure 3 shows the time evolution of the susceptible, infectious, and recovered subpopulations during the first 100 days, while Figure 4 and Figure 5 display, respectively, that of the whole population during the first 100 days and along 100,000 days.

In a longer time performed simulation, one can see that the model tends to a DFE point where the number of individuals is $N_{DFE}=1050$ with $S_{DFE}=182$, $I_{DFE}=0$ and $R_{DFE}=868$. This result is coherent with Theorems 3 and 4, since $R_c=2.4687$, ${\overline{R}}_c=3.9644$ and $F\left(q\right)=-0.0732<0$ for the used value $q=0.1$ so that the DFE point is locally and globally asymptotically stable. Moreover, the choice of the control parameters satisfies the conditions of Theorem 2, which implies the non-existence of the EE point. Then, the DFE point is the unique equilibrium point for the controlled model. In summary, this example exhibits the high mortality reduction due to the application of the proposed feedback control actions in the propagation of an infectious disease of a high basic reproduction number $R_0$. Figure 6 displays the time evolution of the control efforts. One can see that the treatment effort reaches a peak value when the number of infectious individuals reaches its maximum value, and then, such an effort converges to zero as a consequence of the convergence to zero of the infectious subpopulation. On the other hand, the vaccination effort converges to a constant value $v_{DFE}=7$.

The vaccination and treatment efforts can be interpreted as the number of vaccines and antivirals, or some appropriated medicaments, applied to the susceptible and infectious individuals, respectively, during the transition from the initial state to the DFE point. In this context, one can consider that the DFE point is reached when the number of infectious individuals is less than 1. Moreover, although the vaccination effort has a constant value $v_{DFE}=7$ given by (10) at the DFE point, one can stop the vaccination control effort once the number of infectious individuals is smaller than 1. In this way, the vaccination can be set $v\left(t\right)=0$ once the DFE point is reached since the transition rate from the susceptible subpopulation to the infectious one is zero in the absence of infectious individuals due to such a rate depending on $\beta \frac{I\left(t\right)}{N\left(t\right)}$. In other words, there is no propagation of the disease in absence of infectious individuals, and then, the vaccination can be stopped. However, the vaccination with the value $v_{DFE}$ could be maintained as a preventive measure against possible new outbreaks of the disease due to immigration or other causes. This preventive measure gives place to a transition from the susceptible subpopulation to the recovered one with a rate that compensates for the transition rate $\rho$ from the recovered subpopulation to the susceptible one caused by the immunity loss. In case that the vaccination is stopped once the DFE point is reached, which is at the 64th day in this example, the number of vaccines and medicaments applied to control the propagation in the transient from the initial day until achieving the DFE point are, respectively, 578 and 483. Such numbers can be interpreted as the control cost to eradicate the persistence of the disease. Figure 7 displays the number of applied vaccines and medicaments each day during the transition to the DFE point. The numbers of vaccines and medicaments applied each day, respectively $n_v\left(t\right)$ and $n_{tr}\left(t\right)$, are obtained by evaluating $n_v\left(t\right)={{\displaystyle \int }}_{t-1}^{t}v\left(\tau \right)d\tau$ and $n_{tr}\left(t\right)={{\displaystyle \int }}_{t-1}^t{t}_r\left(\tau \right)d\tau$ for all integer $t\in \left[1, 2\dots ,t_f\right]$ with $t_f$ denoting the day at which the DFE point is reached.

3.3. Example 3: SIRS Model with a Defective Vaccination and Treatment

The same values as those used in Example 2 for the parameters of the model are considered except that for $c_2$, which is replaced by $c_2=0.2$. In this way, $R_c=7.2945$, while ${\overline{R}}_c=3.9644$. Then, the conditions of Theorem 2 are not satisfied, which implies that the non-existence of the EE point cannot be guaranteed. In fact, the time evolution of the subpopulations under these conditions is displayed in Figure 8, and one can see that the EE point given by (14) and (15) is reached. Concretely, the whole number of individuals at the EE point is $N_{EE}=1022$ with $S_{EE}=284$, $I_{EE}=21$, and $R_{EE}=717$, and the vaccination and treatment control signals converge, respectively, to the values $v_{EE}=2$ and $t_r{}_{EE}=3$.

The time evolution of the control efforts is displayed in Figure 9.

3.4. Example 4: Study of the Influence of the Control Parameter $c_2$ on the Behaviour of the Controlled Model

Models (1)–(5) with the same values for the parameters as those given in Example 2, except that for $c_2$, are considered. In this way, ${\overline{R}}_c=3.9644$ in all the analysed cases. Concretely, four cases are studied. Each one considers a different value for the parameter $c_2$ in order to change slightly the ratio $\frac{c_2}{c_1}$, which implies a different value for the control reproduction number $R_c$. The condition $R_c<{\overline{R}}_c$ as well as the conditions of the Theorems 1, 2, 3, and 4 are satisfied in the four cases so that the DFE point is the unique equilibrium point, and it is locally and globally asymptotically stable, while the solutions of the controlled model are non-negative for all time under any non-negative initial condition. The interest of this study is to analyse the influence of the parameter $c_2$ on the transition of the controlled model solutions from the initial state to the DFE point by evaluating the number of infectious individuals on each day as well as the cost of both vaccination and treatment control efforts. For such a purpose, the values for the control parameter $c_2$ and then that of $R_c$ in the four studied cases are: (i) $c_2=1$, $R_c=2.4687$, (ii) $c_2=0.9$, $R_c=2.6912$, (iii) $c_2=0.8$, $R_c=2.9579$, and (iv) $c_2=0.7$, $R_c=3.2832$. Note that case (i) is that studied in Example 2, and it is used as the reference one for the current analysis. The initial condition is that used in Example 2 for all the cases. Figure 10 displays the time evolution of the infectious subpopulation in the four cases. Figure 11 and Figure 12 show, respectively, the time evolution of the vaccination and treatment efforts in the four cases.

One can see in these figures that the peaks in the infectious population, the vaccination and treatment efforts depend on the value of the control reproduction number $R_c$, which is established by the control parameters $c_1$ and $c_2$. Concretely, the peak of the infectious subpopulation as well as that of the treatment effort increases as $R_c$ increases, while the peak of the vaccination effort decreases as $R_c$ increases, while maintaining $R_c<{\overline{R}}_c$. Table 1 summarises the most relevant specifications, included the aforementioned ones, about the four cases studied in this subsection.

The results show that the duration of the transient from the initial state to the DFE point increases as $R_c$ increases. However, the number of vaccines applied during the whole transient period decreases when the value of $R_c$ increases. Such a result is due to the fact that the peak of the vaccination effort decreases and, also, the number of vaccines used per day, when $R_c$ increases. In summary, the number of required vaccines decreases while that of the required medicaments increases as $R_c$ increases. Then, a tradeoff between the cost of vaccination and the number of infectious individuals with its associated treatment cost has to be taken into account for adjusting the values for the control parameters $c_1$ and $c_2$. Such a tradeoff is going to depend on constraints about the number of available vaccines for the newborns and the susceptible individuals and/or medicaments for treatment of the infectious subpopulation. Finally, Figure 13 and Figure 14 display the number of applied vaccines and medicaments each day during the transition to the DFE point for the four analysed cases.

One can see by analysing the results in Table 1 and Figure 10, Figure 11, Figure 12, Figure 13 and Figure 14 that the increase of the value of $R_c$, by decreasing the values of $c_2$, and maintaining the rest of the control parameters in the prescribed values, gives place to a decrease in the vaccination effort but both the infectious subpopulation and the applied treatment effort increase at the same time. The best of the studied cases from the viewpoint of the vaccination cost is case (iv), but it is not good enough from the viewpoint of both the treatment cost and the number of individuals who experience the infectious disease during the transition from the initial state until reaching the DFE point. In this context, it also seems interesting to examine the influence of the control parameter $c_6$, which take part in the dynamics of the treatment effort, on the specifications of the transient behaviour.

3.5. Example 5: Study of the Influence of the Control Parameter $c_6$ on the Behaviour of the Controlled Model

In the current example, the same values for the parameters $b$, $\mu$, $\beta$, $\gamma$, $\rho$, $\alpha$ and $q$ and for the initial condition as those given in Example 4 are considered. Furthermore, the values of $c_1=25$, $c_2=0.7$, and then $R_c=3.2832$, $c_3=0.0001$ and $c_4=0.07$ corresponding with the case of Example 4 with the smallest vaccination cost, concretely case (iv), which requires the application of 467 vaccines as one can see in Table 1, are considered. Such a case is the worst scenario in Example 4 from the viewpoint of the time evolution of the infectious subpopulation and the cost in applied medicaments during the transient from the initial state to the DFE point; namely, 570 medicaments are required, as Table 1 points out. The objective is to analyse the influence of the control parameter $c_6$, which acts on the value of ${\overline{R}}_c$, on the specifications of such a transient time interval. Concretely, the value $c_5=2$ is going to be considered in four proposed cases while the value for the parameter $c_6$ will be modified in order to change slightly the ratio $\frac{c_6}{c_5}$. In this way, each case has a different transition rate from the infectious to recovered subpopulation due to the treatment action. Moreover, such a transition rate is added to the transition rate $\gamma =0.0455 d^{-1}$ from the infectious to the recovered subpopulation in the absence of treatment, i.e., that derived from the natural immunity system of the infectious individuals to fight against the disease. Such a fact implies a reduction in average for the recovering time of the infectious individuals. The analysed cases are (i) $c_6=0.24$, ${\overline{R}}_c=3.635$, (ii) $c_6=0.27$, ${\overline{R}}_c=3.9644$, (iii) $c_6=0.3$, ${\overline{R}}_c=4.2937$, and (iv) $c_6=0.33$, ${\overline{R}}_c=4.6231$. Note that $R_c<{\overline{R}}_c$ as well as the conditions of the Theorems 1, 2, 3 and 4 are satisfied in all the cases so that the DFE point is the unique equilibrium point while being locally and globally asymptotically stable. Note also that case (ii) of this current example is the same as case (iv) of Example 4, which is used as the reference one for the current analysis. In case (i), the recovering time is reduced on average from ${\gamma }^{-1}=$ 22 days, in absence of treatment, to ${\left(\gamma +\frac{c_6}{c_5}\right)}^{-1}=$ 6.04 days if a treatment is applied following the rule (5). For cases (ii), (iii), and (iv), the recovering time passes from 22 to 5.54, 5.12, and 4.75 days, respectively. The values of $S_{DFE}=242$, $R_{DFE}=808$, and $v_{DFE}=7$ are reached at the DFE point in the four analysed cases, as it can be deduced from (10), since all cases have the same values for $c_1$ and $c_2$. Figure 15 displays the time evolution of the infectious subpopulation, while Table 2 summarises the specifications for these four studied cases.

One can see that the peak of the infectious subpopulation is decreasing as the value of $c_6$ and then also ${\overline{R}}_c$, increases. In addition, the number of required vaccines and medicaments during the transient decreases as the value for $c_6$, and then also ${\overline{R}}_c$, increases. Such a result is mainly due to the fact that the duration of the transient decreases as the value of $c_6$ increases, since there is a small duration on average for the transition of individuals from infectious to the recovered subpopulation. Note also that the peak in the vaccination effort is equal for the four studied cases, since such a peak depends on the values of the parameters $c_1$ and $c_2$, which are the same for the four cases. Case (iv) is the best of the analysed ones in the current study. It requires a strong treatment so that an infectious individual overcomes the infection after an average time of 4.75 days. In such a case, the transient from the initial state until arriving at the DFE point has the following characteristics:

• The evolution of the infectious population reaches a peak of 215 individuals, i.e., approximately 21% of the initial whole population.

• The vaccination cost during the transient supposes 428 vaccines, i.e., the percentage of the susceptible subpopulation to be vaccinated is around of the 43% of the initial population, assuming one vaccine per individual.

• The treatment cost during the transient is of 554 medicaments, i.e., the percentage of infectious subpopulation to be treated is around 55% of the initial population, assuming one medicament per individual.

If the applied treatment has a lower performance so that the average recovering time is 5.54 days, corresponding to case (ii), instead of the 4.75 days of case (iv), then the peak of the infectious population is of 261 individuals, i.e., around the 26% of the initial population. In such a situation, the costs of vaccination and treatment are, respectively, 467 vaccines and 570 medicaments. In the considered worst case, i.e., case (i), which requires an average time of 6.04 days for the transition from the infectious to the recovered subpopulation, the number of necessary vaccines and medicaments is 525 and 573, respectively. As a conclusion, the transient duration and the cost in vaccines and medicaments are related with the performance of the applied treatment. If the applied medicaments have a poor potency such that the reduction of the recovery time is of a few days, then the duration of the transient from the initial state until reaching the DFE point can be very long. In such a case, the transient can be very expensive in relation with the number of required vaccines and medicaments or, even, the disease evolution can converge to the EE point if the resources in the number of medicaments and/or vaccines is not enough.

3.6. Example 6: Study of the Behaviour of the Controlled Model under Vaccination

The objective of this subsection is to study the dynamics of the controlled model when there are not medicaments for treating the infectious subpopulation, and then the only measure to control the propagation of the disease is the planning of a vaccination campaign to the newborns and the susceptible subpopulation. For this purpose, models (1)–(5) with $c_6=0$ can be used with the initial condition $t_r\left(0\right)=0$. In this way, $t_r\left(t\right)=0$ for all time. Moreover, ${\overline{R}}_c=1$ from (16) such that $R_c<1$ is required to guarantee the non-existence of the EE point, and then, the DFE point is the unique equilibrium point while being locally and globally asymptotically stable. Since $R_c=\frac{\beta \left[\mu \left(1-q\right)+\rho \right]}{\left(\mu +\alpha +\gamma \right)\left[\mu +\rho +\frac{c_2}{c_1}\right]}=\frac{\mu \left(1-q\right)+\rho }{\mu +\rho +\frac{c_2}{c_1}}{R}_0$, it is necessary that $\frac{c_2}{c_1}>{\left(\frac{c_2}{c_1}\right)}_{min}=\left[\mu \left(1-q\right)+\rho \right]R_0-\left(\mu +\rho \right)$ to guarantee that $R_c<1$, and then, the eradication of the infectious disease can be achieved. Table 3 compares the results obtained for seven different cases of values for the pair $c_1$ and $c_2$, satisfying the conditions of the Theorems 1, 2, 3, and 4 so that the solutions of the controlled model are non-negative and the DFE point is the unique equilibrium point while being locally and globally asymptotically stable. The same values for the parameters $b$, $\mu$, $\beta$, $\gamma$, $\rho$, $\alpha$, $q$, $c_3$ and $c_4$ as those given in Example 2 are used, and the same initial condition is used as well.

The results in Table 3 point out that the duration of the transient and also the peak in the infectious subpopulation decreases as the value for the parameter $R_c$ decreases or, equivalently, as the ratio $\frac{c_2}{c_1}$ increases. Such a decrease in the transient duration implies a smaller cost in the number of vaccines as $R_c$ decreases. However, even in the cheapest case shown in Table 3, i.e., case 7, the vaccination cost supposes the use of 1569 vaccines. Such a cost may be decreased even more by considering small values for $R_c$. In any case, the results show that the fight against the disease when there is not treatment and then, only the vaccination is available, is quite expensive. Moreover, the propagation of the disease can reach the EE point if there are not enough resources regarding vaccines. Note that the average rate for the transition from the susceptible subpopulation to the recovered one is approximately $\frac{c_2}{c_1}$ or, equivalently, the average time in the transition of vaccinated individuals from the susceptible subpopulation to the recovered one is $\frac{c_1}{c_2}$ days. In the cases studied in this subsection, such an average time is constrained between the values of 2.5 days, corresponding to case 1, and 1 day for case 7. Such a fact implies that the effect of the vaccination in the population is quite fast, less than 3 days, in any case. In this context, a more realistic scenario can be analysed by considering, for instance, $c_1=25$ and $c_2=5$. In such a case, $R_c=0.5731$, and the average time in the transition from the susceptible to the recovered subpopulation is 5 days. The specifications in the transient for this case are as follows:

• Infectious peak: 92;

• Transient duration: 235 days;

• Number of vaccines: 2559;

• $S_{DFE}=42$, $R_{DFE}=1008$ and $v_{DFE}=8$.

One can see that the transient is very large, as it is the number of vaccines to be applied. In this context, the slower the effect of the vaccination in the population is, the larger the value of the parameter $R_c$ is as well as the duration of the transient and the number of vaccines to guarantee the eradication of the disease. Another alternative to deal with a more realistic scenario in the absence of treatment action could be the inclusion of a delay, corresponding to the reaction time of the vaccines, in the epidemic models (1)–(3).

3.7. Example 7: Study of the Behaviour of the Controlled Model under Treatment

The objective of this subsection is to study the dynamics of the controlled model when there are not vaccines for the newborns and the susceptible subpopulation, and then, the only measure to control the propagation of the disease is the application of some treatment to the infectious subpopulation. For this purpose, models (1)–(5) with $q=c_2=c_3=c_4=0$ can be used with the initial condition $v\left(0\right)=0$. In this way, $v\left(t\right)=0$ for all time. In this situation, $R_c=R_0=\frac{\beta }{\mu +\alpha +\gamma }$ from (13), and such a value cannot be modified by control. However, a treatment procedure on the infectious population based on law (5) with appropriate values for the parameters $c_5$ and $c_6$ so that $R_c<{\overline{R}}_c$ can be designed. Since ${\overline{R}}_c=1+\frac{c_6}{c_5\left(\mu +\alpha +\gamma \right)}$, it is necessary that $\frac{c_6}{c_5}>\beta -\left(\mu +\alpha +\gamma \right)$ to guarantee that $R_c<{\overline{R}}_c$ and then the non-existence of the EE point and the local and global asymptotic stability of the DFE point, which is the unique one in that case. In this way, the eradication of the infectious disease can be achieved. The current subsection compares the results obtained for some different cases satisfying the conditions of Theorems 1, 2, 3, and 4 so that the solutions of the controlled models are non-negative and the DFE point is the unique equilibrium point while being locally and globally asymptotically stable. Note that $\frac{c_6}{c_5}$ is the rate of transition from the infectious subpopulation to the recovered one by means of treatment of infectious individuals and that such a transition rate is added to the average transition rate $\gamma$ between such subpopulations in the absence of treatment. By considering the same values for the parameters $\mu$, $\alpha$, $\gamma$ and $\beta$ as those used in Example 1, then $R_c=R_0=14.2728$, and the transition rate $\frac{c_6}{c_5}$ from treatment actions has to be strictly larger than ${\left(\frac{c_6}{c_5}\right)}_{min}=0.6045$ so that $R_c<{\overline{R}}_c$, and then, the DFE point is the unique equilibrium point while being locally and globally asymptotically stable. This supposes an average time strictly smaller than ${\left(\gamma +{\left(\frac{c_6}{c_5}\right)}_{min}\right)}^{-1}=$ 1.5387 days for recovering from the infectious disease. Such a situation is not realistic, since the existence of medicaments with such an extraordinary performance is improbable. By assuming that the non-existence of the EE point is not secured only with the applied treatment to the infectious subpopulation, at least, such a treatment campaign can reduce the effects of the disease in the whole population. Such a result is illustrated by considering several cases, concretely, those studied in Example 5 but without applying the vaccination of either the susceptible subpopulation, i.e., $c_2=c_3=c_4=0$ with $v\left(0\right)=0$, or the newborns, i.e., $q=0$. The results are presented in Table 4 below.

One can see the decrease of the value of ${\overline{R}}_c$ due to the decrease of the value of $\frac{c_6}{c_5}$, which implies an increase of the average time that an infectious individual stays in the infectious subpopulation before passing to the recovered one, leads to a decrease in the number $N_{EE}$ of individuals in the whole population when the EE point is reached while increasing the percentage of infectious individuals at such an equilibrium point. In any case, the application of treatment reduces considerably the effects of the infectious disease if one compares the results of Table 4 with those of Example 1 where no treatment is applied, and then, the whole population at the EE point is of 877 individuals with 127 infectious, i.e., 14.48% of individuals are infectious. Thus, the application of treatment reduces notably the high mortality and maintains the percentage of infectious subpopulation in reasonably small numbers at the achieved EE point.

3.8. Example 8: Study of the Behaviour of the Controlled Model under Different Rates for the Loss of Immunity

The objective of this subsection is to study the dynamics of the controlled model under both control actions, vaccination and treatment, with different values of the rate $\rho$ for the loss of immunity in the recovered subpopulation or, equivalently, different average periods of immunity for the recovered individuals. For this purpose, case (iv) of Example 4 can be taken as the reference one for this analysis, since it is one of the most convenient from the viewpoint of vaccination and treatment costs. The values of the model parameters are $b=0.0384 d^{-1}$, $\mu =3.653\times {10}^{-5} d^{-1}$, $\beta =0.65 d^{-1}$, $\gamma =\frac{1}{22}=0.0455 d^{-1}$, $\rho =\frac{1}{120}=0.0083 d^{-1}$ and $\alpha =5\times {10}^{-5} d^{-1}$, while those of the control actions are $q=0.1$, $c_1=25$, $c_2=0.7$, $c_3=0.0001$, $c_4=0.07$, $c_5=2$ and $c_6=0.27$ in such a reference case. Other values for the parameter $\rho$ are considered to analyse the influence of such a parameter on the transient specifications. Concretely, the following four cases are treated: (i) $\rho =\frac{1}{100}$, $R_c=3.7647$ (ii) $\rho =\frac{1}{120}$, $R_c=3.2832$, (iii) $\rho =\frac{1}{150}$, $R_c=2.7554$, and (iv) $\rho =0$, $R_c=0.0167$. Note that case (ii) is used as a reference for this study and also note that the value of $\rho$ has an influence on the value for $R_c$. In case (i), the average duration of the immunity is 100 days, and that of case (ii) is 120 days, while that of case (iii) is 150 days. Finally, there is no loss of immunity in case (iv). All the cases have the same values for ${\overline{R}}_c$, namely ${\overline{R}}_c=3.9644$, so that $R_c<{\overline{R}}_c$ and then, the DFE point is locally and globally asymptotically stable for all the analysed cases. The results are presented in Table 5 below.

One can see that the number of necessary vaccines and medicaments in the transient from the initial state to the DFE point decreases when the immunity period increases. Such a fact is mainly due to the transient period decreasing as the immunity period increases. Moreover, the number of recovery individuals at the DFE point increases when the immunity period increases.

4. Conclusions

An SIRS epidemic model with the vaccination of newborns and susceptible individuals and treatment of infectious ones has been investigated. The vaccination and the treatment are governed by a control subsystem containing several free-design parameters. This control subsystem provides two more free-design parameters comparing with those available in the more usual SIRS models with vaccination and treatment. In the proposed model, the intensity of both control actions is not directly proportional to the susceptible and/or infectious subpopulation, as it usually happens in the SIRS models. Here, both actions are provided by the control subsystem, and the parameters defining the dynamics of the controller are also available to shape the vaccination and treatment actions. This control strategy allows the achievement of several relevant results. First, an appropriate adjustment of the control parameters guarantees the positivity of the controller model, as it has been mathematically proved in Theorem 1. Note that such a property is required for coherence in epidemic models where all the subpopulations and the control actions, such as vaccination and treatment, have to be non-negative. Then, the equilibrium points of the proposed controlled epidemic model are calculated as a function of the system parameters. There are two equilibrium points: (i) the DFE point where the infectious subpopulation is zero and then the whole population is composed by susceptible and recovery subpopulations; and (ii) the EE point where the three subpopulations of the model are presented. The DFE one always exists while the existence of the EE point depends on the values of the control parameters. Moreover, the control reproduction number $R_c$ of the controlled epidemic model and a threshold value ${\overline{R}}_c$ are mathematically obtained as functions of the control parameters. This result allows us to analyse the influence of the control parameters on both values $R_c$ and ${\overline{R}}_c$. In summary, the values $R_c$ and ${\overline{R}}_c$ and the existence of the EE point depends on the control parameters. As a consequence, the existence of the EE point can be related with $R_c$ and ${\overline{R}}_c$. In this sense, values of the control parameters doing $R_c<{\overline{R}}_c$ guarantee the non-existence of the EE point, as it has been mathematically proved in Theorem 2. In such a situation, the proposed controlled SIRS model only possesses a unique equilibrium point, namely, the DFE point. On the other hand, the local and global asymptotic stability of the DFE point are mathematically proved in Theorems 3 and 4, respectively. Then, an appropriate adjustment of the control parameters so that $R_c<{\overline{R}}_c$ allows guaranteeing the positivity of the controlled epidemic model while ensuring the non-existence of the EE point. In such a situation, the DFE is the unique equilibrium point of the model, and then, the eradication of the infectious disease can be guaranteed.

Finally, the influence of the control parameters on the time evolution of the disease propagation has been studied by several simulation examples. The first example shows the time evolution of the propagation of the disease without applying control actions. One can see that the model converges to the EE point, and then, the disease is not eradicated. The second example shows the time evolution of the propagation of the disease if the proposed control actions are applied in an effective way, i.e., if the control parameters are chosen such that $R_c<{\overline{R}}_c$. In this case, one can see that the model converges to the DFE point, and the disease is eradicated. The third example shows the time evolution of the propagation of the disease if the proposed control actions are applied in a defective way, i.e., if the control parameters are chosen such that $R_c>{\overline{R}}_c$. One can see that the model converges to the EE point, and the disease is not eradicated. The fourth example analyses the influence of the controller parameter $c_2$ in the evolution of the propagation of the disease. Such a parameter acts in the vaccination control law, and it has an influence on the value of $R_c$. Concretely, the value of $R_c$ is indirectly proportional to $c_2$. This example shows that the model converges to the DFE point whenever $R_c<{\overline{R}}_c$. The number of vaccines and medicaments associated to the control efforts obtained from the control subsystem to achieve the DFE point depends on $c_2$. One can see that by decreasing the value $c_2$, or equivalently increasing the value of $R_c$ but guaranteeing that $R_c<{\overline{R}}_c$, the number of vaccines applied to the susceptible individuals decreases while that of medicaments applied to the infectious individuals increases. The fifth example analyses the influence of the controller parameter $c_6$ in the evolution of the propagation of the disease. Such a parameter acts in the treatment control law, and it has an influence on the value of ${\overline{R}}_c$. Concretely, the value of ${\overline{R}}_c$ is directly proportional to $c_6$. This example shows that the model converges to the DFE point whenever $R_c<{\overline{R}}_c$. The number of vaccines and medicaments associated to the control efforts obtained from the control subsystem to achieve the DFE point depends on $c_6$. One can see that by increasing the value $c_6$, or equivalently increasing the value of ${\overline{R}}_c$ but guaranteeing that $R_c<{\overline{R}}_c$, the number of both vaccines and medicaments applied to the respective subpopulations decreases. Example 6 analyses the behaviour of the model when only the vaccination action is available. Example 7 analyses the behaviour of the model when only the treatment action is available. Example 6 shows that the DFE point is not achieved with an appropriate number of vaccines. Example 7 requires a treatment with a great potency so that the recovery of the infectious individuals after being treated is faster than in a real situation. Finally, Example 8 studies the evolution of the propagation for different rates of loss of immunity. The number of vaccines and medicaments in the convergence of the model to the DFE point decreases as the rate of losing immunity decreases.

Author Contributions

Formal analysis, S.A.-Q. and M.D.l.S.; Funding acquisition, M.D.l.S.; Investigation, S.A.-Q.; Methodology, S.A.-Q.; Resources, R.N.; Software, S.A.-Q. and R.N.; Supervision, M.D.l.S.; Validation, M.D.l.S.; Writing—original draft, S.A.-Q.; Writing—review & editing, M.D.l.S. All authors have read and agreed to the published version of the manuscript.

Funding

This research received support from the Spanish Government through grant RTI2018-094336-B-I00 (MCIU/AEI/FEDER, UE).

Institutional Review Board Statement

Not applicable.

Informed Consent Statement

Not applicable.

Data Availability Statement

All data used are in the references and properly cited within the manuscript.

Acknowledgments

The authors are grateful to the Spanish Government for its support through grant RTI2018-094336-B-I00 (MCIU/AEI/FEDER, UE) and to the Basque Government for its support through grant IT1207-19. They also thank the Spanish Institute of Health Carlos III for its support through Grant COV 20/01213.

Conflicts of Interest

The authors declare no conflict of interest.

Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Figures and Tables

Figure 1. Time evolution of the subpopulations in the SIRS model without control actions.

Figure 2. Time evolution of the whole population in the SIRS model without control actions.

Figure 3. Time evolution of the subpopulations in the SIRS model with vaccination and treatment.

Figure 4. Time evolution of the whole population in the SIRS model with vaccination and treatment within the 100 first days.

Figure 5. Time evolution of the whole population in the SIRS model with vaccination and treatment along 100,000 days.

Figure 6. Time evolution of the control efforts: vaccination and treatment.

Figure 7. Number of vaccines and medicaments applied each day.

Figure 8. Time evolution of the subpopulations in the SIRS model under an insufficient vaccination and treatment.

Figure 9. Time evolution of the control efforts under an insufficient vaccination and treatment.

Figure 10. Time evolution of the infectious subpopulation for the four considered values of the control parameter [Forumla omitted. See PDF.].

Figure 11. Time evolution of the vaccination effort for the four considered values of the control parameter [Forumla omitted. See PDF.].

Figure 12. Time evolution of the treatment effort for the four considered values of the control parameter [Forumla omitted. See PDF.].

Figure 13. Number of applied vaccines each day for the four considered values of the control parameter [Forumla omitted. See PDF.].

Figure 14. Number of applied medicaments each day for the four considered values of the control parameter [Forumla omitted. See PDF.].

Figure 15. Time evolution of the infectious subpopulation for the four considered values of the control parameter [Forumla omitted. See PDF.].

Table 1

Specifications of the model behaviour for the four considered values of the parameter $c_2$.

	${\mathit{c}}_2=1$ ${\mathit{R}}_{\mathit{c}}=2.4687$	${\mathit{c}}_2=0.9$ ${\mathit{R}}_{\mathit{c}}=2.6912$	${\mathit{c}}_2=0.8$ ${\mathit{R}}_{\mathit{c}}=2.9579$	${\mathit{c}}_2=0.7$ ${\mathit{R}}_{\mathit{c}}=3.2832$
Infectious peak	215	229	245	261
Transient duration (days)	64	66	69	72
Vaccination peak	39	35	31	27
Treatment peak	29	31	33	35
SDFE	182	198	218	242
RDFE	868	852	832	808
$v_{DFE}$	7	7	7	7
Vaccines	578	537	502	467
Medicaments	483	515	542	570

Table 2

Specifications of the model behaviour for the four considered values of the parameter $c_6$.

	${\mathit{c}}_6=0.24$ ${\overline{\mathit{R}}}_{\mathit{c}}=3.635$	${\mathit{c}}_6=0.27$ ${\overline{\mathit{R}}}_{\mathit{c}}=3.9644$	${\mathit{c}}_6=0.3$ ${\overline{\mathit{R}}}_{\mathit{c}}=4.2937$	${\mathit{c}}_6=0.33$ ${\overline{\mathit{R}}}_{\mathit{c}}=4.6231$
Infectious peak	286	261	237	215
Transient duration (days)	87	72	63	58
Vaccination peak	27	27	27	27
Treatment peak	34	35	35	35
Vaccines	525	467	439	428
Medicaments	573	570	565	554

Table 3

Specifications of the model behaviour for seven considered pairs of values of the control parameters $c_1$ and $c_2$ when there is not treatment to control the disease propagation.

	Infectious Peak	Transient Duration (Days)	Vaccines	DFE Point
Case 1$c_1=25 ; c_2=10$$R_c=0.2924$	33	120	1869	$S_{DFE}=22$ $R_{DFE}=1028$ $v_{DFE}=9$
Case 2$c_1=22 ; c_2=10$$R_c=0.258$	29	110	1800	$S_{DFE}=19$ $R_{DFE}=1031$ $v_{DFE}=9$
Case 3$c_1=19 ; c_2=10$$R_c=0.2233$	25	101	1740	$S_{DFE}=16$ $R_{DFE}=1034$ $v_{DFE}=9$
Case 4$c_1=25 ; c_2=15$$R_c=0.1963$	22	94	1693	$S_{DFE}=14$ $R_{DFE}=1036$ $v_{DFE}=9$
Case 5$c_1=25 ; c_2=17$$R_c=0.1735$	20	88	1651	$S_{DFE}=13$ $R_{DFE}=1037$ $v_{DFE}=9$
Case 6$c_1=25 ; c_2=20$$R_c=0.1477$	18	82	1611	$S_{DFE}=11$ $R_{DFE}=1039$ $v_{DFE}=9$
Case 7$c_1=25 ; c_2=25$$R_c=0.1184$	16	76	1569	$S_{DFE}=9$ $R_{DFE}=1041$ $v_{DFE}=9$

Table 4

Specifications of the model behaviour for the four considered pairs of values of the control parameters $c_5$ and $c_6$ governing the treatment effort.

	${\mathit{c}}_5=2$ ${\mathit{c}}_6=0.24$ ${\overline{\mathit{R}}}_{\mathit{c}}=3.635$	${\mathit{c}}_5=2$ ${\mathit{c}}_6=0.27$ ${\overline{\mathit{R}}}_{\mathit{c}}=3.9644$	${\mathit{c}}_5=2$ ${\mathit{c}}_6=0.3$ ${\overline{\mathit{R}}}_{\mathit{c}}=4.2937$	${\mathit{c}}_5=2$ ${\mathit{c}}_6=0.33$ ${\overline{\mathit{R}}}_{\mathit{c}}=4.6231$
Infectious peak	424	395	369	342
Treatment peak	51	53	55	56
$S_{EE}$	255	280	304	329
$I_{EE}$	36	32	29	26
$R_{EE}$	710	694	677	659
$N_{EE}$	1001	1006	1010	1014
$\frac{I_{EE}}{N_{EE}}\times 100$	3.59%	3.2%	2.87%	2.59%
$t_r{}_{EE}$	4	4	4	4

Table 5

Specifications of the model behaviour for the four considered values of the parameter $\rho$.

.	$\mathit{\rho }=\frac{1}{100}$	$\mathit{\rho }=\frac{1}{120}$	$\mathit{\rho }=\frac{1}{150}$	$\mathit{\rho }=0$
Infectious peak	262	261	259	253
Transient duration (days)	91	72	64	49
Vaccination peak	27	27	27	27
Treatment peak	35	35	35	34
SDFE	277	242	203	1
RDFE	773	808	847	1049
$v_{DFE}$	8	7	6	0
Vaccines	615	467	400	290
Medicaments	589	570	556	512

Appendix A

First, the non-negativity of $S\left(t\right)$ is proved by contradiction. Suppose that the result is false because there exists some $t_s>0$ such that $S\left(t\right)\ge 0$, $I\left(t\right)\ge 0$, $R\left(t\right)\ge 0$, $v\left(t\right)\ge 0$, $t_r\left(t\right)\ge 0$ for all 0 $\le t<t_s$ and $S(t_s)<0$. Equations (1) and (4) can be written as:(A1)$$\begin{gathered} \left[\begin{array}{c}\dot{S}\left(t\right)\\ \dot{v}\left(t\right)\end{array}\right]=\left[\begin{array}{cc}-\mu -\beta \frac{I\left(t\right)}{N\left(t\right)}& -1\\ c_2+c_4\frac{I\left(t\right)}{N\left(t\right)}& -c_1\end{array}\right]\left[\begin{array}{c}S\left(t\right)\\ v\left(t\right)\end{array}\right]+\left[\begin{array}{c}b\left(1-q\right)+\rho R\left(t\right)\\ c_{3}I\left(t\right)\end{array}\right] \\ =\left[\begin{array}{cc}-\left(\mu +\beta \right)& -1\\ c_2+c_4& -c_1\end{array}\right]\left[\begin{array}{c}S\left(t\right)\\ v\left(t\right)\end{array}\right]+\left[\begin{array}{c}b\left(1-q\right)+\rho R\left(t\right)+\beta \frac{S\left(t\right)\left[S\left(t\right)+R\left(t\right)\right]}{N\left(t\right)}\\ c_{3}I\left(t\right)-c_4\frac{S\left(t\right)\left[S\left(t\right)+R\left(t\right)\right]}{N\left(t\right)}\end{array}\right] \end{gathered}$$ where the fact that $I\left(t\right)=N\left(t\right)-S\left(t\right)-R\left(t\right)$ has been used. Then, the subsystem (A1) may be compactly written as:(A2)$$\left[\begin{array}{c}\dot{S}\left(t\right)\\ \dot{v}\left(t\right)\end{array}\right]=A^{SV}\left[\begin{array}{c}S\left(t\right)\\ v\left(t\right)\end{array}\right]+\left[\begin{array}{c}{u}_1\left(t\right)\\ u_2\left(t\right)\end{array}\right]$$ with:(A3)$$A^{SV}=\left[\begin{array}{cc}-\left(\mu +\beta \right)& -1\\ c_2+c_4& -c_1\end{array}\right]$$ and:(A4)$$u_1\left(t\right)=b\left(1-q\right)+\rho R\left(t\right)+\beta \frac{S\left(t\right)\left[S\left(t\right)+R\left(t\right)\right]}{N\left(t\right)} ; u_2\left(t\right)=c_{3}I\left(t\right)-c_4\frac{S\left(t\right)\left[S\left(t\right)+R\left(t\right)\right]}{N\left(t\right)}.$$

From (A2), it follows that:(A5)$$\left[\begin{array}{c}S\left(t\right)\\ v\left(t\right)\end{array}\right]={\varphi }^{SV}\left(t\right)\left[\begin{array}{c}S\left(0\right)\\ v\left(0\right)\end{array}\right]+{{\displaystyle \int }}_0^t{\varphi }^{SV}\left(t-\tau \right)\left[\begin{array}{c}{u}_1\left(\tau \right)\\ u_2\left(\tau \right)\end{array}\right]d\tau \forall t\ge 0$$ with ${\varphi }^{SV}\left(t\right)=e^{A^{SV}t}=L^{-1}\left\{{\left(s{\mathbb{I}}_2-A^{SV}\right)}^{-1}\right\}$, where $L^{-1}$ denotes the inverse of the Laplace transform, ${\mathbb{I}}_2$ is the second-order identity matrix and $s$ is the Laplace variable in the complex domain. By direct calculations, one obtains that:(A6)$$\begin{gathered} {\varphi }_{11}^{SV}\left(t\right)=\frac{\left(c_1+{\lambda }_1^{SV}\right)e^{{\lambda }_1^{SV}t}-\left(c_1+{\lambda }_2^{SV}\right)e^{{\lambda }_2^{SV}t}}{{\lambda }_1^{SV}-{\lambda }_2^{SV}} ; {\varphi }_{12}^{SV}\left(t\right)=-\frac{e^{{\lambda }_1^{SV}t}-e^{{\lambda }_2^{SV}t}}{{\lambda }_1^{SV}-{\lambda }_2^{SV}} \\ {\varphi }_{21}^{SV}\left(t\right)=\frac{\left(c_2+c_4\right)\left(e^{{\lambda }_1^{SV}t}-e^{{\lambda }_2^{SV}t}\right)}{{\lambda }_1^{SV}-{\lambda }_2^{SV}} ; {\varphi }_{22}^{SV}\left(t\right)=\frac{\left(\mu +\beta +{\lambda }_1^{SV}\right)e^{{\lambda }_1^{SV}t}-\left(\mu +\beta +{\lambda }_2^{SV}\right)e^{{\lambda }_2^{SV}t}}{{\lambda }_1^{SV}-{\lambda }_2^{SV}} \end{gathered}$$ where:(A7)$${\lambda }_{1,2}^{SV}=\frac{1}{2}\left[-\left(\mu +\beta +c_1\right)\pm \sqrt{{\left[c_1-\left(\mu +\beta \right)\right]}^2-4\left(c_2+c_4\right)}\right]$$ are the eigenvalues of $A^{SV}$, with ${\lambda }_1^{SV}$ corresponding to the sign “+” and ${\lambda }_2^{SV}$ corresponding to the sign “$-$”. If the control parameters $c_1$, $c_2$, and $c_4$ fulfil condition (i), then ${\lambda }_1^{SV}$ and ${\lambda }_2^{SV}$ are real, and ${\lambda }_2^{SV}<{\lambda }_1^{SV}<0$. From (A4) and (A5), it follows that:(A8)$$\begin{gathered} S\left(t\right)=\frac{\left(c_1+{\lambda }_1^{SV}\right)e^{{\lambda }_1^{SV}t}-\left(c_1+{\lambda }_2^{SV}\right)e^{{\lambda }_2^{SV}t}}{{\lambda }_1^{SV}-{\lambda }_2^{SV}}S\left(0\right)-\frac{e^{{\lambda }_1^{SV}t}-e^{{\lambda }_2^{SV}t}}{{\lambda }_1^{SV}-{\lambda }_2^{SV}}v\left(0\right) \\ +\frac{1}{{\lambda }_1^{SV}-{\lambda }_2^{SV}}{{\displaystyle \int }}_0^t\left[\left(c_1+{\lambda }_1^{SV}\right)e^{{\lambda }_1^{SV}\left(t-\tau \right)}-\left(c_1+{\lambda }_2^{SV}\right)e^{{\lambda }_2^{SV}\left(t-\tau \right)}\right]u_1\left(\tau \right)d\tau \\ -\frac{1}{{\lambda }_1^{SV}-{\lambda }_2^{SV}}{{\displaystyle \int }}_0^t\left[e^{{\lambda }_1^{SV}\left(t-\tau \right)}-e^{{\lambda }_2^{SV}\left(t-\tau \right)}\right]u_2\left(\tau \right)d\tau . \end{gathered}$$

From (A4), (A7), (A8), and the condition $v\left(0\right)=0$, one obtains:(A9)$$\begin{gathered} S(t_s)\ge \frac{\left(c_1+{\lambda }_1^{SV}\right)e^{{\lambda }_1^{SV}{t}_s}-\left(c_1+{\lambda }_2^{SV}\right)e^{{\lambda }_2^{SV}{t}_s}}{{\lambda }_1^{SV}-{\lambda }_2^{SV}}S\left(0\right) \\ +\frac{b\left(1-q\right)}{{\lambda }_1^{SV}-{\lambda }_2^{SV}}{{\displaystyle \int }}_0^{t_s}\left[\left(c_1+{\lambda }_1^{SV}\right)e^{{\lambda }_1^{SV}\left(t_s-\tau \right)}-\left(c_1+{\lambda }_2^{SV}\right)e^{{\lambda }_2^{SV}\left(t_s-\tau \right)}\right]d\tau \\ -\frac{c_3}{{\lambda }_1^{SV}-{\lambda }_2^{SV}}{{\displaystyle \int }}_0^{t_s}\left[e^{{\lambda }_1^{SV}\left(t_s-\tau \right)}-e^{{\lambda }_2^{SV}\left(t_s-\tau \right)}\right]I\left(\tau \right)d\tau , \end{gathered}$$ where the facts that $S\left(t\right)\ge 0$, $I\left(t\right)\ge 0$, $R\left(t\right)\ge 0$, $N\left(t\right)\ge 0$, and $\frac{S\left(t\right)}{N\left(t\right)}\le 1$ jointly to $e^{{\lambda }_1^{SV}t}>e^{{\lambda }_2^{SV}t}$, for all 0 $\le t<t_s$, provided that condition (i) is satisfied by the control parameters, have been used. Then, from (A9) and by direct calculations, one obtains:(A10)$$\begin{gathered} S(t_s)\ge \frac{\left(c_1+{\lambda }_1^{SV}\right)e^{{\lambda }_1^{SV}{t}_s}-\left(c_1+{\lambda }_2^{SV}\right)e^{{\lambda }_2^{SV}{t}_s}}{{\lambda }_1^{SV}-{\lambda }_2^{SV}}S\left(0\right) \\ +\frac{b\left(1-q\right)}{{\lambda }_1^{SV}-{\lambda }_2^{SV}}\left[-\frac{c_1+{\lambda }_1^{SV}}{{\lambda }_1^{SV}}\left(1-e^{{\lambda }_1^{SV}{t}_s}\right)+\frac{c_1+{\lambda }_2^{SV}}{{\lambda }_2^{SV}}\left(1-e^{{\lambda }_2^{SV}{t}_s}\right)\right] \\ -\frac{c_3{I}_{max}}{{\lambda }_1^{SV}-{\lambda }_2^{SV}}\left[-\frac{1-e^{{\lambda }_1^{SV}{t}_s}}{{\lambda }_1^{SV}}+\frac{1-e^{{\lambda }_2^{SV}{t}_s}}{{\lambda }_2^{SV}}\right], \end{gathered}$$ where $I_{max}=\underset{0\le t<\infty }{\max }\left\{I\left(t\right)\right\}\ge \underset{0\le t<t_s}{\max }\left\{I\left(t\right)\right\}\ge 0$. From (A10), it follows that:(A11)$$\begin{gathered} S(t_s)\ge \frac{\left(c_1+{\lambda }_1^{SV}\right)\left[\left|{\lambda }_1^{SV}\right|S\left(0\right)-b\left(1-q\right)\right]+c_3{I}_{max}}{\left({\lambda }_1^{SV}-{\lambda }_2^{SV}\right)\left|{\lambda }_1^{SV}\right|}{e}^{{\lambda }_1^{SV}{t}_s} \\ -\frac{\left(c_1+{\lambda }_2^{SV}\right)\left[\left|{\lambda }_2^{SV}\right|S\left(0\right)-b\left(1-q\right)\right]+c_3{I}_{max}}{\left({\lambda }_1^{SV}-{\lambda }_2^{SV}\right)\left|{\lambda }_2^{SV}\right|}{e}^{{\lambda }_2^{SV}{t}_s}+\frac{b\left(1-q\right)c_1-c_3{I}_{max}}{{\lambda }_1^{SV}{\lambda }_2^{SV}}. \end{gathered}$$

Now, by introducing (A7) in (A11) and since $c_1>\mu +\beta$ from condition (i), direct calculations lead to:(A12)$$\begin{array}{ll}S(t_s)\ge & \frac{\left[\left(\mu +\beta \right)S\left(0\right)-b\left(1-q\right)\right]\sqrt{{\left[c_1-\left(\mu +\beta \right)\right]}^2-4\left(c_2+c_4\right)}}{2\left({\lambda }_1^{SV}-{\lambda }_2^{SV}\right)}\left(\frac{e^{{\lambda }_1^{SV}{t}_s}}{\left|{\lambda }_1^{SV}\right|}+\frac{e^{{\lambda }_2^{SV}{t}_s}}{\left|{\lambda }_2^{SV}\right|}\right)+\\ & \frac{\left[c_1-\left(\mu +\beta \right)\right]\left[\left(\mu +\beta \right)S\left(0\right)-b\left(1-q\right)\right]+2\left[\left(c_2+c_4\right)S\left(0\right)+c_3{I}_{max}\right]}{2\left({\lambda }_1^{SV}-{\lambda }_2^{SV}\right)}\left(\frac{e^{{\lambda }_1^{SV}{t}_s}}{\left|{\lambda }_1^{SV}\right|}-\frac{e^{{\lambda }_2^{SV}{t}_s}}{\left|{\lambda }_2^{SV}\right|}\right)+\frac{b\left(1-q\right)\left(\mu +\beta \right)-c_3{I}_{max}}{{\lambda }_1^{SV}{\lambda }_2^{SV}}.\end{array}$$

Then, $S(t_s)\ge 0$ from the conditions (i)–(iii) jointly with the facts that ${\lambda }_2^{SV}<{\lambda }_1^{SV}<0$ and $e^{{\lambda }_1^{SV}{t}_s}>e^{{\lambda }_2^{SV}{t}_s}$. Such a result contradicts the existence of a time instant $t_s>0$ such that $S(t_s)<0$. Then, the non-negativity of $S\left(t\right)$ is proved. Now, the non-negativity of $I\left(t\right)$ is proved by contradiction. Suppose that the result is false because there exists some $t_I>0$ such that $S\left(t\right)\ge 0$, $I\left(t\right)\ge 0$, $R\left(t\right)\ge 0$, $v\left(t\right)\ge 0$, $t_r\left(t\right)\ge 0$ for all 0 $\le t<t_I$ and $I(t_I)<0$. Then, Equations (2) and (5) can be compactly written as:(A13)$$\left[\begin{array}{c}\dot{I}\left(t\right)\\ \dot{t_r}\left(t\right)\end{array}\right]=A^{IT}\left[\begin{array}{c}I\left(t\right)\\ t_r\left(t\right)\end{array}\right]+\left[\begin{array}{c}1\\ 0\end{array}\right]u_3\left(t\right)$$ where $u_3\left(t\right)=\beta \frac{S\left(t\right)I\left(t\right)}{N\left(t\right)}$ and:(A14)$$A^{IT}=\left[\begin{array}{cc}-\left(\mu +\alpha +\gamma \right)& -1\\ c_6& -c_5\end{array}\right].$$

From (A14), it follows that:(A15)$$\left[\begin{array}{c}I\left(t\right)\\ t_r\left(t\right)\end{array}\right]={\varphi }^{IT}\left(t\right)\left[\begin{array}{c}I\left(0\right)\\ t_r\left(0\right)\end{array}\right]+{{\displaystyle \int }}_0^t{\varphi }^{IT}\left(t-\tau \right)\left[\begin{array}{c}1\\ 0\end{array}\right]u_3\left(\tau \right)d\tau \forall t\ge 0,$$ with ${\varphi }^{IT}\left(t\right)=e^{A^{IT}t}=L^{-1}\left\{{\left(s{\mathbb{I}}_2-A^{IT}\right)}^{-1}\right\}$. By direct calculations, one obtains that:(A16)$$\begin{gathered} {\varphi }_{11}^{IT}\left(t\right)=\frac{\left(c_5+{\lambda }_1^{IT}\right)e^{{\lambda }_1^{IT}t}-\left(c_5+{\lambda }_2^{IT}\right)e^{{\lambda }_2^{IT}t}}{{\lambda }_1^{IT}-{\lambda }_2^{IT}} ; {\varphi }_{12}^{IT}\left(t\right)=-\frac{e^{{\lambda }_1^{IT}t}-e^{{\lambda }_2^{IT}t}}{{\lambda }_1^{IT}-{\lambda }_2^{IT}} \\ {\varphi }_{21}^{IT}\left(t\right)=\frac{c_6\left(e^{{\lambda }_1^{IT}t}-e^{{\lambda }_2^{IT}t}\right)}{{\lambda }_1^{IT}-{\lambda }_2^{IT}} ; {\varphi }_{22}^{IT}\left(t\right)=\frac{\left(\mu +\alpha +\gamma +{\lambda }_1^{IT}\right)e^{{\lambda }_1^{IT}t}-\left(\mu +\alpha +\gamma +{\lambda }_2^{IT}\right)e^{{\lambda }_2^{IT}t}}{{\lambda }_1^{IT}-{\lambda }_2^{IT}} \end{gathered}$$ where:(A17)$${\lambda }_{1,2}^{IT}=\frac{1}{2}\left[-\left(\mu +\alpha +\gamma +c_5\right)\pm \sqrt{{\left[c_5-\left(\mu +\alpha +\gamma \right)\right]}^2-4c_6}\right]$$ are the eigenvalues of $A^{IT}$, with ${\lambda }_1^{IT}$ corresponding with the sign “+” and ${\lambda }_2^{IT}$ corresponding with the sign “$-$”. If the control parameters $c_5$ and $c_6$ fulfill condition (iv), then both eigenvalues are real, and ${\lambda }_2^{IT}<{\lambda }_1^{IT}<0$. From (A15) and (A16), it follows that:(A18)$$\begin{array}{r}I\left(t_I\right)=\frac{\left(c_5+{\lambda }_1^{IT}\right)e^{{\lambda }_1^{IT}{t}_I}-\left(c_5+{\lambda }_2^{IT}\right)e^{{\lambda }_2^{IT}{t}_I}}{{\lambda }_1^{IT}-{\lambda }_2^{IT}}I\left(0\right)-\frac{e^{{\lambda }_1^{IT}{t}_I}-e^{{\lambda }_2^{IT}{t}_I}}{{\lambda }_1^{IT}-{\lambda }_2^{IT}}{t}_r\left(0\right)\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\\ +\frac{\beta }{{\lambda }_1^{IT}-{\lambda }_2^{IT}}{{\displaystyle \int }}_0^{t_I}\left[\left(c_5+{\lambda }_1^{IT}\right)e^{{\lambda }_1^{IT}(t_I-\tau )}-\left(c_5+{\lambda }_2^{IT}\right)e^{{\lambda }_2^{IT}(t_I-\tau )}\right]\frac{S\left(\tau \right)I\left(\tau \right)}{N\left(\tau \right)}d\tau .\end{array}$$

Condition (iv) about the control parameters $c_5$ and $c_6$ implies that $c_5+{\lambda }_1^{IT}>0$. Under such a condition jointly with $t_r\left(0\right)=0$, one obtains that $I\left(t_I\right)\ge 0$ from (A18) since $e^{{\lambda }_1^{IT}{t}_I}>e^{{\lambda }_2^{IT}{t}_I}>0$ and $\frac{S\left(t\right)I\left(t\right)}{N\left(t\right)}\ge 0$ for all 0 $\le t<t_I$. Then, such a result contradicts the existence of a time instant $t_I>0$ such that $I(t_I)<0$. Then, the non-negativity of $I\left(t\right)$ is proved. Now, the non-negativity of $R\left(t\right)$ is proved by contradiction. Suppose that the result is false because there exists some $t_R>0$ such that $S\left(t\right)\ge 0$, $I\left(t\right)\ge 0$, $R\left(t\right)\ge 0$, $v\left(t\right)\ge 0$, and $t_r\left(t\right)\ge 0$ for all 0 $\le t<t_R$ and $R(t_R)<0$. From (3), it follows that:(A19)$$R\left(t\right)=e^{-\left(\mu +\rho \right)t}R\left(0\right)+{{\displaystyle \int }}_0^t{e}^{-\left(\mu +\rho \right)\left(t-\tau \right)}\left[bq+\gamma I\left(\tau \right)+v\left(\tau \right)+t_r\left(\tau \right)\right]d\tau \forall t\ge 0.$$

Then, it follows that $R(t_R)\ge 0$ for any $R\left(0\right)\ge 0$. Such a result contradicts the existence of a time instant $t_R>0$ such that $R(t_R)<0$. Then, the non-negativity of $R\left(t\right)$ is proved. Now, the non-negativity of $v\left(t\right)$ is proved by contradiction. Suppose that the result is false because there exists some $t_V>0$ such that $S\left(t\right)\ge 0$, $I\left(t\right)\ge 0$, $R\left(t\right)\ge 0$, $v\left(t\right)\ge 0$, $t_r\left(t\right)\ge 0$ for all 0 $\le t<t_V$ and $v(t_V)<0$. From (4), it follows that:(A20)$$v\left(t\right)=e^{-c_{1}t}v\left(0\right)+{{\displaystyle \int }}_0^t{e}^{-c_1\left(t-\tau \right)}\left[c_{2}S\left(\tau \right)+c_{3}I\left(\tau \right)+c_4\frac{S\left(\tau \right)I\left(\tau \right)}{N\left(\tau \right)}\right]d\tau \forall t\ge 0.$$

Then, it follows that $v(t_V)\ge 0$ for any $v\left(0\right)\ge 0$. Such a fact contradicts the existence of a time instant $t_V>0$ such that $v(t_V)<0$. Then, the non-negativity of $v\left(t\right)$ is proved. Finally, the non-negativity of $t_r\left(t\right)$ is proved by contradiction. Suppose that the result is false because there exists some $t_T>0$ such that $S\left(t\right)\ge 0$, $I\left(t\right)\ge 0$, $R\left(t\right)\ge 0$, $v\left(t\right)\ge 0$, $t_r\left(t\right)\ge 0$ for all 0 $\le t<t_T$ and $t_r(t_T)<0$. From (5), it follows that:(A21)$$t_r\left(t\right)=e^{-c_{5}t}{t}_r\left(0\right)+c_6{{\displaystyle \int }}_0^t{e}^{-c_5\left(t-\tau \right)}I\left(\tau \right)d\tau \forall t\ge 0.$$

Then, it follows that $t_r(t_T)\ge 0$ for any $t_r\left(0\right)\ge 0$, since $I\left(t\right)\ge 0$ $\forall t\ge 0$. Such a result contradicts the existence of a time instant $t_T>0$ such that $t_r(t_T)<0$. Then, the non-negativity of $t_r\left(t\right)$ is proved. In summary, none of the controlled epidemic model variables take negative values.

References

1. Kermack, W.O.; McKendrick, A.G. Contributions to the mathematical theory of epidemics, part 1. Proc. R. Soc. A; 1927; 115, pp. 700-721.

2. Li, J.; Ma, Z. Global analysis of SIS epidemic models with variable total population size. Math. Comput. Model.; 2004; 39, pp. 1231-1242.

3. Pang, J.; Cui, J.; Zhou, X. Dynamical behavior of a hepatitis B virus transmission model with vaccination. J. Theor. Biol.; 2010; 265, pp. 572-578. [DOI: https://dx.doi.org/10.1016/j.jtbi.2010.05.038]

4. Meng, X.; Chen, L.; Cheng, H. Two profitless delays for the SEIRS epidemic disease model with nonlinear incidence and pulse vaccination. Appl. Math. Comput.; 2007; 186, pp. 516-529. [DOI: https://dx.doi.org/10.1016/j.amc.2006.07.124]

5. Cai, L.; Li, X.; Ghosh, M.; Guo, B. Stability analysis of an HIV/AIDS epidemic model with treatment. J. Comput. Appl. Math.; 2009; 229, pp. 313-323. [DOI: https://dx.doi.org/10.1016/j.cam.2008.10.067]

6. De la Sen, M.; Alonso-Quesada, S.; Ibeas, A. On the stability of an SEIR epidemic model with distributed time-delay and a general class of feedback vaccination rules. Appl. Math. Comput.; 2015; 270, pp. 953-976. [DOI: https://dx.doi.org/10.1016/j.amc.2015.08.099]

7. Trawicki, M.B. Deterministic SEIRS epidemic model for modeling vital dynamics, vaccinations, and temporary immunity. Mathematics; 2017; 5, 7. [DOI: https://dx.doi.org/10.3390/math5010007]

8. Gbadamosi, B.; Ojo, M.M.; Oke, S.I.; Matadi, M.B. Qualitative analysis of a Dengue fever model. Math. Comput. Appl.; 2018; 23, 33. [DOI: https://dx.doi.org/10.3390/mca23030033]

9. Tulu, T.W.; Tian, B.; Wu, Z. Modeling the effect of quarantine and vaccination on Ebola disease. Adv. Differ. Equ.; 2017; 2017, 178. [DOI: https://dx.doi.org/10.1186/s13662-017-1225-z]

10. Buonomo, B.; Della Marca, R. Oscillations and hysteresis in an epidemic model with information-dependent imperfect vaccination. Math. Comput. Simul.; 2019; 162, pp. 97-114. [DOI: https://dx.doi.org/10.1016/j.matcom.2019.01.005]

11. Ullah, S.; Khan, M.A.; Farooq, M.; Gul, T. Modeling and analysis of tuberculosis (TB) in Khyber Pakhtunkhwa, Pakistan. Math. Comput. Simul.; 2019; 165, pp. 181-199. [DOI: https://dx.doi.org/10.1016/j.matcom.2019.03.012]

12. Liu, Q.; Jiang, D.; Hayat, T.; Alsaedi, A. Dynamics of a stochastic multigroup SIQR epidemic model with standard incidence rates. J. Frankl. Inst.; 2019; 356, pp. 2960-2993. [DOI: https://dx.doi.org/10.1016/j.jfranklin.2019.01.038]

13. Nkamba, L.N.; Ntaganda, J.M.; Abboubakar, H.; Kamgang, J.C.; Castelli, L. Global stability of a SVEIR epidemic model: Application to poliomyelitis transmission dynamics. Open J. Model. Simul.; 2017; 5, pp. 98-112. [DOI: https://dx.doi.org/10.4236/ojmsi.2017.51008]

14. Zaman, G.; Kang, Y.H.; Cho, G.; Jung, I.H. Optimal strategy of vaccination & treatment in an SIR epidemic model. Math. Comput. Simul.; 2017; 136, pp. 63-77.

15. Alonso-Quesada, S.; De la Sen, M.; Ibeas, A. On the discretization and control of an SEIR epidemic model with a periodic impulsive vaccination. Commun. Nonlinear Sci. Numer. Simul.; 2017; 42, pp. 247-274. [DOI: https://dx.doi.org/10.1016/j.cnsns.2016.05.027]

16. Ling, L.; Jiang, G.; Long, T. The dynamics of an SIS epidemic model with fixed-time birth pulses and state feedback pulse treatments. Appl. Math. Model.; 2015; 39, pp. 5579-5591. [DOI: https://dx.doi.org/10.1016/j.apm.2015.01.022]

17. Huang, H.; Chen, Y.; Yan, Z. Impacts of social distancing on the spread of infectious diseases with asymptomatic infection: A mathematical model. Appl. Math. Comput.; 2021; 398, pp. 1-13. [DOI: https://dx.doi.org/10.1016/j.amc.2021.125983]

18. Nadim, S.S.; Ghosh, I.; Chattopadhyay, J. Short-term predictions and prevention strategies for COVID-19: A model based study. Appl. Math. Comput.; 2021; 404, pp. 1-19. [DOI: https://dx.doi.org/10.1016/j.amc.2021.126251]

19. Han, D.; Shao, Q.; Li, D.; Sun, M. How the individuals’ risk aversion affect the epidemic spreading. Appl. Math. Comput.; 2020; 369, pp. 1-10. [DOI: https://dx.doi.org/10.1016/j.amc.2019.124894]

20. Tanimoto, J. Sociophysics Approach to Epidemics; Springer: Singapore, 2021.

21. Kabir, K.M.A.; Tanimoto, J. Modelling and analysing the coexistence of dual dilemmas in the proactive vaccination game and retroactive treatment game in epidemic viral dynamics. Proc. R. Soc. A; 2019; 475, pp. 1-20. [DOI: https://dx.doi.org/10.1098/rspa.2019.0484]

22. Moualeu, D.P.; Weiser, M.; Ehrig, R.; Deuflhard, P. Optimal control for a tuberculosis model with undetected cases in Cameroon. Commun. Nonlinear Sci. Numer. Simul.; 2015; 20, pp. 986-1003. [DOI: https://dx.doi.org/10.1016/j.cnsns.2014.06.037]

23. Sharma, S.; Singh, F. Bifurcation and stability analysis of a cholera model with vaccination and saturated treatment. Chaos Solitons Fractals; 2021; 146, pp. 1-15. [DOI: https://dx.doi.org/10.1016/j.chaos.2021.110912]

24. den Driessche, P.V.; Watmough, J. Reproduction numbers and sub-threshold endemic equilibria for compartmental models of disease transmission. Math. Biosci.; 2002; 180, pp. 29-48. [DOI: https://dx.doi.org/10.1016/S0025-5564(02)00108-6]

25. Zhou, L.; Fan, M. Dynamics of an SIR epidemic model with limited medical resources revisited. Nonlinear Anal. Real World Appl.; 2012; 13, pp. 312-324. [DOI: https://dx.doi.org/10.1016/j.nonrwa.2011.07.036]

26. Dumont, Y.; Chiroleu, F.; Domerg, C. On a temporal model for the Chikungunea disease: Modeling, theory and numerics. Math. Biosci.; 2008; 213, pp. 80-91. [DOI: https://dx.doi.org/10.1016/j.mbs.2008.02.008] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/18394655]

27. Mwasunda, J.A.; Irunde, J.I.; Kajunguri, D.; Kuznetsov, D. Modeling and analysis of taeniasis and cysticercosis transmission dynamics in human, pigs and cattle. Adv. Differ. Equ.; 2021; 2021, 176. [DOI: https://dx.doi.org/10.1186/s13662-021-03341-9]

28. Vidyasagar, M. Nonlinear Systems Analysis; 2nd ed. Prentice Hall: Englewood Cliffs, NJ, USA, 1993.