© Revista de Matemática: Teoría y Aplicaciones 2021 28(1) : 55–78 cimpa – ucr issn: 1409-2433 (Print), 2215-3373 (Online)

doi: https:// doi. org/ 10.15517/ rmta. v28i1. 42077

optimizing the quarantine cost for

suppression of the covid-19 epidemic in mexico

optimización del costo de la cuarentena para la supresión de la

epidemia del covid-19 en méxico

Abdon E. Choque-Rivero^{∗}  Evgenii N. Khailov^{†}

Ellina V. Grigorieva^{‡}

Received: 07/Aug/2020; Revised: 12/Nov/2020; Accepted: 12/Nov/2020

∗Universidad Michoacana de San Nicolás de Hidalgo, Institute of Physics and Mathematics, Morelia, México. E-Mail: abdon.choque@umich.mx

†Lomonosov Moscow State University, Faculty of Computational Mathematics and Cybernetics, Moscow, Russia. E-Mail: khailov@cs.msu.su

‡Texas Woman’s University, Department of Mathematics and Computer Sciences, Denton, United States. E-Mail: egrigorieva@twu.edu

56 a. choque-rivero — e. khailov — e. grigorieva

Abstract

This paper is one of the few attempts to use the optimal control theory to find optimal quarantine strategies for eradication of the spread of the COVID-19 infection in the Mexican human population. This is achieved by introducing into the SEIR model a bounded control function of time that reflects these quarantine measures. The objective function to be minimized is the weighted sum of the total infection level in the population and the total cost of the quarantine. An optimal control problem reflecting the search for an effective quarantine strategy is stated and solved analytically and numerically. The properties of the corresponding optimal control are established analytically by applying the Pontryagin maximum principle. The optimal solution is obtained numerically by solving the two-point boundary value problem for the maximum principle using MATLAB software. A detailed discussion of the results and the corresponding practical conclusions are presented.

Keywords: coronavirus; quarantine cost; Pontryagin maximum principal; optimal control.

Resumen

En este trabajo empleamos la teoría de control óptimo para encontrar una cuarentena óptima y estrategias para la erradicación de la propagación de la infección por COVID-19 en la población humana mexicana. En un modelo SEIR, introducimos un control acotado que es una función respecto del tiempo, la cual refleja las medidas de la cuarentena. La función objetivo a minimizar es la suma ponderada del nivel total de infección en la población y el costo total de la cuarentena. Planteamos un problema de control óptimo que representa la búsqueda de una estrategia eficaz de una cuarentena. Resolvemos este problema analíticamente y numéricamente. Establecemos analíticamente las propiedades del control óptimo correspondiente aplicando el principio del máximo de Pontryagin. La solución óptima se obtiene resolviendo un problema de valor de frontera de dos puntos asociado al principio del máximo. Usamos el software MATLAB. Presentamos una discusión detallada de los resultados y las correspondientes conclusiones prácticas.

Palabras clave: coronavirus; costo de una cuarentena; principio del máximo de Pontryagin; control óptimo.

Mathematics Subject Classification: 49N90, 58E25, 92D25, 92D30.

1 Introduction

The coronavirus (COVID-19) pandemic is the number one topic worldwide. This pandemic has literally affected all countries, and their authorities are trying to take all possible (including emergency) measures to contain and fight this virus: state border closures, quarantine, self-isolation, the termination of work for many businesses and institutions, as well as the transition working and training at home. The regions most affected by the pandemic are the United States, Brazil, Europe (Russia, Italy, Spain, Germany, France, the United Kingdom), China, the Republic of Korea, and Iran ([32]).

In the fight against COVID-19, different countries are helped by various mathematical models that can predict possible variants in the development of the pandemic and the onset of its peak and duration, depending on the initial data and measures to contain the disease. One of the most famous mathematical tools for predicting the development of epidemics and taking appropriate measures for combating them is compartment models: SIR (Susceptible-Infectious-Recovered) and SEIR (Susceptible-Exposed-Infectious-Recovered) ([2, 29]).

In contrast to previous epidemics associated with the spread of coronaviruses (for example, SARS in 2003 to 2004 and MERS in 2012), the COVID-19 pandemic has two distinctive features:

• a long (more than 14 days) incubation period and

• a large number of asymptomatic patients who, having contacted coronavirus, do not demonstrate clinical manifestations, although they could infect others.

These features make it necessary to modify the existing models. Currently, there are many studies published as preprints and articles in which SIR and SEIRtype models are proposed that consider the features of COVID-19. A detailed analysis of their properties is provided, specific parameters of such models are determined to describe the spread of the coronavirus infection in a particular country or group of countries ([1, 4, 16, 21, 22, 30, 33]).

In addition, for each such model, the reproductive ratio R_{0} is calculated, which is the basic characteristic of the coronavirus pandemic. For the COVID-19 pandemic, this ratio takes the value between 2.0 and 6.0 ([3, 4, 17, 21, 24, 34, 35]). For any epidemic, it is known ([2, 29]) that with R_{0} < 1, the epidemic gradually leaves, and with R_{0} > 1 the disease accompanying this epidemic expands exponentially.

In the framework of such SIR and SEIR-type models, attempts are made to find the effective values of their parameters, reflecting, for example, quarantine and self-isolation. At these values the current values of R_{0} decreases. In [26] the Bayesian approach is employed to estimate model parameters and the value R_{0} related with the study of the transmission of the COVID-19 in Mexico.

Usually mathematical modeling helps medical doctors and epidemiologists to answer the following questions:

• How fast will the virus spread?

• How many people will get infected by COVID-19?

• If there are no vaccine and a reliable treatment, then how many will remain uninfected?

• How long the pandemic will last and will there be the second wave of it?

On one hand, most people follow the rules, wear masks and stay at their homes. On the other hand, especially recently, from TV screens and internet we hear stories about huge unemployment and the damage that the virus caused to the economies in all the countries, the government of which imposed “stay at homes” order. We see people going on strikes demanding to reopen businesses, go back to work and to a normal life style, even we all know that the virus has not been eradicated yet. As an example, we hear about Sweden, a small (by the size of the population) European country, that did not impose quarantine restrictions at all and relied on so-called herd immunity of their citizens and their own selfdiscipline and wise decision. While herd immunity is something to think about, with high R_{0} = 4.0 that was reported in some regions of Italy, should we allow 75% of the citizen to get sick, knowing that some of them may never recover from the disease? Is that OK to play Russian roulette with innocent people?

