Modeling and optimal control analysis of COVID-19: Case studies from Italy and Spain.

Coronavirus disease 2019 (COVID-19) is a viral disease which is declared as a pandemic by WHO. This disease is posing a global threat, and almost every country in the world is now affected by this disease. Currently, there is no vaccine for this disease, and because of this, containing COVID-19 is not an easy task. It is noticed that elderly people got severely affected by this disease specially in Europe. In the present paper, we propose and analyze a mathematical model for COVID-19 virus transmission by dividing whole population in old and young groups. We find disease-free equilibrium and the basic reproduction number (R 0). We estimate the parameter corresponding to rate of transmission and rate of detection of COVID-19 using real data from Italy and Spain by least square method. We also perform sensitivity analysis to identify the key parameters which influence the basic reproduction number and hence regulate the transmission dynamics of COVID-19. Finally, we extend our proposed model to optimal control problem to explore the best cost-effective and time-dependent control strategies that can reduce the number of infectives in a specified interval of time.

INTRODUCTION

A novel coronavirus (COVID‐19) originated from Wuhan, China. Now, it has spread to 213 countries worldwide. Up to May 15, 2020, the total number of confirmed cases were 4,525,420 with a death toll of 303,372. Italy and Spain are in the list of top five countries which are severely affected by this disease. The transmission of COVID‐19 virus is primarily through droplets of saliva or discharge from the nose of an infected individual while coughing or sneezing. At present, there is no vaccines or specific treatments for COVID‐19. Trials of vaccines are in progress in different part of the world, and it may take few months for a approved vaccine to come into the market. In this situation, the best way to prevent and retard the transmission of this disease is to be well informed about the COVID‐19 virus, the disease it causes, and modes of its transmission. The protective measures as per WHO guidelines involve washing hands or using an alcohol‐based rub frequently and not touching the face. As there are large number of asymptomatic cases, maintaining a safe distance while roaming out, i.e., physical distancing, is another effective measure to reduce the chance of getting infected by this virus.

There are several studies on COVID‐19 based on real data. , , A simple SEIR type mathematical model for the transmission of COVID‐19 is proposed and analyzed in Lin et al. where the authors estimated the basic reproduction number based on the data from Wuhan, China. A stochastic model combined with data on number of COVID‐19 cases in Wuhan, China, is studied by Kucharski et al. Here, the authors concluded that the there was more than 50% decline in the basic reproduction number (R 0) after the introduction of travel control measures. Singh and Adhikari and Sardar et al. discussed the impact of different scenario of lock‐down to study the transmission dynamics of COVID‐19 in India. As the correct estimation of asymptomatic cases are not easy, the predictions did not go well with the current situation of COVID‐19 spread in India. As India is comparatively young country, the number of deaths in India is much less compared to many developed nations. Italy and Spain are among the top five worst affected countries in the World. It is observed that most of the deaths took place among elderly people and those who had other existing health issues.

In this paper, we have constructed mathematical model for COVID‐19 taking simple mass action type incidence. Here, we formulate our model by keeping in mind that COVID‐19 behaves differently with elderly compared to young people. The remaining of this paper is organized as follows: Section 2 describes the model; Section 3 deals with the basic reproduction number; Section 4 deals with data scenario and parameter estimation; Section 5 describes optimal control problem and the simulation results of the optimal control model; and finally, Section 6 concludes the paper.

THE MODEL

The main route of transmission of COVID‐19 is human to human. Here, we formulate our mathematical model for COVID‐19 by divide the total human population N(t) into two groups, namely, group of elderly individuals and group of young individuals keeping in mind that major death reported in elderly people. Again, we divide these groups into different compartments, namely, susceptible individuals who are young S 1(t), exposed individuals who are young E 1(t), infected individuals who are young I 1(t), susceptible individuals who are old S 2(t), exposed individuals who are old E 2(t), infected individuals who are old I 2(t), home‐isolated/hospitalized individuals of both groups who are identified as COVID‐positive and under medical supervision H(t), and recovered individuals R(t). Here, we assume that I 1 and I 2 are undetected infectives and the rate of transmission due to individuals in these two groups is very high. Here, individuals in H(t) are also infectious, but as they are under medical supervision, transmission due to individuals in H class is very low. It is assumed that the rates of transmission, rates of reinfection, rates of screening/detection, rates of movement of exposed individuals to infected compartment, and disease‐related deaths in groups of elderly people and young people are different. As elderly people can have some existing health issues, so it is assumed that the rate of transmission in elderly individuals will be more compared to rate of transmission in young individuals. Here, the compartments H(t) and R(t) contain both young and elderly individuals. The schematic diagram of our proposed model is shown in Figure 1, and the mathematical model is given as follows: where β 2 > β 1. Here, β 3 and β 4 can be taken equal as these correspond to transmission of COVID‐19 from patients under medical supervision. Additionally, we assume that ν 2 > ν 1 as the rate of detection in elderly people will be more as they will fall sick faster than younger individuals. Here, all the parameters are considered positive, and its description is given in Table 1.

