key: cord-0281670-jl9i75zv
authors: Zaitri, Mohamed Abdelaziz; Bibi, Mohand Ouamer; Torres, Delfim F. M.
title: Transport and Optimal Control of Vaccination Dynamics for COVID-19
date: 2022-02-01
journal: nan
DOI: 10.1016/b978-0-32-390504-6.00007-3
sha: c4090fb7f464546f3a4d7ae19b041b025b59140f
doc_id: 281670
cord_uid: jl9i75zv

We develop a mathematical model for transferring the vaccine BNT162b2 based on the heat diffusion equation. Then, we apply optimal control theory to the proposed generalized SEIR model. We introduce vaccination for the susceptible population to control the spread of the COVID-19 epidemic. For this, we use the Pontryagin minimum principle to find the necessary optimality conditions for the optimal control. The optimal control problem and the heat diffusion equation are solved numerically. Finally, several simulations are done to study and predict the spread of the COVID-19 epidemic in Italy. In particular, we compare the model in the presence and absence of vaccination.

BNT162b2 is an mRNA-based vaccine candidate against SARS-CoV-2, currently being developed by Pfizer and BioTech [9] . As announced on 9 th November 2020, BNT162b2 shows an efficacy against COVID-19 in patients without prior evidence of SARS-CoV-2 infection. A first interim efficacy analysis was conducted by an external, independent Data Monitoring Committee from the Phase 3 clinical study, and the case split, between vaccinated individuals and those who received the placebo, indicates a vaccine efficacy rate above 90%, at seven days after the second dose, of the 94 cases reviewed [20] .

The major obstacle that must be overcome is related to the process of transporting the vaccine, which must be stored at −70 o C [24] . Pfizer indicates that the vaccine will be distributed from its factories in the USA, Belgium and Germany. The American Wall Street Journal revealed that Pfizer has developed a special box packed with dry ice and a GPS tracker, which can hold 5000 doses of the vaccine under the right conditions for 10 days. Moreover, there is another obstacle related to the cost of the transportation boxes, where a similar box of 1200 doses in −8 o C costs 6868 USD, which is very expensive.

The transport of the vaccine must comply with the general standards for drug storage and the recommended conditions. Although many transport vehicles are equipped with refrigeration devices, assuring recommended storage conditions, simple insulated transport boxes are often used. In this study, we use the heat diffusion equation and assume that the shape of the vaccine bottle is cylindrical [18] . We perform the calculations to find out an initial temperature that ensures the arrival of the vaccine while fulfilling the required condition of −70 o C, by using insulated transfer boxes with the internal temperature at 0 o C [22] .

Optimal control is a mathematical theory that consists of finding a control that optimizes a functional on a domain described by a system of differential equations. This theory is applied in various fields of the engineering sciences: aeronautics, physics, biomedicine, etc. The Pontryagin minimum principle is used to find the necessary conditions for optimal controls [21] .

Several models were presented to predict the spread of COVID-19 [14, 16, 17, 23, 26, 28] . These studies used the SIR and SEIR models [11] and the generalized SEIR model [19] . Most of them were implemented to evaluate the strategy of the preventive measures [2, 3, 4, 13, 15, 27] .

In [1] , the authors present a mathematical model to analyze the Ebola epidemic and two optimal control problems related to the transmission of Ebola disease with vaccination. In [19] , the authors present a mathematical model to analyze the COVID-19 epidemic based on a dynamic mechanism that incorporates the intrinsic impact of hidden latent and infectious cases on the entire process of the virus transmission. The authors validate this model by analyzing data correlation and forecasting available general data. Their model reveals the key parameters of the COVID-19 epidemic. Here, we modify the model analyzed in [19] and consider an optimal control problem. More precisely, we introduce an extra variable for the number of vaccines used. Secondly, we study the associated optimal control problem, solving it numerically. Moreover, in order to find out the main parameters, we have performed a numerical simulation of the spread of COVID-19 in Italy from 01 st November 2020 to 31 th January 2021. Finally, we have presented another simulation to find the optimal control, and we have compared the models with and without vaccination.

The paper is organized as follows. We begin by formulating the vaccination transport model in Section 2. In Section 3, we recall the generalized SEIR model. Then, in Section 4, we formulate the generalized SEIR model with vaccination as an optimal control problem. The obtained optimal control problem is solved numerically in Section 5. In Section 6, we present a discussion concerning the spread of COVID-19 in Italy during three months, starting from 1 st November 2020. We end with Section 7 of conclusion, including some future research directions.

In this section, we present a model to maintain the effectiveness of the vaccine while transporting it from the factory storage area to the desired destination. The aim is to know the initial temperature that maintains the effectiveness of the vaccine, less than −70 o , and this by using the available mobile boxes at 0 o C. Thus, we propose the following mathematical model:

