domain: math.math-st math sdv

text:
Matthieu Oliver
email: matthieu.oliver@mines-paristech.fr

Didier Georges
email: didier.georges@gipsa-lab.grenoble-inp.fr

Clémentine Prieur
email: clementine.prieur@univ-grenoble-alpes.fr

Spatialized Epidemiological Forecasting applied to Covid-19 Pandemic at Departmental Scale in France

Keywords: COVID-19, pandemic modeling and simulation, data assimilation, epidemiological forecasting, model reduction

published or not. The documents may come

Introduction

Over the past year, the spread of Covid-19 has had more serious consequences than expected. This pandemic has brought about massive organizational changes in many countries, with repeated lockdowns and the use of remote solutions for work, study and personal life. It has also had massive impacts on hospital care services with the need to expand hospital capacity and postpone operations for non-Covid-19 patients. Finally, the pandemic has killed several million people around the world [see, e.g., 2]. One of the most difficult specifics to model regarding the spread of Covid-19 is the amount of asymptomatic infected individuals, as described, for example, in [START_REF] Al-Tawfiq | Asymptomatic coronavirus infection: MERS-CoV and SARS-CoV-2 (COVID-19)[END_REF], which makes the pandemic difficult to track and has led to massive testing campaigns to assess the evolution of the pandemic and the effectiveness of health measures. The efficiency of different strategies of lockdown and testing has been studied [see, e.g., 5]. The incubation period of more than 10 days resulted in a lag between the implementation of restriction rules and their effect on hospital admissions, and made the effectiveness of restrictions more difficult to measure. Forecasting the number of infections and modeling hospital needs is critical to pandemic management and has been the subject of numerous articles over the past year.

Since the start of the epidemic, many studies have been carried out to model and forecast the spread of Covid-19 [see, e.g., [START_REF] Jordan | Optimization in the context of COVID-19 prediction and control: A literature review[END_REF], for a review on optimization studies]. The asymptomatic infected individuals were typically taken into account by building compartmental models that are more detailed than the usual susceptible-infected-recovered (SIR) model [see, e.g., [START_REF] Liu | Predicting the number of reported and unreported cases for the COVID-19 epidemics in China, South Korea, Italy, France, Germany and United Kingdom[END_REF][START_REF] Khajanchi | Forecasting the daily and cumulative number of cases for the COVID-19 pandemic in India[END_REF][START_REF] Berger | An SEIR Infectious Disease Model with Testing and Conditional Quarantine[END_REF][START_REF] Roques | Impact of lockdown on the epidemic dynamics of COVID-19 in France[END_REF][START_REF] Charpentier | COVID-19 pandemic control: balancing detection policy and lockdown intervention under ICU sustainability[END_REF], and references therein]. These works focus specifically on the effect of lockdown on the propagation of Covid-19 by taking into account detected and undectected infected individuals as well as hospitalizations. Forecasting the number of infections of Covid-19 has been studied in [START_REF] Liu | Predicting the number of reported and unreported cases for the COVID-19 epidemics in China, South Korea, Italy, France, Germany and United Kingdom[END_REF] by using a 3-step approach on the transmission rate during a pandemic outbreak with linear, then exponential increase in infections and finally an exponential decrease of the transmission rate due to major public interventions and social distancing measures.

