Modeling the pandemic trend of 2019 Coronavirus with optimal control analysis.

In this work, we propose a mathematical model to analyze the outbreak of the Coronavirus disease (COVID-19). The proposed model portrays the multiple transmission pathways in the infection dynamics and stresses the role of the environmental reservoir in the transmission of the disease. The basic reproduction number R 0 is calculated from the model to assess the transmissibility of the COVID-19. We discuss sensitivity analysis to clarify the importance of epidemic parameters. The stability theory is used to discuss the local as well as the global properties of the proposed model. The problem is formulated as an optimal control one to minimize the number of infected people and keep the intervention cost as low as possible. Medical mask, isolation, treatment, detergent spray will be involved in the model as time dependent control variables. Finally, we present and discuss results by using numerical simulations.

Introduction

The 2019 novel Coronavirus (COVID-19) is a new virus that causes respiratory illness in people. This virus probably, originally, emerged from an animal source but now seems to be spreading from person-to-person. On 12th December 2019, a pneumonia case of unknown etiology was reported in Wuhan City, Hubei Province, China, and on 31st December 2019, the disease outbreak was reported to World Health Organization (WHO) [1], [2], [3], [4].

Cases have been detected in most countries worldwide and community spread is being detected in a growing number of countries. On March 11, the COVID-19 outbreak was characterized as a pandemic by the WHO [5], [6].

This virus is viewed as having zoonotic origin, which is going to be the third zoonotic human Coronavirus that is emerging in the 21st century, after the severe acute respiratory syndrome Coronavirus (SARS-CoV) that emerged in 2003 and (MERS-CoV) that emerged in 2012. The symptoms of COVID-19 including cough, fever, fatigue, breathing difficulties similar to those caused by (SARS-CoV) and (MERS-CoV) infection [7], [8].

Mathematical models have played an increasingly important role in predicting the behavior of outbreaks, optimizing control strategies, understanding the immune response and so on [9], [10], [11], [12], [13], [14], [15]. At the current stage, there are many unclear aspects about this novel Coronavirus, e.g., the contribution of asymptomatic transmission, the case fatality rate, and the initial source of 2019-nCoV. We believe that modelers can certainly make a substantial contribution to understand the virus and the transmission dynamics of the disease it causes. A number of modeling studies have already been performed for the COVID-19 pandemic. Based on reported data from December 31, 2019 to January 28, 2020, Wu et al. [16] introduced a susceptible-exposed-infectious-recovered (SEIR) model to describe the transmission dynamics, and forecasted the national and global spread of the disease. Read et al. [17] reported a value of 3:1 for the basic reproductive number based on data fitting of a SEIR model. A compartmental model which incorporate the clinical progression of the disease was proposed by Tang et al. [18].

In the current study, we investigate an epidemic model of COVID-19. This model will consist of six epidemiological classes i,e susceptible population , exposed population , infected population , asymptomatic population , hospitalized population , recovered population and reservoir for COVID-19 . In every disease the role of threshold parameter is very important for the transmission potential of a disease. We find the threshold quantity by using next generation method. We discuss sensitivity analysis in order to analyze the important of every epidemic parameter in the disease transmission. We use Routh Hurwitz criteria for the local stability of the proposed model, for global stability we use Lyapunov function theory and geometrical approach. We use optimal control strategy to minimize infected people and maximize the number of recovered people in the population. Medical mask, isolation, treatment and detergent will be involved in the model as time dependent control variable t. Finally, all the theoretical results will be verified with the help of numerical simulation for easy understanding.

This article is arrange as follow: In Section 2, we present the flowchart for the transmission of COVID-19 between reservoir and people. In Section 3, we develop COVID-19 virus transmission model. The basic reproductive number along with sensitivity analysis are presented in Section 4. In Sections 5, 6 we investigate the stability of model (1). Numerical simulation of the stability results are presented in Section 7, to verify our analytical results. In Section 8, we discuss optimal control of the proposed model and its numerical simulation.

Flow chart

The reservoir for nCoronavirus (COVID-19) is denoted as W. In Fig. 1 the population is divided into six compartment: is susceptible people; is the exposed people; is infectious people; infectious but asymptotic class ; hospitalized ; remover or recovery class . b is the birth rate and is the death rate. The susceptible people will be infected through sufficient contact with W and and the transmission rate were defined by . The transmissibility of was times and that of hospitalize was q time .