where T (t, x, y, z) represents the temperature of the vaccine at the point (x, y, z) and the time t; t * is the arrival time of the vaccine; and 0 o C is the temperature inside the box. The sets Ω and ∂Ω represent the interior and the border of the bottle containing the vaccine, respectively, r and h are the radius and height of the bottle, respectively, and α is the thermal diffusivity defined by

where k is the thermal conductivity, c ρ is the specific heat capacity, and ρ is the density.

The generalized SEIR model proposed by Peng et al. [19] is expressed by a sevendimensional dynamical system defined by

where the state variables are subjected to the following initial conditions:

In this model, the population is divided into the following compartments: susceptible individuals S(t), exposed individuals E(t), infected individuals I(t), quarantined individuals Q(t), recovered individuals R(t), death individuals D(t), and insusceptible/protected individuals P (t). These variables, in total, constitute the whole population, denoted by N:

The parameters ω, β, γ, δ, λ(t), and κ(t) represent the protection rate, infection rate, inverse of the average latent time, rate at which infectious people enter in quarantine, time-dependent recovery rate, and the time-dependent mortality rate, respectively. The recovery λ(t) and mortality κ(t) rates are analytical functions of time, defined by

where the parameters λ 1 , λ 2 , λ 3 , κ 1 , κ 2 and κ 3 are empirically determined in Section 6.

We now introduce the vaccine for the susceptible population in order to control the spread of COVID-19. Let us introduce in model (3) 

where t f represents the final time of the vaccination program. Hence, our model with vaccination is given by the following system of eight nonlinear ordinary differential equations:

where the state variables are subject to the initial conditions:

A schematic diagram of model (7) is given in Figure 1 . 

We consider the model with vaccination (7) and formulate an optimal control problem to determine the vaccination strategy u that minimizes the cost of treatment and vaccination:

where w 1 and w 2 represent the weights associated with the cost of treatment and vaccination, respectively. We assume that the control function u takes values between 0 and 1. When u(t) = 0, no susceptible individual is vaccinated at time t and if u(t) = 1, then all susceptible individuals are vaccinated at time t. Let x(t) = (x 1 (t), . . . , x 8 (t)) = (S(t), E(t), I(t), Q(t), R(t), D(t), P (t), W (t)) ∈ R 8 . The optimal control problem consists in finding the controlũ and the associated optimal trajectoryx, satisfying the control system (7) with the given initial conditions

where the controlũ ∈ Γ,

minimizes the objective functional (8) . With the new variables, problem (7)-(10) becomes

where

The existence of the optimal controlũ and the associated optimal trajectoryx comes from the convexity of the integrand of the cost functional (8) with respect to the control u and the Lipschitz property of the state system with respect to the state vector x(t) (see [7] for existence results of optimal solutions). According to the Pontryagin Minimum Principle [21] , ifũ ∈ Γ is optimal for the problem (11) with fixed final time t f , then there exists

, called the adjoint vector, such thaṫ

where the Hamiltonian H is defined by

The adjoint functions satisfyψ

wherē

holds almost everywhere on [0, t f ]. Moreover, the transversality conditions assert that ψ i (t f ) = 0, i = 1, . . . , 8. It follows from the the Pontryagin minimum principle that the extremal controlũ p is given bỹ

The current study aims to find the initial temperature to maintain the effectiveness of the vaccine during the transportation process, as well as determining an optimal vaccination strategy to limit the spread of COVID-19 in Italy. For that we reduce the costs of treatment and vaccination, during the three months starting from 1 st November 2020. We use the MATLAB R2020b program to perform all numerical computations. The initial conditions and real data are taken from the public database Dati COVID-19 Italia, available from https://github.com/pcm-dpc/COVID-19.

We assume r = 3 cm, h = 4 cm, α = 0.0137 W/(m · C), ρ = 2600 kg/m 3 , c ρ = 750 W · s/(kg · C) and t * = 7200 s, with the heat transfer coefficient equal to 1. In Fig. 2 we present the numerical solution of the heat diffusion equation (1), which gives the initial temperature equal to −94.5 o C.

We consider the following initial guesses: ω = 0.06, β = 1, γ = 5, δ = 0.5, (λ 1 , λ 2 , λ 3 ) = (0.01, 0.1, 10) and (κ 1 , κ 2 , κ 3 ) = (0.001, 0.001, 10).

The parameters of the generalized SEIR model are computed simultaneously by a nonlinear least-squares solver [8] . These parameters over the period starting from 1 st November 2020 till 31 th January 2021 are: ω = 0.0547, β = 0.5425, γ = 0.0873, δ = 0.3425, (λ 1 , λ 2 , λ 3 ) = (0.0999, 0.0501, 38.8542) and (κ 1 , κ 2 , κ 3 ) = (0.0021, 0.0125, 66.6652).

In Fig. 3 we show the recovery rate λ(t) and the mortality rate κ(t).

We fixed w 1 = w 2 = 1. The numerical solutions to the non-linear differential equations that represent the generalized SEIR model (3) , the generalized SEIR model with vaccination (7) , and the real data of the quarantined, recovered and death cases, from 1 st November till 6 th December 2021, are shown in Fig. 4 .