So, is there any optimal strategy which would simultaneously save the lives and won’t sink the economy? The answer to this question can be obtained by using optimal control theory, which has proven itself in the search for the optimal intervention strategies for the 2014 Ebola epidemic ([14, 15]).

This paper is one of the few attempts to use the optimal control theory to find optimal quarantine strategies for a SEIR-type model that describes the spread of the coronavirus infection in the Mexican human population. This is achieved by introducing into the model a bounded control function of time that reflects these quarantine measures. The objective function to be minimized is the weighted sum of the total infection level in the population and the total cost of this quarantine (Section 2). In Section 3, the Pontryagin maximum principle is applied and the uniqueness and other important properties of the optimal control are established analytically. In Section 4, the model parameters are calculated with respect to real COVID-19 data reported in Mexico between February 28 and May 21, 2020. These official data are divided into three time periods (phase one, two and three). For each phase, we calculated the basic reproductive ratio with and without optimal control. In this section we also found the maximum strength of the quarantine measures (the threshold value of u_{max}^{new}) that mathematically leads to R_{0} < 1 and would mean the end of the pandemic. It is found that in order to eradicate the virus, there is not need for the optimal control to take its maximum values (for example, 0.95) but its value depends on the current basic reproductive ratio in the region and changes from 0.5 − 0.83 for R_{0} = 2.0−6.0, respectively. The optimal solutions are obtained numerically using a computer program written in MATLAB and their graphs are presented for two months time period and two different basic reproductive ratios (R_{0} = 2.65 and R_{0} = 3.55) for Mexico and Mexico City, respectively. A detailed discussion of the obtained results and the corresponding practical conclusions are presented in Section 5.

2 Statement of the optimal control problem

At a given time interval [0, t_{f}], let us examine the spread of the virus in a human population of size N(t) (not counting the deceased people), which is divided into the five following compartments: susceptible people S(t), exposed people E(t), symptomatic infected people I(t), asymptomatic infected people A(t) and recovered people R(t). Hence, we imply the natural equality:

S(t) + E(t) + I(t) + A(t) + R(t) = N(t), t ∈ [0, t_{f}].  (1)

During the epidemic period, the population’s natural birth and death rate are at a relatively low level. Therefore, we can exclude them from consideration. Furthermore, we can consider the natural situation when the virus is imported from outside. In addition, we will assume that at initial moment t = 0 there are no recovered people.

Thus, the change in the size of the compartments is described by the following system of differential equations:

\begin{cases}S′(t) & = −β̃S(t)(I(t) + νA(t)), t ∈ [0, t_{f}], \\ E′(t) & = β̃S(t)(I(t) + νA(t)) − ωE(t), \\ I′(t) & = (1 − δ)ωE(t) − γI(t), \\ A′(t) & = δωE(t) − ηA(t), \\ R′(t) & = γI(t) + (1 − q)ηA(t) \\ N′(t) & = −qηA(t), \end{cases}  (2)

with the corresponding initial conditions

\begin{aligned}S(0) & = S_{0}, E(0) = E_{0}, I(0) = I_{0}, \\ A(0) & = A_{0}, R(0) = R_{0}, N(0) = N_{0}. \end{aligned}  (3)

We will assume that at the initial moment t = 0 the values S_{0}, E_{0}, I_{0}, A_{0}, N_{0} are positive and R_{0} = 0. Moreover, the equality

S_{0} + E_{0} + I_{0} + A_{0} = N_{0}.  (4)

holds, where N_{0} is the initial population size.

A system similar to (2) was used in [4] to describe the spread of the COVID-19 virus among the Wuhan (China) population.

It is easy to see that firstly the equations of system (2) together with the initial conditions (3) and equality (4) imply relationship (1), and, secondly, the value N(t) varies (the decreases caused by disease-induced mortality).

In system (2), the incubation period of the human infection is defined as ω^{−1}. The infectious periods of I(t) and A(t) are defined as γ^{−1} and η^{−1}, respectively. The proportion of asymptomatic infection is defined as δ. The S(t) will be infected through sufficient contact with I(t) and A(t), and the transmission rates are defined as β̃ and νβ̃, respectively. Here we assume that the transmissibility of A(t) is ν times that of I(t), where ν ∈ (0, 1). Moreover, we consider that the values γI(t) and (1 − q)ηA(t) are the number of recovered individuals in I(t) and A(t), respectively, and the value qηA(t) determines the number of people who die from the disease. Here the value q ∈ (0, 1) sets the proportion of deaths in ηA(t).

Let us perform the normalization of the phase variables for system (2) with the initial conditions (3) using the following formulas:

\begin{aligned}s(t) & = N_{0}^{−1}S(t), e(t) = N_{0}^{−1}E(t), i(t) = N_{0}^{−1}I(t), \\ a(t) & = N_{0}^{−1}A(t), r(t) = N_{0}^{−1}R(t), n(t) = N_{0}^{−1}N(t). \end{aligned}

Moreover, we introduce the new transmission rate β by the formula β = β̃N_{0}.

As a result, we obtain the following system of equations:

\begin{cases}s′(t) & = −βs(t)(i(t) + νa(t)), t ∈ [0, t_{f}], \\ e′(t) & = βs(t)(i(t) + νa(t)) − ωe(t), \\ i′(t) & = (1 − δ)ωe(t) − γi(t), \\ a′(t) & = δωe(t) − ηa(t), \\ r′(t) & = γi(t) + (1 − q)ηa(t) \\ n′(t) & = −qηa(t), \end{cases}  (5)

with the corresponding initial conditions

\begin{aligned}s(0) & = s_{0}, e(0) = e_{0}, i(0) = i_{0}, \\ a(0) & = a_{0}, r(0) = 0, n(0) = 1, \end{aligned}  (6)

where s_{0}, e_{0}, i_{0}, a_{0}. are positive and satisfy the equality

s_{0} + e_{0} + i_{0} + a_{0} = 1.  (7)

following from (4).

Note that with this normalization the relationship (1) is converted to equality:

s(t) + e(t) + i(t) + a(t) + r(t) = n(t), t ∈ [0, t_{f}].

