Epidemiological Forecasting with Model Reduction of Compartmental Models. Application to the COVID-19 Pandemic

Simple Summary

Using tools from the reduced order modeling of parametric ODEs and PDEs, including a new positivity-preserving greedy reduced basis method, we present a novel forecasting method for predicting the propagation of an epidemic. The method takes a collection of highly detailed compartmental models (with different initial conditions, initial times, epidemiological parameters and numerous compartments) and learns a model with few compartments which best fits the available health data and which is used to provide the forecasts. We illustrate the promising potential of the approach to the spread of the current COVID-19 pandemic in the case of the Paris region during the period from March to November 2020, in which two epidemic waves took place.

Abstract

We propose a forecasting method for predicting epidemiological health series on a two-week horizon at regional and interregional resolution. The approach is based on the model order reduction of parametric compartmental models and is designed to accommodate small amounts of sanitary data. The efficiency of the method is shown in the case of the prediction of the number of infected people and people removed from the collected data, either due to death or recovery, during the two pandemic waves of COVID-19 in France, which took place approximately between February and November 2020. Numerical results illustrate the promising potential of the approach.

1. Introduction

Providing reliable epidemiological forecasts during an ongoing pandemic is crucial to mitigate the potentially disastrous consequences for global public health and the economy. As the ongoing pandemic of COVID-19 sadly illustrates, this is a daunting task in the case of new diseases due to the incomplete knowledge of the behavior of the disease and the heterogeneities and uncertainties in the health data count. Despite these difficulties, many forecasting strategies exist, and we can cast them into two main categories: the first type is purely data-based and involves statistical and learning methods such as time series analysis, multivariate linear regression, grey forecasting or neural networks [1,2,3,4,5]; the second approach uses epidemiological models, which are appealing since they provide an interpretable insight of the mechanisms of the outbreak. They also provide high flexibility in the level of detail to describe the evolution of a pandemic, ranging from simple compartmental models that divide the population into a few exclusive categories to highly detailed descriptions involving numerous compartments or even agent-based models (see, e.g., [6,7,8] for general references on mathematical epidemiological models and [9,10,11] for some models focused on COVID-19). One salient drawback of using epidemiological models for forecasting purposes lies in the very high uncertainty in the estimation of the relevant parameters. This is due to the fact that the parameters cannot often be inferred from real observations, and the available data are insufficient or too noisy to provide any reliable estimation. The situation is aggravated by the fact that the number of parameters can quickly become large even in moderately simple compartmental models [10]. As a result, forecasting with these models involves making numerous a priori hypotheses which can sometimes be difficult to justify by data observations.

In this paper, our goal is to forecast the time-series of infected, removed and dead patients with compartmental models that involve as few parameters as possible in order to infer these series solely from the data. The available data are only given for hospitalized people; one can nevertheless estimate the total number of infected people through an adjustment factor taken from the literature. Such a factor takes into account the proportion of asymptomatic people and infected people who do not go to hospital. The model that includes the least number of parameters is probably the susceptible–infected–removed (SIR) model [12], which is based on a partition of the population into the following groups:

• Uninfected people, called susceptible (S);
• Infected and contagious people (I), with more or less marked symptoms;
• People removed (R) from the infectious process, either because they were cured or unfortunately died after being infected.

If N denotes the total population size that we assume to be constant over a certain time interval $[0,T]$, we have

$$N=S\left(t\right)+I\left(t\right)+R\left(t\right),\forall t\in [0,T],$$

and the evolution from S to I and from I to R is given for all $t\in [0,T]$ by

$$\begin{array}{cc}\hfill \frac{dS}{dt}\left(t\right)& =-\frac{\beta I\left(t\right)S\left(t\right)}{N}\hfill \\ \hfill \frac{dI}{dt}\left(t\right)& =\frac{\beta I\left(t\right)S\left(t\right)}{N}-\gamma I\left(t\right)\hfill \\ \hfill \frac{dR}{dt}\left(t\right)& =\gamma I\left(t\right).\hfill \end{array}$$

The SIR model has only two parameters:

• $\gamma >0$ represents the recovery rate. In other words, its inverse ${\gamma }^{-1}$ can be interpreted as the length (in days) of the contagious period;
• $\beta >0$ is the transmission rate of the disease. It essentially depends on two factors: the contagiousness of the disease and the contact rate within the population. The larger this second parameter is, the faster the transition from susceptible to infectious will be. As a consequence, the number of hospitalized patients may increase very quickly and may lead to a collapse of the health system [13]. Strong distancing measures such as confinement can effectively act on this parameter [14,15], helping to keep it low.

Our forecasting strategy is motivated by the following observation: by allowing the parameters $\beta$ and $\gamma$ to be time-dependent, we can find optimal coefficients ${\beta }^{*}\left(t\right)$ and ${\gamma }^{*}\left(t\right)$ that exactly fit any series of infected and removed patients. In other words, we can perfectly fit any observed health series with an SIR model with time-dependent coefficients.

As we explain below, the high fitting power stems from the fact that the parameters $\beta$ and $\gamma$ are searched in $L^{\infty }([0,T],{\mathbb{R}}_{+})$—the space of essentially bounded measurable functions. For our forecasting purposes, however, this space is too large to give any predictive power, and we need to find a smaller manifold that simultaneously has good fitting and forecasting properties. To this end, we developed a method based on model order reduction. The idea of the method was to find a space with a reduced dimension that can host the dynamics of the current epidemic. This reduced space is learnt from a series of detailed compartmental models based on precise underlying mechanisms of the disease. One major difficulty in these models is the fitting of the correct parameters. In our approach, we do not seek to estimate these parameters; instead, we consider a large range of possible parameter configurations with a uniform sampling that allows us to simulate virtual epidemic scenarios in a longer range than the fitting window $[0,T]$. We next cast each virtual epidemic from the family of detailed compartmental models into the family of SIR models with time-dependent coefficients, as explained below. This procedure yields time-dependent parameters $\beta$ and $\gamma$ for each detailed virtual epidemic. The set of all such $\beta$ (or $\gamma$) is then condensed into a reduced basis with a small dimension. We finally use the available health data on the time window $[0,T]$ to find the functions $\beta$ and $\gamma$ from the reduced space that best fit the current epidemic over $[0,T]$. Since the reduced basis functions are defined over a longer time range $[0,T+\tau ]$ with $\tau >0$ (e.g., two weeks), the strategy automatically provides forecasts from T to $T+\tau$. Its accuracy will be related to the pertinence of the mechanistic mathematical models that have been used in the learning process.

We note that an important feature of our approach is that all virtual simulations are considered equally important in the first stage, and the procedure automatically learns what the best scenarios (or linear combinations of scenarios) to describe the available data are. Moreover, the approach even mixes different compartmental models to accommodate these available data. This is in contrast to other existing approaches which introduce a strong a priori belief regarding the quality of a certain particular model. Note also that we can add models related to other illness and use the large manifold to fit a possible new epidemic. It is also possible to mix the current approach with other purely statistical or learning strategies by means of expert aggregation. One salient difference with these approaches which is important to emphasize is that our method hinges on highly detailed compartmental models which have clear epidemiological interpretations. Our collapsing methodology into the time-dependent SIR is a way of “summarizing” the dynamics into a few terms. One may expect that reducing the number of parameters comes at the cost of losing the interpretability of parameters, and this is true in general. Nevertheless, the numerical results of the present work show that a reasonable tradeoff between the “reduction of the number of parameters” and “interpretability of these few parameters” can be achieved.

The paper is organized as follows. In Section 2, we present the forecasting method in the case of a single region with a constant population. For this, in Section 2.1, we briefly introduce the epidemiological models involved in the procedure, namely the SIR model with time-dependent coefficients and more detailed compartmental models used for the training step. In Section 2.2, after proving that the SIR model with time-dependent coefficients in $L^{\infty }\left([0,T]\right)$ is able to fit any admissible epidemiological evolution (as explained below), we present the main steps of the forecasting method. The method involves a collapsing step from detailed models to SIR models with time-dependent coefficients and model reduction techniques. We detail these points in Section 2.3 and Section 2.4. In Section 3, we explain how the method can easily be extended to a multi-regional context involving population mobility and regional health data observations (provided, of course, that mobility data are available). In Section 3.1, we begin by clarifying that the nature of the mobility data will dictate the kind of multi-regional SIR model to use in this context. In Section 3.2, we outline how to adapt the main steps of the method to the multi-regional case. Finally, in Section 4, we present numerical results for the the two pandemic waves of COVID-19 in France in 2020, which took place approximately between February and November 2020. Concluding comments are given in Section 5, followed by two Appendix A and Appendix B that present details about the processing of the data noise and the forecasting error.

2. Methodology for a Single Region

For the sake of clarity, we first consider the case of a single region with a constant population and no population fluxes with other regions. Here, the term region is generic and may be applied to very different geographical scales, ranging from a full country to a department within a country or even smaller partitions of a territory.

2.1. Compartmental Models

The final output of our method is a mono-regional SIR model with time-dependent coefficients as explained below. This SIR model with time-dependent coefficients is evaluated with reduced modeling techniques involving a large family of models with finer compartments proposed in the literature. Before presenting the method in the next section, we here introduce all the models that we use in this paper along with useful notations for the rest of the paper.

2.1.1. SIR Models with Time-Dependent Parameters

We fit and forecast the series of infected and removed patients (dead and recovered) with SIR models where the coefficients $\beta$ and $\gamma$ are time-dependent:

$$\begin{array}{cc}\hfill \frac{dS}{dt}\left(t\right)& =-\frac{\beta \left(t\right)I\left(t\right)S\left(t\right)}{N}\hfill \\ \hfill \frac{dI}{dt}\left(t\right)& =\frac{\beta \left(t\right)I\left(t\right)S\left(t\right)}{N}-\gamma \left(t\right)I\left(t\right)\hfill \\ \hfill \frac{dR}{dt}\left(t\right)& =\gamma \left(t\right)I\left(t\right).\hfill \end{array}$$

In the following, we use bold-faced letters for past-time quantities. For example, $\mathbf{f}:=\left\{f\right(t):0\le t\le T\}$ for any function $\mathbf{f}\in L^{\infty }\left([0,T]\right)$. Using this notation, for any given $\mathbf{\beta }$ and $\mathbf{\gamma }\in L^{\infty }\left([0,T]\right)$ we denote by

$$(\mathbf{S},\mathbf{I},\mathbf{R})=\mathrm{SIR}(\mathbf{\beta },\mathbf{\gamma },[0,T\left]\right)$$

the solution of the associated SIR dynamics in $[0,T]$.

2.1.2. Detailed Compartmental Models

Models involving many compartments offer a detailed description of epidemiological mechanisms at the expense of involving many parameters. In our approach, we use them to generate virtual scenarios. One of the initial motivations behind the present work is to provide forecasts for the COVID-19 pandemic; thus, we have selected the two following models which are specific for this disease, but note that any other compartmental model [9,10,16] or agent-based simulation could also be used.

• First model, SEI5CHRD: This model is inspired by the one proposed in [10]. It involves 11 different compartments and a set of 19 parameters (see

  Table 1. Description of the compartments in Model SEI5CHRD.

  Compartment	Description
  S	Susceptible
  E	Exposed (non infectious)
  $I_p$	Infected and pre-symptomatic (already infectious)
  $I_a$	Infected and a-symptomatic (but infectious)
  $I_{ps}$	Infected and paucisymptomatic
  $I_{ms}$	Infected with mild symptoms
  $I_{ss}$	Infected with severe symptoms
  H	Hospitalized
  C	Intensive Care Unit
  R	Removed
  D	Dead

  ). The dynamics of the model are illustrated in Figure 1, and the system of equations reads as follows:

$$\begin{array}{cc}\hfill \frac{dS}{dt}\left(t\right)& =-\frac{1}{N}S\left(t\right)\left({\beta }_p{I}_p\left(t\right)+{\beta }_a{I}_a\left(t\right)+{\beta }_{ps}{I}_{ps}\left(t\right)+{\beta }_{ms}{I}_{ms}\left(t\right)+{\beta }_{ss}{I}_{ss}\left(t\right)+{\beta }_{H}H\left(t\right)+{\beta }_{C}C\left(t\right)\right)\hfill \\ \hfill \frac{dE}{dt}\left(t\right)& =\frac{1}{N}S\left(t\right)\left({\beta }_p{I}_p\left(t\right)+{\beta }_a{I}_a\left(t\right)+{\beta }_{ps}{I}_{ps}\left(t\right)+{\beta }_{ms}{I}_{ms}\left(t\right)+{\beta }_{ss}{I}_{ss}\left(t\right)+{\beta }_{H}H\left(t\right)+{\beta }_{C}C\left(t\right)\right)-\epsilon E\left(t\right)\hfill \\ \hfill \frac{d{I}_p}{dt}\left(t\right)& =\epsilon E\left(t\right)-{\mu }_p{I}_p\left(t\right)\hfill \\ \hfill \frac{d{I}_a}{dt}\left(t\right)& =p_a{\mu }_p{I}_p\left(t\right)-\mu I_a\left(t\right)\hfill \\ \hfill \frac{d{I}_{ps}}{dt}\left(t\right)& =p_{ps}(1-p_a){\mu }_p{I}_p\left(t\right)-\mu I_{ps}\left(t\right)\hfill \\ \hfill \frac{d{I}_{ms}}{dt}\left(t\right)& =p_{ms}(1-p_a){\mu }_p{I}_p\left(t\right)-\mu I_{ms}\left(t\right)\hfill \\ \hfill \frac{d{I}_{ss}}{dt}\left(t\right)& =p_{ss}(1-p_a){\mu }_p{I}_p\left(t\right)-\mu I_{ss}\left(t\right)\hfill \\ \hfill \frac{dC}{dt}\left(t\right)& =p_c\mu I_{ss}\left(t\right)-({\lambda }_{C,R}+{\lambda }_{C,D})C\left(t\right)\hfill \\ \hfill \frac{dH}{dt}\left(t\right)& =(1-p_c)\mu I_{ss}\left(t\right)-({\lambda }_{H,R}+{\lambda }_{H,D})H\left(t\right)\hfill \\ \hfill \frac{dR}{dt}\left(t\right)& ={\lambda }_{C,R}C\left(t\right)+{\lambda }_{H,R}H\left(t\right)\hfill \\ \hfill \frac{dD}{dt}\left(t\right)& ={\lambda }_{C,D}C\left(t\right)+{\lambda }_{H,D}H\left(t\right)\hfill \end{array}$$

The different parameters involved in the model are described in

Table 2. Description of the parameters involved in Model SEI5CHRD.

Parameter	Description
${\beta }_p$	Relative infectiousness of $I_p$
${\beta }_a$	Relative infectiousness of $I_a$
${\beta }_{ps}$	Relative infectiousness of $I_{ps}$
${\beta }_{ms}$	Relative infectiousness of $I_{ms}$
${\beta }_{ss}$	Relative infectiousness of $I_{ss}$
${\beta }_H$	Relative infectiousness of $I_H$
${\beta }_C$	Relative infectiousness of $I_C$
${\epsilon }^{-1}$	Latency period
${\mu }_p^{-1}$	Duration of prodromal phase
$p_a$	Probability of being asymptomatic
${\mu }^{-1}$	Infectious period of $I_a$, $I_{ps}$, $I_{ms}$, $I_{ss}$
$p_{ps}$	If symptomatic, probability of being paucisymptomatic
$p_{ms}$	If symptomatic, probability of developing mild symptoms
$p_{ss}$	If symptomatic, probability of developing severe symptoms (note that $p_{ps}+p_{ms}+p_{ss}=1$)
$p_C$	If severe symptoms, probability of going in C
${\lambda }_{CR}$	If in C, daily rate entering in R
${\lambda }_{CD}$	If in C, daily rate entering in D
${\lambda }_{HR}$	If hospitalized, daily rate entering in R
${\lambda }_{HD}$	If hospitalized, daily rate entering in D

and detailed in the appendix of [10].

We denote by

$$\begin{array}{cc}\hfill \mathbf{u}=\mathrm{SEI}5\mathrm{CHRD}({\mathbf{u}}_0,{\beta }_p,& {\beta }_a,{\beta }_{ps},{\beta }_{ms},{\beta }_{ss},{\beta }_H,{\beta }_C,\hfill \\ & \epsilon ,{\mu }_p,p_a,\mu ,p_{ps},p_{ms},p_{ss},p_C,\hfill \\ & {\lambda }_{CR},{\lambda }_{CD},{\lambda }_{HR},{\lambda }_{HD},[0,T])\hfill \end{array}$$

the parameter-to-solution map where $\mathbf{u}=(\mathbf{S},\mathbf{E},{\mathbf{I}}_{\mathbf{p}},{\mathbf{I}}_{\mathbf{a}},{\mathbf{I}}_{\mathbf{ps}},{\mathbf{I}}_{\mathbf{ms}},{\mathbf{I}}_{\mathbf{ss}},\mathbf{C},\mathbf{H},\mathbf{R},\mathbf{D})$.

• Second model, SE2IUR: This model is a variant of the one proposed in [9]. It involves five different compartments (see

  Table 3. Description of the compartments in model SE2IUR.

  Compartment	Description
  S	Susceptible
  $E_1$	Exposed (non infectious)
  $E_2$	Infected and pre-symptomatic (already infectious)
  I	Infected and symptomatic
  U	Un-noticed
  R	dead and removed

  ) and a set of six parameters. The dynamics of the model are illustrated in Figure 2 and the set of equations is as follows:

  $$\begin{array}{cc}\hfill \frac{dS}{dt}\left(t\right)& =-\frac{1}{N}\beta S\left(t\right)(E_2\left(t\right)+U\left(t\right)+I\left(t\right))\hfill \\ \hfill \frac{d{E}_1}{dt}\left(t\right)& =\frac{1}{N}\beta S\left(t\right)(E_2\left(t\right)+U\left(t\right)+I\left(t\right))-\delta E_1\left(t\right)\hfill \\ \hfill \frac{d{E}_2}{dt}\left(t\right)& =\delta E_1\left(t\right)-\sigma E_2\left(t\right)\hfill \\ \hfill \frac{dI}{dt}\left(t\right)& =\nu \sigma E_2\left(t\right)-{\gamma }_{1}I\left(t\right)\hfill \\ \hfill \frac{dU}{dt}\left(t\right)& =(1-\nu )\sigma E_2\left(t\right)-{\gamma }_{2}U\left(t\right)\hfill \\ \hfill \frac{dR}{dt}\left(t\right)& ={\gamma }_{1}I\left(t\right)+{\gamma }_{2}U\left(t\right)\hfill \end{array}$$
  We denote by

  $$\mathbf{u}=\mathrm{SE}2\mathrm{IUR}({\mathbf{u}}_0,\beta ,\delta ,\sigma ,\nu ,{\gamma }_1,{\gamma }_2,[0,T])$$

  the parameter-to-solution map where $\mathbf{u}=(\mathbf{S},\mathbf{E}\mathbf{1},\mathbf{E}\mathbf{2},\mathbf{I},\mathbf{U},\mathbf{R})$. The different parameters involved in the model are described in

  Table 4. Description of the parameters involved in model SE2IUR.

  Parameter	Description
  $\beta$	Relative infectiousness of I, U, $E_2$
  ${\delta }^{-1}$	Latency period
  ${\sigma }^{-1}$	Duration of prodromal phase
  $\nu$	Proportion of I among $I+U$
  ${\gamma }_1$	If I, daily rate entering in R
  ${\gamma }_2$	If U, daily rate entering in R

  .

• Generalization: In the following, we abstract the procedure as follows. For any Detailed_Model with d compartments involving a vector $\mu \in {\mathbb{R}}^p$ of p parameters, we denote by

  $$\mathbf{u}\left(\mu \right)=\mathrm{Detailed}\_\mathrm{Model}(\mu ,[0,\tilde{T}]),\mathbf{u}\left(\mu \right)\in L^{\infty }([0,\tilde{T}],{\mathbb{R}}^d)$$

  the parameter-to-solution map, where $\tilde{T}$ is some given time simulation that can be as large as desired because this is a virtual scenario. Note that, because the initial condition of the illness at time 0 is not known, we include the initial condition ${\mathbf{u}}_0$ in the parameter set.

2.2. Forecasting Based on Model Reduction of Detailed Models

We assume that we are given health data in a time window $[0,T]$, where $T>0$ is assumed to be the present time. The observed data are the series of infected people, denoted $I_{\mathrm{obs}}$, and removed people, denoted $R_{\mathrm{obs}}$. They are usually given at a national or a regional scale and on a daily basis. For our discussion, it is useful to work with time-continuous functions, and $t\to I_{\mathrm{obs}}\left(t\right)$ denotes the piecewise constant approximation in $[0,T]$ from the given data (and similarly for $R_{\mathrm{obs}}\left(t\right)$). Our goal is to give short-term forecasts of the series in a time window $\tau >0$ whose size is about two weeks. We denote by $I\left(t\right)$ and $R\left(t\right)$ the approximations to the series $I_{\mathrm{obs}}\left(t\right)$ and $R_{\mathrm{obs}}\left(t\right)$ at any time $t\in [0,T]$.

As mentioned above, we propose to fit the data and provide forecasts with SIR models with time-dependent parameters $\mathbf{\beta }$ and $\mathbf{\gamma }$. The main motivation for using such a simple family is that it possesses optimal fitting and forecasting properties for our purposes, as explained above. We define the cost function

$$\mathcal{J}(\mathbf{\beta },\mathbf{\gamma }|I_{\mathrm{obs}}\left(t\right),R_{\mathrm{obs}}\left(t\right),[0,T]):={\int }_0^T\left(|I\left(t\right)-I_{\mathrm{obs}}{\left(t\right)|}^2+{|R\left(t\right)-R_{\mathrm{obs}}\left(t\right)|}^2\right)\mathrm{d}t$$

(1)

such that

$$(\mathbf{S},\mathbf{I},\mathbf{R})=\mathrm{SIR}(\mathbf{\beta },\mathbf{\gamma },[0,T\left]\right),$$

and the fitting problem can be expressed at the continuous level as the optimal control problem of finding

$$J^{*}=\underset{(\mathbf{\beta },\mathbf{\gamma })\in L^{\infty }\left([0,T]\right)\times L^{\infty }\left([0,T]\right)}{inf}\mathcal{J}(\mathbf{\beta },\mathbf{\gamma }|{\mathbf{I}}_{\mathrm{obs}},{\mathbf{R}}_{\mathrm{obs}},[0,T]).$$

(2)

The following result ensures the existence of a unique minimizer under very mild constraints.

Proposition 1.

Let$N\in {\mathbb{N}}^{*}$and$T>0$. For any real-valued functions$S_{obs},I_{obs},R_{obs}$of class${\mathcal{C}}^1$, defined on$[0,T]$satisfying

(i) 

$S_{obs}\left(t\right)+I_{obs}\left(t\right)+R_{obs}\left(t\right)=N$for every$t\in [0,T]$,

(ii) 

$S_{obs}$in nonincreasing on$[0,T]$,

(iii) 

$R_{obs}$is nondeacreasing on$[0,T]$,

there exists a unique minimizer$({\mathbf{\beta }}_{obs}^{*},{\mathbf{\gamma }}_{obs}^{*})$to Equation (2).

Proof. 

One can set