FIGURE 1

Schematic diagram of the model [Colour figure can be viewed at wileyonlinelibrary.com]

TABLE 1

Description of parameters

Parameter		Description
β 1	:	Transmission rate from I 1 or I 2 to S 1,
β 2	:	Transmission rate from I 1 or I 2 to S 2,
β 3	:	Transmission rate from H to S 1,
β 4	:	Transmission rate from H to S 2
δ 1	:	Disease related death rate in I 1 compartment,
δ 2	:	Disease related death rate in I 2 compartment,
δ 3	:	Disease related death rate in H compartment,
ν 1	:	Rate of detection/isolation in I 1 compartment,
ν 2	:	Rate of detection/isolation in I 2 compartment,
η 1	:	Rate of progression of individuals from E 1 to I 1,
η 2	:	Rate of progression of individuals from E 1 to I 1,
γ 1	:	Rate of reinfection in E 1 compartment,
γ 2	:	Rate of reinfection in E 2 compartment,
α	:	Recovery rate of home isolated/hospitalized people.

THE BASIC REPRODUCTION NUMBER

We consider the system (1) and find the disease‐free equilibrium. For our model, we have disease‐free equilibrium as We find the basic reproduction number R 0 by following the next generation matrix method described in Driessche and Watmough. Following the same notations as in Driessche and Watmough, we find the vector and as follows: F = Jacobian of at

and V = Jacobian of at

and it follows that where Three eigenvalues of the above matrix are zero, and remaining two are the roots of the following quadratic equation: and

The basic reproduction number (R 0) is the largest positive root of the above quadratic and is computed as follows: where

DATA SCENARIO AND PARAMETER ESTIMATION

The total number of cases recorded in Italy as on May 10, 2020, was 219,070, and total deaths was 30,560. In Spain too, the number of cases is increasing day by day. The total number of cases recorded in Spain as on May 10, 2020, was 264,663, and total deaths was 26,621. The high rate of death from COVID‐19 in Italy and Spain may be explained by the country's relatively high proportion of elderly people. Here, we assume that 60% of total population is of young age and 40% of total population is elderly. Research has shown that the death rate is very high in elderly. At the beginning of 2020, Italy had an estimated population of 60.3 million, and at the end of the first decade of the 21st century, one in five Italians was over 65 years old. And the estimated population of Spain was 46.75 million in 2020. Keeping in view of these data, we did parameter estimation by least square method using R software. We calibrated our 2019‐nCoV model (1) to the active COVID‐19 cases for both Italy and Spain. Daily active COVID‐19 cases are collected for the period February 15, 2020, to May 10, 2020 (from https://www.worldometers.info/coronavirus/country/italy/). For the Spain, daily active COVID‐19 cases are collected for the period February 23, 2020, to May 10, 2020 (from https://www.worldometers.info/coronavirus/country/spain/).

We fit the model (1) to active cases of COVID‐19 in the Italy. We estimate the diseases transmission rates β 1,  β 2, and rate of detection of infected individuals ν 1 and ν 2. The other parameter values and the estimated values are listed in Tables 2 and 3, respectively. The observed active cases and fitted one for Italy and Spain can be seen in Figures 2 and 3, respectively. We also perform sensitivity analysis for the parameters involved in reproduction number (R 0), which reflects that increase or decrease in these parameter causes increase or decrease in (R 0). The sensitivity of R 0 to different parameters is shown in Figure 4. It is used to discover the parameters that have a high impact on R 0 and should be targeted by intervention strategies. Sensitivity indices allow to measure the relative change in a variable when parameter changes. For that, we use the forward sensitivity index of a variable, with respect to a given parameter, which is defined as the ratio of the relative change in the variable to the relative change in the parameter. If such variable is differentiable with respect to the parameter, then the sensitivity index is defined using partial derivatives. The normalized forward sensitivity index of R 0, which is differentiable with respect to a given parameter ϵ, is defined by The above formula can be used to compute the analytical expression for the sensitivity of R 0 to each parameter that it includes. From Figure 4, we can conclude that β and ν for are very sensitive parameters as small variation in these parameters can cause large variation in the value of R 0. So correct estimation of these parameters is very important to predict the transmission of this disease.