The important properties of solutions for system (5) are established by the following lemma.

Lemma 2.1 Let system (5) with the initial conditions (6) have the solutions s(t), e(t), i(t), a(t), r(t), n(t). They are then defined in the entire interval [0, t_{f}], and are also positive and bounded on (0, t_{f}].

The Proof of this fact is standard, so we omit it. (For example, similar proofs are given in [6, 19]). Lemma 2.1 implies that all solutions s(t), e(t), i(t), a(t), r(t), n(t) for system (5) with the initial conditions (6) retain their biological meanings for all t ∈ [0, t_{f}].

Let us introduce the control function u(t) into system (5). This control reflects the intensity of the quarantine and all indirect protective measures (mask wearing, as well as various educational and informational campaigns), which are set up in the population to limit the spread of the virus and are aimed at the reduction of its transmission. In addition, we assume that there is no vaccine and no drug available for the disease treatment. This control function satisfies the restrictions:

0 ≤ u(t) ≤ u_{max} < 1.  (8)

This leads to the following control system:

\begin{cases}s′(t) & = −β(1 − u(t))s(t)(i(t) + νa(t)), t ∈ [0, t_{f}], \\ e′(t) & = β(1 − u(t))s(t)(i(t) + νa(t)) − ωe(t), \\ i′(t) & = (1 − δ)ωe(t) − γi(t), \\ a′(t) & = δωe(t) − ηa(t), \\ r′(t) & = γi(t) + (1 − q)ηa(t) \\ n′(t) & = −qηa(t), \end{cases}  (9)

with the corresponding initial conditions (6).

Note that for u(t) = 0 (the absence of the quarantine measures) system (9) becomes system (5) with the rate of virus transmission as β. When u(t) > 0 (the presence of the quarantine), in system (9) such a transmission rate is reduced.

In constructing the control model (9), extensive experience has been used in similar SEIR-type control models of Ebola epidemics ([11, 12, 13, 14, 15]).

Now, for formulating optimal control problem, let us introduce the set Δ(t_{f}) of all admissible controls, which is formed by all possible Lebesgue measurable functions u(t) that for almost all t ∈ [0, t_{f}] satisfy restrictions (8).

Next, for the control system (9) in the set Δ(t_{f}) of all admissible controls, we consider the objective function:

J(u(·)) = (e(t_{f}) + i(t_{f}) + a(t_{f})) + ∫_{0}^{t_{f}}(e(t) + i(t) + a(t))dt + 0.5α ∫_{0}^{t_{f}} u^{2}(t)dt,

(10)

where α is a positive weighting coefficient. The first two terms in (10) reflect the level of disease in the population caused by COVID-19: the level at the end of quarantine period [0, t_{f}] and cumulative level over the entire period. The last term determines the total cost of the quarantine.

Since only phase variables e(t), i(t), a(t) are present in the objective function (10), then by considering these functions together with system (9), the last two differential equations can be excluded from it.

As a result, we state the following optimal control problem (OC problem) consisting of minimizing the objective function (10) in the set Δ(t_{f}) of all admissible controls for the system:

\begin{cases}s′(t) & = −β(1 − u(t))s(t)(i(t) + νa(t)), t ∈ [0, t_{f}], \\ e′(t) & = β(1 − u(t))s(t)(i(t) + νa(t)) − ωe(t), \\ i′(t) & = (1 − δ)ωe(t) − γi(t), \\ a′(t) & = δωe(t) − ηa(t), \end{cases}  (11)

with the corresponding initial conditions

s(0) = s_{0}, e(0) = e_{0}, i(0) = i_{0}, a(0) = a_{0},  (12)

where s_{0}, e_{0}, i_{0}, a_{0} are positive and satisfy equality (7).

Lemma 2.1 and the fulfillment of the easily verified condition of Theorem 4 (Chapter 4, [20]) guarantee for the OC problem the existence of an appropriate optimal solution, which consists of the optimal control u_{∗}(t) and the corresponding optimal solutions s_{∗}(t), e_{∗}(t), i_{∗}(t), a_{∗}(t) to system (11).

3 Pontryagin maximum principle

For the analytical study of the OC problem, we use the Pontryagin maximum principle ([25]). According to it, we first write down the Hamiltonian of this problem:

H(s,e, i, a, ψ_{1}, ψ_{2}, ψ_{3}, ψ_{4}, u) = −β(1 − u)s(i + νa)(ψ_{1} − ψ_{2})

− ωe(ψ_{2} − (1 − δ)ψ_{3} − δψ_{4}) − γiψ_{3} − ηaψ_{4} − (e + i + a) − 0.5αu^{2},

where ψ_{1}, ψ_{2}, ψ_{3}, ψ_{4} are the adjoint variables.

Next, for this Hamiltonian we calculate the required partial derivatives:

\begin{aligned}H′_{s}(s, e, i, a, r, ψ_{1}, ψ_{2}, ψ_{3}, ψ_{4}, u) & = −β(1 − u)(i + νa)(ψ_{1} − ψ_{2}), \\ H′_{e}(s, e, i, a, r, ψ_{1}, ψ_{2}, ψ_{3}, ψ_{4}, u) & = ω(ψ_{2} − (1 − δ)ψ_{3} − δψ_{4}) − 1, \\ H′_{i}(s, e, i, a, r, ψ_{1}, ψ_{2}, ψ_{3}, ψ_{4}, u) & = −β(1 − u)s(ψ_{1} − ψ_{2}) − γψ_{3} − 1, \\ H′_{a}(s, e, i, a, r, ψ_{1}, ψ_{2}, ψ_{3}, ψ_{4}, u) & = −βν(1 − u)s(ψ_{1} − ψ_{2}) − ηψ_{4} − 1, \\ H′_{u}(s, e, i, a, r, ψ_{1}, ψ_{2}, ψ_{3}, ψ_{4}, u) & = βs(i + νa)(ψ_{1} − ψ_{2}) − αu. \end{aligned}

