A profusion of mathematical models of the COVID-19 epidemic propagation has been developed in a very short time lapse. Specific models were designed to address particular objectives, such as forecasting the epidemics, understanding the impact of various interventions such as restraining mobility and contacts between users, or deploying contact tracing and case isolation.
The goal of this study is to propose a single model of COVID-19 epidemics based on a well understood family of probabilistic models -that of multi-type branching processes -and to detail simple extensions of this model that capture the impact of contact tracing procedures as well as user mobility.
Our model can serve several purposes: It can be used to assess from available time series so-called contact rates, and in turn the impact of user mobility on epidemic propagation. It can also be used in surveillance platforms to infer the status of specific sub-populations (in particular fraction of infectious or immunized persons), a prerequisite for the implementation of targeted testing campaigns.
Paper organization We introduce in Section 1 a discrete time epidemiological model of COVID-19 for a homogeneous population. We detail a probabilistic version of the model, that is a multitype branching process. It represents in a precise manner the durations of the principal phases of a patient's condition and the transition rates between such phases. A deterministic version of this model is then described. We then introduce a simple extension of the model to represent the impact of contact tracing and case isolation.
In Section 2, we describe two methods that both use model structure for inference of parameters and prediction. For the deterministic, or mean-field version of the model, we propose an approach to perform inference that relies on Perron's theory. For the probabilistic version of the model, by exploiting the fact that we are dealing with a hidden Markov model, we propose a Kalman filtering approach to perform inference.
Next, we present a numerical illustration in Section 3, where the model is fitted to daily counts of hospitalizations in the Paris area. There we estimate piecewise constant contact rates corresponding to distinct time periods, namely pre-confinement, intermediate, and confinement periods. We discuss the application of this fitting procedure to forecasting the future number of hospitalizations.
We then describe in Section 4 an extension for the contact rates to take daily mobility -or other covariates -into account, making these parameters depend on time.
We finally discuss in Section 5 the extension of the model to multiple sub-populations, corresponding e.g. to geographical regions and age ranges. Key parameters in this extension are contact matrices, describing the frequency of contacts between sub-populations. The routing dimension of this extension is detailed. we revisit the question of estimating model parameters, specifically contact matrices in this multi sub-populations model. We assume potentially timevarying contact matrices and consider estimation of such contact matrices based on user mobility traces.
Related work SIR or SEIR-type models of COVID-19 featuring asymptomatic or undocumented infectious individuals have been considered in [11] , [14] and [2] . Studies [14, 2] consider same contact rates for symptomatic and asymptomatic cases. Some works [13] also consider a split of the infectious phase to explicit the occurrence of aggravated symptoms.
General distributions of times spent in a given state have been studied. Prodromic phase, during which individuals are infectious, yet have not displayed symptoms, is considered in [3] . The model of [3] comprises a detailed breakdown in various phases, yet relies on geometric/exponential durations for time spent in each phase, while the model of [15] describes general discrete time distributions for the time spent in states.
A broad range of studies have addressed the issue of estimating the time spent in a given state. In [12] , the incubation period is modeled with a log-normal distribution with mean 5.8 days. In [9] , the median incubation period was estimated to be 5.1 days, and 97.5% of those who develop symptoms will do so within 11.5 days of infection.
[1] reviews studies of the infectious period's duration, reporting several estimates with substantial variations, giving for asymptomatic cases an estimate of 6.5-9.5 days. Pre-symptomatic infectious period is estimated to 0-4 days across studies. [16] also investigates the time between infection and remission/death. The model of Figure 1 is a simplified version of the model considered in [3] , which features additional states, including a breakdown between hospitalized cases entering Intensive Care Units or not, and a distinction between several severity levels of the infection. It is similar to the model considered in [13] , with a key difference in our model that is that we take the discrete time t into account. Figure 1 : State diagram for individual's condition in model 1. Dotted arrows show the contamination of new cases: the mean number of new E cases per day per A/P (resp. per I 1 /I 2 ) individual is α a (resp. α i ).
We consider the discrete time SEIR-type model depicted in Figure 1 . We classify the states of an individual as E for exposed; P for prodromic; A for infectious asymptomatic; I 1 for phase 1 infectious symptomatic; H for hospitalized; I 2 for (non-hospitalized) phase 2 infectious symptomatic; and R for removed from the contamination chain. After having moved to prodromic phase P , an individual becomes either asymptomatic (state A) or symptomatic, entering state I 1 . After some time in the I 1 state, two situations arise: either the individual is hospitalized, entering state H, either there is no aggravation, leading to state I 2 . Then, the last state is always R, where the individual is either dead or recovered.
We shall in fact adopt a more detailed representation for the state of an individual: its detailed state will be represented by a couple (a, t) where a ∈ {E, A, I 1 , H, I 2 , R} is the stage of the disease, and t the number of days spent in this stage. For instance, an individual in state (E, 3) just completes their third day in the exposed phase. Remark 1.1. (i) In epidemiological models one often tracks individuals in so-called susceptible state S: these are those individuals who may get exposed to infection. In the case of a negligible fraction of infected individuals, in particular in the early stages of the epidemics, on which we focus, it is not necessary to track them, hence our choice not to include state S in the diagram of Figure 1 .
(ii) We distinguish between asymptomatic and symptomatic individuals, as new symptomatic cases can be observed when the corresponding individuals consult their physician or enter hospital after experiencing serious symptoms.
(iii) We distinguish between hospitalized and non-hospitalized cases. The number of new hospital admissions can be measured, for instance, from medical emergency services records.
(iv) We distinguish according to the number of days an individual has spent in its current state for the following reason. Traditional SEIR-type models implicitly assume that the time that an individual spends in a particular state admits an exponential distribution (for continuous time ODE-based models) or a geometric distribution (for discrete time models). However the durations spent in the states of interest for COVID-19 appear to be very far from exponentially distributed. Thus, for an accurate short-term prediction of the epidemics evolution, one needs to take this number of days into account.
We represent the dynamics by means of the following state variables. On each day t, we represent the population by a vector x(t) with entries:
x E,d (t) : number of individuals in phase E, having spent d days in this phase;
x P,d (t) : number of individuals in phase P , having spent d days in this phase;
x I1,d (t) : phase 1 infectious individuals, having spent d days in this phase;
x I2,d (t) : phase 2 infectious individuals having spent d days in this phase;
x A,d (t) : asymptomatic individuals having spent d days in this phase;
x H (t) : number of individuals newly hospitalized on day t.
By considering the state variable x H (t), instead of x H,d (t), we ignore the number of days an individual spends at hospital. This stems from the assumption that such cases are efficiently isolated and do not contribute to propagating the epidemics, an assumption that could be revised.
For each τ ∈ {E, P, I 1 , I 2 , A}, define r τ (d) the probability that corresponding phase will end on the following day, given that it has lasted d days. Denoting by p τ (d) the probability that phase τ lasts d days, r τ (d) is the associated failure rate, defined by
Available statistics suggest that typical durations lie in the ranges given on Figure 1 . A reasonable baseline assumption could consist in taking
• p P (1) = p P (2) = 1/2 and p P (d) = 0 otherwise, Note however that far longer durations have been observed, e.g. incubation lasting 10 days or more; I 1 phase between first symptoms and aggravation lasting as long as 20 days. This baseline is therefore by no means an accepted characterization of the dynamics of an individual's condition.
• p i : probability that previously exposed individual becomes symptomatic at end of incubation. A baseline assumption consists in taking p i = 0.7. Again, no consensus has emerged yet on the fraction of asymptomatic infected users and competing assumptions have been proposed (for baselines, data from APHP https://www.aphp.fr can be used, as done e.g. in [3] ).
• p h : probability that previously phase 1 infectious individual develops aggravated form at end of phase 1. A baseline assumption consists in taking p h = 0.05, again with the same cautionary note.
• α i : average number of new exposed individuals generated by single infectious individual (I 1 , I 2 ) in a single day. Here we do not distinguish between phases 1 and 2, an assumption which could be revisited.
• α a : average number of new exposed individuals generated by single asymptomatic or prodromic individual in a single day.
In the sequel we assume that α a = α i . This assumption is often made, for instance in [14, 2] . In contrast, [11] makes an explicit distinction in contact rates α for documented cases and µα for undocumented ones. Especially, an estimate for µ in China between January 10th and January 23th is given: µ ∼ 0.55 (95% confidence interval: (0.46, 0.62)). Denoting by p d the probability of death conditional on having been hospitalized, and assuming that only hospitalized individuals die, this leads to two other state variables:
The number of deaths as written neglects time between hospitalization and death.
Observables We assume that coordinates x H (t) are observed by the end of day t. This would correspond to the number of new hospitalized patients with COVID-19 pathology. In contrast, the states variables x E (t), x A (t), x P (t) are never observed. They will later be referred to as hidden trajectories. It would be straightforward to adjust the model to incorporate other observables, such as the number of patients consulting their physician on each day.
This is a deterministic, linear model, based on the assumption that the fraction of susceptible individuals in the total population is close to 1. In other words it is meant to represent stages of the epidemics where collective immunization is yet negligible. Simple adjustments can be made to transform this model to account for a non-negligible, time-evolving, fraction of immunized population.
Assume there is a maximal number of days, h, that each of the phases E, P , I 1 , I 2 , A can last. Let x E = (x E,1 , . . . , x E,h ) , x P = (x P,1 , . . . , x P,h ) , and similarly for h-dimensional vectors x A , x I1 , x I2 . The state space on day t, x(t), is given by
The dynamics are then given by
where matrix M is specified by the detailed equations, for d ∈ [h − 1]:
We now transform the previous mean field model into a probabilistic counterpart: a multi-type (or multi-dimensional) branching process. Consideration of probabilistic behaviour is particularly relevant when we deal with not-so-large infected populations, or in sub-critical cases where the largest eigenvalue λ 1 (M ) of matrix M defined in (3) verifies λ 1 (M ) ≤ 1. Indeed, when λ 1 (M ) > 1, by results of Kesten and Stigum [7] , there is a random variable Z such that for large t, the state X(t) renormalized by λ 1 (M ) −t converges to Zu, where u is the eigenvector of M associated with eigenvalue λ 1 (M ). This provides a theoretical justification of the mean field model in the case of large infected populations. Figure 2 gives a numerical illustration of this result.
The number of newly exposed individuals now follows a Poisson distribution with mean as specified in the mean field model. Determination of durations in each state, and choice of next state at branches in the state diagram of Figure 1 , are made with independent random coins with the prescribed probabilities.
For clarity, to distinguish between the two versions, random variables will be denoted in capital letters through the rest of this article. Binomial (resp. Poisson) distributions are denoted by B (resp. P).
Thus, conditioning to state vector X(t) at time t, the state X(t + 1) at time t + 1 verifies: 
One motivation of the present study is the control of the epidemics via case isolation, that is quarantining of identified cases so that they do not contribute to further virus propagation. The main approach to deploy case isolation relies on identification of infectious cases via tests. Thus, when testing capacity is limited, efficient case isolation crucially depends on targeted test campaigns. Contact tracing is of paramount importance for targeting tests: it allows to determine who has been in contact with infectious cases, and hence focus tests on such contacts.
We now describe an extension of our previously introduced model to include contact tracing and case isolation. In a nutshell, the model assumes that those individuals subsequently infected by a specific person (its children in the branching tree associated to the branching process) get tested with some probability once that specific person is positively tested (and subsequently isolated). This formalization of contact tracing has already been considered in [8] for a simple branching process, analyzing the efficiency of contact tracing plus case isolation to control epidemics. We give here a more general version, providing a generic method to expand the state space of a multi-type branching process in order to capture the introduction of contact tracing plus case isolation.
Consider then a multi-type branching process where user types are indexed by j ∈ [J]. For the sake of generality, we shall assume that types j encode the succession of future phases φ that an individual will visit over the coming days. Thus an individual with type j = {E, E, P, I 1 , H} is in state E for the present day and the next, and will then spend a single day in the states P , then I 1 , and then H subsequently. Also, an individual with type j = {φ 1 , . . . , φ } on a given day will become an individual with type j = {φ 2 , . . . , φ } on the next day.
Note that the previous branching process gives rise to an equivalent branching process with this set of types, as can be seen by sampling, upon birth of an individual, the sequence of states φ it will visit after its birth.
We now assume that there is a set Φ 0 of distinguished states φ ∈ Φ 0 such that an individual entering a state φ in Φ 0 automatically incurs a positive test. Concretely, the state φ = H, namely hospitalization of the individual, would automatically trigger a positive test, and we may think of Φ 0 as the singleton {H}.
We then assume that, for some fixed tracing probability p t > 0, the contact between an infectious individual and another individual that became infected is being recorded. The original branching process dynamics are then modified as follows: once an individual is tested positive, it is isolated (removed) and its traced contacts are subsequently tested positive on the following day, and subsequently isolated, and so on.
The key observation is that this modification preserves the branching structure, provided we perform the following extension of the original type space [J]. The new type space becomes
where the second coordinate d represents the number of days until the corresponding individual will be positively tested.
Upon birth of an individual whose parent has type (j, d), the newly born individual will have type (j , d ). We now specify the distribution of the pair (j , d ). First sample a sequence j = {φ 1 , . . . , φ } ∈ [J] conditionally on j with the same conditional distribution as in the original branching process. Let d 1 = inf{k ≥ 1 : φ k ∈ Φ 0 } be the number of days before the newly born individual enters some state in Φ 0 .
We then let
this reflects the fact that if traced, the new individual will become tested either when it enters some state in Φ 0 , or on the day after its parent gets tested, whichever happens first. If untraced, it only becomes tested on the day it enters some state in Φ 0 . Finally, we let j = {φ 1 , . . . , φ d }, i.e. we truncate the sequence j to length d , thereby encoding the impact of case isolation, by which a positively tested individual can no longer infect anyone else.
Note that further extensions can be considered that preserve the branching process structure. Here is an example of such an extension. Assume that on each day, each infected individual may be submitted to a randomly administered test, and trigger a positive test reply, this happening independently with probability ε on any given day. To reflect such random tests, the abovedescribed dynamics would then need to be amended in the following manner. For an individual newly born from a (j, d)-parent, sample a random variable X distributed as
where G follows a geometric distribution of parameter ε.
With j and d 1 as above, let then
The construction then goes through unchanged, letting j = {φ 1 , . . . , φ d } and the newly born individual's type being (j , d ).
The above extensions to the original branching process can be used to assess how the tracing probability p t , or the fraction of random tests, affects the potency of the epidemics. In particular, one may consider the mean progeny matrix M , indexed by extended types j ∈ [J ] and associated with the branching process that captures contact tracing. A crude criterion for success of contact tracing consists in requiring that the spectral radius ρ(M ) be strictly less than 1. Under such conditions, given a vector x(0) capturing the initial number of individuals of each type, then the expectation of X infected , the total number individuals ultimately infected, reads
In the above, e denotes the all-ones vector and the second equality is valid whenever ρ(M ) < 1. This expression for E(X infected ) can be further processed to obtain simple upper bounds. Particularly simple upper bounds have been obtained in [4] in the case where M is symmetric. However symmetry does not hold for the matrices M considered here, so that we cannot rely on the bounds of [4] .
In this section, we propose inference methods for our model. We consider a general deterministic (respectively, probabilistic) setting where the state vector at time t is denoted by x(t) (resp. X(t)).
Individual types are now denoted by j ∈ [J], and M i,j is the mean number of type i−children from a type j−parent. M is known as the mean-progeny matrix in branching process theory. It depends on a parameter of interest θ ∈ Θ, that encapsulates for instance p i , p h , α a , α i and time distributions in Section 1, as well as p t and c max if we consider contact tracing as in Section 1.3.
We assume that the observed trajectories are the first h coordinates of x(t) (resp. X(t)), denoted by x H (t) (resp. X H (t)), all other coordinates being hidden trajectories.
For simplicity, assume that the remaining unknown parameters of the model are the contact rates α a and α i , all other parameters having correctly been estimated, e.g. from patient statistics compiled at hospitals.
In the case where contact rates are constant at θ 0 = (α i,0 , α a,0 ), the state at time t is given by x(t) = M (θ 0 ) t x(0). Its behavior is determined at first order by the largest eigenvalue of matrix M = M (θ 0 ), λ 1 (M ), x(t) being close to cλ 1 (M ) t u, where u is the eigenvector of M associated with eigenvalue λ 1 (M ) and c some constant. The exponential rate of growth λ 1 (M ) of observable x H (t) can be estimated, e.g. by fitting a regression line to the semilog plot up to time t of the corresponding time series, thus producing the estimateλ(t).
For a single population, in the mean field model, estimation of the exponent λ 1 (M ) gives one relation constraining the unknown parameters, namely:
Remark 2.1. Let p be a polynomial such that p(M ) = 0 (p could be the characteristic polynomial of M for instance). Writing
In other words, the observables x H (t) follow an auto-regressive dynamics. Instead of estimating unknown coefficients of matrix M , one could therefore estimate directly the auto-regressive coefficients a i . The latter approach is potentially simpler, and suffices if one is only interested in forecasting x H (t). The coefficients a i however do not have a direct physical interpretation, in contrast to the coefficients of M . Estimation of M is therefore better suited to assess the impact of specific measures.
If on given day t 0 (start of confinement say), the parameters α i and α a abruptly change, the dynamics from t 0 onwards would again be linear, with a modified matrix M (θ 1 ) with θ 1 = (α i,1 , α a,1 ) corresponding to the new contact rates. This suggests the following approach for inferring these parameters. Estimate the state vector x(t 0 ) as Cu, for some constant C = C(θ 0 ) such that Cu H = x H (t 0 ), where u H is the projection of u = u(θ 0 ) on the first h vectors of the canonical basis.
Leverage then the transient behaviour of the post-shock dynamics as follows. Making the dependency of the post-shock matrix M (θ 1 ) on the unknown θ 1 explicit, estimate the unknown (θ 0 , θ 1 ) as the minimizers of the optimization problem:
Remark 2.2. The above optimization criterion is essentially motivated by the assumption of constant contact parameters before and after the shock. This is a strong assumption, verified at best in an approximate sense. See [10] which considers multiple phases in an SIR dynamics, to be fitted to observations. In Section 3 we will consider two shocks. Time-varying contact rates will be discussed below in Section 4.
We just saw how to leverage Perron's theory to perform inference in the mean-field version of the model. Let us now describe an inference method in the probabilistic model which -as already explained -is more relevant when dealing with small populations, e.g. when studying sub-critical evolutions in a decay phase, or spread relaunch phases.
In the following paragraph we describe an inference method in the probabilistic setting that enables to estimate the state X(t) (including hidden trajectories) by updating the estimation at every new time step. The vector of average numbers per type at generation t, conditionally on the initial state to be X(0) ∈ R J , is M t X(0).
First, we describe a model for measures, that is the observed data X * H , as follows:
where P H is the h × J matrix of the projection on the h last components of the canonical basis, and v(t) is a centered noise, with covariance matrix R(t), representing the confidence in measurements.
Model errors Now, we introduce the filtration F t := σ (X * H (0), . . . , X * H (t)), that is the σ−algebra generated by all measurements up to time t. We provide estimates of states X(t), as well as the estimated errors. We will use the following notations:
•X(t): an estimate of state X(t) given F t , and P (t) its covariance matrix,
•X(t|t − 1): an estimate of X(t), given F t−1 , and P (t|t − 1) its covariance matrix.
At time t, we write
with u(t) = X(t) − M X(t − 1) a centered vector. For j ∈ [J], let S i denote the covariance matrix of vector of children from type i-parent. The covariance matrix Q(t) of u(t) can be computed as follows:
The fact that closed-form expressions for second moments are available is a consequence of the Markovian structure of our model. The same Markovian structure also implies that u(t) is decorrelated from any F t−1 −measurable variable.
The strategy is then to computeX(t|t − 1) at each step, and to adjust the prediction with the measurement X * H (t) to buildX(t), assuming thatX(t) is of the form
Optimizing K(t) in the quadratic error sense leads to Kalman filtering (see e.g. [6] ), which is the best linear recursive estimator in our setting, namely a hidden Markov model. With this method, at each time t, we perform the following predictions:
then we adjust our estimates, updating with the new measure X * H (t):
The intermediate matrix K(t) of size J × h is usually called the Kalman gain. A possible initial condition consists for instance in setting X(0) = e i , that is one initial individual in a given state i ∈ [J], with P (0) = 0 J , the J × J null matrix.
The above discussion shows how to leverage the model structure to provide an alternative approach to maximum likelihood estimation approaches.
Such techniques can be used in a single population scenario, as we illustrate in the next section. They could however become prohibitive for the multi-population scenarios we shall consider later, and for which the Kalman filtering approach may be an appealing alternative.
In this section, we illustrate the probabilistic version of our model by fitting it to daily hospitalization data, and show how to use it to forecast future numbers of hospitalizations.
The population we will consider is that of department 75 (Paris, France). We still assume that the only observed trajectory is X * H (t), the number of new hospitalizations, all other trajectories being hidden. A noisy version of this data, namely the number of new COVID-19 hospitalizations every day in Paris, can be found in public data provided by SurSaUD syndromic surveillance system.
Though Kalman filtering approach -as described in Section 2 -would be natural here, we take a more direct, brute-force method, using Monte-Carlo simulation to infer the parameters.
For simulations of our probabilistic model, we take p i = 0.7, p h = 0.05 (see Transition probabilities in Section 1). We initialize the time range on February 8th, which is 18 days before the first local SARS-CoV2 death case in France. We focus on the pre-confinement and confinement phases, so the fitting is made up to May 5th, on a total of 72 days for the train set. The period from May 6th to May 12th is kept as the test set for prediction (see below).
Recall that for contact parameters, we assume α i = α a = α. Furthermore we assume that it is piecewise constant on the time range, as follows:
• between day 0 (February 8th) and day t 2 , α = α 1 .
• between day t 2 and day t 3 , α = α 2 .
• between day t 3 and day T (May 5th), α = α 3 .
The three phases can be identified to the pre-confinement phase, the confinement phase I and confinement phase II, each phase being the result of several changes: mainly political decisions, but also work habits, observance of the rules, etc. In order to get a fitting as realistic as possible, the boundaries of these three phases are also inferred.
In total, six parameters (X E (0), α 1 , t 2 , α 2 , t 3 , α 3 ) are tuned. X E (0) is the initial condition, that is the initial number of exposed individuals. This quantity is then divided and allocated to all variables X E,d (0) according to distribution p E described in Section 1.
Given (X E (0), α 1 , t 2 , α 2 , t 3 , α 3 ), we sample several independent trajectories of X(t). Now we consider two choices for the loss function.
We first define the L 1 loss function as follows:
where X H is the random trajectory under model 1, and X * H is the observed data.
Second, since we expect exponential growth for the variables of interest (see e.g. 1.2 or 2.1), we shall also consider the L 1 -log loss function:
For both choices, we aim to minimize the loss function on our train set. To do so, we perform a naive grid search in six dimensions, computing at each step an empirical expectation of the loss for the chosen parameters over sampled trajectories, in order to compute approximate solutions
for Loss ∈ L, L .
Prediction After having estimated all the parameters, we use our model to predict the evolution of trajectories in a future time range. In our example, this time range runs from T + 1 (May 6th) to T pred (May 12th), that is a period of one week.
The performance of our predicition (over our test set) is evaluated in terms of the mean L 1 norm
where the probability P fit is defined for previously inferred parameters (X E (0),α 1 ,t 2 ,α 2 ,t 3 ,α 3 ) in equation (16). We chose as our test criterion in (17) absolute deviation rather than deviation between logarithms because this corresponds to an error in the absolute numbers of hospital entrances, which has arguably a clearer operational meaning.
The inferred parameters with SurSaUD data, together with the training and test errors, are summed up in Table 1 . Corresponding curves of inferred parameters, predicted hidden variables and hospital entrances are given in Figure 3 and Table 1 : Approximate values of inferred parameters, for department 75, approximated with the probabilistic model, with SurSaUD data. All errors are given with their 95% confidence interval obtained by Monte Carlo simulations.
Several remarks can be made from these results:
• The L pred prediction error is of order 4, which is to be compared to the typical diameter of a 95% confidence interval in the prediction time range in Figure 3 and Figure 4 , which is about 10. We measure a prediction error which is less than half a typical uncertainty on a sampled stochastic trajectory X H (t) with fitted parameters. This supports the claim that a prediction error of order 4 is the best possible with our model.
• For both loss functions, the optimal transition points t 1 and t 2 are identified on March 15th and March 20th. Some other strategies could have been used to identify phases, e.g. by fitting a piecewise linear regression model on the observations log X * H (t), and reporting them earlier on the contact rates, according to the distributions of phase durations, as done in recent work by Gaubert et al. [5] . For the Paris population, the major political decisions and events occurred on March 12th (closure of all schools and universities), March 14th (closure of all non essential public places), March 15th (municipal elections) and March 17th (lockdown). These phases do not necessarily exactly coincide with the transition points, this can be caused e.g. by assumptions made for time parameters (see Section 1), and delays due to observance of the rules. In addition, we fixed three phases, but this number is arbitrary and could be discussed (see e.g. [10] ).
• From the figures, the L 1 -log Loss function L seems to give better results for the fitting part (which is not surprising since we are in log scale), but this is also true for prediction, as presented in Table 1 . This can justify the use of the L 1 -log Loss function L in future work.
• These results, though based on noisy data and perfectible with other datasets or broader inference, validate our model: with few parameters, the evolution of epidemics can be described in a satisfactory way: with the L 1 -log Loss function for instance, the mean relative error is ∼ 12.1% (in logarithmic scale) for the trajectory of X H (t).
As shown in this section, the model of Figure 1 can be useful to describe the evolution of crucial variables such as the number of daily hospitalized individuals, but also show hidden trajectories, and predict hospital load on future days.
We considered piecewise constant contact rates in Section 2 for inference purposes, and performed corresponding numerical illustration in the present Section. We shall now consider the case of contact rates that are no longer piecewise constant, being potentially dependent not just on "shocks" such as the start of confinement, but also on other covariates. For instance, regular day-to-day mobility or commuting flows of population could also impact contact rates, as we now investigate.
Contact rates are now assumed to depend on piece-wise constant rates depending on specific temporal phases, as described before, and daily variables that are responsible for additional variability. Among these daily variables, some are related to user mobility, and in particular on the daily outflow, that is the number of day-return trips outside the region of interest.
The model we advocate is then exactly as in Section 1, except that contact rates α i , α a are now time-dependent.
Assuming that the daily outflows f (t) (assumed to be normalized so as to have zero mean and standard deviation 1) are known, the contact rates α a , α i are now mobility-dependent in the following way:
In the above, φ(t) is a discrete phase (e.g. φ(t) ∈ {1, 2, 3} as in Section 3), and F i , F a are fixed functions. These could be taken as logistic functions in the variable f (t); alternatively, one might take Figure 4 : Best fitting of our probabilistic model with SurSaUD data for Paris with log-L 1 cost function, with prediction, some hidden trajectories, and variation of contact rates α a = α i .
With this formulation, variables α i,φ(t) (resp. α a,φ(t) ) is the value of α i (resp. α a ) when the commuting flow is at its equilibrium, and γ i , γ a are non-negative parameters. Figure 5 below gives illustrates the evolution of the contact rates for Paris based on Equation (19), where the observables f (t) are obtained from SFR mobile operator data. The shape of the resulting contact rates reflects the discrete phases as before, together with weekly periodicity patterns. Contact rates are lower on Saturdays and Sundays, when daily commuting trips are less numerous, as expected.
Remark 4.1. In its mean-field version, the model with variable contact rates is
for fixed matrices M 0 , . . . , M K and scalars α k (t) that depend on mobility covariates relative to time t. Assume for tractability that the covariates α k (t) are small, so that heuristically we have the first order approximation, for τ ∈ [d],
Following on Remark 2.1, let p(z) = z d − d i=1 a i z d−1 be a degree-d polynomial such that p(M 0 ) = 0. We then have the first-order expansion in the α(t):
Provided that for all i = 0, . . . , d − 1 and all k ∈ [K], there exist matrices N k,i and polynomials
, then up to terms of first order one has
Plugged into the previous display, this entails that the observations x H (t) satisfy an auto-regressive relation with coefficients involving the covariates α, of general form:
x
where P i,j are constant matrices. A possible approach could then consist, for prediction purposes, in fitting such an auto-regressive relation while learning suitable functions α(t) of available covariates.
Besides commuting outflows mentioned above, another aspect of mobility that plays an important role is that of routing, that is change of residence location: this happens during holiday periods, but has also happened at the start of lockdown, where a significant exodus from dense urban areas towards countryside has been observed.
We now describe an extension of our model that captures such routing.
To introduce routing, we must first consider an extension of our basic model to cover multiple cohorts, or sub-populations. For the sake of readability we only describe the generalized model in its mean-field version. However a probabilistic version could easily be described along the same lines as in the single population case. We thus introduce a subdivision of the global population into cohorts c ∈ C. These could for instance be given by a couple (r, a) characterizing geographical regions r and age ranges a.
A potential alternative to the above definition of cohort c is to define sub-population according to a given age range a, the region r in which the corresponding persons usually live, and the region r in which they slept the previous night. In that case, we have c = (r, r , a). This notion of cohort is particularly important to describe our routing model, and account for viral dissemination through regions.
We denote by N c (t) the size of cohort c at time t.
The model is as follows. Each cohort c ∈ C has its own biological parameters θ c (encapsulating e.g. contact rates, distributions of phase durations, etc.). A natural assumption is to take these to depend only on age range. Per age parameters could be obtained from hospital-collected statistics.
Routing dynamics of N c (t) For a given cohort c ∈ C, the time evolution of N c (t) is the result of population flows. We represent such flows in the following manner: we have
where R c ,c (t) is the fraction of individuals counted in N c (t) on day t, that have migrated to be counted in N c (t + 1) on day t + 1. The case c = c of sedentary individuals (no migration) is also taken into account. New external arrivals in c (say, from a foreign country) are accounted for by E c (t + 1). We also assume that for all cohort c , c∈C R c ,c (t) ≤ 1,
allowing for departures outside the considered population (say, abroad) when the sum is strictly less than 1.
Daily activity, contact between sub-populations On day t, a contact intensity factor n c,c (t) gives the average number of individuals from cohort c that a typical individual from cohort c will encounter. It is natural to assume that N c (t)n c,c (t) = N c (t)n c ,c (t)
Parameter q c,c represents the probability that upon a contact between individuals from population c and population c , infection gets propagated from that in c to that in c . It is natural to assume that q c,c only depends on the age ranges of c and c . We then let α c,c (t) := q c,c n c,c (t),
the infection rate from c to c on day t. For simplicity we do not distinguish between asymptomatic/prodromic (previously α a ) and symptomatic individuals (previously α i ) here. Note that the above dynamics, and in particular Equations (23), (24) generalize the mobilitydependent contact rates of Section 4 to the multi-subpopulation case.
Each sub-population c on day t is characterized by state vector x c (t), given by concatenation of vectors x c E , x c P , x c I1 , x c A , x c I2 , x c H . As before this is a vector of size 5h + 1, h being the maximal number of days each phase lasts.
To specify the transition between day t and day t + 1 we consider the following two-stage mechanism: a first stage consists in an operation similar to the previously considered dynamics, extended to account for multiple cohorts c, giving vector y(t) as a linear transform of vector x(t). A second stage representing the routing of individuals between cohorts gives x(t + 1) as a linear transform of y(t). Specifically we have 
where the parameters β z are to be inferred, and capture the density of contacts to be expected in (time-space) location z. In this formulation, the parameters to be inferred are: parameters β z , and infection probabilities q c,c . This model is over-parameterized, since multiplying the β z by η and dividing the q c,c by η leaves the model unchanged for all η = 0. This can be easily solved, forcing for instance the parameters β z to have mean 1. While this over-parameterization is easily circumvented, it remains a challenge to effectively fit all the free parameters of the above model, as well as to define the proper spatio-temporal granularity associated with the time-space locations z.
In this study we laid the foundations of a general model of COVID-19 epidemics which captures user mobility and contact tracing. We described it as a multi-type branching process, proposed inference methods to perform short-term prediction, that we illustrated on the Paris hospitalization data. We also introduced extensions to capture contract tracing and case isolation, and finally proposed extensions to capture user mobility, distinguishing between commuting and routing mobility. Future work will further exploit this model together with hospital incidence and mobility data to forecast epidemic progress and assess impact of mobility on infectious contacts.