Fig. 1

Flow chart for the transmission of COVID-19 between reservoir and people.

Model formulation

This section, describe the (COVID-19) virus transmission model between reservoir and people and from people to people. This model contains a composition of differential equations. The compartmental deterministic mathematical model can be represented by nonlinear system of ordinary differential equations: with initial condition:

shows transmission per unit time, shows the approximate transmissibility of hospitalized patient, is the progression at which individuals go to infectious class, is the moving rate from exposed class to infectious class , is that of transmission to asymptotic class . The rate at which infected individuals are hospitalize is and is the recovery rate beyond hospitalization. The recovery rate of hospitalized patient is , is the life time of virus reservoir.

Equilibria and basic reproductive number

Equilibria

We discuss qualitative study of the proposed model. For this we find equilibria of the model (1). In order to find the disease free equilibrium of the proposed model (1), we set the right hand side of all equations equal to zero and set , we get is given by and the endemic equilibrium point is represented by , and it occur when the disease present in the population where: are given in Box I.

Basic reproductive number

A simple but effective measure of the transmissibility of an infectious disease is given by the basic reproduction number , defined as the total number of secondary infections produced by introducing a single infective in a completely susceptible population. In general, for simple epidemic models, if is greater than unity, an epidemic will occur while if is less than unity, an outbreak will most likely not occur.

To find for our proposed model (1) we use the method of Driessche and Watmough [19], we have

is therefore the spectral radius of next generation matrix .

Thus for our proposed model (1) becomes where The basic reproduction number of the proposed model consists of four parts, which represent the four different rout of transmission one from the exposed population, second is from the infected population, third from hospitalized population and fourth from the environmental reservoir to the susceptible population. The transmission from this four population classes shape the overall disease risk of COVID-19 pandemic.

Sensitivity analysis

In this section, we present sensitivity analysis of a few parameters which are used in the proposed model (1). This will makes it easier for us to know that parameters that have highly effect on the reproductive number. For this analysis we apply the technic given in [20]. Sensitivity index of basic reproductive number , is given by where h is parameter. For our model (1) sensitivity analysis is given by:

and summarized by:

Fig. 2 show sensitivity analysis of the basic reproductive number . These indices allow us the importance of different factor involved in the disease transmission.

Fig. 2

The graphs show the variation of different parameters and its effect on the basic reproductive number.

Local stability analysis

We show the local asymptotic stability of disease free equilibrium point and endemic equilibria of the system (1) in the following theorems.

The DFE point , is locally asymptotically stable if , otherwise unstable for .

The Jacobian matrix of the system at DFE point , is given by

Hence, the characteristic equation of matrix (3) is given by if . According to Routh-Herwitz criteria, all the roots of the characteristic polynomial have negative real parts, which complete the proof [21], [22].  □

If , then the endemic equilibrium point is locally asymptotically stable, unstable for .

Linearization of the model (1) around endemic equilibrium point is given by where Using elementary row transformation we get the following matrix : where . It is clear that all of the eigenvalues , for all , of have negative real parts for which completes the proof.  □

Global stability analysis

For the disease free equilibrium of the system is stable globally asymptotically, unstable for .

We define the following Lyapunov function, and show that this function satisfy the condition of Lyapunov function that is function is positive definite and its derivative is negative definite, Here where are arbitrary constant, which are determined later by differentiating equation (7), and using the system (1) we have By choosing the positive parameter , , and after interpretation we get, where if and and , if . By Lasala inverience principle [23], [24], and .

All the condition of Lyapunov function are satisfied that function is positive definite and its derivative is negative definite. Thus the disease free equilibrium is globally asymptotically stable in .  □

If , then the endemic equilibrium point is globally asymptotically stable and unstable otherwise.

To prove the global asymptotic stability of the proposed model (1) at endemic equilibrium , we use castilo chevez method [25], [26] let us consider the subsystem of (1): Taking the Jacobean as well as the additive compound matrix of order 2 for the above system (8), which may take the form is given by: Consider the function , then , the time derivative of the function, , implies that Now and .

, which can be written as where Let be a vector in and the of given by, Now we take the Lozinski measure described by [27], , where for and , which implies that where , , and . Therefore and becomes, such that, and , which show that . Hence . Taking integral of , we get Hence model (1) is globally asymptotically stable.  □

Numerical simulation of stability results

