key: cord-1047362-k18wa1r7
authors: Antonelli, Erminia; Piccolomini, Elena Loli; Zama, Fabiana
title: Switched forced SEIRDV compartmental models to monitor COVID-19 spread and immunization in Italy
date: 2021-11-12
journal: Infect Dis Model
DOI: 10.1016/j.idm.2021.11.001
sha: 7c100057b51320dcaa101c1345152c860567d709
doc_id: 1047362
cord_uid: k18wa1r7

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 in the first six months. 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. We tested the model on italian data and on local data from Emilia-Romagna region. The calibration of the model shows an excellent representation of the epidemic behaviour in both cases. Thirty days forecasts have proven to well reproduce the infection spread, better for regional than for national data. Both models produce accurate predictions of infected, but the exponential-based one perform better in most of the cases. 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 period, lethality rate. Moreover, through such models, it is possible to estimate the number of 18 secondary cases produced by a single infected person at start time (basic reproduction number 19 R 0 ) and during the epidemic evolution (effective time-dependent reproduction number R t ). In 20 particular, the trend of R t is of great importance to check the epidemic evolution over time. 21 The COVID-19 pandemic, caused by the Sars-CoV-2 virus, has renewed interest in studying 22 these models and a significant number of papers appeared on this subject since the beginning of 23 2020 (refer to LitCovid database for up to date literature [1] ). They differ each other for the 24 type of model proposed, the external events considered, such as movement restrictions imposed by 25 governments or quarantines, and the regions where models are applied. 26 Starting from the first SIR (susceptible (S), infected (I), and recovered (R)) model, proposed 27 in 1927 by Kermack and McKendrick [2] , several generalizations have been formulated over the 28 years by increasing the number of compartments, such as, for example, the susceptible -exposed - 29 infectious -recovered (SEIR) and the susceptible -exposed -infected -recovered -dead (SEIRD) 30 schemes. Further extensions have been proposed to model the COVID-19 outbreak considering the 31 different social distancing policies and control measures applied in the various geographic areas to 32 contain the epidemic spread. More compartments have been added, making the models more and 33 more complex (see [3] , [4] , [5] , [6] , to mention only a few of the most recent). 34 In this paper we intend to consider the effects of the vaccine on the epidemic spread during the 35 first months since the vaccination campaign started. We introduce a new scheme, named SEIRDV, such as surveillance, social distancing, social relaxation, quarantining, patient treatment/isolation 41 (see [8, 9] and references therein). Following the approach in [10], we introduce a switching rule 42 that governs the SEIRDV model state at any given time. Besides producing optimal fit to epidemic 43 data, introducing such a hybrid approach allows us to represent disease evolution when restric-44 tion policies and virus variants cause changes in fundamental parameters such as infection rate, 45 recovery periods, and death rates. Although switched models are widespread in various engineer- 46 ing applications, studies about epidemic models are less common; see, for example, [11] (SIRV), 47 [12] (SIR) and [13] (SEIRD). In particular the authors in [13] propose a hybrid SEIRD model with 

The model calibration is carried out by solving a sequence of constrained minimizations of the 59 weighted least-squares residuals between the measured epidemic data and the value of the state 60 variables, which satisfy the initial value ordinary differential system representing the SEIRDV 61 model.

This paper is an extension of our previous work [14] , where we proposed a SEIRD model (before 63 the availability of vaccines), with two different forcing functions, to monitor the first phase of the 64 evolution of COVID-19 in Italy (2020/02/24-2020/05/24). Compared to [14] we modify the model 65 as follows: we include the vaccinated compartment , we represent the proposed scheme into the 66 hybrid models theoretical frame, and finally we change the expression of the forcing functions.

Moreover in the calibration phase, we add weights in the fitting objective function and bounded 68 constraints, thus improving the model computational effectiveness and accuracy.

In the experimental section, we report the results obtained by our simulations on the Italian 70 national and regional epidemic data in the period 2020/12/27-2021/06/17. in the prediction phase since the vaccination campaign is at its first beginning. Nevertheless, we 77 believe that such an extension and the distinction in vaccination doses and ages will significantly 78 improve the monitoring of the autumn-winter epidemic spread.

Contributions. We summarize here the main contributions of this paper. Limitations. The present study is concerned with modelling and prediction COVID-19 evolution 90 during the first six month of vaccination campaign in Italy. We consider vaccination data regarding 91 a single vaccine dose without distinguishing the vaccine type.

The rest of the paper is organized as follows. In Section 2, we present the proposed switched 

Then the hybrid model is represented as [10]:

with state variable (u, Θ) T and F Θ (t, u) is the right-hand-side of the system (A.1) , and with model 101 parameters represented by the piecewise constant function Θ(t).

However, using a constant value for the infection rate β k does not represent the epidemic be-103 haviour in a sufficiently flexible way [14] ; therefore, we introduce a continuous time-dependent 104 infection rate β(t). In this case, the epidemic model is known in the literature as forced model (see 105 for example [17] chp 6). In this paper, we represent the infection rate as piecewise linear and expo-106 nential 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 108 4 J o u r n a l P r e -p r o o f β k ≡ β(t k ), k = 0, . . . , p . The SEIRDV_pwl defines the infection rate as:

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 We describe the calibration procedure in the interval [t 0 , T ] supposed of n days. We calibrate 121 the parameters on the sub-intervals ∆ k , k = 1, . . . p with uniform size L = n/p and n = L · p. In 122 the case n = L · p, the length of the last interval ∆ p is set as

In the case n − L · p > L/2 , the number of the intervals is increased by one (p = p + 1), and the 124 length of the last interval is n − (p − 1)L. In this way we avoid that the last interval is has too 125 small with respect to the previous ones. 

where the function U(t) is the linear interpolation of the finite difference approximation of V at each day. Let V 1 , V 2 , . . . , V n be the values obtained from the vaccination database on days 1, . . . , n,

we consider the finite difference approximation of V :

We describe now the parameter estimation process in a single sub-interval ∆ k . We collect the 133 observed data about infected, recovered, dead compartments in vectors I, R and D of size L and we

in ∆ k , we calculate the model parameters Θ k = (α k , β k , γ k , δ k ) solving a weighted constrained 137 nonlinear least squares problem of the form:

where the positive weights µ j are introduced to compensate different data scales. The bounded set B is defined as To suitably choose the first starting guess Θ 0 , which has a fundamental role in the quality of 145 the final solution, we compute the solutionΘ of problem (5) on a unique short time interval of 146 about ten days using as starting guess (1/6, 0.02, 0.05, 0.001) as discussed in appendix AppendixA. 147 We set Θ 0 =Θ. The forward differential problem (4) Concerning the starting value E 0 of the exposed compartment which is not available from data, 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) [19] , defined as follows:

where N θ is the number of the estimated parameters, Xd i represents the acquired compartment 182 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 184 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, in the hypothesis that previously acquired data is not available, and using the new infected value (I new ), available in the epidemic 194 data repository, which represents the number of daily new cases. We test, for t d = 1, . . . , 10, In Figure 5 reproduction number function R t (t) (on the right) computed as follows:

obtained by extending (A.2) to time dependent parameters β and γ. 217 We observe in Figure 7 that the red line relative to the exponential model changes more rapidly Moreover, the decreasing values of β(t) and R t in the interval ∆ 4 (Figure 7) suggest that the end 225 of epidemic growth could be close, as soon as R t becomes less than one. 226 We now discuss the parameters computed in the calibration step and reported in Table 2 .

Concerning the incubation rate α, both models report a decreasing behaviour in the period ∆ 1 −∆ 4 .

It corresponds to an incubation period between 1.1 and 4.6 days.

The removal rate γ is very similar for both methods and gives the following removal periods: 

We use SEIRDV to predict the future behaviour of the disease evolution in short-medium 

we obtain the vaccination rate ν p as the ratio between U and S.

In this paper we have adopted the following two strategies for prediction. 2. We set in (A.1) the parameters Θ σ = (α p , β p , γ p , η p , σ · ν p ), with σ > 1, to simulate an 243 increased vaccination rate. We compute also in this case the prediction using both the linear 244 and exponential β functions. We highlight that the reported forecast refers to the vaccination rate ν = 0.0018. From In the second and third lines of the table, we report the number of infected in the peak day 263 together with the date of peak if we suppose to increase the vaccination rate. Graphically, the 264 behaviour of infected is plotted in Figure 9 , where we represent the predictions given by the two 265 models in 40 days using the vaccination rates in Table 3 . Comparing the two models, we see that 266 the SEIRDV_pwe gives the more realistic prediction. 

268 Finally, we present the prediction obtained using more homogeneous and smaller-scale data 269 acquired in the Emilia-Romagna region after performing the calibration on the same sub-intervals 270 as in (6). In Figure 10 we plot the prediction obtained with the same procedure as in Figure 8 .

Differently to what happens for the Italian case, in the linear model (Figure 10 (a) ) the prediction 272 obtained with the red curve is entirely inaccurate, whereas for the exponential model (Figure 10 (b) ) 273 the red and blue lines define a region containing the infected data. Therefore, SEIRDV_pwe can 274 be used to make reliable predictions with both strategies. Finally, in Figure 11 we show the results 275 for increasing vaccination rates, as done in Figure 9 for the national data. The SEIRDV_pwl In Figure 12 

