A Mathematical Modelling and Analysis of COVID-19 Transmission Dynamics with Optimal Control Strategy

,


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 phys-ical 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][3][4][5][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][8][9] assumed a deterministic autonomous dynamical system and others considered fractional-order deterministic system [10][11][12][13].
There are also more studies that have been carried out on the COVID-19 dynamics system using optimal control strategies [13][14][15][16][17][18][19]. In [13], the authors proposed a mathematical model involving lock-down, quarantine, and selfisolating 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][21][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 diseasefree 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.

Model Formulation
The human population is divided into four subclasses. These subclasses of human populations are as follows: susceptible population ðSÞ, infected population ðIÞ, infected individuals under intensive care unit (ICU) (I u ), and recovered population ðRÞ. The total population N at a time t is given by Nð tÞ = SðtÞ + IðtÞ + I u ðtÞ + RðtÞ. The susceptible population S 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 I is generated proportional to the susceptible and infected individuals at a rate of β. However, the infected compartment I is decreased by the rate of progress α from I to I u , rate of recovery γ 1 , and natural death rate μ. The infected population I u is increased by the transmission rate α but declined by the recovery rate γ 2 and natural death rate μ. Moreover, the recovery compartment R gains population from vaccination of susceptible population S at the rate ν, recovery of infected individual I at the rate γ 1 , and recovery of infected individual I u at the rate γ 2 . The class R 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 I is the source of infection, but infected individuals under intensive care unit (IUC) I u 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 R 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: Computational and Mathematical Methods in Medicine Descriptions of the state variables and parameters are illustrated in Abbreviations, and the schematic diagram for the system is illustrated in Figure 1.

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 t. Theorem 1. The solution S, I, I u , and R of system Equation (1) with initial conditions S 0 > 0, I 0 ≥ 0, I u0 ≥ 0, and R 0 ≥ 0 remain positive for all t > 0. Furthermore, Proof. Using a similar approach to [23], we obtain One can observe that the above rates are all nonnegative on the bounding plane R + 4 . 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 R 0 , 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 E 0 = ðΛ/μ + ν, 0, 0Þ.
Let the state variables of system Equation (3) denoted by x ′ = ðS, I, I u Þ t . Thus, system Equation (2) can be rewritten as dx i /dt = F i ðxÞ − V i ðxÞ, where F i ðxÞ is the rate of new infection and V i ðxÞ is the rate of transfer infection compartment into and out of compartment i, and are, respectively, given by The Jacobian matrices of F i ðxÞ and V i ðxÞ at E 0 are, respectively, It follows that the spectral radius of the next generation matrix FV − is R 0 and given by 3.1. Existence of Equilibria and Stability. By straightforward computation, system Equation (3) has the disease-free equilibrium E 0 = ðΛ/μ + ν, 0, 0Þ and a unique endemic equilibrium, E 1 = ðS * , I * , I * u Þ, where In this case, Λβ > ðμ + νÞðγ 1 + μ + αÞ to maintain the feasibility state solutions of system Equation (3).

Local Stability of Disease-Free Equilibrium
Theorem 2. If R 0 < 1, then the disease-free equilibrium E 0 of system Equation (3) is locally asymptotically stable and is unstable otherwise.

Global Stability of Disease-Free Equilibrium
Theorem 3. The disease-free equilibrium, E 0 , of system Equation (3) is globally asymptotically stable whenever R 0 ≤ 1 and unstable otherwise.
Proof. We define the following Lyapunov function to investigate global stability of E 0 of system Equation (3).
The time derivatives of Lyapunov is defined by Substituting dI/dt and dI u /dt from system Equation (3) into Equation (13) gives Since all the model parameters and variables are nonnegative, it follows from Equation (14) that _ V ≤ 0 for R 0 ≤ 1 with _ V = 0 holds if and only if I = 0 or R 0 = 1, and I u = 0. Therefore, V is Lyapunov function of Π. Thus, the state variables I, I u ⟶ 0 as t ⟶ ∞. Putting I = I u = 0 in Equation (3) shows that S ⟶ Λ/μ + ν as t ⟶ ∞. It therefore follows from the LaSalle's invariance principle [26] that, every solution of system Equation (3) with initial conditions in Π approaches E 0 as t ⟶ ∞. Therefore, this completes the proof.

Local Stability of the Endemic Equilibrium
Theorem 4. The endemic equilibrium, E 1 , of system Equation (3) is locally asymptotically stable whenever R 0 > 1, and unstable if R 0 < 1.
Proof. The following variational matrix is computed at E 1 to determine its local stability. That is, Following Equation (16), we found the first eigenvalue λ 1 = −ðγ 2 + δ + μÞ < 0. One can find the other eigenvalues from the following reduced variational matrix: The reduced variational matrix has negative eigenvalues if traceðMÞ < 0 and det ðMÞ > 0 [27]. Consequently, we then have and the endemic equilibrium point E 1 is locally asymptotically stable.

Global Stability of Endemic Equilibrium
Theorem 5. If R 0 > 1, the endemic equilibrium E 1 is locally asymptotically stable in the interior of Π.
Proof. Construct the Lyapunov function, V, in SI-plane as At the endemic equilibrium, from system Equation (3), we have Computing the time derivative of VðS, IÞ along the solutions of system Equation (3), we obtain Computational and Mathematical Methods in Medicine Using Equation (19), we obtain Let x = S/S * . Since the arithmetic mean is greater than or equal to the geometric mean, the function x − 2 + 1/x is nonnegative for all x ≥ 0, and S, I, I u > 0 ensures that _ V ≤ 0: Moreover, the equality _ V = 0 holds if and only if S = S * , and I = I * . Hence, the largest invariant compact set in fðS, I, I u Þ ∈ Π : _ V = 0g is the singleton fE 1 g. By the LaSalle's invariant principle [26], the endemic equilibrium E 1 , if it exists, is globally asymptotically stable inside Π.
3.3. Sensitivity Analysis of R 0 . 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, R 0 , in relation to the various model parameters to quantify the impact of each parameter on decreasing or increasing R 0 .
Following a similar approach [24,28], the sensitivity analysis of the basic reproduction number, R 0 , with respect to the parameter Φ, whose sensitivity is to be determined, is computed by The normalized forward sensitivity index of R 0 with respect to the parameters is given by Sensitivity indices of the basic reproduction number, R 0 , are illustrated in Table 1.
One can observe that the sensitivity indices μ, α, and γ 1 are all negative, while the other sensitivity indices Λ and β remain positive. The reproduction number, R 0 , 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 R 0 . Moreover, an increase (decrease) in the parameter α causes to a corresponding decrease (increase) in R 0 . On the other hand, our sensitivity analysis shows that R 0 is less sensitive to γ 1 and μ as presented in Table 1.

Model with Optimal Control
In this section, we will formulate an optimal control problem and analyze it using four time variable control measures, u 1 ðtÞ, u 2 ðtÞ, u 3 ðtÞ, and u 4 ðtÞ, to reduce the disease burden and associated costs. The control u 1 ðtÞ represents vaccination of the susceptible population to prevent new infection in society, while the controls u 2 ðtÞ and u 3 ðtÞ represent efforts to provide appropriate treatments to the infected groups, I and I u , respectively, to reduce the death rate of individuals due to the disease and the associated treatment costs. And the control u 4 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 γ 1 and γ 2 were considered as constant parameters. However, in this section, we consider them as time variable and control variables u 1 ðtÞ, u 2 ðtÞ, and u 3 ðtÞ, 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 Table 1: Sensitivity indices of R 0 calculated at the parameter values presented in Table 2.
We introduce the set of all admissible measurable controls defined by The time T 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 xðtÞ = ðS, I, I u , RÞ solves the system Equation (24) for the specified control uðtÞ. Moreover, the constants A and B, respectively, represent the positive weights of the infected individuals I and I u , while the constant a i > 0 denotes the weight of costs. Further, a 1 u 2 1 /2 represents the cost of vaccination for the susceptible humans, a 1 u 2 2 /2 denotes the cost of treatment for the infectious human individuals I, a 3 u 2 3 /2 represents the cost of infected people under intensive care unit I u , and a 4 u 2 4 /2 denotes the cost of preventive measures. The aim of the optimal control problem is determining the control u * , such that subject to Equation (24). If such a control u * exists, it is, clearly, known as an optimal control. The optimal control, u * , associated to the corresponding solution, x * = ðS * , I * , I * u , R * Þ, gives the optimal control pair ðx * , u * Þ.
To achieve this, we need to check the existence of an optimal control for the system Equation (24) and the optimality system.
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 U 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 U and is bounded below by C 1 ðju 1 j 2 + ju 2 j 2 + ju 3 j 2 + ju 4 j 2 Þ μ/2 − C 2 , where C 1 > 0, C 2 > 0 and μ > 1: Proof.
It follows that Computational and Mathematical Methods in Medicine where C 1 = Minða 1 , a 2 , a 3 , a 4 Þ, μ = 2 and C 2 > 0: Thus, the condition is justified.

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 u. Then, the Hamiltonian (H) is given, for all t ∈ ½0, T, by Following the approach [30], if the control u * and corresponding state x * are the optimal control pair, then necessarily there exists a nonzero adjoint vector function λðtÞ fulfilling the following equality.

Computational and Mathematical Methods in Medicine
More precisely, Equation (42) can be written in the form of Equation (40). The second derivative of the Hamiltonian (H) with respect to u is positive definite. That is, ∂ 2 Hðx, u, λÞ/∂u 2 i > 0. Therefore, the control u * is then a minimizer.

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, I, reach its maximum value 2127 at the approximate period of time t = 5:28 days, while the infectious population under intensive care unit, I u , attains its maximum value 1641 at an approximate time t = 15 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 β = 0:0000218&ν = 0:1 (hence, R 0 = 0:1530 < 1). 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) Table 2 and, hence, R 0 = 1:5296 > 1. 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.

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 fourthorder 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: Λ = 100, β = 0:000218, μ = 4:7788 × 10 −5 , α = 0:01, δ = 0:0212: Moreover, the assumed weight factors are A = 10, B = 5, a 1 = 5, a 2 = 10, a 3 = 10, a 4 = 5: Among the I (infected) and I u (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 ðA = 10Þ is assumed to be to I and relatively less (B = 5) to I u . 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 a 1 = a 4 = 5, whereas treatment is relatively more expensive; hence, set a 2 = a 3 = 10. 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 u 1 max = u 2 max = u 3 max = u 4 max = 1: All the optimal control simulations are carried out by using the initial state variables given in Table 2.   Computational and Mathematical Methods in Medicine 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 u 1 , u 2 , u 3 , and u 4 (ii) Strategy B. Implementing controls u 1 , u 3 , and u 4 (iii) Strategy C. Implementing controls u 1 , and u 3 (iv) Strategy D. Implementing controls u 1 and u 4 (v) Strategy A. Implementing controls u 2 and u 3 5.1.1. Strategy A: Implementing All Controls u 1 , u 2 , u 3 , and u 4 . 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 I and I u quickly decreased and lead to zero, respectively, at times t = 6:5 and t = 8 days. Specifically, the total number of infectious, I and I u , 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. From the control profile Figure 5(a), the control effort u 3 maintained its maximum level till time t = 16:7 days and gradually dropped to its lower bound in the remaining time.
To minimize the burden of the disease, the optimal control measures u 1 and u 2 should be kept at their maximum level for about the first 18.5 and 19 days, respectively, while the control effort u 4 remains at its maximum rate for the whole period. Controls u 1 , u 3 , and u 4 . This is the case where the vaccination control u 1 , the treatment control u 3 , and the public health education control u 4 are applied to minimize the infected population while we set the treatment control u 2 to zero. One can observe from Figure 6(a) that due to the implementation of control  Computational and Mathematical Methods in Medicine measures, an increased number is observed for the first one day in the infectious individuals I and then slowly decreased to 135 for the remaining periods, while the infected population under intensive care unit I u gradually declined and reaches to its lower value 16 at time t = 20 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 u 1 , treatment u 3 , and public health education u 4 are illustrated in Figure 5(b). 5.1.3. Strategy C: Implementing Controls u 1 and u 3 . To optimize our objective function, we used the optimal control vaccination u 1 and the treatment control u 3 , while the treatment control u 2 and the public health education control u 4 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 I with controls shows an initial increment and gradually drops to 208 at the time t = 20 days, while the infected population under the intensive care unit I u revealed quickly declined at the initial time and tends to 23 at the end of the period, whereas the number of cases I and I u 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 I + I u is 2308 without controls and 231 with implementing the controls.

Strategy B: Implementing
To reduce the burden of infections in the community, we implemented the controls u 1 and u 3 illustrated in Figure 5(c). The results in Figure 5(c) show that the vaccination control u 1 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 u 3 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. . From this, we notice that the infected people I increases and reaches its maximum peak 716 at time of one day, while the infectious individuals under the intensive care unit I u attain its maximum peak 832 during the whole periods. However, the maximum peak values of the cases I and I u without implemented controls become 2127 at an approximate time t = 5:4 days and 1641 at an approximate time t = 15 days, respectively. This shows that there is a significant difference in the number of cases with and without controls. Figure 5(d) shows the control profiles of the vaccination u 1 and the public health education u 4 in which the optimal controls u 1 and u 4 are maintained at their upper bound for about the first 19.4 and the whole periods, respectively. Moreover, the control u 1 drops rapidly to its lower bound for the remaining time.
5.1.5. Strategy E: Implementing Controls u 2 and u 3 . We use the two treatment controls, u 2 and u 3 with u 1 = 0 and u 4 = 0, to optimize the objective functional J found in Equation (26). We observed from Figure 9(a) that the number of infectious individuals I shows a sharp declined and then tends to zero near to the time t = 16 days, while the number of infectious population with the intensive care unit I u leads to zero in the approximate time t = 20 days as can be seen in Figure 9(b). Figure 5(a) illustrates the controls u 2 and u 3 against time. To minimize the number of cases, the optimal control u 2 is kept at its upper bound for about 18 days and then gradually dropped to its lower bound, whereas the optimal control u 3 is kept at its upper bound for about 19.8 days and then sharply dropped to its lower bound. The treatment controls u 2 and u 3 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 (u 2 and u 3 ) show the most significant

Conclusion
In this study, we formulated and analyzed a mathematical model with vaccination (v) as a constant control that describes the dynamics of COVID-19 disease. We analyzed the model for the existence and global stability of diseasefree 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][37][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 R 0 is independent of the parameter recovery rate of the ICU γ 2 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.

13
Computational and Mathematical Methods in Medicine α: Rate of transmission from I to I u γ 1 : Rate of recovery of the infected population I γ 2 : 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.