In this paper, we focus on deterministic models which can be interpreted as limits for a large population size of stochastic models [see, e.g. , 8]. Most of deterministic compartmental models encountered in the literature on COVID-19 pandemic consist of systems of ordinary differential equations (ODEs). Then a common way to tackle spatial variation in such systems is to define regional compartments on homogeneous geographical sub-areas and to add coupling terms [see, e.g., 13, 14, 15, and references therein]. However, the increase in the number of parameters to describe such coupled systems of ODEs increase rapidly with the number of subareas. Then, PDE-based compartmental models is a natural alternative to describe more finely spatial heterogeneity. There exist a few papers analyzing COVID-19 pandemic with such models. Let us mention [START_REF] Viguerie | Diffusion-reaction compartmental models formulated in a continuum mechanics framework: application to COVID-19, mathematical analysis, and numerical study[END_REF][START_REF] Viguerie | Simulating the spread of COVID-19 via a spatially-resolved susceptible-exposed-infectedrecovered-deceased (SEIRD) model with heterogeneous diffusion[END_REF][START_REF] Guglielmi | Delay differential equations for the spatially resolved simulation of epidemics with specific application to COVID-19[END_REF] based on diffusion-reaction compartmental models, or [START_REF] Bertaglia | Spatial spread of COVID-19 outbreak in Italy using multiscale kinetic transport equations with uncertainty[END_REF] [see also [START_REF] Li | Prediction of Epidemics Dynamics on Networks with Partial Differential Equations: A Case Study for COVID-19 in China[END_REF] in which the authors couple a transport dynamic of the commuter population at large spatial scales, based on kinetic equations, with a diffusion model for non commuters at the urban scale.

In the present paper, our aim is to propose a spatially distributed compartmental model, accounting for spatial heterogeneity of the dynamics while including several features of the recent COVID-19 outbreak. We also propose to exploit meshless numerical scheme properties to simulate our model on a complex spatial domain. Finally we rely on very recent time-dependent model reduction techniques to lower short-term forecast complexity. More precisely, we propose a spatialized version of a detailed compartmental model inspired by [START_REF] Charpentier | COVID-19 pandemic control: balancing detection policy and lockdown intervention under ICU sustainability[END_REF], we then restrict the spatial domain to the French department of Isère and propose a calibration procedure over the period of the first French confinement, extending the whole methodology exposed in [START_REF] Guan | Transport effect of COVID-19 pandemic in France[END_REF]. This procedure can be adapted to any scale, provided that corresponding data is available. It allows to reproduce the macroscopical behaviour of Covid-19 pandemic on different geographical areas. A preliminary step to the calibration and a second contribution of this paper is the simulation of the model. Our model is defined as a set of spatio-temporal partial differential equations over a complex domain, namely the French department of Isère in our study. We propose a meshless method based on a finite difference scheme relying on radial basis function approximation of the spatial operators (the RBF-FD method, see [START_REF] Fornberg | Solving PDEs with radial basis functions[END_REF]) to solve it, as such an approach is known to be flexible to solve problems on complex domains. To the best of our knowledge, it is the first time RBF-FD is applied in the framework of epidemiology to study the spatio-temporal spread of a disease. Our methodology allows to simulate the spread of Covid-19 pandemic at a departmental level, and for each compartment.

However the simulation cost prevents from online short-term forecast. Therefore we propose as a final contribution of the paper to rely on reduced order modeling tools introduced in [START_REF] Bakhta | Epidemiological forecasting with model reduction of compartmental models. application to the covid-19 pandemic[END_REF] to compute short-term forecasts of infection number. The strategy consists in learning a timedependent reduced order model with few compartments from a collection of evaluations of our spatialized detailed model, by varying initial conditions, initial times and parameter values. A set of reduced bases is learnt in an offline phase while the projection on each reduced basis and the selection of the best projection is performed online, allowing short-term forecast of the global number of infections in the department.

Finally, we would like to stress that the original approach we introduce in this paper is generic and not limited to the analysis of COVID-19 dynamics. Indeed, it could be adapted to model and simulate other dynamics described by a model with spatially distributed parameters of the type diffusion-reaction on complex domains. Also, the time-dependent model reduction techniques we introduced could be leveraged to compute control strategies related to such dynamics.

The paper is organized as follows. In Section 2, we present an extension of the classical SIR epidemiological model that models the infections, recoveries and hospitalizations spatially over a given geographical area of interest. In Section 3, we present the meshless RBF-FD scheme for solving the set of differential equations from model proposed in Section 2. In Section 4, the model is calibrated using hospital data from Covid-19 outbreak in the department of Isère, France. In Section 5, we use model order reduction techniques to build a surrogate model and extrapolate the number of infected individuals on a 10-day horizon.

Spatialized SIRHU D +/-model

In this section, we present the usual susceptible-infected-recovered (SIR) compartmental model, then we add compartments to adapt to Covid-19 specificities. Finally, we propose an extension of this detailed model to account for the spatial spread of Covid-19 pandemic on a geographical setting, in our framework the Isère department in France.

Usual SIR model

The model presented in this paper is an extension of the usual SIR epidemiological model from [START_REF] Brauer | Mathematical Models in Population Biology and Epidemiology[END_REF]. The typical SIR model describes a population of constant size, where each individual can be either Susceptible S to a disease, Infected I with the disease or Recovered R from the disease. The transmission of the disease is ruled by the following equations:

dS (t) dt = -β(t)S (t)I(t), dI(t) dt = β(t)S (t)I(t) -γ(t)I(t), (1)
dR(t) dt = γ(t)I(t), with
S (t) + I(t) + R(t) = N
where β is the transmission rate of the disease and γ is the recovery rate.

The S IRHUD +/-model

More detailed models such as the S IRHUD +/-model in [START_REF] Charpentier | COVID-19 pandemic control: balancing detection policy and lockdown intervention under ICU sustainability[END_REF] can be better suited to model Covid-19 pandemic as they account for detected and undetected infected individuals I + and I -, detected and undetected recovered individuals R + and R -as well as hospitalized individuals H, individuals receiving intensive care U and deceased individuals D. This approach is particularly useful in the case of Covid-19 because undetected infections and hospital saturation both play a key role in the handling of the pandemic.

The flow diagram of this model is sketched out in Figure 1. The dynamics of each compartment is modelled by the following equations:

dS (t) dt = -β(t)S (t)I -(t), dI -(t) dt = β(t)S (t)I -(t) -λ 1 I -(t) -(a IR + a IH + a IU )I -(t), dI + (t) dt = λ 1 I -(t) -(a IR + a IH + a IU )I + (t), dR -(t) dt = a IR I -(t) -λ 2 R -(t), (2)
dR

+ (t) dt = a IR I + (t) + λ 2 R -(
S (t) + I -(t) + I + (t) + R -(t) + R + (t) + H(t) + U(t) + D(t) = N.
This model was calibrated in [START_REF] Guan | Transport effect of COVID-19 pandemic in France[END_REF] to fit the Covid-19 outbreak at the regional scale in France. In this paper, we aim to enrich this model to better account for spatial non-homogeneity of Covid-19 propagation.

Spatialization of the S IRHUD +/-model

In our model, we keep the same variables as in the S IRHUD +/-model from [START_REF] Charpentier | COVID-19 pandemic control: balancing detection policy and lockdown intervention under ICU sustainability[END_REF], but each state variable now depends on location x and time t. The governing equations of our model are then:

∂S ∂t (x, t) = -β(t)S (x, t)I -(x, t) + ∇(k(x, t)∇S (x, t)), ∂I - ∂t (x, t) = β(t)S (x, t)I -(x, t) -(λ 1 (t) + a IR + a IH + a IU )I -(x, t) +∇(k(x, t)∇I -(x, t)), ∂I + ∂t (x, t) = λ 1 (t)I -(x, t) -(a IR + a IH + a IU )I + (x, t), ∂R - ∂t (x, t) = a IR I -(x, t) -λ 2 (t)R -(x, t) + ∇(k(x, t)∇R -(x, t)), (3)
∂R + ∂t (x, t) = a IR I + (x, t) + λ 2 (t)R -(x, t) + ∇(k(x, t)∇R + (x, t)), ∂H ∂t (x, t) = a IH (I -(x, t) + I + (x, t)) -(a HR + a HU + a HD )H(x, t), ∂U ∂t (x, t) = a IU (I -(x, t) + I + (x, t)) + a HU H(x, t) -(a UR + a UD )U(x, t), ∂D ∂t (x, t) = a UD U(x, t) + a HD H(x, t),
where ∇ denotes the spatial gradient operator. Note that each of the compartments S , I -, R -and R + now feature a diffusion term k(x, t) which models local displacement of individuals inside geographical area of interest. On the contrary, we assume that individuals from compartments I + , H, U and D are isolated and do not transmit the disease. It means that we neglect, e.g., the spread of the virus by individuals tested positive (I + ) during eventual transportation to hospital. The diffusion coefficient k(x, t), that may depend on the position x and vary in time due to changes in sanitary measures, is modeled as follows:

k(x, t) = a di f f β(t)ρ(x),
where ρ(x) is the local population density, β(t) is the transmission rate and a di f f is a positive constant. It seems indeed reasonable (see, e.g., arguments in the french mobility report [START_REF] Atif | Report 1: Feedback on mobility during the COVID-19 epidemic[END_REF]) to assume that local movement is more intensive in high density regions with low sanitary restrictions. As a first approximation, we assume that parameter a di f f is constant, independent from compartment S , I -, R -or R + . This parameter will be calibrated with other uncertain parameters (see Section 4). Parameters a AB are derived from the probability of moving from compartment A to B and the average duration in compartment A. These variables are more precisely described in [START_REF] Guan | Transport effect of COVID-19 pandemic in France[END_REF]. Parameter λ 1 (t) corresponds to the proportion of undetected infected I -that are being tested positive, hence moving to compartment I + . This models virological testing such as PCR tests. λ 2 (t) corresponds to serological testing of recovered people from R -that are moved to R+ when tested. These parameters are time-dependent because the testing policies and capabilities have changed during the pandemic. Moreover they can also be used for optimal control, because individuals from I + isolate themselves and do not spread the virus. In this work, we considered that the serological testing was negligible: λ 2 (t) = 0. For the calibration presented in Section 4, virological testing was set to λ 1 (t) = 0.01 because there was very few testing during the period on which we calibrate the model. However, in the forecasting part of the study that is described in Section 5, this parameter is varied over a wide range of values to replicate different phases of the testing policies and capabilities. In this model, we do not account for the loss of immunity (which could be done by adding a transition from the recovered compartment to the susceptible one), this is because our modeling is done at scale of days to a couple of months, at which scale the loss of immunity is neglected.

A meshless approach for solving partial differential equations

In this work, we used the radial basis function finite difference (RBF-FD) method from [START_REF] Fornberg | Solving PDEs with radial basis functions[END_REF] to solve the set of partial differential equations defined by (3) in Section 2. To the best of our knowledge, this is the first implementation of the RBF-FD approach to simulate an epidemiological model. Meshless approaches are particularly well-suited to our framework because the domain on which we solve PDE's is a geographical area that has a complex boundary. In the following, we take N d points denoted (x i ) 1≤i≤N d in the interior of the domain. We also take N b points denoted (x i ) N d +1≤i≤N d +N b over the boundary. In the following of this section, we focus on compartment of susceptible individuals S . Note that all the compartments involved in (3) are treated in a similar way.

Radial Basis Functions

Radial basis functions (RBF) are functions ϕ(x k , x l ) that depend on the distance r = ∥x kx l ∥ between x k and x l . We will use the set of points (x i ) 1≤i≤N b +N d to get an approximation of variable S . For each point x i we denote I n (x i ) the indices of the n nearest neighbors of x i in (x j ) 1≤ j≤N d +N b . The set of points indexed by I n (x i ) defines the RBF-FD stencil of x i . A local approximation of S around x i can then be obtained from the radial basis functions ϕ k (x) = ϕ(x, x k ) f or k ∈ I n (x i ). In this paper, to compute the approximation S a of S , we used the inverse multiquadric functions:

ϕ(x, x k ) = 1 √ 1 + ϵ 2 r 2 with r = ∥x -x k ∥ , ϵ = 3 and stencil size n = 25.
Equation ( 4) below gives a local approximation of S around x i :

S a (x, t) = k∈I n (x i ) S i k (t)ϕ k (x). (4)
Similarly, we denote X a the approximation of X ∈ {I -, R -, R + } whose dynamic is governed by reaction-diffusion PDEs of system (3).

Radial Basis Function Finite Differences

The approximation of the variables around stencil points allows us to locally discretize the diffusion operators that appears in the set of equations (3).

Operator discretization on the domain

Using (4), we can approximate the diffusion operator for variable S appearing in (3) as follows:

∀ i ∈ [1, N d ] , ∇(k(x i , t)∇S a (x i , t)) = k∈I n (x i ) S i k (t)∇(k(x i , t)∇ϕ k (x i )).
The sum can be written in matrix form

∇(k(x i , t)∇S a (x i , t)) = B(x i )K(x i ) -1 S i (t), (5)
with B(x i ) containing the operator values at points in the stencil of x i , indexed by

I n (x i ) = { j 1 , . . . , j n }, B(x i ) = (∇(k(x i , t)∇ϕ j 1 (x i )), . . . , ∇(k(x i , t)∇ϕ j n (x i ))),
K is a n × n symmetric matrix containing the values of the radial basis functions at each couple of points in the stencil of x i ,

K(x i ) =             ϕ(x j 1 , x j 1 ) . . . ϕ(x j 1 , x j n ) . . . . . . . . . ϕ(x j n , x j 1 ) . . . ϕ(x j n , x j n )             , (6)
S i (t) is a vector containing the values of S a at each point of the stencil of x i :

S i (t) = (S a (x j 1 , t), . . . , S a (x j n , t)) T . (7)
The same diffusion operator approximations as given by ( 5) hold for the compartment variables X ∈ {I -, R -, R + } governed by reaction-diffusion PDEs of system (3).

Discretization of the boundary conditions

Boundary conditions can be defined by operator B c , leading for variable S to:

B c S (x, t) = S b (x, t), (8)
where S b (x, t) is the boundary term that is imposed on point x. Applying the same method as in the previous section to the boundary points, we obtain for points (x i ) N d +1≤i≤N d +N b of the boundary:

B c Φ(x i )K(x i ) -1 S i (t) = S b (x i , t), (9)
where K(x i ) and S i (t) are defined by ( 6) and ( 7) respectively, and Φ(x i ) is the vector containing the values of the radial basis functions between x i and the n nearest neighbors of

x i : Φ(x i ) = (ϕ j 1 (x i ), . . . , ϕ j n (x i ))
. Similar boundary operator approximations as given by ( 9) hold for other compartment variables X ∈ {I -, R -, R + } governed by reaction-diffusion PDEs of system (3).

Reducing the equations to a system of ordinary differential equations (ODEs)

We finally apply (5) to the diffusion terms in (3) which leads, for points (x i ) 1≤i≤N d in the interior of the domain, to:

Ṡ a (x i , t) = F S (S a (x i , t), I - a (x i , t)) + B(x i )K(x i ) -1 S i (t), İ- a (x i , t) = F I -(S a (x i , t), I - a (x i , t)) + B(x i )K(x i ) -1 I - i (t), İ+ a (x i , t) = F I + (I - a (x i , t), I + a (x i , t)), Ṙ- a (x i , t) = F R -(I - a (x i , t), R - a (x i , t)) + B(x i )K(x i ) -1 R - i (t), (10)
Ṙ+ a (x i , t) = F R + (I + a (x i , t), R + a (x i , t)) + B(x i )K(x i ) -1 R + i (t), Ḣa (x i , t) = F H (I - a (x i , t), I + a (x i , t), H a (x i , t)), Ua (x i , t) = F U (I - a (x i , t), I + a (x i , t), H a (x i , t), U a (x i , t)), Ḋa (x i , t) = F D (H a (x i , t), U a (x i , t)),
where F α is the source term of compartment α in (3). Additionally, for points x i on the boundary, the boundary condition for each compartment variable X ∈ {I -, R -, R + } can be written as

B c Φ(x i )K(x i ) -1 X i (t) = X b (x i , t)
. This set of equations with the boundary conditions can be summarized as the following algebraic-differential equations:

Ż1 d = F 1 (Z 1 d , Z 2 ) + D 1 Z 1 d + D 2 Z 1 b , Ż2 = F 2 (Z 1 d , Z 2 ), G = D 3 Z 1 d + D 4 Z 1 b ,
where G is the vector of boundary conditions for the compartment variables at each boundary point:

            S b (x N d +1 ) I - b (x N d +1 ) R - b (x N d +1 ) R + b (x N d +1 ) . . . . . . . . . . . . S b (x N d +N b ) I - b (x N d +N b ) R - b (x N d +N b ) R + b (x N d +N b )             , Z 1 d and Z 1
b are the vectors of variables S , I -, R -, R + evaluated at each point of the interior and boundary of the domain respectively, Z 2 is the vector of variables I + , H, U, D evaluated at each of the points (interior and on the boundary) of the domain, and matrices D 1 , D 2 , D 3 , D 4 contain the coefficients of the discretized diffusion and boundary operators using ( 5) and [START_REF] Charpentier | COVID-19 pandemic control: balancing detection policy and lockdown intervention under ICU sustainability[END_REF].

When applying Dirichlet conditions as proposed in what follows, we have D 3 = 0 and D 4 = I d , leading to the following system of nonlinear ODEs:

Ż1 d = F 1 (Z 1 d , Z 2 ) + D 1 Z 1 d + D 2 G, Ż2 = F 2 (Z 1 d , Z 2 ). ( 11
)

Application to the epidemiological model

The methodology described in the previous section allows us to solve the set of differential equations (3) on a domain given a set of collocation points in the interior of the domain and over its boundary. In our application we want to model the spread of Covid-19 in the department of Isère, France. However our methodology is generic and could be applied to any geographical area of interest, provided the relevant data are available.

Collocation points

For generating the collocation points in the interior of the domain, we used a two-dimensional Halton sequence [see, e.g., 20] defined on a square and only selected points in the interior of the region, as shown in Figure 2. We thus get a low discrepancy sequence of N d = 1000 points in the interior of the domain. Each collocation point was given the population density value of the closest city taken from [START_REF]Population des villes de France (Population of cities in France)[END_REF]. A polygon describing the boundary of the department was obtained from [START_REF] Lexman | Carte des départements (map of the departments)[END_REF] then interpolated to obtain N b = 800 points along the boundary.

Boundary conditions

As shown in Figure 3, the population density is very heterogeneous over the domain, with the south-east being very rural with densities under 100 hab/km 2 and the north-east and center of the department being a lot more urbanised. We chose to apply null boundary condition on the boundary points x i with density ρ(x i ) < d min = 200 hab/km 2 . On the rest of the boundary we assumed that the population of each compartment was equal to the average population of that compartment over the domain. More precisely, for

N d + 1 ≤ i ≤ N d + N b with ρ(x i ) ≥ d min , we assumed X b (x i , t) = 1 N d N d j=1
X(x j , t).

Initial Conditions

At t = 0, the proportion of individuals in each compartment should be defined in every point of the domain, which represents 8 × N d = 8000 parameters. These proportions are not known everywhere nor for each compartment and thus have to be calibrated. To reduce the number of parameters to be calibrated, the proportions of infected and hospitalized individuals are modelled by a uniform distribution over the domain. For i ≤ N d , we assume:

I - 0 (x i ) ∼ U(0, i 0 ), H 0 (x i ) ∼ U(0, h 0 ), I + 0 (x i ) = R + 0 (x i ) = R - 0 (x i ) = U 0 (x i ) = 0, S 0 (x i ) = 1 -I - 0 (x i ) -H 0 (x i ), (12)
where i 0 /2 and h 0 /2 are the initial proportions of infected and hospitalized individuals in the domain. Then we used a numerical integrator based on backward differentiation formulas to solve the stiff system of ordinary differential equations expressed in [START_REF] Di Domenico | Impact of lockdown on COVID-19 epidemic in Île-de-France and possible exit strategies[END_REF]. The time step was set equal to one day.

Model Calibration

In this section, we describe the calibration procedure we propose for our spatialized detailed model. It extends the procedure in [START_REF] Guan | Transport effect of COVID-19 pandemic in France[END_REF], based on hospital data from [START_REF]Données hospitalières relatives a l'épidémie de COVID-19[END_REF].

Initial calibration

The calibration procedure we propose for the spatialized model described by ( 3) is initialized with the calibration results of its non-spatialized version performed in [START_REF] Guan | Transport effect of COVID-19 pandemic in France[END_REF].

Variables to fit

We use data from [START_REF]Données hospitalières relatives a l'épidémie de COVID-19[END_REF], namely the number of hospitalizations H obs , of intensive care hospitalizations U obs , the death toll D obs and the number of recoveries from hospitalization R hobs . Variables H, U and D are already involved in the model description [START_REF] Atif | Report 1: Feedback on mobility during the COVID-19 epidemic[END_REF]. Time evolution of the number of recoveries R h from hospitalization is described by:

∂R h ∂t (x, t) = a HR H(x, t) + a UR U(x, t).
Because of the weekly fluctuations of the data, a rolling average of 7 days is applied to the data.

Calibration period

Following [START_REF] Guan | Transport effect of COVID-19 pandemic in France[END_REF], we calibrate our model over the first wave of Covid-19 in France which took place between March 17 th and May 11 th , as sanitary restrictions remained stable during this period. Moreover, we use hospital data which are reliable data, as opposed to test data that are subject to changes in testing capabilities and policies.

Calibration parameters

Our model depends on numerous uncertain parameters that are presented in this section. Because of the amount of uncertain parameters , we only calibrate a selection of them, fixing remaining parameters at the value calibrated in [START_REF] Guan | Transport effect of COVID-19 pandemic in France[END_REF] as detailed hereafter. Except parameters related to the transmission rate β(t), the parameters we calibrate in this section are assumed to be constant over time. These parameters are disease-related. These parameters will not be re-calibrated for forecast. Only those related to time dependent transmission rate will be recalibrated in Section 5.

Transmission rate β(t)

The transmission rate is one of the most sensitive parameters to tune the model [see, e.g., 10], hence its parametrization has to be chosen carefully. We used the one proposed in [START_REF] Liu | Predicting the number of reported and unreported cases for the COVID-19 epidemics in China, South Korea, Italy, France, Germany and United Kingdom[END_REF] with an additional non-zero asymptotic value β ∞ : β(t) = (β 0 -β ∞ ) exp -µ(t-τ) + β ∞ . Note that under this form, the transmission rate only depends on time. This assumption is motivated by the fact that sanitary measures were the same throughout the region of interest during the first wave of Covid-19 in France, which corresponds to the calibration period we have chosen.

This parametrization of β(t) is well-suited for our calibration period: after time τ, the transmission rate decreases exponentially at speed µ from it's original value β 0 to an asymptotic value β ∞ . An example can be seen on Figure 4. In the current form, the calibration of the transmission rate depends on 4 parameters. Note that for a larger scale implementation of our model, the transmission rate could benefit from a region-wise parametrization, depending on the scale at which sanitary measures are enforced.

Diffusion coefficient

The diffusion coefficient k(x, t) appearing in (3) accounts for the non-homogeneity of the transmission of the disease. It depends on a single parameter a di f f to be calibrated.

Compartment transition parameters

The parameters a XY appearing in (3) depend on transition times and probabilities. The detailed formulas to obtain a XY are written in [START_REF] Guan | Transport effect of COVID-19 pandemic in France[END_REF]. The transition times and probabilities coming into play are:

p a , p h and p u which denote the respective probabilities of needing no hospital care, hospitalization and intensive care, p hu , p hd , p ud which denote the respective probabilities of needing intensive care or dying while hospitalized or in intensive care, -N a and N s which denote the recovery times of asymptomatic and lightly symptomatic individuals with no hospitalization, -N ih , N hd , N hu , N ud , N hr and N ur which denote the average transition times between compartments I ± , H, U and D.

This totals 14 constants to be calibrated to account for recoveries.

Initial Conditions

As detailed in Section 3.4.3, the initial conditions are reduced to the initial proportions i 0 /2 and h 0 /2 of I -and H. If the calibration was done in a later phase of the pandemic, we could add initial proportions of recovered individuals (R ± ) as well as intensive care occupation (U).

Calibration process

In the above, we presented the variables to fit and the corresponding French hospital data, the set of parameter values for the initialization, the time period over which we perform the calibration and finally all the parameters to tune, including initial conditions. In this section, we detail the optimization procedure and we present the results of the calibration.

Loss function

As we use hospital data aggregated on the department of Isère to fit our model, we first compute H sim , U sim , D sim and R h sim by aggregating on the department of Isère simulated data obtained from our spatialized model. We then define the loss function to be minimized over the set of parameters p as:

L(p) = t max 0 X∈{H,U,D,R h } (X obs (t) -X sim (p, t)) 2 X obs (t) dt, (13)
with t max = 55 days, which is length of the calibration period, which corresponds to the first lockdown in France.

Minimization method

For minimizing the calibration loss function [START_REF] Gatto | Spread and dynamics of the COVID-19 epidemic in Italy: Effects of emergency containment measures[END_REF], we use a stochastic optimization Python package noisyopt proposed by [START_REF] Mayer | Noisyopt: A python library for optimizing noisy functions[END_REF]. We also give the solver bounds for each parameter equal to [v 0 /5, 5v 0 ] where v 0 is the nominal value obtained from Section 4.1.

Results

Only a subset of parameters is calibrated with the aforementioned procedure. Parameters p a , p h , p u , p hu , p hd , p ud , N a and N s are kept to the value obtained in the calibration of the nonspatialized version of our model performed in [START_REF] Guan | Transport effect of COVID-19 pandemic in France[END_REF]: p a = 0.9, p h = 0.9, p u = 0.2, p hu = 0.001, p hd = 0.176, p ud = 0.3, N a = 7.82 and N s = 15. This last point is commented in Section 4.5.4. The optimal parameters resulting from the data fitting as obtained on Figure 5 are given in Table 1.

β 0 µ τ β ∞ i 0 h 0 N ih N hd N hu N

Calibration limitations

The calibration was performed under lockdown condition thus the transmission rate β(t) was assumed exponentially decreasing. This assumption does not hold in the general case; the evolution of the transmission rate can follow more complex variations as will be shown in the forecasting of β that is carried out in Section 5. We also assumed that the probabilities of moving from one compartment to the other and the average duration in each compartment remained constant regardless of the sanitary measures, as being more disease-related parameters. Note however that this assumption may not hold when studying subsequent variants of Covid-19 and their propagation in a partially vaccinated population, hence these parameters should be included in the calibration procedure for studying this later phase of the pandemic.

Forecasting using reduced order modeling

As shown in Section 4, our spatialized model is able to reproduce the outbreak of Covid-19 at the scale of the department while giving a very localized information on hospital needs and infection number. However this model is computationally expensive compared to a nonspatialized model, this is why we propose to rely on reduced order modeling to fasten the online computations required for short-term forecasts of the infection number. In this section, we extend the model reduction tools introduced in [START_REF] Bakhta | Epidemiological forecasting with model reduction of compartmental models. application to the covid-19 pandemic[END_REF] to the treatment of our spatialized model. In the following of this section, we explain how to build a reduced basis for extrapolating variables β * obs and γ * obs defined in [START_REF] Giordano | Modelling the COVID-19 epidemic and implementation of population-wide interventions in Italy[END_REF]. We then detail how it is possible to use the reduced basis to compute online the n f orecast -day forecasts of the global infection number in the department of Isère from an observation period of n calib = 10 days.

Fitting epidemiological data with a time dependent SIR model

In this section, we project our spatialized 8-compartment model into a time-dependent SIR model described by [START_REF] Al-Tawfiq | Asymptomatic coronavirus infection: MERS-CoV and SARS-CoV-2 (COVID-19)[END_REF]. This is done by aggregating the spatial output of each compartment over the domain and regrouping separate compartments in either S , I or R as follows:

S (t) = N d i=1 S (x i , t), I(t) = N d i=1 I + (x i , t) + I -(x i , t) + H(x i , t) + U(x i , t), R(t) = N d i=1 R + (x i , t) + R -(x i , t) + D(x i , t).
Starting from an initial epidemiological state (I(t = 0) = i 0 , R(t = 0) = r 0 ), the evolution of the SIR model is ruled by parameters β(t) and γ(t). As a consequence, forecasting the total number of infections and recoveries over the domain can be done by extrapolating β(t) and γ(t). Let S obs (t), I obs (t) and R obs (t) denote the evolution of S (t), I(t), R(t) observed over a time period. We then define:

β * obs (t) := - N I obs (t)S obs (t) dS obs (t) dt , γ * obs (t) := N I obs (t) dR obs (t) dt . (14)
Under mild assumptions that are detailed in [START_REF] Bakhta | Epidemiological forecasting with model reduction of compartmental models. application to the covid-19 pandemic[END_REF], Equations ( 14) allow a simple SIR modeling of the observed data starting from the initial state (I 0 obs , R 0 obs ). Note that in [START_REF] Giordano | Modelling the COVID-19 epidemic and implementation of population-wide interventions in Italy[END_REF], one can replace observation data by simulation data from a model that outputs S , I and R over time, leading to the definition of β * sim and γ * sim that we will use for forecasting as detailed in Section 5.5.2. It is possible to infer I obs and R obs from the data provided by [START_REF]Données hospitalières relatives a l'épidémie de COVID-19[END_REF], applying an adjustment factor (see [START_REF] Bakhta | Epidemiological forecasting with model reduction of compartmental models. application to the covid-19 pandemic[END_REF][Section 4.1] for more details): I obs = f H obs , R obs = f R hobs + D obs , with f = 15 computed from the literature [START_REF] Mizumoto | Estimating the asymptomatic proportion of coronavirus disease 2019 (COVID-19) cases on board the Diamond Princess cruise ship[END_REF][START_REF] Di Domenico | Impact of lockdown on COVID-19 epidemic in Île-de-France and possible exit strategies[END_REF] and assessing the proportion of individuals needing hospitalization when infected with Covid-19.

The evolution of β * obs is shown on Figure 6, on which one observes the influence of lockdown periods and sanitary measures. Note that β * obs takes negative values around July, 2020. This reflects the very low tendency of hospital data in Isère during the summer of 2020. During that period, the SIR dynamics did not fit well the hospital data, leading to negative values for β * obs . However this only happened in a short interval of time during which the pandemic was stable.

Detailed model evaluations

In order to fasten the online computations we built a reduced order model. A key stone in the reduction of the model is the construction of a reduced basis from a set of evaluations of the detailed model. We built K different reduced bases (B k ) 1≤k≤K for extrapolating the short-term evolution of β obs and γ obs , by moving the initial time. We then selected the most suited reduced basis and used it for actual forecasting. Note that for the construction of the reduced bases, the detailed model described in (3) was slightly modified by considering that β and λ 1 do not depend on time. Then the model was run with different set of initial conditions and parameter values. The range of values we considered for initial conditions and for parameters β and λ 1 is provided in Table 2. Even if the bounds for each parameter were chosen to be consistent with the sanitary condition at the start of the second wave of Covid-19 in France (September, 2020), the idea is that the range of admissible values is chosen wide enough to match different calibration periods. Each parameter was sampled from a Halton sequence of length n runs = 500. The use of a Halton sequence allows to better explore the interval in which each parameter is defined. Other parameters were fixed from Table 1 in Section 4.5.3. Then, the detailed model was run with these n runs combinations of parameter values over a period of n eval = 200 days, producing a set of detailed outputs -the number of individuals in each of the 8 compartments over time -that are collapsed into SIR outputs using the following rules:

Parameter β 0 λ 1 i - 0 i + 0 r - 0 r + 0 h 0 u 0 k di f f Min 0.05
S (t) = S sim (t), I(t) = I + sim (t) + I - sim (t) + H sim (t) + U sim (t), (15)
R(t) = R + sim (t) + R - sim (t) + D sim (t).
For each run, we computed a realization of (β * sim ,γ * sim ) using ( 14) on the S IR outputs. The set of realizations of β * sim and γ * sim are respectively denoted B and G in the following sections. The infection number over time is shown in Figure 7 for a selection of 50 sets of parameter values. In Sections 5.3 and 5.4 we detail the construction of the reduced order model from the set of realizations of (β * sim , γ * sim ). On most of the simulations in Figure 7, we observe a wave growing over 50 days after the start of the simulation, as the bounds in Table 2 were chosen to be coherent with the sanitary conditions in France at the start of the second wave in France (September, 2020). However, recall that we widen the bounds in order to build a more flexible set of reduced bases, which explains why we observe other scenarios, such as simulations with a wave starting around 110 days after the starting time, or even simulations without any wave.

Even though Figure 8 was obtained from detailed model evaluations with β not depending on time, the projection to a simple SIR model of the output as described in [START_REF] Guan | Transport effect of COVID-19 pandemic in France[END_REF] led to a set of realizations of time-dependent β * sim .

Simulation output reduction

From the set of realizations of (β * sim , γ * sim ) shown in Figure 8 we built respective sets of functions B and G defined on n eval = 200 days for β and γ that can reproduce the typical variations of the detailed model. We used two different methods for building B and G: one based on Singular Value Decomposition and another one on a greedy algorithm. In the following, we detail the construction of B from B, the same steps were applied to construct G from G.

Singular Value Decomposition

This approach consists in discretizing the functions in B at a 1-day timestep. We then obtain a n eval × n runs matrix denoted M β containing all the β * sim in B. We then compute eigenvalues and eigenvectors of the matrix A := M T β M β . The n eval positive and ordered eigenvalues of A are shown in Figure 9. Then B S V D is composed with the n span eigenvectors with the largest eigenvalues. The basis we obtained for n span = 4 is shown on Figure 10. We notice on Figure 9 that the magnitude of the 10 largest eigenvalues is significantly larger than the rest, which is why with this approach we never used more than 10 functions to build B S V D afterwards. As can be seen on Figure 10, the elements in B S V D are not constrained to remain positive. Section 5.3.2 tackles this issue by constructing in a greedy way a basis of positive functions.

Greedy algorithm

The method is based on the algorithm called Enlarged Nonnegative Greedy (ENG) presented in [START_REF] Bakhta | Epidemiological forecasting with model reduction of compartmental models. application to the covid-19 pandemic[END_REF]. We first use a greedy algorithm to build a set A of nonnegative functions selected from B. This set is initialized as follows: a 0 = argmax a∈B (∥a∥ 2 ). Then, while A has k ≤ n span elements (A = {a 1 , . . . , a k }), we add the function a new ∈ B that has the maximal distance from S pan (A):

a new = argmax a∈B min c 1 ,...,c k ≥0 a - k i=1 c i a i 2 .
We thus obtain a set of n span nonnegative elements of B that will be used to build the set B ENG using the enlarge cone algorithm from [START_REF] Bakhta | Epidemiological forecasting with model reduction of compartmental models. application to the covid-19 pandemic[END_REF]:

k i * = arg max (α > 0 | [a i -αa k -j i σ j i a j ](t) > 0 ∀t ∈ [0, n eval ]) σ k i = σ k i + α k i * 2 while α k i * ≥ ϵ do α k i * = arg max (α > 0 | [a i -αa k -j i σ j i a j ](t) > 0 ∀t ∈ [0, n eval ]) σ k i = σ k i + α k i * 2 end end ψ i = a i -j i σ j i a j end
Note that, due to numerical integration errors, some of the curves in Figure 8 take a few values below zero. However in practice, applying Algorithm 1 to these functions, we obtained positive functions in B ENG as can be seen on Figure 11. The same happened for G ENG (defined from G as B ENG from B).

Building a set of reduced bases for short-term forecasting by moving the initial time

In this section, we detail the way B (obtained from SVD or ENG) was used to build a set of reduced bases for extrapolating β * obs .

Enrichment of the reduced bases

In Section 5.4.1 we built a set of reduced bases (B k ) 1≤k≤K for extrapolating β * obs from the reduction of our detailed model evaluations. Additionally, we can add well-chosen functions to these reduced bases thanks to prior knowledge of the evolution of β * obs . As can be seen on Figure 6, the variations of β * obs are mostly exponential. Starting at time t 0 and given a n calib -day evolution of β * obs , we can fit the following exponential function:

E f it (t, α 0 , α ∞ , α ′ ) = (α 0 -α ∞ ) exp(-α ′ t) + α ∞ ,
to β * obs ([t 0 , t 0 + n calib ]) using least squares and add the function E f it ([t 0 , t 0 + n calib + n f orecast ]) to the reduced basis. This is illustrated in Figure 13. In the following, this raw exponential extrapolation will be used as benchmark extrapolation of β * obs and γ * obs or in combination with the reduced bases obtained from B S V D and B ENG .

Short-term forecasts of the infection number

We now have K potential candidates for the reduced basis of β * obs extrapolation. In this section, we present the way we selected the best reduced basis, and then how we used it to get a forecast on the infection number.

Selection of the best forecast for beta

In order to compute a forecast on [t 0 +n calib +1, t 0 +n calib +n f orecast ], we divided the calibration period in T 1 = [t 0 , t 0 + n 1 -1] and T 2 = [t 0 + n 1 , t 0 + n calib ]. Then for each reduced basis B k = (b k i ) 1≤i≤n span , we defined the loss function:

L 1 k (a 1 , . . . , a n span ) = T 1 (β * obs (t) - n span i=1 a i b k i (t)) 2 dt. (16)
By optimizing it, we computed a function β k calib (t) = n span i=1 a * k i b k i (t) defined on [t 0 , t 0 + n calib + n f orecast ] fitting at best β * obs on T 1 . We then computed:

L 2 k (k) = T 2 (β * obs (t) -β k calib ) dt, (17)
and selected the function β k calib with the lowest evaluation of the loss L 2 k for extrapolating β * obs on [t 0 + n calib + 1, t 0 + n calib + n f orecast ].

Short-term forecast of the infection number

By applying the methodology presented in Section 5.5.1 to the parameter β, we obtained a calibration-forecast function β calib as presented on Figure 14. The methodology was applied in parallel to the parameter γ. From β calib and γ calib , we could run the time-dependent SIR model defined by (1) with initial conditions (I obs (t 0 ), R obs (t 0 )) and parameters β = β calib , γ = γ calib . The model output on the infection number is shown on Figure 15.

Comparison of the SVD and greedy algorithms

To evaluate the error of our models, we generated 10-day forecasts every 5 days between March 2020 and May 2021. In this section and in the following, the forecast error for a calibration starting at t 0 is:

e(t 0 ) = t 0 +n calib +n f orecast t 0 +n calib +1 (I sim (t) -I obs (t)) 2 I obs (t) dt.
We compared the SVD and Greedy algorithms with and without the addition of an exponential extrapolation of β * obs , the number of basis elements was set to n span = 5. The average error on the full timeframe are shown on Figure 16. The ENG and SVD algorithms performed very similarly, and the exponential enrichment of the reduced bases did not improve significantly the forecasting score. Therefore in the following we computed 10 day forecasts with the SVD algorithm, as it is much faster. Also, in the following sections, we did not use the exponential enrichment. 16: 10-day forecast error for the greedy algorithm (ENG), the enriched greedy algorithm (ENG exp ), the SVD algorithm (S V D), the enriched SVD algorithm (S V D exp ) and the exponential benchmark algorithm (exp). We can see that the enrichment does not improve significantly the forecasting score, and that the SVD and ENG algorithms perform very similarly.

Building an aggregated model

In this section, we chose to focus on the SVD algorithm that we evaluated on 10-day forecasts. The choice of n span is crucial because keeping fewer basis functions may not be sufficient to capture the complexity of the problem while taking a too large span of basis functions could lead to overfitting over the calibration and deteriorating the forecast accuracy of the reduced basis. We computed the forecast error on the infection number during a full year of 10-day forecasts for models built with 1 to 8 SVD basis functions in the reduced bases for both β and γ. The average output of these 8 models was also used as a forecast model, the so-called aggregated model, denoted by Agg hereafter. The error of each of these 9 models was computed using the forecast error introduced in Section 5.6. The results are shown on Figure 17 for the SVD algorithm. We can see that our aggregated model performed better than any individual model, and that this model reached a 5.5% error on 10 day forecasts. The visualization of some individual model predictions can be seen on Figure 19. Moreover, the aggregated model seems more robust than the exponential extrapolation, preventing from very large errors as can be seen, e.g., on Figure 18.

Epidemiological forecast using the detailed model

To evaluate the efficiency of the reduced modeling approach for forecasting, we compared its cost with the one consisting in forecasting the infection number directly from the detailed model. All the computation times mentioned below were obtained on a regular laptop. In the reduced basis approach, the offline phase consists of n runs = 500 evaluations of the detailed model over n days = 200 days. The detailed model runs in τ = 0.15 s.day -1 , which totals C o f f = n runs n days τ = 15000 seconds. The online phase has a cost of C on = 26 seconds, which is the time for projecting β * obs and γ * obs onto the reduced bases (B k ) 1≤k≤K and (G k ) 1≤k≤K , as described in 5.5.1. The cost of running the reduced model itself is negligible because it is run only once, with the optimal extrapolations of β and γ. Hence, the total cost of the reduced model approach for κ forecasts is κ C on + C o f f .

In the direct approach, the detailed model is evaluated over n calib + n f orecast = 20 days and the optimization of the 9 parameters presented in Section 5.2 required n eval = 600 evaluations, totalling a direct time of C = (n calib + n f orecast ) × n eval × τ = 1800 seconds per forecast.

For κ forecasts, the gain is

G(κ) = κ C
C o f f +κ C on . With κ ≈ 10 forecasts, the cost of the offline phase is compensated by the gain in the online phase.

Discussion

This study has shown that extrapolating β * obs and γ * obs using the exponential function as presented in 5.4.2 is a reasonable solution for forecasting the infection number at a 10-day timescale. However our approach led to a slightly better forecasting score with a model that is more robust, hence more reliable. We would like to point out that the forecasting score has an underlying error due to the changes in dynamics of the epidemic with sanitary measures that are taken. In fact, the model has large overestimations of the infection number every time a new measure is taken, as can be seen on the bottom plot of Figure 18. As a consequence, the forecasting error could not be reduced to a very low value without the introduction of the sanitary measures in the model. The model would perform better with no changes in the dynamics, however this model is still very useful because it gives the evolution of the pandemic if no measures were taken, which can help to measure the efficiency of the restriction rules.

Conclusion

In this paper, a spatialized extension of a detailed multi-compartmental epidemiological model has been proposed, allowing to reproduce the evolution of hospital needs at the scale of the department of Isère in France during the Covid-19 outbreak, while giving a detailed information about the geographical heterogeneity of the sanitary conditions. The partial differential equations that define the model were solved with a very efficient meshless scheme on a very irregular domain. From this model we built a reduced basis for extrapolating the transmission rate β and recovery rate γ obtained from hospital data and we used them in a surrogate model to output forecasts of the global infection number in the department of Isère. A simple exponential extrapolation proved to be efficient for the extrapolation of the transmission and recovery rates but the aggregation of surrogate models using different reduced bases gave a more robust forecast. The work presented in the present paper could be used to evaluate the efficiency of restriction rules that are taken by comparing the forecasts as a reference for the evolution of the pandemic with the actual evolution of the infection number when restriction rules are taken. Also, testing policies could be evaluated by controlling parameters λ 1 and λ 2 and by evaluating the dynamics of the pandemic while more people are being isolated after being tested. Note that the spatialized detailed model we introduced could be completed to account, e.g., for vaccination or emergence of a variant. We could add compartments for the vaccinated population by calibrating their own probability of transmission and severity of symptoms. In particular, this could allow to compare different vaccination policies. The spatialized model could take into account age categories using contact, testing and vaccination rates as well as symptom severity in each age category. This could be solved with the same approach, but would benefit from more precise hospital data involving the age of patients for calibration.

Figure 1 :

 1 Figure 1: Flow diagram of the pandemic model described by the set of equations in (2)

Figure 2 :

 2 Figure 2: 1000 collocation points in the interior of the domain obtained by Halton sequence.

Figure 3 :

 3 Figure 3: Population density attributed to each collocation point.

Figure 4 :

 4 Figure 4: Transmission rate β over time with coefficients β 0 = 0.25, µ = 0.2, τ = 5 and β ∞ = 0.05.

Figure 5 :

 5 Figure 5: Fitting of compartments H, U, D and R h in Isère during the first lockdown in France (March-May) using the spatialized model described by (3).

Figure 6 :

 6 Figure 6: Evolution of β * obs obtained from hospital data between March 2020 and July 2021. We can see that the dates of the sanitary measures match the change of variations of β * obs . At the very end of the profile we can see the start of the "delta wave" that occured during the redaction of this work.

Figure 7 :

 7 Figure 7: Evolution of the infection number over n eval days whith 50 sets of parameter values from the training set.

Figure 8 :

 8 Figure 8: Subsets of B (left) and G (right) obtained by applying (14) to the SIR trajectories outputted by the detailed model.

Figure 9 :

 9 Figure 9: Eigenvalues of matrix A.

Figure 10 :

 10 Figure 10: The set B S V D of functions obtained through the Spectral Values Decomposition-based algorithm with n span = 4.

Figure 11 :

 11 Figure 11: The set B ENG of functions obtained through the greedy algorithm combined with Algorithm 1 for n span = 4.

Figure 12 :

 12 Figure 12: Extraction of a reduced basis from the functions obtained with the detailed model evaluations and the reduced order modeling algorithms presented in Sections 5.3.1 and 5.3.2.

Figure 13 :

 13 Figure 13: Exponential fit and forecast of the observed β in Isère department over a 20-day window starting on March 17 th , 2020. Here t 0 = 0, n calib = 10, n f orecast = 10.

Figure 14 :

 14 Figure 14: 10-day fit followed by a 10-day forecast of β * obs using the SVD basis presented in 5.3.1.

Figure 15 :

 15 Figure 15: Forecast of the infection number in Isère department using the SVD reduced basis with n span = 4.

Figure

  Figure16: 10-day forecast error for the greedy algorithm (ENG), the enriched greedy algorithm (ENG exp ), the SVD algorithm (S V D), the enriched SVD algorithm (S V D exp ) and the exponential benchmark algorithm (exp). We can see that the enrichment does not improve significantly the forecasting score, and that the SVD and ENG algorithms perform very similarly.

Figure 17 :

 17 Figure 17: Forecast error of the aggregated model and of each individual model obtained with n span ∈ [1, 8]. The aggregated model is the average of models 3, 4 and 5. We can see that the aggregated model performs slightly better than any of these individual models.

Figure 18 :

 18 Figure 18: Comparison of forecasts obtained with the aggregated and exponential models at periods where the dynamics of Covid-19 changed quickly. We can see that the exponential extrapolation of β *obs is more prone to exploding forecasts while the aggregated model has a more robust behaviour. Additionnaly, the third plot shows an example of a change in the dynamics that creates a large error as the forecasts can not anticipate a brutal change, e.g., in sanitary restrictions.

Figure 19 :

 19 Figure 19: Visualization of the individual model forecasts in different phases of the pandemic. The bottom left plot illustrates the fact that the aggregated model performs better than any individual model.

Table 1 :

 1 Optimal calibration parameters.

Table 2 :

 2 Bounds for the training set of parameters.

		0.1 0.03 0.0015 0.03 0.01	9e-5	8e-6	1e-8	.
	Max	0.15 0.5 0.06 0.0025 0.06 0.1 1.1e-4 1.2e-5 5e-7

Acknowledgements

We would like to thank the reviewers for their valuable comments and suggestions. We also thank Olga Mula, Associate Professor at Eindhoven University of Technology, and the members of the SEEPIA research group and the MODCOV19 platform for the many enriching discussions.