Abstract
We proposed a deterministic compartmental model for the transmission dynamics of COVID-19 disease. We performed qualitative and quantitative analysis of the deterministic model concerning the local and global stability of the disease-free and endemic equilibrium points. We found that the disease-free equilibrium is locally asymptotically stable when the basic reproduction number is less than unity, while the endemic equilibrium point becomes locally asymptotically stable if the basic reproduction number is above unity. Furthermore, we derived the global stability of both the disease-free and endemic equilibriums of the system by constructing some Lyapunov functions. If , it is found that the disease-free equilibrium is globally asymptotically stable, while the endemic equilibrium point is globally asymptotically stable when . The numerical results of the general dynamics are in agreement with the theoretical solutions. We established the optimal control strategy by using Pontryagin’s maximum principle. We performed numerical simulations of the optimal control system to investigate the impact of implementing different combinations of optimal controls in controlling and eradicating COVID-19 disease. From this, a significant difference in the number of cases with and without controls was observed. We observed that the implementation of the combination of the control treatment rate, , and the control treatment rate, , has shown effective and efficient results in eradicating COVID-19 disease in the community relative to the other strategies.
1. Introduction
Since the emergence of the disease in December 2019, the coronavirus disease 2019 has been one of the greatest global threats, causing diverse political, economic, and social crises. It is now officially known as COVID-19 by the World Health Organization (WHO). The disease was announced as a pandemic on March 11, 2020, by WHO [1]. The burden of the disease in developing nations like Africa is very serious. The main causes of transmission of COVID-19 are sneezing, coughing, contacting infected people, touching items used daily, and so on [2].
Many countries are involved in taking various measures to prevent the pandemic, such as quarantining of suspected individuals, isolation of the infected individuals, lockdown of the community, contact tracing, mask wearing, and physical distancing in collaboration with the WHO. Even though the trial of the production of the medication was not successful, the production of the vaccination seems effective. Vaccinating the susceptible population has officially started globally, including for developing nations.
Mathematical modelling of infectious diseases plays a vital role in understanding the underlying mechanisms by which an infectious disease spreads, in forecasting the long-term effect of an epidemic, and in suggesting intervention methods for controlling an existing disease [2–6]. Currently, several authors have studied a mathematical model of COVID-19 dynamics. The majority of the studies on COVID-19 assumed the proposed model as an autonomous or nonautonomous system. Modelling studies such as by [7–9] assumed a deterministic autonomous dynamical system and others considered fractional-order deterministic system [10–13].
There are also more studies that have been carried out on the COVID-19 dynamics system using optimal control strategies [13–19]. In [13], the authors proposed a mathematical model involving lock-down, quarantine, and self-isolating control functions to minimize the burden of the disease on society and the associated costs. The authors in [14] developed a deterministic mathematical model with time variable controls. The results of their study revealed that an optimal implementation of personal protective measures reduced the transmission of COVID-19. Quarantine and treatment of the infected individuals are considered optimal controls in the model stated by Abbasi et al. [15]. In addition, the authors employed an impulsive epidemic model to show the abrupt growth in the population. The model proposed by Kada et al. [16] considered two control strategies in their dynamical system. That is, the treatment of COVID-19-infected individuals and using masks for the sensitive body parts, respectively. They found that implementation of the combined control strategies exhibited an effective strategy in reducing COVID-19 epidemic disease in society compared to implementing the individual strategy.
A few others proposed a nonautonomous optimal control system [20–22]. In addition to the prevention and treatment measures being implemented so far, many countries have launched a nationwide vaccination campaign against COVID-19, which is a huge step towards reducing the spread of the disease. Therefore, it is vital to study COVID-19 disease dynamics using mathematical modelling, considering it as an autonomous and nonautonomous system with the implementation of vaccination, treatment, and public health education as constant and continuous controls.
Therefore, in this paper, we propose a deterministic compartmental model for the dynamics of COVID-19 in a single host population and investigate the mathematical and numerical analysis of the model. More importantly, it establishes the global stability of the disease-free and endemic equilibria and the implication of vaccination, treatment, and public health education on the disease dynamics. The system is further extended into an optimal control problem with the implementation of continuous controls: vaccination, treatment, and public health education. And finally, the best control measure for controlling the disease is recommended.
The next part of the paper is organized as follows. Section 2 presents the new dynamics model for COVID-19 disease. In Section 3, we study the qualitative analysis including the basic reproduction number, the existence of the disease-free and endemic equilibria, their local and global stability, and the sensitivity of the parameters. The proposed model is further extended into an optimal control problem and is investigated using optimal control theory in Section 4. The numerical methods employed and the numerical simulation of the system with parameter values are performed in Section 5. Finally, the conclusion and recommendation are briefly presented in Section 6.
2. Model Formulation
The human population is divided into four subclasses. These subclasses of human populations are as follows: susceptible population , infected population , infected individuals under intensive care unit (ICU) (), and recovered population . The total population at a time is given by . The susceptible population is increased by the constant recruitment rate . On the other hand, it is reduced by the rate of transmission , vaccination rate , and natural death rate . The infected population is generated proportional to the susceptible and infected individuals at a rate of . However, the infected compartment is decreased by the rate of progress from to , rate of recovery , and natural death rate . The infected population is increased by the transmission rate but declined by the recovery rate and natural death rate . Moreover, the recovery compartment gains population from vaccination of susceptible population at the rate , recovery of infected individual at the rate , and recovery of infected individual at the rate . The class declined through a natural death rate .
The following assumptions are considered in the formulation of the model: the population is homogeneously mixed, the exposed compartment is not considered, the infected population is the source of infection, but infected individuals under intensive care unit (IUC) are assumed to be hospitalized and restricted, and people contracted with these groups are assumed to take all forms of protection and are thus not assumed to be the source of infection.
Following the population scheme in Figure 1, the differential equations that describe the transmission of COVID-19 become
with initial conditions
The variable is not involved in the first three equations of the system Equation (1) and does not influence the dynamics of the first three equations. Hence, we only consider the following reduced system in the analysis that follows:
Descriptions of the state variables and parameters are illustrated in Abbreviations, and the schematic diagram for the system is illustrated in Figure 1.
3. Mathematical Analysis of the Model
In this section, we must demonstrate that the system Equation (1) is biologically feasible and realistic whenever all of the system Equation (1) state variables are nonnegative for all time .
Theorem 1. The solution , , and of system Equation (1) with initial conditions and remain positive for all . Furthermore,
Proof. Using a similar approach to [23], we obtain
One can observe that the above rates are all nonnegative on the bounding plane . Then, it is easy to show that the region is positively invariant and attracting [24]. The region of attracting for system Equation (1) can be given as
Thus, the closed set is positively invariant and attracting region with respect to system Equation (1). The basic reproduction number, denoted by , is the very important quantity for qualitative analysis of the dynamics model. It is stated as the average number of secondary cases produced by a single infected individual introduced into a completely susceptible population [7, 24].
Using the next-generation matrix approach, as stated in [25], it is easy to determine the basic reproduction number of system Equation (3). One can easily find the disease-free equilibrium of system Equation (3) as .
Let the state variables of system Equation (3) denoted by . Thus, system Equation (2) can be rewritten as , where is the rate of new infection and is the rate of transfer infection compartment into and out of compartment , and are, respectively, given by
The Jacobian matrices of and at are, respectively,
It follows that the spectral radius of the next generation matrix is and given by
3.1. Existence of Equilibria and Stability
By straightforward computation, system Equation (3) has the disease-free equilibrium and a unique endemic equilibrium, , where
In this case, to maintain the feasibility state solutions of system Equation (3).
3.1.1. Local Stability of Disease-Free Equilibrium
Theorem 2. If , then the disease-free equilibrium of system Equation (3) is locally asymptotically stable and is unstable otherwise.
Proof. Linearizing system Equation (3) about , we get the characteristic equation at . It follows that and . This returns out that the disease-free equilibrium, , is locally asymptotically stable.
3.1.2. Global Stability of Disease-Free Equilibrium
Theorem 3. The disease-free equilibrium, , of system Equation (3) is globally asymptotically stable whenever and unstable otherwise.
Proof. We define the following Lyapunov function to investigate global stability of of system Equation (3). The time derivatives of Lyapunov is defined by Substituting and from system Equation (3) into Equation (13) gives
Since all the model parameters and variables are nonnegative, it follows from Equation (14) that for with holds if and only if or , and . Therefore, is Lyapunov function of . Thus, the state variables , as . Putting in Equation (3) shows that as . It therefore follows from the LaSalle’s invariance principle [26] that, every solution of system Equation (3) with initial conditions in approaches as . Therefore, this completes the proof.
3.2. Local Stability of the Endemic Equilibrium
Theorem 4. The endemic equilibrium, , of system Equation (3) is locally asymptotically stable whenever , and unstable if .
Proof. The following variational matrix is computed at to determine its local stability. That is, Following Equation (16), we found the first eigenvalue . One can find the other eigenvalues from the following reduced variational matrix:
The reduced variational matrix has negative eigenvalues if and [27]. Consequently, we then have
and the endemic equilibrium point is locally asymptotically stable.
3.2.1. Global Stability of Endemic Equilibrium
Theorem 5. If , the endemic equilibrium is locally asymptotically stable in the interior of .
Proof. Construct the Lyapunov function, , in -plane as At the endemic equilibrium, from system Equation (3), we have Computing the time derivative of along the solutions of system Equation (3), we obtain Using Equation (19), we obtain Let . Since the arithmetic mean is greater than or equal to the geometric mean, the function is nonnegative for all , and ensures that Moreover, the equality holds if and only if , and . Hence, the largest invariant compact set in is the singleton . By the LaSalle’s invariant principle [26], the endemic equilibrium , if it exists, is globally asymptotically stable inside .
3.3. Sensitivity Analysis of
It would be of interest to understand the relative importance of model parameters for the transmission and prevalence of communicable and noncommunicable diseases. This will contribute to identifying the critical model parameters that should be taken into account when considering an intervention strategy. In this section, we use a sensitivity analysis of the basic reproduction number, , in relation to the various model parameters to quantify the impact of each parameter on decreasing or increasing .
Following a similar approach [24, 28], the sensitivity analysis of the basic reproduction number, , with respect to the parameter , whose sensitivity is to be determined, is computed by
The normalized forward sensitivity index of with respect to the parameters is given by
Sensitivity indices of the basic reproduction number, , are illustrated in Table 1.
One can observe that the sensitivity indices , and are all negative, while the other sensitivity indices and remain positive. The reproduction number, , is the most sensitive to the parameters , and . The implication of this is that an increase (decrease) in the parameters and causes a corresponding increase (decrease) in . Moreover, an increase (decrease) in the parameter causes to a corresponding decrease (increase) in . On the other hand, our sensitivity analysis shows that is less sensitive to and as presented in Table 1.
4. Model with Optimal Control
In this section, we will formulate an optimal control problem and analyze it using four time variable control measures, , , , and , to reduce the disease burden and associated costs. The control represents vaccination of the susceptible population to prevent new infection in society, while the controls and represent efforts to provide appropriate treatments to the infected groups, and , respectively, to reduce the death rate of individuals due to the disease and the associated treatment costs. And the control represents the public health educational rate. That is, educating the community about the importance of handwashing, social distancing, using face masks, staying at home, etc., through different media outlets, to suppress the spread of COVID-19.
In system Equation (1), the vaccination rate and the treatment rates and were considered as constant parameters. However, in this section, we consider them as time variable and control variables , and , respectively, as stated above. This is a similar approach in [29, 30]. By imposing the time-dependent control variables in system Equation (1), we formulate the optimal control problem and becomes
We introduce the set of all admissible measurable controls defined by
The time is the terminal time. For disease control models, it is important to find the control that minimizes the prevalence and cost of controlling the disease. More specifically, we define the objective functional: where solves the system Equation (24) for the specified control . Moreover, the constants and , respectively, represent the positive weights of the infected individuals and , while the constant denotes the weight of costs. Further, represents the cost of vaccination for the susceptible humans, denotes the cost of treatment for the infectious human individuals , represents the cost of infected people under intensive care unit , and denotes the cost of preventive measures.
The aim of the optimal control problem is determining the control , such that subject to Equation (24). If such a control exists, it is, clearly, known as an optimal control. The optimal control, , associated to the corresponding solution, , gives the optimal control pair .
To achieve this, we need to check the existence of an optimal control for the system Equation (24) and the optimality system.
4.1. Existence of the Optimal Control
Theorem 6. Consider subject to system Equation (24) with , then there exists an optimal control and corresponding that minimizes .
We need to check the following assumptions based on [31] to prove the existence of an optimal control. (i)The set of controls and corresponding state variables are nonempty(ii)The control set is convex and closed(iii)The right hand side of the system Equation (24) is bounded by a linear function in the state and control
Integrand of the system Equation (26) is convex on and is bounded below by , where and
Proof. (i)The control variable is nonempty set of measurable function on the finite interval time in the real numbers(ii)Clearly, the control set is convex and closed by definition(iii)The right hand side of system Equation (24) is, by definition, continuous on the finite time interval , and can be written as a linear function of the control with coefficients depending on time and state. In addition, all the state and control variables are bounded on To verify this condition, let , and It follows that where and Thus, the condition is justified.
4.2. Characterization of the Optimal Controls
The characterization of the optimal control is based on the Pontryagin’s maximum principle [32]. This principle converts the optimal control problem into the problem of minimizing pointwise a Hamiltonian with respect to . Then, the Hamiltonian () is given, for all , by
Following the approach [30], if the control and corresponding state are the optimal control pair, then necessarily there exists a nonzero adjoint vector function fulfilling the following equality.
Finally, it follows that
Theorem 7 (see [29–31]). Consider the optimal control and corresponding trajectory with for the system Equation (24), then it is necessary that there exist nontrivial adjoint vector function satisfying Hence, Furthermore, the compact form of the optimal control, , is given by
Proof. We find the adjoint system Equation (37) by differentiating the Hamiltonian () in Equation (31) with respect to . More precisely, We assume that the final state is free. Hence, the transversality condition holds that . Moreover, the optimal controls, , are obtained based on the optimality condition
Thus, the controls are become as follows:
More precisely, Equation (42) can be written in the form of Equation (40). The second derivative of the Hamiltonian () with respect to is positive definite. That is, . Therefore, the control is then a minimizer.
5. Numerical Results and Discussions
In this case, we carried out numerical simulation of the model to support our analytical results using reasonable parameter values of the system illustrated in Table 2. In Figure 2, numerical simulation of cases without the implementation of any optimal control strategy against time can be seen. In this regard, the number of cases increases initially and then slowly drops to its endemic equilibrium. From Figure 2, it is evident that the infectious individuals, , reach its maximum value at the approximate period of time days, while the infectious population under intensive care unit, , attains its maximum value at an approximate time days from 2200 population.
As can be observed in Figures 3(a) and 3(b), numerical simulation of the infected population has been carried out with different initial values for (hence, ). The other parameters are as indicated in Table 2. These numerical results illustrate the global stability of the disease-free equilibrium of the system Equation (3). Moreover, Figures 3(c) and 3(d) present the time series of system Equation (3) for parameters illustrated in Table 2 and, hence, . From this numerical simulation, one can observe that the solution of system Equation (3) ultimately converges to its endemic equilibrium. This is an evident that our numerical results are in a good agreement with the analytical results.
(a) The infected against time for
(b) The infected against time for
(c) The infected against time for
(d) The infected against time for
5.1. The Effect of Optimal Control Strategies
In this section, we examine the results of numerical simulations with the implementation of optimal control strategies for systems Equation (24).
We employ the forward-backward sweep method to solve the optimality system. The procedure can be summarized as follows: first, we solve the control system Equation (32) with an initial guess value using the forward fourth-order Runge-Kutta method. By employing the backward fourth-order Runge-Kutta method, we therefore solved the adjoint system Equation (37) using the current iteration solution of the state variables. The controls are updated by using a convex combination of the previous controls and the values from Equation (42). The process continues till the results at the consecutive iterations are too enough.
The following are chosen as parameter values for the optimal control simulation: Moreover, the assumed weight factors are Among the (infected) and (ICU), the number of infected humans is considered the most important ones because it is assumed in the model formulation that disease transmission is due to the contact between the susceptible and the infected groups. Therefore, the highest attention is assumed to be to and relatively less () to . The cost of controlling the disease is relatively smaller but still important. For example, health education through creating awareness and vaccination is comparatively cheaper and assumed to be , whereas treatment is relatively more expensive; hence, set . These values may differ depending on different regions and nations. On the assumption that the maximum rate of intervention per day is 100%, we set All the optimal control simulations are carried out by using the initial state variables given in Table 2.
Accordingly, we study numerically to show the impact of different control strategies on the control of the spread of the disease. (i)Strategy A. Implementing all controls , and (ii)Strategy B. Implementing controls , and (iii)Strategy C. Implementing controls and (iv)Strategy D. Implementing controls and (v)Strategy A. Implementing controls and
5.1.1. Strategy A: Implementing All Controls , , , and
In this case, we use all the control strategies to optimize the objective function in Equation (26). There is a significant difference in number of cases in simulation of the system with and without controls. It can be observed from Figure 2 that the number of cases is plotted without any controls. In this case, the number of infections showed an increasing alarm and slowly adjusted to their endemic equilibrium. On the other hand, Figures 4(a) and 4(b) indicate that the number of cases and quickly decreased and lead to zero, respectively, at times and days. Specifically, the total number of infectious, and , at the end of 120 days become 743 and 1565, respectively, without controls, and zero with implementing all control strategies. This result shows the significant effect of the optimal control measures on the dynamics of the disease.
(a) with and without controls
(b) with and without controls
From the control profile Figure 5(a), the control effort maintained its maximum level till time days and gradually dropped to its lower bound in the remaining time. To minimize the burden of the disease, the optimal control measures and should be kept at their maximum level for about the first 18.5 and 19 days, respectively, while the control effort remains at its maximum rate for the whole period.
(a) Optimal control profiles when
(b) Optimal control profiles when
(c) Optimal control profiles when
(d) Optimal control profiles when
(e) Optimal control profiles when
5.1.2. Strategy B: Implementing Controls , , and
This is the case where the vaccination control , the treatment control , and the public health education control are applied to minimize the infected population while we set the treatment control to zero. One can observe from Figure 6(a) that due to the implementation of control measures, an increased number is observed for the first one day in the infectious individuals and then slowly decreased to 135 for the remaining periods, while the infected population under intensive care unit gradually declined and reaches to its lower value at time days Figure 6(b). In Figures 6(a) and 6(b), it can be observed that the number of cases is reduced when optimal controls are implemented compared to uncontrolled ones. The optimal controls for vaccination , treatment , and public health education are illustrated in Figure 5(b).
(a) with and without interventions
(b) with and without controls
5.1.3. Strategy C: Implementing Controls and
To optimize our objective function, we used the optimal control vaccination and the treatment control , while the treatment control and the public health education control are set to zero. It can be observed from Figures 7(a) and 7(b) that there is a significant difference in the number of cases with and without controls. That means, the number of infectious individuals with controls shows an initial increment and gradually drops to 208 at the time days, while the infected population under the intensive care unit revealed quickly declined at the initial time and tends to 23 at the end of the period, whereas the number of cases and without controls shows a rapid increased at the initial time and slowly decreased until they adjusted to their stable endemic equilibrium. More importantly, at the end of 120 days, the total number of cases is 2308 without controls and 231 with implementing the controls.
(a) The number of cases, , with and without controls
(b) The number of cases, , with and without controls
To reduce the burden of infections in the community, we implemented the controls and illustrated in Figure 5(c). The results in Figure 5(c) show that the vaccination control stayed in its maximum level for the first 19 days and afterwards dropped to its lower bound for the remaining periods, while the treatment rate is at its maximum rate until the period of 19.5 days and be declined to its lower bound for about the period of half a day.
5.1.4. Strategy D: Implementing Controls and
In this case, we implemented the vaccination control and the public health education control in order to minimize the number of infected people in the absence of treatment controls and . The corresponding simulation results are presented in Figures 8(a) and 8(b). From this, we notice that the infected people increases and reaches its maximum peak 716 at time of one day, while the infectious individuals under the intensive care unit attain its maximum peak 832 during the whole periods. However, the maximum peak values of the cases and without implemented controls become 2127 at an approximate time days and at an approximate time days, respectively. This shows that there is a significant difference in the number of cases with and without controls.
(a) Changes in number of cases, , with and without controls
(b) Changes in number of cases, , with and without controls
Figure 5(d) shows the control profiles of the vaccination and the public health education in which the optimal controls and are maintained at their upper bound for about the first 19.4 and the whole periods, respectively. Moreover, the control drops rapidly to its lower bound for the remaining time.
5.1.5. Strategy E: Implementing Controls and
We use the two treatment controls, and with and , to optimize the objective functional found in Equation (26). We observed from Figure 9(a) that the number of infectious individuals shows a sharp declined and then tends to zero near to the time days, while the number of infectious population with the intensive care unit leads to zero in the approximate time days as can be seen in Figure 9(b).
(a) Time series for the number of cases, , with and without controls
(b) Time series for the number of cases, , with and without controls
Figure 5(a) illustrates the controls and against time. To minimize the number of cases, the optimal control is kept at its upper bound for about 18 days and then gradually dropped to its lower bound, whereas the optimal control is kept at its upper bound for about 19.8 days and then sharply dropped to its lower bound. The treatment controls and in this strategy have to be prolonged as compared to in strategy A before dropping to their lower bounds (see Figures 5(a) and 5(e)).
In conclusion, we observed from the optimal control strategies that the combined all controls and the combined treatment controls ( and ) show the most significant effect in eradicating the burden of COVID-19 in society relative to the other combinations. The possible implication is that the disease can be eradicated by implementing the combination of all controls and or the combination of the treatments and . However, implementation of the intervention measures and or their combinations can provide an important contribution in reducing the peak values in a number of cases, as can be seen in Figures 8(a) and 8(b).
6. Conclusion
In this study, we formulated and analyzed a mathematical model with vaccination () as a constant control that describes the dynamics of COVID-19 disease. We analyzed the model for the existence and global stability of disease-free and endemic equilibria. For global stability, we used Lyapunov function methods. In particular, the global stability analysis of the endemic equilibrium is done by constructing a logarithmic Lyapunov function in the SI plane, which is similar to the works of Beretta et al., Korobeinikov et al., and Vargas-De Leon [36–38]. The system formulated is further extended into an optimal control problem by introducing continuous controls: vaccination, treatments, and public health education into the system to minimize the burden of the disease in the community and the associated costs. We verified analytically the existence of the optimal controls and characterized them by employing Pontryagin’s maximum principle. We also numerically conducted comparative studies by designing different combinations of control strategies to propose better intervention measures in the reduction and eradication of the disease.
The following can be highlighted from the analysis in this study: (i)It is found from the mathematical analysis that both the disease-free and endemic equilibria are globally asymptotically stable by constructing suitable Lyapunov functions(ii)We found also from the local sensitivity analysis that the most sensitive parameters are the recruitment rate , the transmission rate , and the vaccination rate . It is moreover found out that the reproduction number is independent of the parameter recovery rate of the ICU and the death rate of the ICU. Thus, the spread of the disease does not depend on the increase or decrease of these parameters
This work has limitations. The model formulated is deterministic and based on homogeneous mixing with the assumption that all susceptible individuals have the same chance of being infected with the disease if they come into contact with the infected ones [2, 5, 7, 30, 39]. Another limitation in this study is the assumption that the ICU is not the source of infection. But, it is vital to consider an epidemic model with the ICU as the source of infection in addition to the outpatients. It is also important that the model proposed be validated with real data. It is assumed in all the above cases that the disease emerges in a large population. However, if the population is small, susceptible individuals can possibly contact the infectious at random, and reporting the number of contacts as a function of the population density is difficult [40]. In this case, a stochastic model is preferred. A mathematical model for COVID-19 disease that fills these gaps will be considered part of our future work.
Abbreviations
| : | Susceptible human individuals at time |
| : | Infected human individuals at time |
| : | Infected human individuals under ICU at time |
| : | Recovered human individuals at a time |
| : | The total size of the human population at time |
| : | Recruitment rate of the population due to birth/immigration |
| : | The rate of transmission from susceptible to infectious |
| : | Vaccination rate of individuals |
| : | Rate of transmission from to |
| : | Rate of recovery of the infected population |
| : | Rate of recovery of the ICU |
| : | The rate of death rate of the ICU |
| : | Natural death rate. |
Data Availability
All supporting data are included within the article.
Conflicts of Interest
The authors declare that they have no any conflict of interest.
Authors’ Contributions
All authors contributed equally to this work.