key: cord-1056388-uindxuj6
authors: Antonelli, E.; Loli Piccolomini, E.; Zama, F.
title: Switched forced SEIRDV compartmental models to monitor COVID-19 spread and immunization in Italy
date: 2021-06-25
journal: nan
DOI: 10.1101/2021.06.21.21259230
sha: 2f6495c616341c4777884e56edb8799389738e58
doc_id: 1056388
cord_uid: uindxuj6

This paper presents a new hybrid compartmental model for studying the COVID-19 epidemic evolution in Italy since the beginning of the vaccination campaign started on 2020/12/27 and shows forecasts of the epidemic evolution in Italy. The proposed compartmental model subdivides the population into six compartments and extends the SEIRD model proposed in [E.L.Piccolomini and F.Zama, PLOS ONE, 15(8):1, 17, 08 2020] by adding the Vaccinated population and framing the global model as a hybrid-switched dynamical system. Aiming to represent the quantities that characterize the epidemic behaviour from an accurate fit to the observed data, we partition the observation time interval into sub-intervals. The model parameters change according to a switching rule depending on the data behaviour and the infection rate continuity condition. In particular, we study the representation of the infection rate both as linear and exponential piecewise continuous functions. We choose the length of sub-intervals balancing the data fit with the model complexity through the Bayesian Information Criterion. The calibration of the model shows an excellent representation of the epidemic behaviour, and thirty days forecasts have proven to reproduce the infection spread reliably. Finally, we discuss different possible forecast scenarios obtained by simulating an increased vaccination rate.

Compartmental models are essential mathematical tools in the analysis of the evolution of 14 epidemics, for prediction and simulation of future strategies which can be used by governments 15 and policymakers to allocate sanitary and economic resources. The parameters of such models (paragraph 2.2). In paragraph 2.3, we discuss the algorithmic details of the calibration proce-91 dure consisting of a sequence of bound-constrained optimization problems defined by the model 92 switches. Finally, in paragraph 2.4 we briefly discuss how we use the model to make predictions. The SEIRDV model characterized by constant parameters is obtained from the SEIRD model [14] by adding the new compartment V representing the vaccinated population. The following differential system represents the populations' dynamics:

where the total population, assumed of constant size N , is subdivided into six compartments: Sus- Infected. Standard models, as well as our SEIRDV, assume this relationship to be linear.

The parameter α > 0 represents the incubation rate for the transition from Exposed to Infected 102 states. Such value relates to the incubation period A I as follows: A I = 1/α. The average 103 incubation ranges from 2 to 14 days (d) (see https://www.worldometers.info/coronavirus/ 104 coronavirus-incubation-period/). According to [15] , more than 97 percent of people who con-105 tract SARS-CoV-2 show symptoms within 11.5 days of exposure. Recently a comparative study 106 assesses the incubation period around 6.5 days [16] . The cited studies assess the value of α in the 107 interval [0.14, 0.5].

The parameter γ > 0 representing the removal rate relates to the average infectious period T I as 109 γ = 1/T I . At the beginning of the outbreak, an average value T I 20d has been measured [17], population age, the virus spread, medical care availability and treatments. Finally, the parameter 115 ν > 0 represents the vaccination rate. Its value is particularly useful in the prediction phase to 116 obtain different scenarios.

117 Important information about the epidemic development is obtained from the number of infection 118 cases generated from a single infectious individual, i.e. the basic reproduction number R 0 , defined 119 as follows (see details in appendix Appendix B):

(2)

It is well known that the epidemic occurs when R 0 > 1; however, this information refers to 121 the initial stage, assuming that the entire population is Susceptible. 

Then the hybrid model is represented as [10]:

with state variable (u, Θ) T and F Θ (t, u) as in (1) [19] chp 6). In this paper, we represent the infection rate as piecewise linear and exponential 131 interpolating functions, yielding to SEIRDV pwl and SEIRDV pwe models, respectively.

Let us define β k (t) the restriction of β(t) to the interval ∆ k , k = 1, . . . , p, and set the values 133 5 . CC-BY-NC-ND 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. (which was not certified by peer review) 

whereas SEIRDV pwe defines the infection rate as follows:

We observe that for both models it holds : interval ∆ k , is represented in Figure 2 , where the model populations, for t ∈ ∆ k , are given by . CC-BY-NC-ND 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. (which was not certified by peer review)

The copyright holder for this preprint this version posted June 25, 2021. ; https://doi.org/10.1101/2021.06.21.21259230 doi: medRxiv preprint (special care is taken in the case n = L · p to avoid a too small length of ∆ p ). 145 Generally, we are interested in keeping L as large as possible to guarantee a proper balancing 146 between the data fit and the model complexity, evaluated in terms of the number of parameters to 147 be calibrated. In section 3 we discuss the details of the choice of a proper value for L. 148 We describe now the parameter estimation process in a single sub-interval ∆ k . We collect the 

where the positive weights µ j are introduced to compensate different data scales. The bound set B is defined as contain the values estimated in the literature.

To solve the minimization problem (6) numerically, we use iterative solvers as discussed in To suitably choose the first starting guess Θ 0 , which has a fundamental role in the quality of 160 the final solution, we solve problem (6) on a unique short time interval T i of about ten days with 161 7 . CC-BY-NC-ND 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. 

Compute Θ k as the solution of (6) in the sub-interval ∆ k with starting guess Θ k−1 5. end

The forward differential problem (1) is solved by a fourth order variable step Runge-Kutta 165 method. The initial conditions in each sub-interval ∆ k are given by the observed values of the 166 compartments I, R, D, V in the initial day of ∆ k . Concerning the starting value E 0 of the Exposed 167 compartment, which is not available from data, in the starting interval ∆ 1 , we relate it with the 168 delay time t d between the contact with the infectious agent and the onset of symptoms or signs of 169 infection, as follows:

(see Section 3 for more details on the values of t d ). In the successive intervals ∆ k , k > 1, we set 171 E 0 as the last value of Exposed computed in the interval ∆ k−1 . The starting value of Susceptible 172 is the difference between the total population N and the sum of all the other compartments. 2. We set in (1) the parameters Θ σ = (α p , β p , γ p , η p , σ · ν p ), with σ > 1, to simulate an increased 181 vaccination rate. We compute also in this case the prediction using both the linear and 182 exponential β functions. To analyse the numerical solution of the calibration in the interval [t 0 , T ] constituted of n days, we compute, for any considered population, a Relative RESidual defined as:

Xd 2 i and the Bayesian Information Criterion (BIC) [21], defined as follows:

where N θ is the number of the estimated parameters, Xd i represents the acquired compartment 199 data and X i is the corresponding value computed by the calibrated model at day i, i = 1, . . . , n.

The BIC takes into account the number of model estimated parameters and tends to penalize the 201 inclusion of additional parameters. The lower this quantity, the better the model will be.

To set a convenient initial value E 0 for the Exposed compartment we run SEIRDV pwl and 203 SEIRDV pwe on the first 15 days, from 2020/12/27 to 2021/01/10, choosing E 0 as in (7) . CC-BY-NC-ND 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

The copyright holder for this preprint this version posted June 25, 2021. ; https://doi.org/10.1101/2021.06.21.21259230 doi: medRxiv preprint cinated populations, respectively. We can appreciate the good quality of data-fit of SEIRDV pwl 230 and SEIRDV pwe. In Figure 6 (b) we represent the difference between the infected population ob-231 tained by SEIRDV exp and SEIRDV pwl algorithms. We observe that the main differences occur 232 in the ∆ 3 and ∆ 4 intervals.

(a) (b) 

. CC-BY-NC-ND 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. In Figure 8 we plot the values of both the infection rate function β(t) (on the left) and the 233 reproduction number function R t (t) (on the right) computed as follows

The red line relative to the exponential model changes more rapidly when the epidemic spread 12 . CC-BY-NC-ND 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

The copyright holder for this preprint this version posted June 25, 2021. ; https://doi.org/10.1101/2021.06.21.21259230 doi: medRxiv preprint the incubation rate α, both models report a decreasing behaviour in the period ∆ 1 − ∆ 4 . It 243 corresponds to an incubation period between 2 and 4 days. The removal rate γ is very similar 244 for both methods and gives the following removal periods: 34 d(∆ 1 ), 25.9 d(∆ 2 ), 26.9 d(∆ 3 ), 245 34.2 d(∆ 4 ) producing the average removal period of 30.3 d.