TABLE 2

Values of parameters

Parameter		value
β 3	:	0.000513	Assumed
β 4	:	0.000672	Assumed
γ 1	:	0.14	assumed
η 1	:	0.08	1–14 days 18
η 2	:	0.1	1–14 days 18
γ 2	:	0.2	Assumed
δ 1	:	0.013	Assumed
δ 2	:	0.014	Assumed
δ 3	:	0.015	0.001–0.1 19
α	:	0.071	14–28 days 20 , 21

TABLE 3

Values of parameters

Country	Estimated values	Value of R 0
Italy	β_1=0.0028
	β_2=0.0086
	ν_1=0.031	2.644
	ν_2=0.058
Spain	β_1=0.0024
	β_2=0.0085
	ν_1=0.043	2.137
	ν_2=0.053

FIGURE 2

Plots of the output of the fitted model (1) and the observed Corona active cases for Italy. Dotted line shows data points and line showed model solution [Colour figure can be viewed at wileyonlinelibrary.com]

FIGURE 3

Plots of the output of the fitted model (1) and the observed Corona active cases for Spain. Dotted line shows data points, and line showed model solution [Colour figure can be viewed at wileyonlinelibrary.com]

FIGURE 4

Forward sensitivity analysis of the parameters on R 0 (a) Italy and (b) Spain [Colour figure can be viewed at wileyonlinelibrary.com]

THE OPTIMAL CONTROL MODEL

Here, the mathematical model (1) is extended to formulate optimal control problem. Generally, control policies depend upon the severity of the epidemic in the area under investigation. It is clear from the sensitivity analysis of our proposed model that the parameters related to transmission of disease, i.e., β ′ s, and screening/detection, i.e., ν ′ s, are very important and it can have great impact in reducing the infection prevalence. Keeping this in view, we incorporate optimal control in our proposed model by considering two types of control parameters, namely, u 1(t) and u 2(t). Here, the control variable u 1(t) represents the reduction in the transmission between human to human via social distancing, awareness of transmission of disease, and sanitization.

The control variable u 2(t) corresponds to the increase in testing facility which can lead to fast detection of COVID‐19 positive cases, and we add additional time‐dependent parameter a u 2(t) in the rates of detection ν for . Keeping in view of the above assumptions, the optimal control model is formulated as follows:

The optimal control problem

In this section, we study the behavior of the proposed model by using optimal control theory. The objective functional for fixed time t is given by

Here, the parameter C 1 ≥ 0, C 2 ≥ 0, C 3 ≥ 0, C 4 ≥ 0, and they represent the weight constants. Our objective is to find the control and , such that where Ω is the control set and is defined as : measurable and 0 ≤ u 1, u 2 ≤ 1} and t ∈ [0, t ].

The Lagrangian of this problem is defined as For our problem, the associated Hamiltonian is given by where λ for is the adjoint variables. Now, adjoint variables in the form of differential equation can be written as follows: Let , , , , , , , and be the optimum values of S 1, E 1, I 1, S 2, E 2, I 2, H, and R, respectively, and , , , , , , , and be the solution of the above system of differential equations.

By using Lenhart and Workman and Pontryagin et al, we state and prove the following theorem:

There exist optimal controls such that subject to system (2).

Proof: To prove this theorem, we use Pontryagin et al. Here, all the state variables and the controls are taken as positive. For this minimizing problem, the necessary convexity of the objective functional in (u 1, u 2) is satisfied. The control variable set u 1, u 2 ∈ Ω is also convex and closed by the definition. The integrand of the functional is convex on the control set Ω, and the state variables are bounded.

Since there exist optimal controls for minimizing the functional subject to Equation 4, we use Pontryagin's maximum principle to derive the necessary conditions to find the optimal solutions as follows:

If (x, u) is an optimal solution of an optimal control problem, then there exist a non‐trivial vector function satisfying the following equalities. With the help of Pontryagin's maximum principle and Theorem 5.1, we proved the following theorem:

The optimal controls , which minimizes J over the region Ω given by where

Using optimally condition, we get This implies and gives Again, upper and lower bounds for these control are 0 and 1, respectively, i.e., if u 1 < 0 and u 2 < 0, and if and ; otherwise, and . Hence, for these controls , we get optimum value of the function J. □

NUMERICAL SIMULATION