$$\begin{array}{c}\hfill \left\{\begin{array}{cc}{\beta }_{\mathrm{obs}}^{*}\left(t\right)\hfill & :=-{\displaystyle \frac{N}{I_{\mathrm{obs}}\left(t\right)S_{\mathrm{obs}}\left(t\right)}}{\displaystyle \frac{d{S}_{\mathrm{obs}}}{dt}}\left(t\right)\hfill \\ {\gamma }_{\mathrm{obs}}^{*}\left(t\right)\hfill & :={\displaystyle \frac{1}{I_{\mathrm{obs}}\left(t\right)}}{\displaystyle \frac{d{R}_{\mathrm{obs}}}{dt}}\left(t\right)\hfill \end{array}\right.\end{array}$$

(3)

so that

$$({\mathbf{S}}_{\mathrm{obs}},{\mathbf{I}}_{\mathrm{obs}},{\mathbf{R}}_{\mathrm{obs}})=\mathrm{SIR}({\mathbf{\beta }}^{*},{\mathbf{\gamma }}^{*},[0,T])$$

and

$$\mathcal{J}({\mathbf{\beta }}_{\mathrm{obs}}^{*},{\mathbf{\gamma }}_{\mathrm{obs}}^{*},[0,T])=0$$

which obviously implies that $J^{*}=0$.    □

Note that one can equivalently set

$$\begin{array}{c}\hfill \left\{\begin{array}{cc}{\beta }_{\mathrm{obs}}^{*}\left(t\right)\hfill & :=-{\displaystyle \frac{N}{I_{\mathrm{obs}}\left(t\right)S_{\mathrm{obs}}\left(t\right)}}{\displaystyle \frac{d{S}_{\mathrm{obs}}}{dt}}\left(t\right)\hfill \\ {\gamma }_{\mathrm{obs}}^{*}\left(t\right)\hfill & :={\displaystyle \frac{1}{I_{\mathrm{obs}}\left(t\right)}}\left[{\displaystyle \frac{d{I}_{\mathrm{obs}}}{dt}}\left(t\right)-{\displaystyle \frac{{\beta }_{\mathrm{obs}}^{*}\left(t\right)I_{\mathrm{obs}}\left(t\right)S_{\mathrm{obs}}\left(t\right)}{N}}\right]\hfill \end{array}\right.\end{array}$$

or again

$$\begin{array}{c}\hfill \left\{\begin{array}{cc}{\gamma }_{\mathrm{obs}}^{*}\left(t\right)\hfill & :={\displaystyle \frac{1}{I_{\mathrm{obs}}\left(t\right)}}{\displaystyle \frac{d{R}_{\mathrm{obs}}}{dt}}\left(t\right)\hfill \\ {\beta }_{\mathrm{obs}}^{*}\left(t\right)\hfill & :={\displaystyle \frac{N}{I_{\mathrm{obs}}\left(t\right)S_{\mathrm{obs}}\left(t\right)}}\left[{\displaystyle \frac{d{I}_{\mathrm{obs}}}{dt}}\left(t\right)-{\gamma }_{\mathrm{obs}}^{*}\left(t\right)I_{\mathrm{obs}}\left(t\right)\right]\hfill \end{array}\right.\end{array}$$

This simple observation means that there exists a time-dependent SIR model which can perfectly fit the data of any (epidemiological) evolution that satisfies properties (i), (ii), and (iii). In particular, we can perfectly fit the series of sick people with a time-dependent SIR model (with a smoothing of the local oscillations due to noise). Since the health data are usually given on a daily basis, we can approximate ${\mathbf{\beta }}_{\mathrm{obs}}^{*},{\mathbf{\gamma }}_{\mathrm{obs}}^{*}$ by approximating the derivatives by classical finite differences in Equation (3).

The fact that we can build ${\mathbf{\beta }}_{\mathrm{obs}}^{*}$ and ${\mathbf{\gamma }}_{\mathrm{obs}}^{*}$ such that $\mathcal{J}({\mathbf{\beta }}_{\mathrm{obs}}^{*},{\mathbf{\gamma }}_{\mathrm{obs}}^{*})=J^{*}=0$ implies that the family of time-dependent SIR models is rich enough not only to fit the evolution of any epidemiological series but also to deliver perfect predictions of the health data. It is however important to note that since the ${\mathbf{\beta }}_{\mathrm{obs}}^{*},{\mathbf{\gamma }}_{\mathrm{obs}}^{*}$ are derived exclusively from the data and depend on time, we lose the direct interpretations of these coefficients in terms of the length of the contagious period or the transmission rate that these coefficients present when they are considered constant. The great approximation power comes also at the cost of defining the parameters $\mathbf{\beta }$ and $\mathbf{\gamma }$ in $L^{\infty }\left([0,T]\right)$ which is a space that is too large to be able to define any feasible prediction strategy.

In order to pin down a smaller manifold where these parameters may vary without sacrificing much in terms of the fitting and forecasting power, our strategy is the following:

• Learning phase: The fundamental hypothesis of our approach is that the specialists of epidemiology have understood the mechanisms of infection transmission sufficiently well. The second hypothesis is that these accurate models suffer from the proper determination of the parameters they contain. We thus propose to generate a large number of virtual epidemics, following these mechanistic models, with parameters that can be chosen in the vicinity of the suggested parameter values, including the different initial conditions.

  (a)

  Generate virtual scenarios using detailed models with constant coefficients:

  • Define the notion of a Detailed_Model which is most appropriate for the epidemiological study. Several models could be considered. For our numerical application, the detailed models are defined in Section 2.1.
  • Define an interval range $\mathcal{P}\subset {\mathbb{R}}^p$ where the parameters $\mu$ of etailed_Model vary. We call the solution manifold $\mathcal{U}$ the set of virtual dynamics over $[0,T+\tau ]$, namely

    $$\mathcal{U}:=\left\{\mathbf{u}\right(\mu )=\mathrm{Detailed}\_\mathrm{Model}(\mu ,[0,T+\tau ]):\mu \in \mathcal{P}\}.$$
  • Draw a finite training set

    $${\mathcal{P}}_{\mathrm{tr}}=\{{\mu }_1,\dots ,{\mu }_K\}\subseteq \mathcal{P}$$

    of $K\gg 1$ parameter instances and compute $\mathbf{u}\left(\mu \right)=\mathrm{Detailed}\_\mathrm{Model}$$(\mu ,[0,T+\tau \left]\right)$ for $\mu \in {\mathcal{P}}_{\mathrm{tr}}$. Each $\mathbf{u}\left(\mu \right)$ is a virtual epidemiological scenario. An important detail for our prediction purposes is that the simulations are done in $[0,T+\tau ]$; that is, we simulate not only in the fitting time interval but also in the prediction time interval. We call

    $${\mathcal{U}}_{\mathrm{tr}}=\{\mathbf{u}\left(\mu \right):\mu \in {\mathcal{P}}_{\mathrm{tr}}\}$$

    the training set of all virtual scenarios.

  (b)

  Collapse every detailed model $\mathbf{u}\left(\mu \right)\in {\mathcal{U}}_{\mathrm{tr}}$ into an SIR model following the ideas explained in Section 2.3. For every $\mathbf{u}\left(\mu \right)$, the procedure gives time-dependent parameters $\mathbf{\beta }\left(\mu \right)$ and $\mathbf{\gamma }\left(\mu \right)$ and associated SIR solutions $(\mathbf{S},\mathbf{I},\mathbf{R})\left(\mu \right)$ in $[0,T+\tau ]$. This yields the sets

  $${\mathcal{B}}_{\mathrm{tr}}:=\{\mathbf{\beta }\left(\mu \right):\mu \in {\mathcal{P}}_{\mathrm{tr}}\}\mathrm{and}{\mathcal{G}}_{\mathrm{tr}}\{\mathbf{\gamma }\left(\mu \right):\mu \in {\mathcal{P}}_{\mathrm{tr}}\}.$$

  (4)

  (c)

  Compute reduced models:

  We apply model reduction techniques using ${\mathcal{B}}_{\mathrm{tr}}$ and ${\mathcal{G}}_{\mathrm{tr}}$ as training sets in order to build two bases

  $${\mathrm{B}}_n=\mathrm{span}\{b_1,\dots ,b_n\},{\mathrm{G}}_n=\mathrm{span}\{g_1,\dots ,g_n\}\subset L^{\dots }([0,T+\tau ],\mathbb{R}),$$

  which are defined over $[0,T+\tau ]$. The space ${\mathrm{B}}_n$ is such that it approximates all functions $\mathbf{\beta }\left(\mu \right)\in {\mathcal{B}}_{\mathrm{tr}}$ well (or all $\mathbf{\gamma }\left(\mu \right)\in {\mathcal{G}}_{\mathrm{tr}}$ can be well approximated by elements of ${\mathrm{G}}_n$). In Section 2.4, we outline the methods we have explored in our numerical tests.

• Fitting on the reduced spaces: We next solve the fitting problem (2) in the interval $[0,T]$ by searching $\mathbf{\beta }$ (or $\mathbf{\gamma }$) in ${\mathrm{B}}_n$ (or in ${\mathrm{G}}_n$) instead of in $L^{\infty }\left([0,T]\right)$; that is,

  $$J_{({\mathrm{B}}_n,{\mathrm{G}}_n)}^{*}=\underset{(\mathbf{\beta },\mathbf{\gamma })\in {\mathrm{B}}_n\times {\mathrm{G}}_n}{min}\mathcal{J}(\mathbf{\beta },\mathbf{\gamma }|{\mathbf{I}}_{\mathrm{obs}},{\mathbf{R}}_{\mathrm{obs}},[0,T]).$$

  (5)
  Note that the respective dimensions of ${\mathrm{B}}_n$ and ${\mathrm{G}}_n$ can be different; for simplicity, we consider them to be equal in the following. Obviously, since ${\mathrm{B}}_n$ and ${\mathrm{G}}_n\subset L^{\infty }\left([0,T]\right)$, we obtain

  $$J^{*}\le J_{({\mathrm{B}}_n,{\mathrm{G}}_n)}^{*},$$

  but we numerically observe that the function $n\mapsto J_{({\mathrm{B}}_n,{\mathrm{G}}_n)}^{*}$ decreases very rapidly as n increases, which indirectly confirms the fact that the manifold generated by the two above models accommodates the COVID-19 epidemic well.
  The solution of problem (5) gives us the coefficients ${\left(c_i^{*}\right)}_{i=1}^n$ and ${\left({\tilde{c}}_i^{*}\right)}_{i=1}^n\in {\mathbb{R}}^n$ such that the time-dependent parameters

  $$\begin{array}{cc}\hfill {\beta }_n^{*}\left(t\right)& =\sum _{i=1}^n{c}_i^{*}{b}_i\left(t\right),\forall t\in [0,T+\tau ],\hfill \\ \hfill {\gamma }_n^{*}\left(t\right)& =\sum _{i=1}^n{\tilde{c}}_i^{*}{g}_i\left(t\right).\hfill \end{array}$$

  achieve the minimum (5).
• Forecast: For a given dimension n of the reduced spaces, we propagate in $[0,T+\tau ]$ the associated SIR model, as follows:

  $$({\mathbf{S}}_n^{*},{\mathbf{I}}_n^{*},{\mathbf{R}}_n^{*})=\mathrm{SIR}({\mathbf{\beta }}_n^{*},{\mathbf{\gamma }}_n^{*},[0,T+\tau ])$$
  The values $I_n^{*}\left(t\right)$ and $R_n^{*}\left(t\right)$ for $t\in [0,T[$ are by construction close to the observed data ${\mathbf{I}}_{\mathrm{obs}},{\mathbf{R}}_{\mathrm{obs}}$ (up to some numerical optimization error). The values $I_n^{*}\left(t\right)$ and $R_n^{*}\left(t\right)$ for $t\in [T,T+\tau ]$ are then used for prediction.
• Forecast combination/aggregation of experts (optional step): By varying the dimension n and using different model reduction approaches, we can easily produce a collection of different forecasts, and thus the question of how to select the best predictive model arises. Alternatively, we can also resort to forecast combination techniques [17]: denoting $(I_1,R_1),\dots ,(I_P,R_P)$ as the different forecasts, the idea is to search for an appropriate linear combination

  $$I^{\mathrm{FC}}\left(t\right)=\sum _{p=1}^P{w}_p{I}_p\left(t\right)$$

  and perform a similar operation for R. Note that these combinations do not need to involve forecasts from our methodology only; other approaches such as time series forecasts could also be included. One simple forecast combination is the average in which all alternative forecasts are given the same weight $w_p=1/P,p=1,\dots ,P$. More elaborate approaches consist in estimating the weights that minimize a loss function involving the forecast error.

Before going into detail regarding some of the steps, three points should be highlighted:

• To bring out the essential mechanisms, we have idealized some elements in the above discussion by omitting certain unavoidable discretization aspects. To start with, the ODE solutions cannot be computed exactly but only up to a certain level of accuracy given by a numerical integration scheme. In addition, the optimal control problems (2) and (5) are non-convex. As a result, in practice, we can only find a local minimum. Note, however, that modern solvers find solutions which are very satisfactory for all practical purposes. In addition, note that solving the control problem in a reduced space as in (5) could be interpreted as introducing a regularizing effect with respect to the control problem (2) in the full $L^{\infty }\left([0,T]\right)$ space. It is to be expected that the search of global minimizers is facilitated in the reduced landscape.
• routine-IR and routine-$\mathbf{\beta }$$\mathbf{\gamma }$: A variant for the fitting problem (5) as studied in our numerical experiments is to replace the cost function $\mathcal{J}(\mathbf{\beta },\mathbf{\gamma }|{\mathbf{I}}_{\mathrm{obs}},{\mathbf{R}}_{\mathrm{obs}},[0,T])$ by the cost function

  $$\tilde{\mathcal{J}}(\mathbf{\beta },\mathbf{\gamma }|{\mathbf{\beta }}_{\mathrm{obs}}^{*},{\mathbf{\gamma }}_{\mathrm{obs}}^{*},[0,T]):={\int }_0^T\left(|\mathbf{\beta }-{\mathbf{\beta }}_{\mathrm{obs}}^{*}{|}^2+|\mathbf{\gamma }-{\mathbf{\gamma }}_{\mathrm{obs}}^{*}{\left)\right|}^2\right)\mathrm{d}t.$$

  (6)
  In other words, we use the variables ${\mathbf{\beta }}_{\mathrm{obs}}^{*}$ and ${\mathbf{\gamma }}_{\mathrm{obs}}^{*}$ from (3) as observed data instead of working with the observed health series ${\mathbf{I}}_{\mathrm{obs}}$, ${\mathbf{R}}_{\mathrm{obs}}$. In Section 4, we refer to the standard fitting method as routine-IR and to this variant as routine-$\mathbf{\beta }$$\mathbf{\gamma }$.
• The fitting procedure works both on the components of the reduced basis and the initial time of the epidemics to minimize the loss function; however, for simplicity, this last optimization is not reported here.

2.3. Details on Step 1-(b): Collapsing the Detailed Models into SIR Dynamics

Let

$$\mathbf{u}\left(\mu \right)=\mathrm{Detailed}\_\mathrm{Model}(\mu ,[0,T+\tau ])\in L^{\infty }([0,T+\tau ],{\mathbb{R}}^d)$$

be the solution in $[0,T+\tau ]$ to a detailed model for a given vector of parameters $\mu \in {\mathbb{R}}^d$. Here, d is possibly large ($d=11$ in the case of the SEI5CHRD model and $d=5$ in the case of SE2IUR’s model). The goal of this section is to explain how to collapse the detailed dynamics $\mathbf{u}\left(\mu \right)$ into SIRdynamics with time-dependent parameters. The procedure can be understood as a projection of a detailed dynamics into the family of dynamics given by SIR models with time-dependent parameters.

For the SEI5CHRD model, we collapse the variables by setting

$$\begin{array}{cc}\hfill {\mathbf{S}}^{\mathrm{col}}& =\mathbf{S}+\mathbf{E}\hfill \\ \hfill {\mathbf{I}}^{\mathrm{col}}& ={\mathbf{I}}_{\mathbf{p}}+{\mathbf{I}}_{\mathbf{a}}+{\mathbf{I}}_{\mathbf{ps}}+{\mathbf{I}}_{\mathbf{ms}}+{\mathbf{I}}_{\mathbf{ss}}+\mathbf{C}+\mathbf{H}\hfill \\ \hfill {\mathbf{R}}^{\mathrm{col}}& =\mathbf{R}+\mathbf{D}\hfill \end{array}$$

Similarly, for the SE2IUR model, we set

$$\begin{array}{cc}\hfill {\mathbf{S}}^{\mathrm{col}}& =\mathbf{S}+{\mathbf{E}}_{1i}\hfill \\ \hfill {\mathbf{I}}^{\mathrm{col}}& ={\mathbf{E}}_{2i}+{\mathbf{I}}_i+{\mathbf{U}}_i\hfill \\ \hfill {\mathbf{R}}^{\mathrm{col}}& =\mathbf{R}\hfill \end{array}$$

Note that ${\mathbf{S}}^{\mathrm{col}}$, ${\mathbf{I}}^{\mathrm{col}}$ and ${\mathbf{R}}^{\mathrm{col}}$ depend on $\mu$, but we have omitted this dependence to simplify the notation.

Once the collapsed variables are obtained, we identify the time-dependent parameters $\mathbf{\beta }$ and $\mathbf{\gamma }$ of the SIR model by solving the fitting problem

$$(\mathbf{\beta },\mathbf{\gamma })\in \underset{(\mathbf{\beta },\mathbf{\gamma })\in L^{\infty }([0,T+\tau ],\mathbb{R})\times L^{\infty }([0,T+\tau ],\mathbb{R})}{arg\; inf}\mathcal{J}(\mathbf{\beta },\mathbf{\gamma }|{\mathbf{I}}^{\mathrm{col}},{\mathbf{R}}^{\mathrm{col}},[0,T+\tau ])$$

(7)

where

$$(\mathbf{S},\mathbf{I},\mathbf{R})=\mathrm{SIR}(\mathbf{\beta },\mathbf{\gamma },[0,T+\tau \left]\right).$$

Note that problem (7) has the same structure as problem (2), with the difference arising from the fact that the collapsed variables ${\mathbf{I}}^{\mathrm{col}}$, ${\mathbf{R}}^{\mathrm{col}}$ in (7) play the role of the health data ${\mathbf{I}}_{\mathrm{obs}}$, ${\mathbf{R}}_{\mathrm{obs}}$ in (2). Therefore, it follows from Proposition 1 that problem (7) has a very simple solution as it suffices to apply formula (3) to solve it. Note here that the exact derivatives of ${\mathbf{S}}^{\mathrm{col}}$, ${\mathbf{I}}^{\mathrm{col}}$, and ${\mathbf{R}}^{\mathrm{col}}$ can be obtained from the Detailed_Model.

Since the solution $(\mathbf{\beta },\mathbf{\gamma })$ to (7) depends on the parameter $\mu$ of the detailed model, repeating the above procedure for every detailed scenario $\mathbf{u}\left(\mu \right)$ for any $\mu \in {\mathcal{P}}_{\mathrm{tr}}$ yields the two families of time-dependent functions ${\mathcal{B}}_{\mathrm{tr}}=\{\mathbf{\beta }\left(\mu \right):\mu \in {\mathcal{P}}_{\mathrm{tr}}\}$ and ${\mathcal{G}}_{\mathrm{tr}}=\{\mathbf{\gamma }\left(\mu \right):\mu \in {\mathcal{P}}_{\mathrm{tr}}\}$ defined in the interval $[0,T+\tau ]$ as introduced in Section (4).

2.4. Details of Model Order Reduction

Model order reduction is a family of methods aiming at the approximation of a set of solutions of parametrized PDEs or ODEs (or related quantities) with linear spaces, which are called reduced models or reduced spaces. In our case, the sets to approximate are

$$\mathcal{B}=\left\{\mathbf{\beta }\right(\mu ):\mu \in \mathcal{P}\}\mathrm{and}\mathcal{G}=\left\{\mathbf{\gamma }\right(\mu ):\mu \in \mathcal{P}\},$$

where each $\mu$ is the vector of parameters of the detailed model which take values over $\mathcal{P}$, and $\mathbf{\beta }\left(\mu \right)$ and $\mathbf{\gamma }\left(\mu \right)$ are the associated time-dependent coefficients of the collapsed SIR evolution. In the following, we view $\mathcal{B}$ and $\mathcal{G}$ as subsets of $L^2\left([0,T]\right)$, and we denote by $\parallel ·\parallel$ and $〈·,·〉$ its norm and inner product. Indeed, in view of Proposition 1, $\mathcal{B}$ and $\mathcal{G}\subset L^{\infty }\left([0,T]\right)\subset L^2\left([0,T]\right)$.

Continuing the discussion if $\mathcal{B}$ (the same will hold for $\mathcal{G}$), of we measure performance in terms of the worst error in the set $\mathcal{B}$, the best possible performance that reduced models of dimension n can achieve is given by the Kolmogorov n-width:

$$d_n{\left(\mathcal{B}\right)}_{L^2\left([0,T]\right)}:=\underset{\begin{array}{c}Y\in L^2\left([0,T]\right)\\ dim\left(Y\right)\le n\end{array}}{inf}\underset{u\in \mathcal{B}}{max}\parallel u-P_{Y}u\parallel$$

where $P_Y$ is the orthogonal projection onto the subspace Y. In the case of measuring errors in an average sense, the benchmark is given by

$${\delta }_n^2{(\mathcal{B},\nu )}_{L^2\left([0,T]\right)}:=\underset{\begin{array}{c}Y\in L^2\left([0,T]\right)\\ dim\left(Y\right)=n\end{array}}{inf}{\int }_{\mathcal{P}}\parallel u\left(y\right)-P_Y{u\left(y\right)\parallel }^2\mathrm{d}\nu \left(y\right)$$

where $\nu$ is some given measure on $\mathcal{P}$.

In practice, building spaces that meet these benchmarks is generally not possible. However, it is possible to build sequences of spaces for which the error decay is comparable to that given by ${\left(d_n{\left(\mathcal{B}\right)}_{L^2\left([0,T]\right)}\right)}_n$ or ${\left({\delta }_n{\left(\mathcal{B}\right)}_{L^2\left([0,T]\right)}\right)}_n$. As a result, when the Kolmogorov width decays quickly, the constructed reduced spaces will deliver a very good approximation of the set $\mathcal{B}$ with few modes (see [18,19,20,21]).

We next present the reduced models used in our numerical experiments. Other methods could, of course, be considered, and we refer readers to [22,23,24,25] for general references on model reduction. We continue the discussion in a fully discrete setting in order to simplify the presentation and keep it as close to the practical implementation as possible. All the claims below could be written in a fully continuous sense at the expense of introducing additional mathematical objects such as certain Hilbert–Schmidt operators to define the continuous version of Singular Value Decomposition (SVD).

We build the reduced models using the two discrete training sets of functions ${\mathcal{B}}_{\mathrm{tr}}={\left\{{\mathbf{\beta }}_i\right\}}_{i=1}^K$ and ${\mathcal{G}}_{\mathrm{tr}}={\left\{{\mathbf{\gamma }}_i\right\}}_{i=1}^K$ from $\mathcal{B}$ and $\mathcal{G}$, where K denotes the number of virtual scenarios considered. The sets have been generated in step 1-(b) of our general pipeline (see Section 2.2).

We consider a discretization of the time interval $[0,T+\tau ]$ into a set of $Q\in {\mathbb{N}}^{*}$ points as follows: $\{t_1=0,\cdots ,t_P=T,\cdots ,t_Q=T+\tau \}$ where $P<Q$. Thus, we can represent each function ${\mathbf{\beta }}_i$ as a vector of Q values

$${\mathbf{\beta }}_i={({\beta }_i\left(t_1\right),\cdots ,{\beta }_i\left(t_Q\right))}^T\in {\mathbb{R}}_{+}^Q.$$

and thus assemble all the functions of the family ${\left\{{\mathbf{\beta }}_i\right\}}_{i=1}^K$ into a matrix $M_{\mathcal{B}}\in {\mathbb{R}}_{+}^{Q\times K}$. The same remark applies for the family ${\left\{{\mathbf{\gamma }}_i\right\}}_{i=1}^K$ which gives a matrix $M_{\mathcal{G}}\in {\mathbb{R}}_{+}^{Q\times K}$.

• SVD: The eigenvalue decomposition of the correlation matrix $M_{\mathcal{B}}^T{M}_{\mathcal{B}}\in {\mathbb{R}}^{K\times K}$ gives

  $$M_{\mathcal{B}}^T{M}_{\mathcal{B}}=V\mathsf{\Lambda }{V}^T,$$

  where $V=\left(v_{i,j}\right)\in {\mathbb{R}}^{K\times K}$ is an orthogonal matrix and $\mathsf{\Lambda }\in {\mathbb{R}}^{K\times K}$ is a diagonal matrix with non-negative entries, which we denote as ${\lambda }_i$ and present in decreasing order. The ${\ell }^2\left({\mathbb{R}}^Q\right)$-orthogonal basis functions $\{{\mathbf{b}}_1,\dots ,{\mathbf{b}}_K\}$ are then given by the linear combinations

  $${\mathbf{b}}_i=\sum _{j=1}^K{v}_{j,i}{\mathbf{\beta }}_j,1\le i\le K.$$
  For $n\le K$, the space

  $${\mathrm{B}}_n=\mathrm{span}\{{\mathbf{b}}_1,\dots {\mathbf{b}}_n\}$$

  is the best n-dimensional space to approximate the set ${\left\{{\mathbf{\beta }}_i\right\}}_{i=1}^K$ in the average sense. We have

  $${\delta }_n{\left({\left\{{\mathbf{\beta }}_i\right\}}_{i=1}^K\right)}_{{\ell }^2\left({\mathbb{R}}^{Q+1}\right)}={\left(\frac{1}{K}\sum _{i=1}^K\parallel {\mathbf{\beta }}_i-P_{{\mathrm{B}}_n}{\mathbf{\beta }}_i{\parallel }_{{\ell }^2\left({\mathbb{R}}^{Q+1}\right)}^2\right)}^{1/2}={\left(\sum _{i>n}^K{\lambda }_i\right)}^{1/2}$$

  and the average approximation error is given by the sum of the tail of the eigenvalues.
  Therefore, the SVD method is particularly efficient if there is a fast decay of the eigenvalues, meaning that the set ${\mathcal{B}}_{\mathrm{tr}}={\left\{{\mathbf{\beta }}_i\right\}}_{i=1}^K$ can be approximated by only few modes. However, note that, by construction, this method does not ensure positivity in the sense that $P_{{\mathrm{B}}_n}{\beta }_i\left(t\right)$ may become negative for some $t\in [0,T]$ although the original function ${\mathbf{\beta }}_i\left(t\right)\ge 0$ for all $t\in [0,T]$. This is due to the fact that the vectors ${\mathbf{b}}_i$ are not necessarily nonnegative. As we will see later, in our study, ensuring positivity especially for extrapolation (i.e., forecasting) is particularly important and motivates the next methods.
• Nonnegative Matrix Factorization (NMF, see [26,27]): NMF is a variant of SVD involving nonnegative modes and expansion coefficients. In this approach, we build a family of non-negative functions ${\left\{{\mathbf{b}}_i\right\}}_{i=1}^n$ and we approximate each ${\mathbf{\beta }}_i$ with a linear combination

  $${\mathbf{\beta }}_i^{\mathrm{NMF}}=\sum _{j=1}^n{a}_{i,j}{\mathbf{b}}_j,1\le i\le K,$$

  (8)

  where for every $1\le i\le K$ and $1\le j\le n$, the coefficients $a_{i,j}\ge 0$ and the basis function ${\mathbf{b}}_j\ge 0$. In other words, we solve the following constrained optimization problem:

  $$(W^{*},B^{*})\in \underset{(W,B)\in {\mathbb{R}}_{+}^{K\times n}\times {\mathbb{R}}_{+}^{n\times Q}}{arg\; min}{\parallel M_{\mathcal{B}}-WB\parallel }_F^2.$$
  We refer readers to [27] for further details on the NMF and its numerical aspects.
• The Enlarged Nonnegative Greedy (ENG) algorithm with projection on an extended cone of positive functions: We now present our novel model reduction method, which is of interest in itself as it allows reduced models that preserve positivity and even other types of bounds to be built. The method stems from the observation that NMF approximates functions in the cone of positive functions of $\mathrm{span}{\{{\mathbf{b}}_i\ge 0\}}_{i=1}^n$ since it imposes $a_{i,j}\ge 0$ in the linear combination (8). However, note that the positivity of the linear combination is not equivalent to the positivity of the coefficients $a_{i,j}$ since there are obviously linear combinations involving very small $a_{i,j}<0$ for some j which may still deliver a nonnegative linear combination ${\sum }_{j=1}^n{a}_{i,j}{\mathbf{b}}_j$. We can thus widen the cone of linear combinations yielding positive values by carefully including these negative coefficients $a_{i,j}$. One possible strategy for this is proposed in Algorithm 1, which describes a routine that we call Enlarge_Cone. The routine

  $$\{{\mathbf{\psi }}_1,\dots ,{\mathbf{\psi }}_n\}=\mathrm{Enlarge}\_\mathrm{Cone}[\{{\mathbf{b}}_1,\dots ,{\mathbf{b}}_n\},\epsilon ]$$

  takes a set of nonnegative functions $\{{\mathbf{b}}_1,\dots ,{\mathbf{b}}_n\}$ as an input and modifies each function ${\mathbf{b}}_i$ by iteratively adding negative linear combinations of the other basis functions ${\mathbf{b}}_j$ for $j\ne i$ (see line 11 of the routine). The coefficients are chosen in an optimal way so as to maintain the positivity of the final linear combination while minimizing the $L^{\infty }$-norm. The algorithm returns a set of functions

  $${\mathbf{\psi }}_i={\mathbf{b}}_i-\sum _{j\ne i}{\sigma }_j^i{\mathbf{b}}_j,\forall i\in \{1,\dots ,n\}$$

  with ${\sigma }_j^i\ge 0$. Note that the algorithm requires the setting of a tolerance parameter $\epsilon >0$ for the computation of the ${\sigma }_j^i$.
  Once we have run Enlarge_Cone, the approximation of any function $\mathbf{\beta }$ is then sought as

  $${\mathbf{\beta }}^{\left(EC\right)}=\underset{c_1,\dots ,c_n\ge 0}{arg\; min}{\parallel \mathbf{\beta }-\sum _{j=1}^n{c}_j{\mathbf{\psi }}_j\parallel }_{L^2\left([0,T+\tau ]\right)}^2$$
  We emphasize that the routine is valid for any set of nonnegative input functions. We can thus apply Enlarge_Cone to the functions ${\{{\mathbf{b}}_i\ge 0\}}_{i=1}^n$ from NMF but also to the functions selected by a greedy algorithm such as the following:
  • For $n=1$, find

    $${\mathbf{b}}_1=\underset{\mathbf{\beta }\in {\left\{{\mathbf{\beta }}_i\right\}}_{i=1}^K}{arg\; max}{\parallel \mathbf{\beta }\parallel }_{L^2\left([0,T+\tau ]\right)}^2$$
  • At step $n>1$, we have selected the set of functions $\{{\mathbf{b}}_1,\dots ,{\mathbf{b}}_{n-1}\}$. We next find

    $${\mathbf{b}}_n=\underset{\mathbf{\beta }\in {\left\{{\mathbf{\beta }}_i\right\}}_{i=1}^K}{arg\; max}\underset{c_1,\dots ,c_n\ge 0}{min}{\parallel \mathbf{\beta }-\sum _{j=1}^{n-1}{c}_j{\mathbf{b}}_j\parallel }_{L^2\left([0,T+\tau ]\right)}^2$$
  In our numerical tests, we call the Enlarged Nonnegative Greedy (ENG) method the routine involving the above greedy algorithm combined with our routine Enlarge_Cone.

  Algorithm 1 Enlarge_Cone$[\{{\mathbf{b}}_1,\dots ,{\mathbf{b}}_n\},\epsilon ]\to \{{\mathbf{\psi }}_1,\dots ,{\mathbf{\psi }}_n\}$.

  Input: Set of nonnegative functions $\{{\mathbf{b}}_1,\dots ,{\mathbf{b}}_n\}$. Tolerance $\epsilon >0$.   for $i\in \{1,\dots ,n\}$ do     Set ${\sigma }_j^i=0,\forall j\ne i.$     for $\ell \in \{1,\dots ,n\}$ do       ${\alpha }_{\ell }^{i,*}=arg\; max\{\alpha \ge 0:\left[{\mathbf{b}}_i-{\sum }_{j\ne i}{\sigma }_j^i{\mathbf{b}}_j-\alpha {\mathbf{b}}_{\ell }\right]\left(t\right)>0,\forall t\in [0,T+\tau ]\}$       ${\sigma }_{\ell }^i={\sigma }_{\ell }^i+\frac{{\alpha }_{\ell }^{i,*}}{2}$        while ${\alpha }_{\ell }^{i,*}\ge \epsilon$ do         ${\alpha }_{\ell }^{i,*}=arg\; max\{\alpha \ge 0:\left[{\mathbf{b}}_i-{\sum }_{j\ne i}{\sigma }_j^i{\mathbf{b}}_j-\alpha {\mathbf{b}}_{\ell }\right]\left(t\right)>0,\forall t\in [0,T+\tau ]\}$         ${\sigma }_{\ell }^i={\sigma }_{\ell }^i+\frac{{\alpha }_{\ell }^{i,*}}{2}$        end while      end for     ${\mathbf{\psi }}_i={\mathbf{b}}_i-{\sum }_{j\ne i}{\sigma }_j^i{\mathbf{b}}_j$   end for Output:$\{{\mathbf{\psi }}_1,\dots ,{\mathbf{\psi }}_n\}$

  We remark that, if we work with positive functions that are upper bounded by a constant $L>0$, we can ensure that the approximations, denoted as $\Psi$, and written as a linear combination of basis functions, will also be between these bounds 0 and L by defining on the one hand, and as we have just done, a cone of positive functions generated by the above family ${\left\{{\psi }_i\right\}}_i$, and on the other hand, considering the base of the functions $L-\phi$, $\phi$ to be above the set all greedy elements of the reduced basis, to which we also apply the enlargement of these positive functions. We then require the approximation to be written as a positive combination of the first (positive) functions and for $L-\mathsf{\Psi }$ to also be written as a combination with positive components in the second basis.
  In this frame, the approximation appears under the form of a least-square approximation with $2n$ linear constraints on the n coefficients, expressing the fact that the coefficients are nonnegative in the two above transformed bases.

In addition to the previous basis functions, it is possible to include more general/generic basis functions such as polynomial, radial, or wavelet functions in order to guarantee simple forecasting trends. For instance, one can add affine functions in order to include the possibility of predicting with a simple linear extrapolation to the range of possible forecasts offered by the reduced model. Given the overall exponential behavior of the health data, we have added an exponential function of the form $b_0\left(t\right)=\xi exp(-{\xi }^{\prime }t)$ (or $g_0\left(t\right)=\psi exp(-{\psi }^{\prime }t)$) to the reduced basis functions $\{b_1,\dots ,b_n\}$ (or $\{g_1,\dots ,g_n\}$) where the real-valued nonnegative parameters $\xi ,{\xi }^{\prime },\psi ,{\psi }^{\prime }$ are obtained through a standard exponential regression from ${\mathbf{\beta }}_{\mathrm{obs}}^{*}$ (or ${\mathbf{\gamma }}_{\mathrm{obs}}^{*}$) associated with the targeted profiles of infectious people; that is, the profiles defined in the time interval $[0,T]$ that should be extrapolated to [T, $\tau$]. In other words, the final reduced models that we use for forecasting are

$${\mathrm{B}}_{n+1}=\mathrm{span}\{b_0,b_1,\dots ,b_n\},{\mathrm{G}}_{n+1}=\mathrm{span}\{g_0,g_1,\dots ,g_n\}\subset L^{\infty }([0,T+\tau ],\mathbb{R}),$$

Indeed, including the exponential functions in the reduced models gives easy access to the overall behavior of the parameters $\beta$ and $\gamma$; the rest of the basis functions generated from the training sets catch the higher-order approximations and allow then to improve the extrapolation.

Remark 1.

Reduced models on $\mathcal{I}=\left\{\mathit{I}\right(\mu ):\mu \in \mathcal{P}\}$ and $\mathcal{R}=\left\{\mathit{R}\right(\mu ):\mu \in \mathcal{P}\}$ Instead of applying model reduction to the sets $\mathcal{B}$ and $\mathcal{G}$, as we do in our approach, we could apply the above techniques directly to the sets of solutions $\mathcal{I}$ and $\mathcal{R}$ of the SIR models with time-dependent coefficients in $\mathcal{B}$ and $\mathcal{G}$. In this case, the resulting approximation would however not follow SIR dynamics.

3. Methodology for Multiple Regions Including Population Mobility Data

The forecasting method of Section 2.2 for a single region can be extended to the treatment of multiple regions involving population mobility. The prediction scheme is based on a multi-regional SIR with time-dependent coefficients. Compared to other more detailed models, the main advantage of our approach is that it drastically reduces the number of parameters to be estimated. Indeed, detailed multi-regional models such as multi-regional extensions of the above SEI5CHRD and SE2IUR models from Section 2.3 require a number of parameters that quickly grows with the number P of regions involved. Their calibration thus requires large amounts of data which, in addition, may be unknown, very uncertain, or not available. In a forthcoming paper, we will apply the fully multi-regional setting for the post-lockdown period.

The structure of this section is the same as above for the case of a single region. In Section 3.1, we begin by introducing the multi-regional SIR model with time-dependent coefficients and associated detailed models. As with any multi-regional model, mobility data are required as input data, and the nature and level of detail of the available data motivates certain choices regarding the modeling of the multi-regional SIR (as well as the other detailed models). We then present in Section 3.2 the general pipeline, in which we emphasize the high modularity of the approach.

3.1. Multi-Regional Compartmental Models

In the spirit of fluid flow modeling, there are essentially two ways of describing mobility between regions:

• In a Eulerian description, we take the regions as fixed references for which we record incoming and outgoing travels;
• In a Lagrangian description, we follow the motion of people living in a certain region and record their travels in the territory. We can expect this modeling to be more informative regarding the geographical spread of the disease, but it comes at the cost of additional details regarding the home region of the population.

Note that both descriptions hold at any coarse or fine geographical level, in the sense that what we call the regions could be taken to be full countries, departments within a country, or very small geographical partitions of a territory. We next describe the multi-regional SIR models with the Eulerian and Lagrangian descriptions of population fluxes, which form- the output of our methodology.

3.1.1. Multi-Regional SIR Models with Time-Dependent Parameters

Eulerian description of population flux: Assume that we have P regions and the number of people in region i is $N_i$ for $i=1,\dots ,P$. Due to mobility, the population in each region varies, so $N_i$ depends on t. However, the total population is assumed to be constant and equal to N; that is,

$$N=\sum _{i=1}^P{N}_i\left(t\right),\forall t\ge 0.$$

For any $t\ge 0$, let ${\lambda }_{i\to j}\left(t\right)\in [0,1]$ be the probability that people from i travel to j at time t. In other words, ${\lambda }_{i\to j}\left(t\right)N_i\left(t\right)\delta t$ is the number of people from region i that have travelled to region j between time t and $t+\delta t$. Note that we have

$$\sum _{j=1}^P{\lambda }_{i\to j}\left(t\right)=1,\forall t\ge 0.$$

Since, for any $\delta t\ge 0$,

$$N_i(t+\delta t)=N_i\left(t\right)-\sum _{j\ne i}{\lambda }_{i\to j}\left(t\right)N_i\left(t\right)\delta t+\sum _{j\ne i}{\lambda }_{j\to i}\left(t\right)N_j\left(t\right)\delta t$$

dividing by $\delta t$ and taking the limit $\delta t\to 0$ yields

$$\frac{\mathrm{d}{N}_i}{\mathrm{d}t}\left(t\right)=-\sum _{j\ne i}{\lambda }_{i\to j}\left(t\right)N_i\left(t\right)+\sum _{j\ne i}{\lambda }_{j\to i}\left(t\right)N_j\left(t\right).$$

Note that we have

$$\sum _{i=1}^P\frac{\mathrm{d}{N}_i}{\mathrm{d}t}\left(t\right)=0,\forall t\ge 0.$$

Thus, ${\sum }_i{N}_i\left(t\right)={\sum }_i{N}_i\left(0\right)=N$, which is consistent with our assumption that the total population is constant.

The time evolution of the $N_i$ is known in this case if we are given the ${\lambda }_{i\to j}\left(t\right)$ from Eulerian mobility data. In addition to these mobility data, we also have the data of the evolution of infected and removed people, and our goal is to fit a multi-regional SIR model that is in accordance with this data. Thus, we propose the following model.

Denoting $S_i$, $I_i$ and $R_i$ as the number of susceptible, infectious and removed people in region i at time t, we first have the relation

$$N_i\left(t\right)=S_i\left(t\right)+I_i\left(t\right)+R_i\left(t\right)\iff 1=\frac{S_i\left(t\right)}{N_i\left(t\right)}+\frac{I_i\left(t\right)}{N_i\left(t\right)}+\frac{R_i\left(t\right)}{N_i\left(t\right)}.$$

Note that from the second relation, it follows that

$$0=\frac{\mathrm{d}}{\mathrm{d}t}\frac{S_i}{N_i}+\frac{\mathrm{d}}{\mathrm{d}t}\frac{I_i}{N_i}+\frac{\mathrm{d}}{\mathrm{d}t}\frac{R_i}{N_i}.$$

(9)

To model the evolution between compartments, one possibility is the following SIR model:

$$\begin{array}{c}\frac{\mathrm{d}}{\mathrm{d}t}\frac{S_i}{N_i}=-\left({\beta }_i{\lambda }_{i\to i}\frac{I_i}{N_i}+\sum _{j\ne i}{\beta }_j{\lambda }_{j\to i}\frac{I_j}{N_j}\right)\frac{S_i}{N_i}\hfill \\ \frac{\mathrm{d}}{\mathrm{d}t}\frac{I_i}{N_i}=-\frac{\mathrm{d}}{\mathrm{d}t}\frac{S_i}{N_i}-{\gamma }_i\frac{I_i}{N_i}\hfill \\ \frac{\mathrm{d}}{\mathrm{d}t}\frac{R_i}{N_i}={\gamma }_i\frac{I_i}{N_i},\hfill \end{array}$$

(10)

The parameters ${\beta }_i$, ${\gamma }_i$, $N_i$ depend on t, but we have omitted this dependence for ease of reading. Introducing the compartmental densities

$$s_i=\frac{S_i}{N_i},i_i=\frac{I_i}{N_i},r_i=\frac{R_i}{N_i},$$

the system equivalently reads

$$\begin{array}{c}\frac{\mathrm{d}}{\mathrm{d}t}{s}_i=-\left({\beta }_i{\lambda }_{i\to i}{i}_i+\sum _{j\ne i}{\beta }_j{\lambda }_{j\to i}{i}_j\right)s_i\hfill \\ \frac{\mathrm{d}}{\mathrm{d}t}{i}_i=-\frac{\mathrm{d}}{\mathrm{d}t}{s}_i-{\gamma }_i{i}_i\hfill \\ \frac{\mathrm{d}}{\mathrm{d}t}{r}_i={\gamma }_i{i}_i,\hfill \end{array}$$

(11)

Before going further, some points should be highlighted:

• The model is consistent in the sense that it satisfies (9), and when $P=1$, we recover the traditional SIR model;
• Under lockdown measures, ${\lambda }_{i\to j}\approx {\delta }_{i,j}$ and the population $N_i\left(t\right)$ remains practically constant. As a result, the evolution of each region is decoupled from the others, and each region can be addressed with a mono-regional approach;
• The use of ${\beta }_j$ in Equation (11) is debatable. When people from region j arrive in region i, it may be reasonable to assume that the contact rate is ${\beta }_i$;
• The use of ${\lambda }_{j\to i}$ in Equation (11) is also very debatable. The probability ${\lambda }_{j\to i}$ was originally defined to account for the mobility of people from region j to region i without specifying the compartment. However, in Equation (11), we need the probability of mobility of infectious people from region j to region i, which we denote by ${\mu }_{j\to i}$ in the following. It seems reasonable to think that ${\mu }_{j\to i}$ may be smaller than ${\lambda }_{j\to i}$, because as soon as people become symptomatic and suspect an illness, they will probably not move. Two possible options would be as follows:

  -

  We could try to estimate ${\mu }_{j\to i}$. If symptoms arise, for example, 2 days after infection and if people recover in 15 days on average, then we could say that ${\mu }_{j\to i}=2/15{\lambda }_{j\to i}$.

  -

  As the above seems to be quite empirical, another option would be to use ${\lambda }_{j\to i}$ and absorb the uncertainty in the values of the ${\beta }_j$ that can be fitted.

• We choose not to add mobility in the R compartment as this does not modify the dynamics of the epidemic spread; only adjustments in the population sizes are needed.

Lagrangian description of population flux: We call the above description Eulerian because we have fixed the regions as a fixed reference. Another possibility is to follow the trajectories of inhabitants of each region, in the same spirit as when we follow the trajectories of fluid particles.

Let $S_i$, $I_i$, and $R_i$ now be the number of susceptible, infectious and removed people who are resident in region $i,i=1,\dots ,P$. It is reasonable to assume that $S_i\left(t\right)+I_i\left(t\right)+R_i\left(t\right)$ is constant in time. However, not all the residents of region i may be in that region at time t. Let ${\lambda }_{j\to k}^{\left(i\right)}\left(t\right)$ be the probability that susceptible people resident in i travel from region j to region k at time t. With this notation, ${\lambda }_{i\to i}^{\left(i\right)}\left(t\right)$ is the probability that susceptible people resident at i remain in region i at time t. Similarly, let ${\mu }_{j\to k}^{\left(i\right)}\left(t\right)$ be the probability that infectious people resident in i travel from region j to k at time t. Thus, the total number of susceptible and infectious people that are in region i at time t is

$$\begin{array}{cc}\hfill {\mathcal{S}}_i\left(t\right)& =\sum _{k=1}^P\sum _{j=1}^P\left({\lambda }_{j\to i}^{\left(k\right)}\left(t\right)-{\lambda }_{i\to j}^{\left(k\right)}\left(t\right)\right)S_k\left(t\right)\hfill \end{array}$$

$$\begin{array}{cc}\hfill {\mathcal{I}}_i\left(t\right)& =\sum _{k=1}^P\sum _{j=1}^P\left({\mu }_{j\to i}^{\left(k\right)}\left(t\right)-{\mu }_{i\to j}^{\left(k\right)}\left(t\right)\right)S_k\left(t\right)\hfill \end{array}$$

We can thus write the evolution over $S_i$, $I_i$, $R_i$ as

$$\begin{array}{cc}\hfill \frac{\mathrm{d}{S}_i}{\mathrm{d}t}& =-\sum _{j=1}^P\sum _{k=1}^P{\beta }_k\left(t\right){\lambda }_{j\to k}^{\left(i\right)}\left(t\right)S_i\left(t\right){\mathcal{I}}_k\left(t\right)\hfill \\ \hfill \frac{\mathrm{d}{I}_i}{\mathrm{d}t}& =-\frac{\mathrm{d}{S}_i}{\mathrm{d}t}-{\gamma }_i\left(t\right)I_i\left(t\right)\hfill \\ \hfill \frac{\mathrm{d}{R}_i}{\mathrm{d}t}& ={\gamma }_i\left(t\right)I_i\left(t\right)\hfill \end{array}$$

(12)

Note that $S_i\left(t\right)+I_i\left(t\right)+R_i\left(t\right)$ is constant, which is consistent with the fact that, in our model,

$$\frac{\mathrm{d}}{\mathrm{d}t}(S_i+I_i+R_i)=0.$$

We emphasize that, to implement this model, the Lagrangian mobility data ${\lambda }_{j\to k}^{\left(i\right)}$ are required for all $(i,j,k)\in {\{1,\dots ,P\}}^3$.

Notation: In the following, we gather the compartmental variables in vectors

$$\overrightarrow{\mathbf{S}}:={\left(\mathbf{S}\right)}_{i=1}^P,\overrightarrow{\mathbf{I}}:={\left(\mathbf{I}\right)}_{i=1}^P,\overrightarrow{\mathbf{R}}:={\left(\mathbf{R}\right)}_{i=1}^P$$

as well as the time-dependent coefficients

$$\overrightarrow{\mathbf{\beta }}={\left(\mathbf{\beta }\right)}_{i=1}^P,\overrightarrow{\mathbf{\gamma }}={\left(\mathbf{\gamma }\right)}_{i=1}^P.$$

For any $\overrightarrow{\mathbf{\beta }}$ and $\overrightarrow{\mathbf{\gamma }}\in {\left(L^{\infty }\left([0,T]\right)\right)}^P$, we denote by

$$(\overrightarrow{\mathbf{S}},\overrightarrow{\mathbf{I}},\overrightarrow{\mathbf{R}})=\mathrm{Multiregional}\_\mathrm{SIR}(\overrightarrow{\mathbf{S}}\left(0\right),\overrightarrow{\mathbf{I}}\left(0\right),\overrightarrow{\mathbf{R}}\left(0\right),\overrightarrow{\mathbf{\beta }},\overrightarrow{\mathbf{\gamma }},[0,T])$$

the output of any of the above multi-regional SIR models. For simplicity, in what follows, we omit the initial condition in the notation.

3.1.2. Detailed Multi-Regional Models with Constant Coefficients

In the spirit of the multi-regional SIR, one can formulate detailed multi-regional versions of more detailed models such as those introduced in Section 2.1. We omit the details for the sake of brevity.

3.2. Forecasting for Multiple Regions with Population Mobility

Similar to the mono-regional case, we assume that we are given health data in $[0,T]$ in all regions. The observed data in region i are the series of infected people, denoted as $I_i^{\mathrm{obs}}$, and recovered people, denoted as $R_i^{\mathrm{obs}}$. They are usually given at a national or a regional scale and on a daily basis.

We propose to fit the data and provide forecasts with SIR models with time-dependent parameters ${\mathbf{\beta }}_i$ and ${\mathbf{\gamma }}_i$ for each region i. As in the mono-regional case, we can prove that such a simple family possesses optimal fitting properties for our purposes. In the current case, the cost function reads

$$\begin{array}{cc}\hfill & \mathcal{J}(\overrightarrow{\mathbf{\beta }},\overrightarrow{\mathbf{\gamma }}|{\overrightarrow{\mathbf{I}}}_{\mathrm{obs}},{\overrightarrow{\mathbf{R}}}_{\mathrm{obs}},[0,T]):=\sum _{i=1}^P{\int }_0^T\left(|I_i\left(t\right)-I_i^{\mathrm{obs}}{\left(t\right)|}^2+{|R_i\left(t\right)-R_i^{\mathrm{obs}}\left(t\right)|}^2\right)\mathrm{d}t\hfill \\ \hfill & \mathrm{such}\mathrm{that}\left(\overrightarrow{\mathbf{S}},\overrightarrow{\mathbf{I}},\overrightarrow{\mathbf{R}}\right)=\mathrm{Multiregional}\_\mathrm{SIR}\left(\overrightarrow{\mathbf{\beta }},\overrightarrow{\mathbf{\gamma }},[0,T]\right),\hfill \end{array}$$

and the fitting problem is the optimal control problem of finding

$$J^{*}=\underset{\overrightarrow{\mathbf{\beta }},\overrightarrow{\mathbf{\gamma }}\in {\left(L^{\infty }\left([0,T]\right)\right)}^P\times {\left(L^{\infty }\left([0,T]\right)\right)}^P}{inf}\mathcal{J}(\overrightarrow{\mathbf{\beta }},\overrightarrow{\mathbf{\gamma }}|{\overrightarrow{\mathbf{I}}}_{\mathrm{obs}},{\overrightarrow{\mathbf{R}}}_{\mathrm{obs}},[0,T]).$$

(13)

The following proposition ensures the existence of a unique minimizer under certain conditions. To prove this, it is useful to remark that any of the above multi-regional SIR models (see (11) and (12)) can be written in the general form

$$\begin{array}{cc}\hfill \frac{\mathrm{d}\overrightarrow{\mathbf{S}}}{\mathrm{d}t}& =\mathrm{M}(\mathsf{\Lambda }\left(t\right),\overrightarrow{\mathbf{S}}\left(t\right),\overrightarrow{\mathbf{I}}\left(t\right))\overrightarrow{\mathbf{\beta }}\hfill \\ \hfill \frac{\mathrm{d}\overrightarrow{\mathbf{I}}}{\mathrm{d}t}& =-\frac{\mathrm{d}\overrightarrow{\mathbf{S}}}{\mathrm{d}t}-diag\left(\mathbf{I}\left(t\right)\right)\overrightarrow{\mathbf{\gamma }}\hfill \\ \hfill \frac{\mathrm{d}\overrightarrow{\mathbf{R}}}{\mathrm{d}t}& =diag\left(\mathbf{I}\left(t\right)\right)\overrightarrow{\mathbf{\gamma }},\hfill \end{array}$$

where, by a slight change of notation, the vectors $\overrightarrow{\mathbf{S}}$, $\overrightarrow{\mathbf{I}}$ and $\overrightarrow{\mathbf{R}}$ are the densities of population in the case of the Eulerian approach (see Equation (11)). They are classical population numbers in the case of the Lagrangian approach (see Equation (12)). $diag\left(\mathbf{I}\right(t\left)\right)$ is the $P\times P$ diagonal matrix with diagonal entries given by the vector $\mathbf{I}\left(t\right)$. $\mathrm{M}(\mathsf{\Lambda }\left(t\right),\overrightarrow{\mathbf{S}}\left(t\right),\overrightarrow{\mathbf{I}}\left(t\right))$ is a matrix of size $P\times P$ that depends on the vectors of susceptible and infectious people $\overrightarrow{\mathbf{S}}\left(t\right)$, $\overrightarrow{\mathbf{I}}\left(t\right)$ and on the mobility matrix $\mathsf{\Lambda }$. In the case of the Eulerian description, $\mathsf{\Lambda }\left(t\right)={\left({\lambda }_{i,j}\left(t\right)\right)}_{1\le i,j\le P}$ and in the case of the Lagrangian approach $\mathsf{\Lambda }\left(t\right)$ is the $P\times P\times P$ tensor $\mathsf{\Lambda }\left(t\right)={\left({\lambda }_{j,k}^{\left(i\right)}\left(t\right)\right)}_{1\le i,j,k\le P}$. For example, in the case of the Eulerian model (12), the matrix M reads

$$\mathrm{M}(\mathsf{\Lambda }\left(t\right),\overrightarrow{\mathbf{S}}\left(t\right),\overrightarrow{\mathbf{I}}\left(t\right))=-diag\left(\overrightarrow{\mathbf{S}}\left(t\right)\right){\mathsf{\Lambda }}^{T}diag\left(\overrightarrow{\mathbf{I}}\left(t\right)\right)=-{\left(S_i{\lambda }_{i\to j}{I}_j\right)}_{1\le i,j\le P}$$

(14)

Proposition 2.

If $\mathrm{M}(\mathsf{\Lambda }\left(t\right),\overrightarrow{\mathbf{S}}\left(t\right),\overrightarrow{\mathbf{I}}\left(t\right))$ is invertible for all $t\in [0,T]$, then there exists a unique minimizer $({\overrightarrow{\mathbf{\beta }}}^{*},{\overrightarrow{\mathbf{\gamma }}}^{*})$ to problem (13).

Proof. 

Since we assume that $\mathrm{M}(\mathsf{\Lambda }\left(t\right),\overrightarrow{\mathbf{S}}\left(t\right),\overrightarrow{\mathbf{I}}\left(t\right))$ is invertible for every $t\in [0,T]$, we can set

$$\begin{array}{c}\hfill \left\{\begin{array}{cc}{\overrightarrow{\mathbf{\beta }}}^{*}\left(t\right)\hfill & :=M^{-1}\left(t\right)\frac{\mathrm{d}\overrightarrow{\mathbf{S}}}{\mathrm{d}t}\hfill \\ {\overrightarrow{\mathbf{\gamma }}}^{*}\left(t\right)\hfill & :={diag}^{-1}\left(I\left(t\right)\right)\frac{\mathrm{d}\overrightarrow{\mathbf{R}}}{\mathrm{d}t}\hfill \end{array}\right.\end{array}$$

or equivalently

$$\begin{array}{c}\hfill \left\{\begin{array}{cc}{\overrightarrow{\mathbf{\beta }}}^{*}\left(t\right)\hfill & :=M^{-1}\left(t\right)\frac{\mathrm{d}\overrightarrow{\mathbf{S}}}{\mathrm{d}t}\hfill \\ {\overrightarrow{\mathbf{\gamma }}}^{*}\left(t\right)\hfill & :=-{diag}^{-1}\left(I\left(t\right)\right)\left(\frac{\mathrm{d}\overrightarrow{\mathbf{I}}}{\mathrm{d}t}+\mathrm{M}(\mathsf{\Lambda }\left(t\right),\overrightarrow{\mathbf{S}}\left(t\right),\overrightarrow{\mathbf{I}}\left(t\right)){\overrightarrow{\mathbf{\beta }}}^{*}\right)\hfill \end{array}\right.\end{array}$$

so that

$$({\overrightarrow{\mathbf{S}}}_{\mathrm{obs}},{\overrightarrow{\mathbf{I}}}_{\mathrm{obs}},{\overrightarrow{\mathbf{R}}}_{\mathrm{obs}})=\mathrm{Multiregional}\_\mathrm{SIR}\left({\overrightarrow{\mathbf{\beta }}}^{*},{\overrightarrow{\mathbf{\gamma }}}^{*},[0,T]\right)$$

and

$$\mathcal{J}({\overrightarrow{\mathbf{\beta }}}^{*},{\overrightarrow{\mathbf{\gamma }}}^{*}|{\overrightarrow{\mathbf{I}}}_{\mathrm{obs}},{\overrightarrow{\mathbf{R}}}_{\mathrm{obs}},[0,T])=0$$

which implies that $J^{*}=0$. □

Before continuing, let us comment on the invertibility of $\mathrm{M}(\mathsf{\Lambda }(t),\overrightarrow{\mathbf{S}}\left(t\right),\overrightarrow{\mathbf{I}}\left(t\right))$ which is necessary in Proposition 2. A sufficient condition to ensure this is if the matrix is diagonally dominant row-wise or column-wise. This yields certain conditions on the mobility matrix $\mathsf{\Lambda }\left(t\right)$ with respect to the values of $\overrightarrow{\mathbf{S}}\left(t\right)$, $\overrightarrow{\mathbf{I}}\left(t\right)$. For example, if M is defined as in Equation (14), the matrix is diagonally dominant in each row if for every $1\le i\le P$,

$${\lambda }_{i\to i}>\sum _{j\ne i}{\lambda }_{i\to j}\frac{I_j}{I_i}.$$

Similarly, if for every $1\le j\le P$,

$${\lambda }_{j\to j}>\sum _{i\ne j}{\lambda }_{i\to j}\frac{S_i}{S_j},$$

then the matrix is diagonally dominant for each column and guarantees invertibility. Note that any of the above conditions is satisfied in situations with little or no mobility where ${\lambda }_{i\to i}\approx {\delta }_{i,j}$.

Now that we have exactly defined the set-up for the multi-regional case, we can follow the same steps in Section 2.2 to derive forecasts involving model reduction for the time-dependent variables $\overrightarrow{\mathbf{\beta }}$ and $\overrightarrow{\mathbf{\gamma }}$.

4. Numerical Results

In this section, we apply our forecasting method to the ongoing COVID-19 pandemic, which spread in the year 2020 in France and started approximately in February. Particular emphasis is placed on the first pandemic wave, for which we consider the period from 19 March to 20 May 2020. Due to the lockdown imposed between 17 March and 11 May, inter-regional population mobility was drastically reduced in that period. Studies using anonymized Facebook data have estimated the reduction to be 80% (see [28]). As a result, it is reasonable to treat each region independently from the rest, and we apply the mono-regional setting in Section 2. Here, we focus on the case of the Paris region, and we report different forecasting errors obtained using the methods described in Section 2. Some forecasts are also shown for the second wave for the Paris region between 24 September and 25 November.

The numerical results are presented as follows. Section 4.1 explains the sources of health data. Section 4.2.1 and Section 4.2.2 explore the results in detail and present a study of the forecasting power of the methods in a two-week horizon. Section 4.2.3 displays forecasts for the second wave. Section 4.2.4 aims to illustrate the robustness of the forecasting over longer periods of time. A discussion of the fitting errors of the methods is given in Appendix A. Additional results highlighting the accuracy of the forecasts are shown in Appendix B.

4.1. Data

We use public data from Santé Publique France (https://www.data.gouv.fr/en/datasets/donnees-hospitalieres-relatives-a-lepidemie-de-covid-19/) to get the numbers $I_{\mathrm{obs}}\left(t\right)$ of infected and $R_{\mathrm{obs}}\left(t\right)$ of removed people. As shown in Figure 3, the raw data present some oscillations at the scale of a week, which are due to administrative delays for the cases to be officially reported by hospitals. For our methodology, we have smoothed the data by applying a 7 day moving average filter. In order to account for the total number of infected people, we also multiply the data by an adjustment factor $f_{\mathrm{adj}}=15$ as stated in the literature (indeed, it is said in [29] that “of the 634 confirmed cases, a total of 306 and 328 were reported to be symptomatic and asymptomatic”, and in [10], it is claimed that the probability of developing severe symptoms for a symptomatic patient is 0 for children, 0.1 for adults and 0.2 for seniors; thus, if one takes $p=0.13$ as an approximate value of these probabilities, one may estimate the adjustment factor as $f_{\mathrm{adj}}=\frac{634}{328}\times \frac{1}{p}\approx 15$). Obviously, this factor is uncertain and could be improved in the light of further retrospective studies of the outbreak. However, note that when $S\simeq N$, which is the case at the start of the epidemic, the impact of this factor is negligible in the dynamics as can be understood from (3). In addition, since we use the same factor to provide a forecast of hospitalized people, the influence on the choice is minor.

4.2. Results

Using the observations $I_{\mathrm{obs}}\left(t\right)$ and $R_{\mathrm{obs}}\left(t\right)$, we apply a finite difference scheme in Formula (3) to derive ${\beta }_{\mathrm{obs}}^{*}\left(t\right)$ and ${\gamma }_{\mathrm{obs}}^{*}\left(t\right)$ for $t\in [0,T]$. Figure 4 shows the values of these parameters as well as the basic reproduction number $R_{0,\mathrm{obs}}^{*}={\beta }_{\mathrm{obs}}^{*}/{\gamma }_{\mathrm{obs}}^{*}$ for the first pandemic wave in Paris.

We next follow the steps presented in Section 2.2 to obtain the forecasts. In the learning phase (step 1), we use two parametric detailed models of SE2IUR and SEI5CHRD types to generate training sets ${\mathcal{B}}_{\mathrm{tr}}$ and ${\mathcal{G}}_{\mathrm{tr}}$ composed of $K=2618$ training functions $\mathbf{\beta }\left(\mu \right)$ and $\mathbf{\gamma }\left(\mu \right)$ where $\mu$ are uniformly sampled in the set of parameters ${\mathcal{P}}_{\mathrm{tr}}$ in the vicinity of the parameter values suggested in the literature [9,10]. Based on these training sets, we finish step 1 by building three types of reduced models: SVD, NMF and ENG (see Section 2.4).

Given the reduced bases ${\mathrm{B}}_n$ and ${\mathrm{G}}_n$, we next search for the optimal ${\mathbf{\beta }}_n^{*}\in {\mathrm{B}}_n$ and ${\mathbf{\gamma }}_n^{*}\in {\mathrm{G}}_n$ that best fit the observations (step 2 of our procedure). For this fitting step, we consider two loss functions:

• $\mathrm{routine}-\mathrm{IR}$: loss function $\mathcal{J}(\mathbf{\beta },\mathbf{\gamma }|{\mathbf{I}}_{\mathrm{obs}},{\mathbf{R}}_{\mathrm{obs}},[0,T])$ from (1),
• $\mathrm{routine}-$$\mathbf{\beta }$$\mathbf{\gamma }$: loss function $\tilde{\mathcal{J}}(\mathbf{\beta },\mathbf{\gamma }|{\mathbf{\beta }}_{\mathrm{obs}}^{*},{\mathbf{\gamma }}_{\mathrm{obs}}^{*},[0,T])$ from (6)

We study the performance of each of the three reduced models and the impact of the dimension n of the reduced model in terms of the fitting error. The presentation of these results is presented in Appendix A in order not to overload the main discussion. The main conclusion is that the fitting strategy using SVD-reduced bases provides smaller errors than NMF and ENG, especially when we increase the number of modes n. This is illustrated in Figure 5 where we show the fittings obtained with routine-$\mathbf{\beta }\mathbf{\gamma }$ and $n=10$ for the first wave (from $t_0=19/03/2020$ to $T=20/05/2020$). We observe that SVD is the best at fitting ${\beta }_{obs}^{*}$ and ${\gamma }_{obs}^{*}$, while ENG produces a smoother fitting of the data. Although the smoother fitting with ENG yields larger fitting errors than SVD, we see in the next section that it yields better forecasts.

4.2.1. Forecasting for the First Pandemic Wave with a 14 Day Horizon

In this section, we illustrate the short-term forecasting behavior of our method. We consider a forecasting window of $\tau =14$ days and we examine several different starting days in the course of the first pandemic wave. The results are shown in Figure 6, Figure 7, Figure 8, Figure 9, Figure 10, Figure 11, Figure 12, Figure 13 and Figure 14. Recall that that the forecasting uses the coefficients of the reduced bases obtained by the fitting procedure but also the optimal initial condition of the forecast that minimizes the error on the three days prior to the start of the forecast. For each given fitting strategy (routine-IR, routine-$\mathbf{\beta }$$\mathbf{\gamma }$) and each given type of reduced model (SVD, NMF, ENG), we have chosen to plot an average forecast computed with predictions using reduced dimensions $n\in \{5,6,7,8,9,10\}$. This choice is a simple type of forecast combination, but of course other more elaborate aggregation options could be considered. The labels of the plots correspond to the following:

• $I_{SVD}$, $I_{NMF}$, $I_{ENG}$, $R_{SVD}$, $R_{NMF}$, $R_{ENG}$ are the average forecasts obtained using routine-$\mathbf{\beta }$$\mathbf{\gamma }$.
• $I_{SVD}^{*}$, $I_{NMF}^{*}$, $I_{ENG}^{*}$, $R_{SVD}^{*}$, $R_{NMF}^{*}$, $R_{ENG}^{*}$ are the average forecasts obtained using routine-IR.

The main observation from Figure 6, Figure 7, Figure 8, Figure 9, Figure 10, Figure 11, Figure 12, Figure 13 and Figure 14 is that the ENG-reduced model is the most robust and accurate forecasting method. Fitting ENG with routine-IR or routine-$\mathbf{\beta }$$\mathbf{\gamma }$ does not seem to lead to large differences in the quality of the forecasts, but routine-$\mathbf{\beta }$$\mathbf{\gamma }$ seems to provide slightly better results. This claim is further confirmed by the study of the numerical forecasting errors of the different methods shown in Appendix B.

Figure 6, Figure 7, Figure 8, Figure 9, Figure 10, Figure 11, Figure 12, Figure 13 and Figure 14 also show that the SVD-reduced model is very unstable and provides forecasts that blow up. This behavior illustrates the dangers of overfitting, namely that a method with high fitting power may present poor forecasting power due to instabilities. In the case of SVD, the instabilities stem from the fact that approximations are allowed to take negative values. This is the reason why NMF, which incorporates the nonnegative constraint, performs better than SVD. One of the reasons why ENG outperforms NMF is the enlargement of the cone of nonnegative functions (see Section 2.4). It is important to note that, with ENG, the reduced bases are directly related to well-chosen virtual scenarios, whereas SVD and NMF rely on matrix factorization techniques that provide purely artificial bases. This makes forecasts from ENG more realistic and therefore more reliable.

4.2.2. Focus on the Forecasting with ENG

For our best forecasting method (routine-$\mathbf{\beta }$$\mathbf{\gamma }$ using ENG), we plot in Figure 15, Figure 16, Figure 17, Figure 18, Figure 19, Figure 20, Figure 21, Figure 22 and Figure 23 the forecasts for each dimension $n=5$ to 10. The plots give the forecasts on a 14 day-ahead window for $\beta$, $\gamma$, and the resulting evolution of the infected I and removed R. We see that the method performs reasonably well for all values of n, proving that the results of the previous section with the averaged forecast are not compensating for spurious effects which could occur for certain values of n. We have chosen to display the inaccurate forecasts from 3 April, 7 April, and 11 April as they are among the worst predictions obtained using this method; however, it is important to mention that, despite the lack of accuracy in these cases, plausible epidemic behaviors remain, with different but realistic evolutions for $\beta$ and $\gamma$ compared to the actual evolution. Note that the method was able to predict the peak of the epidemic several days in advance (see Figure 15). We also observe that the prediction on $\gamma$ is difficult at all times due to the fact that ${\gamma }_{\mathrm{obs}}^{*}$ presents an oscillatory behavior. Despite this difficulty, the resulting forecasts for I and R are very satisfactory in general.

4.2.3. Forecasting of the Second Wave with ENG

The review took place during the month of November 2020 as the second COVID-19 pandemic wave hit France. We took advantage of this to enlarge the body of numerical results, and we provide some example forecasts with ENG for this wave in Figure 24, Figure 25 and Figure 26. As the figures illustrate, the method provides very satisfactory forecasts in a 14 day-ahead window. We again observe a satisfactory prediction of the second peak (Figure 24, Figure 25 and Figure 26) and the same difficulty in forecasting $\gamma$ due to the oscillations in ${\gamma }_{\mathrm{obs}}$, but this has not greatly impacted the quality of the forecasts for I and R.

4.2.4. Forecasts with ENG with a 28 Day-Ahead Window

To conclude this section, we extend the forecasting window to 28 days instead of 14 and study whether the introduced ENG method still provides satisfactory forecasts. As shown in Figure 27, Figure 28, Figure 29, Figure 30, Figure 31 and Figure 32, the results of the methods are quite stable for large windows. This shows that, in contrast to standard extrapolation methods using classical linear or affine regressions, the reduced basis catches the dynamics of $\mathbf{\beta }$ and $\mathbf{\gamma }$ not only locally but also at extended time intervals.

5. Conclusions

We have developed an epidemiological forecasting method based on reduced modeling techniques. Of the different strategies that have been explored, the one that outperforms the rest in terms of robustness and forecasting power involves reduced models that are built with an Enlarged Nonnegative Greedy (ENG) strategy. This method is novel and of interest in itself as it allows reduced models that preserve positivity and even other types of bounds to be built. Despite the fact that ENG does not have optimal fitting properties (i.e., interpolation properties), it is well suited for forecasting since, due to the preservation of the natural constraints of the coefficients, it generates realistic dynamics with few modes. The results have been presented for a mono-regional test case, and we plan to extend the present methodology to a multi-regional setting using mobility data.

Last but not least, we would like to emphasize that the developed strategy is very general and could be applied to the forecasting of other types of dynamics. The specificity of each problem may, however, require adjustments in the reduced models. This is exemplified in the present problem through the construction of reduced models that preserve positivity.