We solved the proposed deterministic model using Runge-e-Kutta method of order four [28]. This verify our analytical results. The variable and parameter value in Table 1 were used for simulation. For the purpose of illustrations, we assumed some parameters values. The choice of parameters are taken in the way as to be biologically feasible. The time interval is taken units with initial population for susceptible people , exposed people , infected with COVID-19 , asymptomatic people , hospitalize people , and recovered people , reservoir for COVID-19 .

Table 1

Parameters and its values.

Notation	Value	Source	Parameter	Value	Source
μ_0	0.09	[2]	κ	0.022	Estimated
β_p	0.026	[5]	δ	0.0002	Estimated
β_w	0.05	Estimated	ρ	0.065	Estimated
q	0.023	Estimated	α	0.04	[18]
γ_1	0.004	Estimated	γ_a	0.014	[17]
ϵ	0.01	[16]	γ_2	0.008	Estimated
ϕ_1	0.01	[6]	ϕ_2	0.008	[1]
b	0.00181	[2]	δ	0.008	[1]

The application of Runge-e-Kutta method of order 4th on the proposed model leads to the following system:

Algorithm

Step : .

Step : for . Step 3: for , write , , , , , , .

When we run the above algorithm by using Matlab software, we get the graphs presented in Fig. 3, Fig. 4, which represent the dynamics of susceptible population (); Exposed population (); infected with COVID-19 (); asymptomatic population (); hospitalized population (); recovered population (); and reservoir (). The biological interpretation of this results show that if , then the susceptible population decreases, then become stable and shows that there will be always stable susceptible population. The dynamics of exposed, infected, asymptomatic, hospitalize, recover and reservoir for COVID-19 conclude that the number of these populations will decrease and reach to zero, which show the stability of the proposed model.

Fig. 3

The plots demonstrate the time dynamics of different compartmental population (Susceptible, Exposed, Symptomatic and Infected, Infected but Asymptomatic).

Fig. 4

The plots demonstrate the time dynamics of different compartmental population (Hospitalized, Recovered or Removed and Reservoir for COVID-19).

Optimal control strategy for COVID-19