Here, we use MATLAB to simulate our optimal control. All the parameter values are kept same as described in Tables 2, 3. The weight constants for the optimal control problem are taken as . We solve the optimality system by iterative method by using forward and backward difference approximations. We consider the final time as 120 days, i.e., the time interval as [0, 120]. First, we solve the state equations by using forward difference approximation method; then, we solve the adjoint equation by using the backward difference approximation method. We explore different types of control strategies to visualize the impact of optimal control in the total number of infected human.

Strategy I: when only one type of control is used at a time. Here, we try to find which type of optimal control is more effective in reducing the infective population. So we apply each type of control one by one. We simulate our model first for Italy and then for Spain.

Italy: In Figures 5 and 6, the control profiles of different types of optimal control when they are applied alone are shown and corresponding effects on total number of infectives (I 1 + I 2), and home‐isolated/hospitalized people (H) are shown in Figures 7A, 8A and 7B, 8B, respectively.

FIGURE 5

Control Profile (u 1) when [Colour figure can be viewed at wileyonlinelibrary.com]

FIGURE 6

Control profile (u 2) when [Colour figure can be viewed at wileyonlinelibrary.com]

FIGURE 7

A, Plot of (I 1 + I 2) verses time with and without control. B, Plot of H verses time with and without control [Colour figure can be viewed at wileyonlinelibrary.com]

FIGURE 8

A, Plot of (I 1 + I 2) verses time with and without control. B, Plot of H verses time with and without control [Colour figure can be viewed at wileyonlinelibrary.com]

Spain: In Figures 9 and 10, the control profiles of different types of optimal control when they are applied alone are shown, and corresponding effects on total number of infectives (I 1 + I 2) and home‐isolated/hospitalized people (H) are shown in Figures 11A, 12A and 11B, 12B, respectively. From these figures, it is clear that the optimal control u 1(t) is little more effective compared to other type of control but we need to maintain it to 1 for a longer duration which is not easy to achieve. This is the control through social distancing, awareness about modes of transmission of disease, and sanitization.

FIGURE 9

Control profile (u 1) when [Colour figure can be viewed at wileyonlinelibrary.com]

FIGURE 10

Control profile (u 2) when [Colour figure can be viewed at wileyonlinelibrary.com]

FIGURE 11

A, Plot of (I 1 + I 2) versus time with and without control. B, Plot of H versus time with and without control [Colour figure can be viewed at wileyonlinelibrary.com]

FIGURE 12

A, Plot of (I 1 + I 2) versus time with and without control. B, Plot of H versus time with and without control [Colour figure can be viewed at wileyonlinelibrary.com]

Strategy II: When both the controls are used. Here, both the control mechanisms (u 1, u 2) are used to optimize the objective functional J.

Italy: The variation of total infected human and home‐isolated/hospitalized people (H) with time are shown in Figure 13A,B. Here, it is observed that there is a reasonable decrease in the total number of infectives when both controls are used simultaneously. Figure 14A,B shows the control profiles of u 1 and u 2, respectively.

FIGURE 13

A, Plot of (I 1 + I 2) versus time with and without control. B, Plot of H versus time with and without control [Colour figure can be viewed at wileyonlinelibrary.com]

FIGURE 14

A, Control Profile (u 1). B, Control profile (u 2) [Colour figure can be viewed at wileyonlinelibrary.com]

Spain: The variation of total infected human and isolated/hospitalized people (H) with time are shown in Figure 15A,B. Here, it is easy to observe that there is a reasonable decrease in the total number of infectives when both controls are used simultaneously. Figure 16A,B shows the control profiles of u 1 and u 2, respectively.

FIGURE 15

A, Plot of (I 1 + I 2) versus time with and without control. B, Plot of H versus time with and without control [Colour figure can be viewed at wileyonlinelibrary.com]

FIGURE 16

A, Control profile (u 1). B, Control profile (u 2) [Colour figure can be viewed at wileyonlinelibrary.com]

The simulation result demonstrates the effectiveness of optimal control strategies in reducing the number of infectives. It is observed that combined controls are more useful in reducing the number of infected cases significantly.

CONCLUSION

Here, a mathematical model for COVID‐19 virus disease is formulated and analyzed.We computed disease‐free equilibrium and basic reproduction number R 0. We estimated the key parameters using least square estimation method using real life data fitted with mathematical model for Italy and Spain. Sensitivity analysis is performed to find the key parameters that are very sensitive to basic reproduction number R 0. Further, the proposed model is extended to optimal control problem by incorporating two types of controls. Then, Pontryagin's maximum principle is used to analyze the optimal control problem. The numerical simulation is explored by considering different combinations of optimal controls. Simulation results indicate that optimal control strategy is in fact effective in reducing the total number of infectives if both the controls are applied simultaneously.