Regarding the parameter η, we observe that the average value 2.4% obtained by SEIRDV pwl 

In this paragraph, we apply the calibrated SEIRDV pwl and SEIRDV pwe to make predictions. blue and a red star, respectively). Comparing the prediction curves relative to SEIRD pwl ( Figure   261 9 (a)) and SEIRD pwe (Figure 9 (b)) with the epidemic data represented by magenta empty circles 262 we can see that the exponential model is more accurate. 263 We highlight that the reported forecast refers to the vaccination rate ν = 0.0022, computed in 264 the calibration interval ∆ 4 corresponding to 133572 administration per day. From Table 3 13 . CC-BY-NC-ND 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. increased vaccination rate, we plot in Figure 10 the results given by the two models in 40 days 271 using the vaccination rates in Table 3 . Comparing the two models, we see that the SEIRDV pwe 272 gives the more realistic prediction. Finally, we present the prediction obtained using more homogeneous and smaller-scale data 275 acquired in the Emilia-Romagna region after performing the calibration on the same sub-intervals 276 as in (8). In Figure 11 we plot the prediction obtained with the same procedure as in Figure 9 .

Differently to what happens for the Italian case, in the linear model (Figure 11 (a) ) the prediction 278 obtained with the red curve is entirely inaccurate, whereas for the exponential model ( Figure 11 279 (b)) the red and blue lines define a region containing the Infected data. Therefore, SEIRDV pwe 280 can be used to make reliable predictions with both strategies. Finally, in Figure 12 we show the 281 results for increasing vaccination rates, as done in Figure 10 for the national data. The SEIRDV pwl 282 forecasts are now closer to the Infected data compared to SEIRDV pwe, confirming the importance 283 of having both methods available to make different predictions. The differences between the two models can be observed in 

. CC-BY-NC-ND 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

The copyright holder for this preprint this version posted June 25, 2021. ; https://doi.org/10.1101/2021.06.21.21259230 doi: medRxiv preprint target and both methods. For each target and vaccination rate, the Infected and Dead individuals 305 computed by SEIRDV pwe are smaller than SEIRDV pwl.

The pie graph of Figure 14 shows the percentage of people of all the compartments corre-307 sponding to 70% and 80% vaccination targets. In particular Figure 14 (a) represents SEIRDV pwl 308 reaching 70% vaccination target on 2021/08/29, at the calibrated vaccination rate (ν = 0.0095) 309 whereas Figure 14 

We have proposed two SEIRDV compartmental models, each involving six populations (Suscep- . CC-BY-NC-ND 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

The copyright holder for this preprint this version posted June 25, 2021. ; https://doi.org/10.1101/2021.06.21.21259230 doi: medRxiv preprint curve with a different, increasing, vaccination rate in the same period give the idea of how the 330 epidemic would evolve in these cases. 331 Finally, we applied the model calibrated up to 2021/06/12 to forecast the epidemic behaviour 332 when 70% − 80% of population has received one vaccine dose. At the present vaccination rate the 333 80% immunization is reached on 2021/10/28 whereas with an increase of 60% (best scenario) it is 334 reached on 2021/09/05.

Future studies and extensions of the proposed models will consider the limited duration of Although some vaccines are administered in two different doses we consider, in our model, as 345 vaccinated the people who received at least the first dose. Since some studies report that the 346 positive immunity effects obtained from a single dose are evident in the immediate days after the 347 vaccine, we can use this information for our short-term prevision.

18 . CC-BY-NC-ND 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. In the assumption that at disease free state S * = N we obtain (2). We note that it coincides with the R 0 value of a standard SEIR model [19] .

LitCovid: an open database of COVID-19 353 literature

A contribution to the mathematical theory of epidemics

Suihter: A new mathematical

This paragraph describes the calibration procedure in the interval [t 0 , T ] supposed of n days.Appendix B. Reproduction Number R 0