We formulate control strategies on the basis of sensitivity analysis and dynamic of the proposed model. The maximum sensitivity index parameter is () whose value is (0.9398437, 0.601562) increase in this parameter by 10 percent would increase the threshold quantity by (9.939 and 6. 01562. Therefore to control the spread of the disease we need to minimize this parameters by taking the control variable and representing (awareness about medical mask, hand washing and isolation of infected and non infected people). Moreover the parameters decrease the threshold quantity by 10 percent by increasing this parameter, to increase this we use the control variables representing oxygen therapy, mechanical ventilation and detergent spray.

Our goal here are to reduce COVID-19 in the population through increasing the number of recovered person R(t) and decreasing the number of infectious I(t), and hospitalized H(t), environmental reservoir W(t) by applying the time dependent control variables .

. is the time dependent control variable representing the awareness about medical mask hand washing.

. is the time dependent control variable representing isolation of infected people.

. is the time dependent control variable representing oxygen therapy mechanical ventilation.

. represent the time dependent control variable for environmental reservoir i,e detergent spray.

By using this control variables our optimal control problem which is modified version of (1) become with initial condition The goal here is to show that it is possible to implement time dependent control measures while minimizing the cost of implementation of those techniques [29]. We choose the objective (cost) function by

In Eq. (16) , , , , , , , , represent weight constant. The weight constant , , , , represent relative cost of infectious person , asymptomatic person , hospitalized person and reservoir while , , , , represents the associated cost of control variables. , , , describes self care, isolation, medical treatment, and detergent spray.

Our purpose is to find an optimal control pair , , , such that dependent on system (3), we define the control set,

Existence of optimal control problem

Let us take the control system (15) along initial condition at time and reveal the presence of the control problem. Where as bounded Lebesgue measurable controls, positive initial conditions and positive bounded solutions to the state system occur [30]. To asset the optimal solution, we go back to the optimal control problem (15), (16). First we use the Lagrangian and Hamiltonian considering the optimal control problem (15), (16). We define the Lagrangian in the following equation, For the smallest value of the Lagrangian, we determine Hamiltonian as, Hence in this way presence of control problem, we consider the following consequent.

To proceed further, first we show that these control are exists. We follow the result of, which demonstrate

There exist an optimal control , to the control problem as stated in Eqs. (15) – (16) .

In order to prove the presence of an optimal control, using the result in [31]. Since the control variables and the state variables are positive values. It is also noted that the control variables set U is convex and closed by statement. Moreover the control system is bounded which state the compactness of the problem. The integrand in, is also convex with respect to control set . Which guarantee about the existence of the optimal control .  □

Optimality condition

In order to characterize an optimal solution to (15), (16). First we use the Lagrangian and Hamiltonian considering the optimal control problem (15), (16). Indeed the Lagrangian illustrate the optimal control problem is presented by the following equation, We define the associated Hamiltonian (H), therefore using the notion and then, For the smallest value of the Lagrangian, we determine Hamiltonian for the optimal control problem as, where, and .

Following the Pontryagin’n Maximum Principle [31], [32] for finding the optimal solution to the proposed control problem (15). Using as a notation for the optimal solution then, The maximality condition define the transversally condition as

Let the optimal state variables and control variables are denoted by , , , , , be optimal state for the optimal control problem (15) and (16) . Then the set of adjoint variables satisfying

the transversality conditions (Boundary conditions) is define as, More over, the controls variables , , , are obtained as:

The adjoint system (22) comes from the direct application of the Pontryagin Maximum Principle (20), while the transversal conditions are the direct consequences of . For the set of optimal functions and , we used . We solve the optimality system numerically in the subsequent section. Because it would be easy in understanding for the reader rather then analytical results. The optimality system are characterized by the control system (15), the adjoint system (22), boundary (terminal) conditions, together with the optimal control functions. Clearly, the simulation carried out justified our control strategies to minimize the infected population, asymptomatic population, hospitalize and reservoir, and to maximize the susceptible and recovered population as shown in Fig. 5, Fig. 6.  □

Fig. 5

The graphical results show the dynamics of the compartmental population susceptible, exposed, infected, with and without controls.

Fig. 6

The graphical results show the dynamics of the compartmental population hospitalized, recovered, reservoir with and without controls.

Numerical simulation of optimal control analysis

Here we solve the optimal control system (15) to see the impact of medical mask, isolation, treatment and detergent spray by using the Runge–Kutta method of order four. We use forward Runge-Kutta procedure to solve the state system (16) with initial condition in time [0, 50]. Now to solve the adjoint system (22) we use backward Runge-Kutta procedure with transversality condition and the solution of system (16). We use the following parameters for the simulation purposes:

. This parameters are chosen in such away that are biologically more feasible. Furthermore, the weight constants are assumed to be . The obtained results are presented from Fig. 5, Fig. 6.

Fig. 5, Fig. 6 shows the dynamic of susceptible, exposed individuals, infected individuals, asymptomatic individuals, hospitalize individuals, recovered individuals, and reservoir with and with out control.

Conclusion

In this paper we have developed a mathematical model to investigate the current outbreak of Corona virus disease (COVID-19). The proposed model described the multiple transmission pathways in the infection dynamics, and stressed the role of the environmental reservoir in the transmission of the disease. This model consists of six epidemiological classes i,e susceptible population , exposed population , infected population , asymptomatic population , hospitalized population , recovered population and reservoir for COVID-19 W(t). In every disease the role of threshold parameter is very important for transmission potential of diseases. We have found threshold quantity by using the next generation matrix method. We have used Routh Hurwitz criteria for the local stability of the proposed model, while for the global stability we have used the Lyapunov function theory and geometrical approach. We further, have used optimal control strategy to minimize infected people and maximize the number of recovered people in the population. Medical mask, isolation, treatment and detergent spray have been involved in the model as time dependent control variable t. Finally, all the theoretical results have been verified with the help of numerical simulation for easy understanding.

CRediT authorship contribution statement

BiBi Fatima: Conceptualization, Methodology, Investigation, Visualization, Writing - original draft, Review and editing. Gul Zaman: Conceptualization, Methodology, Investigation, Visualization, Supervision, Review and editing. Manar A. Alqudah: Methodology, Investigation, Visualization, Review, Editing. Thabet Abdeljawad: Methodology, Investigation, Visualization, Review, Editing.

Declaration of Competing Interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Parameter	Sensitivity indices	Parameter	Sensitivity indices
β_p	+	ϵ	–
β_w	+	γ_a	–
δ	–	μ_0	–
ϕ_1	+	ϕ_2	+
b	+	q	–
α	+	κ	+
γ_a	–	γ_2	–
ρ	–