We have proposed two switched SEIRDV compartmental models, each involving six populations 294 (susceptibles, exposed, infected, recovered, dead and vaccinated) to analyse COVID-19 spread from Future studies and extensions of the proposed model will also include the possibility of re-infection 317 and the distinction between one and two vaccination doses. We believe that such an extension, 318 together with a more detailed population age structure, will give its decisive contribution in the 319 autumn-winter epidemic spread.

AppendixA. The SEIRDV model with constant parameters 321

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: The parameter α > 0 represents the incubation rate for the transition from exposed to infected 331 states. Such value relates to the incubation period A I as follows: A I = 1/α. parative study assesses the incubation period of COVID-19 around 6.5 days [21] . The cited studies 336 assess the value of α in the interval [0.14, 0.5]. Hence we consider α 0 = 1/6 as starting guess in 337 our model.

The parameter γ > 0 representing the removal rate relates to the average infectious period T I 339 as γ = 1/T I . At the beginning of the COVID-19 outbreak, an average value T I 20d has been 340 measured [22] , hence γ ∈ [0.03, 0.1]. The starting guess used in this case is γ 0 = 1/20.

After the period T I , the infected people split into recovered and dead with weights 1 − η and η, 342 respectively (0 ≤ η ≤ 1). Hence the parameter η represents the fraction of the removed individuals 343 who die and its value depends on environmental situations that change over time, such as the 344 population age, the virus spread, medical care availability and treatments. We use as starting 345 guess for this value is η 0 = 0.001. Finally, the parameter ν > 0 represents the vaccination rate. 

It is well known that the epidemic occurs when R 0 > 1; however, this information refers to the 351 initial stage, assuming that the entire population is susceptible. In the case of COVID-19, estimates 352 of R 0 in the interval [1.5, 6.68] were obtained during the first months of 2020 [23] . In the assumption that at disease free state S * = N we obtain (A.2). We note that it coincides 355 with the R 0 value of a standard SEIR model [17] . 4. We observe that the decreased vaccine efficacy does not affect the infected curve behaviour and 366 the peak day. On the contrary, the reduced efficacy of the vaccine causes the infected curves at 367 different vaccination rates to be closer each other (compared to Figure 9 ). Comparing the two 368 models, we notice that the exponential model is closer to the infected data and the number of 369 infected people is significantly lower. 

A contribution to the mathematical theory of epidemics

Suihter: A new mathematical 376 model for covid-19. application to the analysis of the second epidemic outbreak in italy

Evolution of disease transmission during the covid-19 pandemic: 378 patterns and determinants

A gener-380 alized mechanistic model for assessing and forecasting the spread of the covid-19 pandemic

Modeling vacci-384 nation rollouts, sars-cov-2 variants and the requirement for non-pharmaceutical interventions 385 in italy

Optimal time-varying vaccine allocation 387 amid pandemics with uncertain immunity ratios

Covid-19 optimal vaccination policies: a modeling study on efficacy, natural and 390 vaccine-induced immunity responses

A two-392 phase stochastic dynamic model for covid-19 mid-term policy recommendations in greece: a 393 pathway towards mass vaccination

Hybrid dynamical systems

Infectious disease models with time-varying parameters and general 398 nonlinear incidence rate

Dynamics behaviors of a hybrid 400 switching epidemic model with levy noise

Modeling and forecasting of 402 covid-19 using a hybrid dynamic model based on seird with arima corrections

Monitoring italian covid-19 spread by a forced 405 seird model

Effectiveness of 408 covid-19 vaccines against the b. 1.617. 2 (delta) variant

Global analysis of an epidemic model with 411 vaccination

Modeling infectious diseases in humans and animals

Covid-19 in 415 italy: Dataset of the italian civil protection department. Data in brief

Multimodel inference: understanding aic and 417 bic in model selection

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

Serial interval and incubation period of covid-19: a 423 systematic review and meta-analysis

Clinical course and risk factors for mortality of adult 427 inpatients with covid-19 in wuhan, china: a retrospective cohort study

R0 and re of covid-19: 429 Can we predict when the pandemic outbreak will be contained? Indian journal of critical 430 care medicine: peer-reviewed, official publication of Indian Society of Critical Care Medicine

On an seir 433 epidemic model with vaccination of newborns and periodic impulsive vaccination with eventual 434 on-line adapted vaccination strategies to the varying levels of the susceptible subpopulation