Therefore, by the Pontryagin maximum principle, for the optimal control u_{∗}(t) and the corresponding optimal solutions s_{∗}(t), e_{∗}(t), i_{∗}(t), a_{∗}(t) to system (11), there exists the vector-function ψ_{∗}(t) = (ψ_{1}^{∗}(t), ψ_{2}^{∗}(t), ψ_{3}^{∗}(t), ψ_{4}^{∗}(t), such that

• ψ_{∗}(t) is the nontrivial solution of the adjoint system

\begin{cases}ψ_{1}^{∗}′(t) & = −H′_{s}(s_{∗}(t), e_{∗}(t), i_{∗}(t), a_{∗}(t), r_{∗}(t), ψ_{1}^{∗}(t), ψ_{2}^{∗}(t), ψ_{3}^{∗}(t), ψ_{4}^{∗}(t), u_{∗}(t)) \\ & = β(1 − u_{∗}(t))(i_{∗}(t) + νa_{∗}(t))(ψ_{1}^{∗}(t) − ψ_{2}^{∗}(t)), \\ ψ_{2}^{∗}′(t) & = −H′_{e}(s_{∗}(t), e_{∗}(t), i_{∗}(t), a_{∗}(t), r_{∗}(t), ψ_{1}^{∗}(t), ψ_{2}^{∗}(t), ψ_{3}^{∗}(t), ψ_{4}^{∗}(t), u_{∗}(t)) \\ & = ω(ψ_{2}^{∗}(t) − (1 − δ)ψ_{3}^{∗}(t) − δψ_{4}^{∗}(t)) + 1, \\ ψ_{3}^{∗}′(t) & = −H′_{i}(s_{∗}(t), e_{∗}(t), i_{∗}(t), a_{∗}(t), r_{∗}(t), ψ_{1}^{∗}(t), ψ_{2}^{∗}(t), ψ_{3}^{∗}(t), ψ_{4}^{∗}(t), u_{∗}(t)) \\ & = β(1 − u_{∗}(t))s_{∗}(t)(ψ_{1}^{∗}(t) − ψ_{2}^{∗}(t)) + γψ_{3}^{∗}(t) + 1, \\ ψ_{4}^{∗}′(t) & = −H′_{a}(s_{∗}(t), e_{∗}(t), i_{∗}(t), a_{∗}(t), r_{∗}(t), ψ_{1}^{∗}(t), ψ_{2}^{∗}(t), ψ_{3}^{∗}(t), ψ_{4}^{∗}(t), u_{∗}(t)) \\ & = βν(1 − u_{∗}(t))s_{∗}(t)(ψ_{1}^{∗}(t) − ψ_{2}^{∗}(t)) + ηψ_{4}^{∗}(t) + 1, \end{cases}

(13)

satisfying the corresponding initial conditions

\begin{aligned}ψ_{1}^{∗}(t_{f}) & = −J′_{s(t_{f})} = 0, ψ_{2}^{∗}(t_{f}) = −J′_{e(t_{f})} = −1, \\ ψ_{3}^{∗}(t_{f}) & = −J′_{i(t_{f})} = −1, ψ_{4}^{∗}(t_{f}) = −J′_{a(t_{f})} = −1. \end{aligned}  (14)

• the control u_{∗}(t) maximizes the Hamiltonian

H(s_{∗}(t), e_{∗}(t), i_{∗}(t), a_{∗}(t), ψ_{1}^{∗}(t), ψ_{2}^{∗}(t), ψ_{3}^{∗}(t), ψ_{4}^{∗}(t), u)  (15)

with respect to u ∈ [0, u_{max}] for almost all t ∈ [0, T], and, therefore the following relationship holds:

u_{∗}(t) = \begin{cases}u_{max} & , if ϕ(t) > u_{max}, \\ ϕ(t) & , if 0 ≤ ϕ(t) ≤ u_{max}, \\ 0 & , if ϕ(t) < 0. \end{cases}  (16)

Here the function ϕ(t) is the so-called indicator function ([27]), which is defined from the formula

H′_{u}(s_{∗}(t), e_{∗}(t), i_{∗}(t), a_{∗}(t), ψ_{1}^{∗}(t), ψ_{2}^{∗}(t), ψ_{3}^{∗}(t), ψ_{4}^{∗}(t), u) = 0

and is written as

ϕ(t) = α^{−1}βs_{∗}(t)(i_{∗}(t) + νa_{∗}(t))(ψ_{1}^{∗}(t) − ψ_{2}^{∗}(t)).  (17)

It determines the behavior of the optimal control u_{∗}(t) according to formula (16).

Now by formulas (14) and (17), we find the relationship

ϕ(t_{f}) = α^{−1}βs_{∗}(t_{f})(i_{∗}(t_{f}) + νa_{∗}(t_{f})),

which because of Lemma 2.1 implies the inequality ϕ(t_{f}) > 0. According to formula (16), this means that the following lemma is true.

Lemma 3.1 The optimal control u_{∗}(t) is positive for t = t_{f} and takes either the value of ϕ(t_{f}) or the value of u_{max}.

Thus, based on this lemma, we can state that for any finite time interval, the optimal protective measures are such that even at the very end the optimal control is not zero. Hence, in absence of vaccine or effective medical treatment for COVID-19, some protective measures (i.e. wearing masks and gloves, escaping crowd) must remain in order or become people habits.

Next formula (16) shows that for all values of t ∈ [0, t_{f}], the maximum of the Hamiltonian (15) is reached with a unique value u = u_{∗}(t). Thus, the next lemma immediately follows from Theorem 6.1 ([5]).

Lemma 3.2 The optimal control u_{∗}(t) is a continuous function on the interval [0, t_{f}].

We can rewrite formula (16) in a more suitable form:

u_{∗}(t) = min{u_{max}; max{0; ϕ(t)}}.  (18)

As a result, systems (11) and (13) with the corresponding initial conditions (12) and (14) together with relationship (18) form the two-point boundary value problem for the maximum principle:

s′(t) = −β(1 − u(t))s(t)(i(t) + νa(t)),
e′(t) = β(1 − u(t))s(t)(i(t) + νa(t)) − ωe(t),
i′(t) = (1 − δ)ωe(t) − γi(t),
a′(t) = δωe(t) − ηa(t),
ψ_{1}′(t) = β(1 − u(t))(i(t) + νa(t))(ψ_{1}(t) − ψ_{2}(t)),
ψ_{2}′(t) = ω(ψ_{2}(t) − (1 − δ)ψ_{3}(t) − δψ_{4}(t)) + 1,
ψ_{3}′(t) = β(1 − u(t))s(t)(ψ_{1}(t) − ψ_{2}(t)) + γψ_{3}(t) + 1,
ψ_{4}′(t) = βν(1 − u(t))s(t)(ψ_{1}(t) − ψ_{2}(t)) + ηψ_{4}(t) + 1,
s(0) = s_{0}, e(0) = e_{0}, i(0) = i_{0}, a(0) = a_{0},
ψ_{1}(t_{f}) = 0, ψ_{2}(t_{f}) = −1, ψ_{3}(t_{f}) = −1, ψ_{4}(t_{f}) = −1,
u(t) = min {u_{max}; max {0; α^{−1}βs(t)(i(t) + νa(t))(ψ_{1}(t) − ψ_{2}(t))}} .

(19)

The optimal control u_{∗}(t) satisfies this boundary value problem together with the corresponding optimal solutions s_{∗}(t), e_{∗}(t), i_{∗}(t), a_{∗}(t) for system (11) and the solutions ψ_{1}^{∗}(t), ψ_{2}^{∗}(t), ψ_{3}^{∗}(t) and ψ_{4}^{∗}(t) to the adjoint system (13). Moreover, arguing as in [18, 23, 28], it is possible to establish the uniqueness of this control because of the boundedness of the state and adjoint variables, as well as, the Lipschitz properties of systems (11) and (13) defining these variables and relationships (18) that establishes such a control.

4 Numerical results and their discussion

For our numerical simulation, we employ the following values for the parameters

\begin{matrix}α = 0.0001 & γ = 1/14 & η = 1/21 & δ = 0.2 \\ ω = 0.18 & ν = 0.1 & q = 0.15. \end{matrix}  (20)

These parameters were adjusted to fit the number of infected individuals in Mexico during the period from February 28 to May 21, 2020.

In our model, February 28 is taken as the start of the pandemic. This initial time is also applied for our simulation with u = 0. For the simulation with u = u_{∗}, the initial time is taken as March 24. On this date, phase two of the government’s contingency plan was announced, suspending all non-essential activities. On April 21 the Mexican authorities implemented phase three of the pandemic in order to reduce the movement of people in public spaces, implementing social distancing.

The graphs in Figure 1 show the relation between the official reported data, marked with dashed line, and the trajectory of our simulation, marked with solid line of the infected population in Mexico. The graph on the left corresponds to phase one (from February 28 till March 24), while the graph on the right corresponds to phase three (from March 24 till May 21). We do not provide the graph for phase two (from March 24 till April 21) because in this period the official data was focused on the cumulated infected people and not on the active infected people.

We use the basic reproductive ratio R_{0} to assess the transmissibility of the virus. It is known ([29]) that this ratio shows the average number of people that one contagious person can infect during the period of the disease, that is, until complete recovery. If R_{0} < 1, then the epidemic will spread slowly, and it will fade. Finally, it will die out. If R_{0} > 1, then the epidemic will occur, and the disease will spread exponentially. What the specific value of R_{0} will be, depends on the properties of the virus, the percentage of the population that becomes immune to it (for example, by testing for antibodies to the virus), and on the measures taken by the population to suppress the epidemic (for example, social distancing and quarantine).

Figure 1: Officially reported data and simulated behavior of the normalized number of infected population during the phase one and phase three in Mexico.

We find the value of R_{0} for system (5) by the next-generation matrix approach ([31]). This value is defined by the following formula:

R_{0} = β (\frac{1 − δ}{γ} + \frac{νδ}{η} ) .  (21)

By using the required values from (20) in (21), we find the equality

R_{0} = 11.62 · β.  (22)

Table 1 shows the values of the parameter β that depends on the value of R_{0} in accordance with formula (22). We select the value of R_{0} from {2.0; 2.5; 3.0; 3.5; 4.0; 6.0} because this ratio took the value between 2.0 − 4.0 in Mexico ([7, 26]).

Table 1: Parameter β and threshold value u_{max}^{new} depending on R_{0}.

\begin{matrix}R_{0} & β & u_{max}^{new} \\ 2.0 & 2.0/11.62 = 0.172117 & 0.500000 \\ 2.5 & 2.5/11.62 = 0.215146 & 0.600000 \\ 3.0 & 3.0/11.62 = 0.258176 & 0.666667 \\ 3.5 & 3.5/11.62 = 0.318182 & 0.714286 \\ 4.0 & 4.0/11.62 = 0.363636 & 0.750000 \\ 6.0 & 6.0/11.62 = 0.545455 & 0.833333 \end{matrix}

Rev.Mate.Teor.Aplic. (ISSN print: 1409-2433; online: 2215-3373) Vol. 28(1): 55–78, Jan–Jun 2021

Now let us obtain the appropriate basic reproductive ratios R_{0}(u) for the control system (9), or, what is the same, for system (11) under the assumption of constancy of the control. By using (21), we find the formula:

R_{0}(u) = β(1 − u) (\frac{1 − δ}{γ} + \frac{νδ}{η} ) = (1 − u)R_{0}.  (23)

Next, we substitute the required values from (20) into the expression (23) and then perform the necessary calculations for u = u_{max} because we want to see that with the maximum intensity of quarantine the epidemic will come to an end. As a result, the following relationship is valid:

R_{0}(u_{max}) = 0.05 · 11.62 · β = 0.581 · β.  (24)

It is easy to see that for all values of β given in Table 1 the inequality

R_{0}(u_{max}) < 1

holds. It means that by using the quarantine with the maximum intensity u_{max} = 0.95, the epidemic will definitely end. In addition, it is clear that the value of u_{max} can be reduced, and an important problem arises of finding a threshold value for u_{max}.

The equality R_{0}(u_{max}^{new}) = 1 and the last expression of (23) lead us to the formula:

u_{max}^{new} = 1 − R_{0}^{−1}.  (25)

By substituting the values of R_{0} from Table 1 for (25), we calculate the corresponding threshold values u_{max}^{new}, which are placed in the last column of this table.

Next we provide numerical calculations for the boundary value problem (19), which were performed using MATLAB software. For these numerical calculations, the values of the parameters for system (11), the weighting coefficient of the objective function (10), and the control restriction from (8) given in (20) and also provided in Table 1 were used.

In our calculations we focus only on Mexico and Mexico City.

4.1 Optimal control of COVID-19 related to Mexico

According to [9], in 2020 the population of Mexico is estimated to be N_{0} = 127792286. We use the following values for the initial conditions s_{0}, e_{0}, i_{0}, a_{0} from (19):

\begin{matrix}s_{0} = 0.9999998122 (S_{0} = 127792262) & e_{0} = 1.56504 × 10^{−7} (E_{0} = 20) \\ i_{0} = 7.8252 × 10^{−9} (I_{0} = 1) & a_{0} = 2.34756 × 10^{−8} (A_{0} = 3). \end{matrix}

We assume that during the first 26 days of the epidemic, no quarantine was enacted. In our simulation, the initial conditions for the control system are taken from the values for s(t), e(t), i(t) and a(t) at t = 26. In Table 2 the official data related to Mexico, the results of the simulation in the period of time from February 28 till May 21 are provided.

Table 2: Confirmed active cases in Mexico between February 28 to May 21, 2020.

Data  March 24  April 21  May 21
Official  405  3185  12905
Simulation with no control \begin{matrix}403 \\ R_{0} = 5.8 \end{matrix}  \begin{matrix}3191 \\ R_{0} = 2.65 \end{matrix}  \begin{matrix}12571 \\ R_{0} = 2.0 \end{matrix}
Simulation with u = u_{∗}  403  1802  1904

In the simulation, the values are the sum of symptomatic and asymptomatic cases. The official data was adopted from [10]. The value of t_{f} was taken from {26; 54; 84}, which models twenty-six days, fifty-four and eighty-one days of protective measures. On the second line, the basic reproductive ratio R_{0} = 5.8 corresponds to the period from February 28 till March 24. Similarly, R_{0} = 2.65 and R_{0} = 2.0 are related to the period from March 24 till April 21 and that from April 21 till May 21, respectively. For the simulation with u = u_{∗}, the basic reproductive ratio R_{0} = 2.65 is employed.

4.2 Optimal control of COVID-19 in Mexico City

The population of Mexico City in 2020 is estimated by [8] to be N_{0} = 9018645. We use the following initial conditions s_{0}, e_{0}, i_{0}, a_{0} from (19):

\begin{matrix}s_{0} = 0.999997 (S_{0} = 9018621) & e_{0} = 2.21763 × 10^{−6} (E_{0} = 20) \\ i_{0} = 1.10881 × 10^{−7} (I_{0} = 1) & a_{0} = 3.32644 × 10^{−7} (A_{0} = 3). \end{matrix}

Table 3 shows the official data concerning Mexico City, the results of the simulation related to equation (11) in the period of time {26; 54; 84}.

In the simulation, the values are the sum of symptomatic and asymptomatic cases. The official data was adopted from [10]. The value of t_{f} was taken from {26; 54; 84}, which models twenty-six days, fifty-four and eighty-one days of protective measures. On the second line, the basic reproductive ratio R_{0} = 2.75 corresponds to the period from February 28 till March 24. Similarly, R_{0} = 3.55 and R_{0} = 1.9 are related to the period from March 24 till April 21 and that from April 21 till May 21, respectively. For the simulation with u = u_{∗}, the basic reproductive ratio R_{0} = 3.55 is employed.

Table 3: Confirmed active cases in Mexico City between February 28 to May 21, 2020.

Data  March 24  April 21  May 21
Official  66  883  3339
Simulation with no control \begin{matrix}65 \\ R_{0} = 2.75 \end{matrix}  \begin{matrix}854 \\ R_{0} = 3.55 \end{matrix}  \begin{matrix}3294 \\ R_{0} = 1.9 \end{matrix}
Simulation with u = u_{∗}  66  179  179

4.3 Graphs of s_{∗}(t), e_{∗}(t), i_{∗}(t), a_{∗}(t) and u_{∗}(t)

The graphs in Figure 2 correspond to the optimal solutions s_{∗}(t), e_{∗}(t), i_{∗}(t), a_{∗}(t) and optimal control u_{∗}(t) at the basic reproductive ratio R_{0} = 2.65 related to Mexico.

Figure 3 depicts the graphs of the optimal solutions s_{∗}(t), e_{∗}(t), i_{∗}(t), a_{∗}(t) and optimal control u_{∗}(t) at the basic reproductive ratio R_{0} = 3.55 related to Mexico City.

We can see that the graphs presented in Figures 2 and 3 support our analytical investigation. In this section we found that the strength of the quarantine measures can be reduced to the threshold value, u_{max}^{new} that leads to R_{0} < 1 and to eradication of pandemic. Thus, the optimal control does not take the values exceeding the calculated threshold values for each current ratio (of 0.62 and 0.718 for R_{0} = 2.65 and R_{0} = 3.55, respectively).

The optimal quarantine strategies for Mexico and Mexico City are quite similar: at first, the optimal control takes its maximal constant threshold value in the 3-6 weeks (depending on the reported basic reproductive ratio) and then begins to decrease according to almost linear law, remaining nonzero at the end of the time interval.

It follows from the graphs of the optimal solutions, that there is also a significant slowdown in the growth of infectious populations compared to those without control. Under optimal control, the graphs of symptomatic and asymptomatic infectious people are delayed in their exponential growth and instead are getting flatten on the first half of the time interval. Moreover, the optimal solutions are characterized by shifting in their future pick to the right and by decreasing in its value. optimizing the quarantine cost for suppression of the... 71

Figure 2: Mexico OC Problem: optimal solutions and optimal control for R_{0} = 2.65 and t_{f} = 58 days: upper row: s_{∗}(t), e_{∗}(t), i_{∗}(t); middle row: a_{∗}(t), r_{∗}(t), n_{∗}(t); lower row: i_{∗}(t) + a_{∗}(t), e_{∗}(t) + i_{∗}(t) + a_{∗}(t), u_{∗}(t).

Rev.Mate.Teor.Aplic. (ISSN print: 1409-2433; online: 2215-3373) Vol. 28(1): 55–78, Jan–Jun 2021

72 a. choque-rivero — e. khailov — e. grigorieva

Figure 3: Mexico City OC Problem: optimal solutions and optimal control for R_{0} = 3.55 and t_{f} = 58 days: upper row: s_{∗}(t), e_{∗}(t), i_{∗}(t); middle row: a_{∗}(t), r_{∗}(t), n_{∗}(t); lower row: i_{∗}(t) + a_{∗}(t), e_{∗}(t) + i_{∗}(t) + a_{∗}(t), u_{∗}(t).

Rev.Mate.Teor.Aplic. (ISSN print: 1409-2433; online: 2215-3373) Vol. 28(1): 55–78, Jan–Jun 2021

5 Conclusions

In this paper, at a given time interval, a SEIR-type model that describes the spread of the COVID-19 virus in a human population of variable size is considered. A bounded control function of time was introduced into the model, which reflected the intensity of quarantine measures conducted in the population. This control reflects all sorts of indirect measures (quarantine, mask wearing, and various educational and informational campaigns) aimed at reducing the possibility of transmission of the virus from infected to healthy individuals. For the resulting control model, the optimal control problem was stated, which consisted of minimizing the Bolza-type objective function.

The terminal part of this function determines the level of disease in the population caused by COVID-19 at the end of the quarantine period, and its integral part is a weighted sum of the cumulative level of disease over the entire quarantine period with the total cost of this quarantine.

A detailed analysis of the optimal solutions to the optimal control problem was made using the Pontryagin maximum principle. The properties of the optimal control and its uniqueness were established. Thus, we proved that the optimal control is nonzero even at the end of the time interval, which means that some quarantine measures should never vanish. The values for the control model parameters, based on the knowledge of the basic reproductive ratios, were then taken.

The results of numerical calculations for Mexico and Mexico City performed using MATLAB software support our analytical investigation. Thus, in this paper we evaluate the maximum strength of the quarantine measures (the threshold value, u_{max}^{new}) that mathematically leads to R_{0} < 1 and to the end of the pandemic. It is found that in order to eradicate the virus, there is not need for the optimal control to take its maximum value (there is no need for the strongest quarantine) but its threshold value depends on the current basic reproductive ratio in the region and changes from 0.5 − 0.83 for R_{0} = 2.0 − 6.0, respectively. Therefore in each case the best results (optimal solutions) in the fight against the pandemic can be achieved with less restrictions and costs.

By comparing the optimal solutions to those without control, it is clear that the optimal strategy does “flattening the curve”. Under optimal control policy, the curve of total infections is getting flatten, shifts to the right and the value of its predicted maximum decreases significantly, which means that hospitals will not be crowded to the maximum capacity with COVID-19 patients needed intensive care. Delaying the pick of the infection and reducing its maximum would help the government to escape catastrophe, prepare, to save the lives and to reduce the costs.

Additionally, based on our mathematical study, we can make the following conclusions (recommendations).

• It is necessary to keep the strongest protective quarantine measures (selfisolations, stay at home order) at the beginning of the planned period.

• The strength of the quarantine and of other protective measures must depend on the reported official basic reproductive ratio and can differ from country to country and from region to region.

• Any restrictions can be slowly reduced toward the end of the optimization period of time t_{f}.

• Some protective measures (wearing face masks or gloves, escaping close contacts, etc.) should never vanish with the end of the time interval and must become hygienic habits for the society.

It is our hope that by using optimal strategy many lives can be saved while keeping the economy running.

Funding

The research of the first author is supported by CONACYT Project A1-S-31524 and CIC-UMSNH, Mexico.

References

[1] I.H. Aslan, M. Demir, M.M. Wise, S. Lenhart, Modeling COVID-19: forecasting and analyzing the dynamics of the outbreak in Hubei and Turkey,

MedRxiv (2020), 1–17. Doi: 10.1101/ 2020. 04.11.20061952

[2] F. Brauer, C. Castillo-Chavez, Mathematical Models in Population Biology and Epidemiology, Springer, New York-Heidelberg-Dordrecht-London,

2012. Doi: 10.1007/ 978-1-4614-1686-9

[3] Z. Cao, Q. Zhang, X. Lu, D. Pfeiffer, Z. Jia, H. Song, D.D. Zeng, Estimating the effective reproduction number of the 2019-nCoV in China,

MedRxiv, (2020), 1–8. Doi: 10.1101/ 2020. 01.27.20018952

[4] T.-M. Chen, J. Rui, Q.-P. Wang, Z.-Y. Zhao, J.-A. Cui, L. Yin, A mathematical model for simulating the phase-based transmissibility of a novel coronavirus, Infect. Dis. Poverty. 9(2020), no. 24, 1–8.

Doi: 10.1186/ s40249-020-00640-3

[5] W.H. Fleming, R.W. Rishel, Deterministic and Stochastic Optimal Control, Springer-Verlag, Berlin, 1975. Doi: 10.1007/

978-1-4612-6380-7

[6] H. Gaff, E. Schaefer, Optimal control applied to vaccination and treatment strategies for various epidemiological models, Math. Biosci. Eng. 6(2009),

no. 3, 469–492. Doi: 10.3934/ mbe. 2009. 6.469

[7] Gobierno de la Ciudad de México, El modelo epidemiológico del Gobierno

de la Ciudad de México, 2020. https:// modelo. covid19. cdmx.

gob. mx/modelo-epidemico

[8] Gobierno de México, Proyecciones de la población de México y de las entidades federativas, 2016-2050, CONAPO, 2019. https:// www. gob.

mx/cms/ uploads/ attachment/  file/ 487395/ 09_CMX. pdf

[9] Gobierno de México, Indicadores demográficos de México de 1950 a

2050, CONAPO, http:// www. conapo. gob. mx/work/ models/

CONAPO/ Mapa_Ind_Dem18/ index_2.html

[10] Gobierno de México, Aviso epidemiológico – Casos de infección respiratoria asociados a Coronavirus (COVID-19), Dirección General de Epide-

miología, 2020. https:// www. gob. mx/salud/ es

[11] E.V. Grigorieva, E.N. Khailov, Optimal intervention strategies for a SEIR control model of Ebola epidemics, Mathematics 3(2015), no. 4,

961–983. Doi 10.3390/ math3040961

[12] E.V. Grigorieva, E.N. Khailov, Estimating the number of switchings of the optimal interventions strategies for SEIR control models of Ebola epidemics, Pure and Applied Functional Analysis 1(2016), no. 4, 541–572.

[13] E.V. Grigorieva, P.B. Deignan, E.N. Khailov, Optimal control problem for a SEIR type model of Ebola epidemics, Revista de Matemática: Teoría y Aplicaciones 24(2017), no. 1, 79–96. Doi: 10.15517/ RMTA. V24I1.

27771

[14] E. Grigorieva, E. Khailov, Optimal preventive strategies for SEIR type model of the 2014 Ebola epidemics, Dyn. Contin. Discrete Impuls. Syst. Ser. B Appl. Algorithms. 24(2017), 155–182.

[15] E. Grigorieva, E. Khailov, Determination of the optimal controls for an Ebola epidemic model, Discrete Cont. Dyn. Syst. Ser. S. 2018, vol. 11 (6),

1071–1101. Doi: 10.3934/ dcdss. 2018062

[16] J. Jia, J. Ding, S. Liu, G. Liao, J. Li, B. Duan, G. Wang, R. Zhang, Modeling the control of COVID-19: impact of policy interventions and meteorological factors, Electron J. Differ. Eq. 23(2020), 1–24.

https:// ejde. math. txstate. edu/ Volumes/ 2020/ 23/
jia. pdf

[17] Y. Jing, L. Minghui, L. Gang, Z.K. Lu, Monitoring transmissibility and mortality of COVID-19 in Europe, Int. J. Infect. Dis. (2020), 1–16.

Doi: 10.1016/ j.ijid. 2020. 03.050

[18] E. Jung, S. Lenhart, Z. Feng, Optimal control of treatments in a twostrain tuberculosis model, Discrete Cont. Dyn. Syst. Ser. B. 2(2002), no.4,

473–482. Doi: 10.3934/ dcdsb. 2002. 2.473

[19] U. Ledzewicz, H. Schättler, On optimal singular controls for a general SIR-model with vaccination and treatment, Discrete Cont. Dyn. Syst. 2011, supplement, 981–990.

[20] E.B. Lee, L. Marcus, Foundations of Optimal Control Theory, John Wiley&Sons, New York NY, USA, 1967. Doi: 10.2307/ 2343766

[21] Y. Liu, A.A. Gayle, A. Wilder-Smith, J. Rocklöv, The reproductive number of COVID-19 is higher compared to SARS coronavirus, J. Travel Med.

2020, 1–4. Doi: 10.1093/ jtm/ taaa021 .

[22] Z. Liu, P. Magal, O. Seydi, G. Webb, Understanding unreported cases in the COVID-19 epidemic outbreak in Wuhan, China, and the importance of major public health interventions, Biology 9(2020), no. 50, 1–12.

Doi: 10.3390/ biology9030050

[23] J.P. Mateus, P. Rebelo, S. Rosa, C.M. Silva, D.F.M. Torres, Optimal control of non-autonomous SEIRS models with vaccination and treatment, Discrete Cont. Dyn. Syst. Ser. S. 11(2018), no. 6, 1179–1199. Doi:10.3934/

dcdss. 2018067

[24] M. Park, A.R. Cook, J.T. Lim, Y. Sun, B.L. Dickens, A systematic review of COVID-19 epidemiology based on current evidence, J. Clin. Med. 2020, 9, 967, 1–13. Doi: 10.3390/ jcm9040967 [25] L.S. Pontryagin, V.G. Boltyanskii, R.V. Gamkrelidze, and

E.F. Mishchenko, Mathematical Theory of Optimal Processes, John Wiley & Sons, New York NY, USA, 1962.

[26] F. Saldaña, H. Flores-Arguedas, J.A. Camacho-Gutiérrez, I. Barrandas, Modeling the transmission dynamics and the impact of the control interventions for the COVID-19 epidemic outbreak, Math. Biosci.

Eng., 2020. https:// www. researchgate.  net/ publication/

340902930_Modeling_the_transmission_dynamics_and_
the_impact_of_the_control_interventions_for_the_
COVID-19_epidemic_outbreak

[27] H. Schättler, U. Ledzewicz, Optimal Control for Mathematical Models of Cancer Therapies: An Application of Geometric Methods, Springer, New York-Heidelberg-Dordrecht-London, 2015.

[28] C.J. Silva, D.F.M. Torres, Optimal control for a tuberculosis model with reinfection and post-exposure interventions, Math. Biosci. 244 (2013),

153 –164. Doi: 10.1016/ j.mbs. 2013. 05.005

[29] R. Smith?, Modelling Disease Ecology with Mathematics, AIMS, Springfield MO, USA, 2008.

[30] B. Tang, X. Wang, Q. Li, N.L. Bragazzi, S. Tang, Y. Xiao, J. Wu, Estimation of the transmission risk of the 2019-nCoV and its implication for public health interventions, J. Clin. Med. 9(2020), no. 462, 1–13.

Doi: 10.3390/ jcm9020462

[31] P. van den Driessche, J. Watmough, Reproduction numbers and subthreshold endemic equilibria for compartmental models of disease transmission, Math. Biosci. 180(2002), 29–48. Doi: 10.1016/

s0025-5564(02)00108-6

[32] World Health Organization, Coronavirus disease (COVID-2019) situa-

tion reports. https:// www. who. int/ emergencies/  diseases/

novel-coronavirus-2019/   situation-reports

[33] J.T. Wu, K. Leung, G.M. Leung, Nowcasting and forecasting the potential domestic and international spread of the 2019-nCoV outbreak originating in Wuhan, China: a modelling study, Lancet, 395(2020), no. 10225 , 689–

697. Doi: 10.1016/ S0140-6736(20)30260-9

[34] S. Zhao, Q. Lin, J. Ran, S.S. Musa, G. Yang, W. Wang,..., M.H. Wang, Preliminary estimation of the basic reproduction number of novel coronavirus (2019-nCoV) in China, from 2019 to 2020: A data-driven analysis in the early phase of the outbreak, Int. J. Infect. Dis. 92(2020), 214–217. Doi:

10.1016/ j.ijid. 2020. 01.050

[35] Z. Zhuang, S. Zhao, Q. Lin, P. Cao, Y. Lou, L. Yang, D. He, Preliminary estimation of the novel coronavirus disease (COVID-19) cases in Iran: A modelling analysis based on overseas cases and air travel data, Int. J. Infect. Dis. 94(2020), 29–31. Doi: /10.1016/ j.ijid. 2020. 01.

050