In Fig. 5 we present the optimal control (15)- (16) and the number of vaccines used starting from 1 st November 2020 till 31 th January 2021.

The orange curves in Fig. 4 represent the real data for the number of the quarantine, recovery and death cases in Italy starting from 1 st November till 6 th December 2020. The (7), respectively without and with vaccination, and real data of Italy from 1 st November till 6 th December 2021 with total population of N = 60.480.000. red curves in Fig. 4 represent the solutions of the generalized SEIR model (3) without vaccination, and they simulate what happen from the beginning of November to the end of January. There is an increase in the number of the recovered, death and insusceptible cases that reach, respectively, 1.830.000, 74.050 and 58.130.000 cases. The red curves for both the number of infected and quarantined individuals have their higher limit values of 103500 cases on 11 th November and 798.500 on 25 th November, respectively, reaching the values 614 and 22.640 cases on 31 th January 2021, respectively. We note that the number of susceptible individuals gradually decrease, reaching 416.600 cases at the end of January 2021.

The green curves in Fig. 4 represent the solutions of the generalized SEIR model (7) with vaccination, and they simulate what happened from the beginning of November to the end of January. There is an increase in the number of recovered, death and insusceptible cases that reach, respectively, 1.135.000, 60.560 and 3.076.000 cases. The green curves for both the number of infected and quarantined individuals have their higher limit values of 84.800 cases on 4 th November and 577.600 cases on 15 th November, respectively, reaching 55 and 7.237 cases on 31 th January 2021, respectively. We note that the number of susceptible individuals decrease rapidly reaching 0 cases on 19 th November 2020.

The red curve in Fig. 5 shows that the optimal vaccination of 100 percent of the susceptible individuals takes 19 days, followed by a rapid decrease in the number of susceptible individuals, which means they move to the class of vaccinated. The green curve in Fig. 5 shows the necessary number of vaccines to eliminate COVID-19, which is estimated at 56.200.100 doses. The total number of vaccinated and insusceptible individuals equal to 59.276.100 of the total Italian population of 60.480.000.

Our results show the importance of the vaccine for COVID-19 control and also the best result that could be obtained if the number of available vaccines satisfies the needs of the population and are distributed according with the theory of optimal control.

Here our optimal control problem has only one control: the vaccine. In reality, there are several other factors to take into account and other variables to control. In a future work, we would like to use the support maximum principle [5, 6, 10] , as well as the hybrid direction method [25] , to elaborate a primal-dual method for solving a more realistic optimal control problem, in presence of multiple inputs [12] .

Ebola model and optimal control with vaccination constraints

Estimating the parameters of susceptible-infected-recovered model of COVID-19 cases in India during lockdown periods

Hamiltonian Structure of Compartmental Epidemiological models

Optimal Control and Cost-Effectiveness Analysis for Dysentery Epidemic Model

Methods for solving linear-quadratic problems of optimal control

Optimization of a linear dynamic system with double terminal constraint on the trajectories, Optimization

Optimization-theory and applications, Applications of Mathematics

Generalized SEIR Epidemic Model (fitting and computation)

Safety and Immunogenicity of Two RNA-Based Covid-19 Vaccine Candidates

Optimality and suboptimality criteria in a quadratic problem of optimal control with a piecewise linear entry

A contribution to the mathematical theory of epidemics

Primal-dual method for solving a linear-quadratic multi-input optimal control problem

The incubation period of coronavirus disease 2019 (COVID-19) from publicly reported confirmed cases: Estimation and application

A new compartmental epidemiological model for COVID-19 with a case study of Portugal

A conceptual model for the coronavirus disease 2019 (COVID-19) outbreak in Wuhan, China with individual reaction and governmental action

Genomic characterisation and epidemiology of 2019 novel coronavirus: Implications for virus origins and receptor binding

Fractional model of COVID-19 applied to Galicia

Heat Conduction

Epidemic analysis of COVID-19 in China by dynamical modeling

Safety and efficacy of the BNT162b2 mRNA Covid-19 vaccine

The mathematical theory of optimal processes

Solving Diffusion Equations with Rough Coefficients in Rough Grids

Optimal control of the COVID-19 pandemic: controlled sanitary deconfinement in Portugal

A prefusion SARS-CoV-2 spike RNA vaccine is highly immunogenic and prevents lung infection in non-human primates

A hybrid direction algorithm for solving optimal control problems

An optimal control analysis of a COVID-19 model

Prediction of the COVID-19 spread in African countries and implications for prevention and control: A case study in South Africa

A stochastic time-delayed model for the effectiveness of Moroccan COVID-19 deconfinement strategy

This research is part of first author's Ph.D. project. Zaitri is grateful to the financial support from the Ministry of Higher Education and Scientific Research of Algeria; Torres acknowledges the financial support from CIDMA through project UIDB/04106/2020.