A parsimonious model for spatial transmission and heterogeneity in the COVID-19 propagation
Lionel Roques, Olivier Bonnefon, Virgile Baudrot, Samuel Soubeyrand, Henri Berestycki
AA parsimonious model for spatial transmission andheterogeneity in the COVID-19 propagation
Lionel Roques a,* , Olivier Bonnefon a , Virgile Baudrot a ,Samuel Soubeyrand a and Henri Berestycki b a INRAE, BioSP, 84914, Avignon, Franceb EHESS, CNRS, CAMS, Francec Senior Visiting fellow, HKUST Jockey Club Institute for Advanced Study, Hong Kong University of Science and Technology* Contact : [email protected] July 21, 2020
Abstract
Raw data on the cumulative number of deaths at a country level generally indicate a spa-tially variable distribution of the incidence of COVID-19 disease. An important issue is todetermine whether this spatial pattern is a consequence of environmental heterogeneities, suchas the climatic conditions, during the course of the outbreak. Another fundamental issue isto understand the spatial spreading of COVID-19. To address these questions, we considerfour candidate epidemiological models with varying complexity in terms of initial conditions,contact rates and non-local transmissions, and we fit them to French mortality data with amixed probabilistic-ODE approach. Using standard statistical criteria, we select the modelwith non-local transmission corresponding to a diffusion on the graph of counties that dependson the geographic proximity, with time-dependent contact rate and spatially constant param-eters. This original spatially parsimonious model suggests that in a geographically middlesize centralized country such as France, once the epidemic is established, the effect of globalprocesses such as restriction policies, sanitary measures and social distancing overwhelms theeffect of local factors. Additionally, this modeling approach reveals the latent epidemiologicaldynamics including the local level of immunity, and allows us to evaluate the role of non-localinteractions on the future spread of the disease. In view of its theoretical and numerical sim-plicity and its ability to accurately track the COVID-19 epidemic curves, the framework wedevelop here, in particular the non-local model and the associated estimation procedure, is ofgeneral interest in studying spatial dynamics of epidemics.
Keywords.
COVID-19 | Spatial diffusion | Mechanistic-statistical model | Non-local transmission | Immunity rate
In France, the first cases of COVID-19 epidemic have been reported on 24 January 2020, althoughit appeared latter that some cases were already present in December 2019 [1]. Since then, the first1 a r X i v : . [ q - b i o . P E ] J u l mportant clusters were observed in February in the Grand Est region and the Paris region. A fewmonths later, at the beginning of June, the spatial pattern of the disease spread seems to have kepttrack of these first introductions. As this spatial pattern may also be correlated with covariatessuch as climate [2] (see also SI 1), a fundamental question is to assess whether this pattern is theconsequence of a heterogeneous distribution of some covariates or if it can be explained by theheterogeneity of the initial introduction points. In the latter case, we want to know if the epidemicdynamics simply reflects this initial heterogeneity and can be modeled without taking into accountany spatial heterogeneity in the local conditions.SIR epidemiological models and their extensions have been proposed to study the spread of theCOVID-19 epidemic at the country or state scale (e.g., [3] in France and [4] in three US states),and at the regional scale ([5, 6] in China, [7] in France and [8] in Italy). In all cases, a differentset of parameters has been estimated for each considered region/province. One of the main goalsof our study is to check whether the local mortality data at a thin spatial scale can still be wellexplained by a single set of parameters, at the country scale, but spatially varying initial conditions.In this study we consider SIR models with time-dependent contact rate, to track changes over timein the dynamic reproductive number as in the branching process considered in [4]. Using standardstatistical criteria, we compare four types of models, whose parameters are either global at thecountry scale or spatially heterogeneous, and with non-local transmission or not. We work withmortality data only, as these data appear to be more reliable and less dependent on local testingstrategies than confirmed cases, in settings where cause of death is accurately determined [9].Additionally, it was already shown that for this type of data, SIR models outperform other classesof models (an SEIR models and a branching process) in the three US states considered in [4].The approach we develop here is in line with the general principle of using parsimonious modelseloquently emphasized in [4].Another important issue that such models may help to address is the quantification of therelative effects of restrictions on inter-regional travels, vs reductions in the probability of infectionper contact at the local scale. France went into lockdown on March 17, 2020, which was foundto be very effective in reducing the spread of the disease. It divided the effective reproductionnumber (number R t of secondary cases generated by an infectious individual) by a factor 5 to7 at the country scale by May 11 [10, 7]. This is to be compared with the estimate of the basicreproduction number R carried out in France at the early onset of the epidemic, before the countrywent into lockdown (with values R = 2 .
96 in [7] and R = 3 . vs a healthy individual, leading to a factor x8 between the number of confirmed cases and2he actual number of cases before the lockdown [3], and x20 at the end of the lockdown period [10].This type of framework – often referred to as ‘mechanistic-statistical modeling’ – aims at connectingthe solution of continuous state models such as differential equations with complex data, such asnoisy discrete data, and identifying latent processes such as the epidemiological process underconsideration here. Initially introduced for physical models and data [12], it is becoming standardin ecology [13].Using this framework the objectives of our study are (i) to assess whether the spatial patternobserved in France is due do some local covariates or is simply the consequence of the heterogeneityin the initial conditions together with global processes at the country scale; (ii) to evaluate the roleof non-local interactions on the spread of the disease; (iii) to propose a tool for real-time monitoringthe main components of the disease in France, with a particular focus on the local level of immunity. Materials and Methods
Data
Mainland France (excluding Corsica island) is made of 94 counties called ’d´epartements’. The dailynumber of hospital deaths – excluding nursing homes – at the county scale are available from SantPublique France since 18 March 2020 (and available as Supplementary Material). The daily numberof observed deaths (still excluding nursing homes) in county k during day t is denoted by ˆ µ k,t .We denote by [ t i , t f ] the observation period and by n d the number of considered counties. Toavoid a too large number of counties with 0 deaths at initial time, the observation period rangesfrom t i = March 30 to t f = June 11, corresponding to n t = 74 days of observation. All the countiesof mainland France (excluding Corsica) are taken into account, but the Ile-de-France region, whichis made of 8 counties with a small area is considered as a single geographic entity. This leads to n d = 87. Mechanistic-statistical framework
The mechanistic-statistical framework is a combination of a mechanistic model that describes theepidemiological process, a probabilistic observation model and a statistical inference procedure.We begin with the description of four mechanistic models, whose main characteristics are given inTable 1.
Mechanistic modelsModel M : SIR model for the whole country. The first model is the standard mean fieldSIRD model that was used in [3, 10]: S (cid:48) ( t ) = − α ( t ) N S I,I (cid:48) ( t ) = α ( t ) N S I − ( β + γ ) I,R (cid:48) ( t ) = β I,D (cid:48) ( t ) = γ I, (1)3eterog. Heterog. Intercounty Nb. parametersinitial data contact rate transmission M no no no n t M yes no no n t M yes yes no n d × n t M yes no yes n t + 2Table 1: Main characteristics of the four models. The quantity n t = 74 corresponds to the numberof days of the observation period and n d = 87 corresponds to the number of administrative units.with S the susceptible population, I the infectious population, R the recovered population, D the number of deaths due to the epidemic and N the total population, in the whole country.For simplicity, we assume that N is constant, equal to the current population in France, therebyneglecting the effect of the small variations of the population on the coefficient α ( t ) /N . Theparameter α ( t ) is the contact rate (to be estimated) and 1 /β is the mean time until an infectiousbecomes recovered. The results in [14] show that the median period of viral shedding is 20 days,but the infectiousness tends to decay before the end of this period: the results in [15] indicate thatinfectiousness starts from 2.5 days before symptom onset and declines within 7 days of illness onset.Based on these observations we assume here that 1 /β = 10 days. The parameter γ corresponds tothe death rate of the infectious. It was estimated independently in [3, 10] and [7], leading to a value γ = 5 · − . This value only takes into account the deaths at hospital, and is therefore consistentwith the data that we used here. Model M : SIR model at the county scale with globally constant contact rate andno spatial transmission. The model M is applied at the scale of each county k , leading tocompartments S k , I k , R k , D k that satisfy an equation of the form (1), with N replaced by N k , thetotal population in the county k . In this approach, the contact rate α ( t ) is assumed to be the samein all of the counties. Model M : SIR model at the county scale with spatially heterogeneous contact rateand no spatial transmission. With this approach, the model M is extended by assuming thatthe contact rate α k ( t ) depends on the considered county. Model M : County scale model with globally constant contact rate and spatial trans-mission. The model M is extended to take into account disease transmission events betweenthe counties: S (cid:48) k ( t ) = − ρ ( t ) N k S k n d (cid:88) j =1 w j,k I j ,I (cid:48) k ( t ) = ρ ( t ) N k S k n d (cid:88) j =1 w j,k I j − ( β + γ ) I k ,R (cid:48) k ( t ) = β I k ,D (cid:48) k ( t ) = γ I k . (2)4he weights w j,k describe the dependence of the contagion rates with respect to the distance betweencounties. We assume a power law decay with the distance: w j,k = 11 + (dist( j, k ) /d ) δ , (3)with dist( j, k ) the geographic distance (in km) between the centroids of counties j and k , d >
0a proximity scale, and δ >
0. Thus, the model involves two new global parameters, d and δ ,compared to model M . This model extends the Kermack-McKendrick SIR model to take intoaccount non-local spatial interactions. It was introduced by Kendall [16] in continuous variables.The model we adopt here is inspired from the study of [17] in a different context where the sametypes of weights have been used. We thus take into account diffusion on the weighted graph ofcounties in France. This amounts to considering that individuals in a given county are infectedby individuals form other counties with a probability that decreases with distance as a power law,in addition to contagious individuals from their own county. This dependence of social spatialinteractions with respect to the distance is supported, notably, by [18] that analyzed the short-timedispersal of bank notes in the US. We also refer to [19] for a thorough discussion on the variousapplications of power law dispersal kernels since they were introduced by Pareto [20].With this non-local contagion model, in contradistinction to epidemiological models with dis-persion such as reaction-diffusion epidemiological models [21], the movements of the individuals arenot modeled explicitly. The model implicitly assumes that infectious individuals may transmit thedisease to susceptible individuals in other counties, but eventually return to their county of origin.This has the advantage of avoiding unrealistic changes in the global population density. Observation model
We denote by D k ( t ) the expected cumulative number of deaths given by the model, in county k . With the mean-field model M , we assume that it is proportional to the population size: D k ( t ) = D ( t ) N k /N, with N k the population in county k and N the total French population. Withmodels M , M and M , we simply have D k ( t ) = D k ( t ). The expected daily increment in thenumber of deaths given by the models in a county k is D k ( t ) − D k ( t − µ k,t in county k follows a Poisson distribution with mean value D k ( t ) − D k ( t − µ k,t ∼ Poisson( D k ( t ) − D k ( t − . (4)Note that the time t in the mechanistic models is a continuous variable, while the observations ˆ µ k,t are reported at discrete times. For the sake of simplicity, we used the same notation t for the daysin both the discrete and continuous cases. In the formula (4) D k ( t ) (resp. D k ( t − t (resp. t − Initial conditions
In models M , M , M , at initial time t i , we assume that the number of susceptible cases is equalto the number of inhabitants in county k : S k ( t i ) = N k , the number of recovered is R ( t i ) = 0 andthe number of deaths is given by the data: D k ( t i ) = ˆ µ k,t i . To initialise the number of infectious,we use the equation D (cid:48) ( t ) = γ I ( t ), and we define I ( t i ) as 1 /γ × (mean number of deaths over the5eriod ranging from t i to 29 days after t i ): I k ( t i ) = 1 γ (cid:88) s = t i ,...,t i +29 ˆ µ k,s . (5)The 30-days window was chosen such that there was at least one infectious case in each county.In model M , the initial conditions are obtained by adding the initial conditions of model M (orequivalently, M ) over all the counties. Statistical inference
Real-time monitoring of the parameters and data assimilation procedure.
To smoothout the effect of small variations in the data, and to avoid identifiability issues due to the largenumber of parameters, while keeping the temporal dependence of the parameters, the parameters α ( t ) and α k ( t ) of the ODE models M , M , M are estimated by fitting auxiliary problems withtime-constant parameters over moving windows ( t − τ / , t + τ /
2) of fixed duration equal to τ days.These auxiliary problems are denoted respectively by ˜ M ,t , ˜ M ,t , and ˜ M ,t (see SI 2 for a preciseformulation of these problems). The initial conditions associated with this system, at the date t − τ / M , M and M , respectively. Inference procedure.
For simplicity, in all cases, we denote by f D k , ˆ µ k ( s ) := ( D k ( s ) − D k ( s − ˆ µ k,s ˆ µ k,s ! e − ( D k ( s ) − D k ( s − the probability mass function associated with the observation process (4) at date s in county k ,given the expected cumulative number of deaths D k given by the considered model in county k .In models ˜ M ,t and ˜ M ,t , the estimated parameter is ˜ α . The likelihood L is defined as theprobability of the observations (here, the increments { ˆ µ k,s } ) conditionally on the parameter. Usingthe assumption that the increments ˆ µ k,s are independent conditionally on the underlying SIRDprocess ˜ M ,t (resp. ˜ M ,t ), we get: L (˜ α ) := P ( { ˆ µ k,s , k = 1 , . . . , n d , s = t − τ / , . . . , t + τ / }| ˜ α )= n d (cid:89) k =1 t + τ/ (cid:89) s = t − τ/ f D k , ˆ µ k ( s ) . We denote by ˜ α ∗ t the corresponding maximum likelihood estimator, and we set α ( t ) = ˜ α ∗ t in model M (resp. M ) . For model ˜ M ,t , the inference of the parameters ˜ α k is carried out independently in each county,leading to the likelihoods: L k (˜ α k ) := P ( { ˆ µ k,s , s = t − τ / , . . . , t + τ / }| ˜ α k )= t + τ/ (cid:89) s = t − τ/ f D k , ˆ µ k ( s ) . We denote by ˜ α ∗ k,t the corresponding maximum likelihood estimator, and we set α k ( t ) = ˜ α ∗ k,t inmodel M . M , we apply a two-stage estimation approach. We first use the estimate obtainedwith model M by setting ρ ( t ) = C α ( t ) , where α ( t ) is the estimated contact rate of model M and C is a constant (to be estimated; note that estimating ρ ( t ) means estimating n t parameters).Thus, given α ( t ), the only parameters to be estimated are the constant C , the proximity scale d and the exponent δ . They are estimated by maximizing: L ( C, d , δ ) := P ( { ˆ µ k,s , k = 1 , . . . , n d , s = t i , t f }| C, d , δ )= n d (cid:89) k =1 t f (cid:89) s = t i f D k , ˆ µ k ( s ) . Model selection.
We use the Akaike information criterion (AIC) [22] and the Bayesian infor-mation criterion (BIC) [23] to compare the models. For both criteria, we need to compute thelikelihood function associated with the model, with the parameters determined with the aboveinference procedure: (cid:32)L( M m ) = n d (cid:89) k =1 t f (cid:89) s = t i f D k , ˆ µ k ( s ) , (6)with m = 0 , , , D k the expected (cumulative) number of deaths in county k given by model( M m ). Given the number of parameters (cid:93) M m estimated in model M m , the AIC score is definedas follows: AIC ( M m ) = 2 (cid:93) M m − M m )) , and the BIC score: BIC ( M m ) = (cid:93) M m ln( K ) − M m )) , with K = n d × n t = 6 438 the number of data points. Numerical methods . To find the maximum likelihood estimator, we used a BFGS constrainedminimization algorithm, applied to − ln((cid:32)L) via the Matlab ® function fmincon . The ODEs weresolved thanks to a standard numerical algorithm, using Matlab ® ode45 solver. The Matlab ® codesare available as Supplementary Material. Results
Model fit and model selection.
To assess model fit, we compared the daily increments in the numberof deaths given by each of the four models with the data. For each model, we used the parameterscorresponding to the maximum likelihood estimators. The comparisons are carried out at theregional scale: mainland France is divided into 12 administrative regions, each of which is made ofseveral counties. The results are presented in Fig. 1. In all the regions, the models M , M and M lead to a satisfactory visual fit of the data, whereas the mean field model M does not manageto reproduce the variability of the dynamics among the regions.The log-likelihood, AIC and BIC values are given in Table 2. Models M , M , M lead tosignificantly higher likelihood values than model M . This reflects the better fit obtained with thesethree models, compared to model M and shows the importance of taking into account the spatial7odel AIC BIC Log-likelihood ∆AIC M . · . · − . · − . · M . · . · − . · − M . · . · − . · − . · M . · . · − . · M ).heterogeneities in the initial densities of infectious cases. On the other hand, the log-likelihood,though higher with model M is close to that obtained with M , and the model selection criteriaare both strongly in favor of model M . This shows that the spatial heterogeneity in the contactrate does not have a significant effect on the epidemic dynamics within mainland France.Model M with spatial transmission leads to an intermediate likelihood value, between thoseof models M and M , with only 2 additional parameters with respect to model M . As aconsequence, the model selection criteria exhibit a strong evidence in favor of the selection of model M . This means a large part of the difference between models M and M can be captured bytaking into account the spatial transmission, which therefore seems to have a significant effect onthe epidemic dynamics.As a byproduct of the estimation of the parameter α ( t ) (resp. α k ( t )) of model M (resp. M ),we get an estimate of the effective reproduction number in each county, which is given by theformula [24]: R kt = α ( t ) β + γ S k ( t ) N k , so that I (cid:48) k ( t ) < R kt <
1, whereas I (cid:48) k ( t ) > R kt > R kt obtained with model M are depicted in Fig. 2, clearly showing a decline in R kt , as alreadyobserved in [10, 7], but there the computation was at a fixed date.For model M , the maximum likelihood estimation gives C = 0 . d = 2 .
16 km and δ = 1 . d , indicatesthat non-local contagion plays a secondary role compared to within-county contagion: the minimumdistance between two counties is 36 km, leading to a weight of 5 . / C is significantlysmaller than 1 (recall that the contact rate in model M is ρ ( t ) = C α ( t ) with α ( t ) the contact ratein model M ) shows that the non-local contagion term plays an important role on the spreadingof the epidemic. Immunity rate.
Using model M , which leads to the best fit, we derive the number of recoveredindividuals (considered here as immune) at a date t , in each county, and the immunity rate R k /N k .It is presented at time t f in Fig. 3. The full timeline of the dynamics of immunity obtained withmodel M since the beginning of April is available as Supplementary Material (see SI 1). Limiting movement vs limiting the probability of transmission per contact.
Before the lockdown,the basic reproduction number R in France was about 3 [7, 3], and was then reduced by a factor8igure 1: Model fit at the regional scale. The data have been smoothed (moving average over 15days), for easier graphical comparison with the models.9igure 2: Dynamics of the effective reproduction rate R kt given by model M over the 87 consideredcounties.Figure 3: Estimated immunity rate, at the county scale, by June 11, 2020, using model M .10 to 7, leading to values around 0.5 (see [10, 7] and Fig. 2). This corresponds to a contact rate α ( t ) ≈ β R t ≈ . α ( t ) ≈ .
05 after the lockdown in model M . Let usnow consider a hypothetical scenario of a new outbreak with an effective reproduction number thatrises again to reach values above 1. Due to the higher awareness of the population with respectto epidemic diseases and the new sanitary behaviors induced by the first COVID-19 wave, thereproduction number will probably not reach again values as high as 3.The new outbreak scenario is described as follows: we start from the state of the epidemicat June 11, and we assume a ’local’ contact rate ρ ( t ) jumping to the value 0 .
11 in model M (corresponding to a 2-fold increase compared to the previous 30 days). In parallel, to describe thelifting of restrictions on individual movements, we set d = 20 km for the proximity scale in model M . This new outbreak runs during 10 days, and then, we test four strategies:- Strategy 1: no restriction. The parameters remain unchanged: ρ ( t ) = 0 .
11 and d = 20 km;- Strategy 2: restriction on intercounty movement. The parameter ρ ( t ) = 0 .
11 is unchanged,but d = 2 .
16 km, corresponding to its estimated value during the period ( t i , t f );- Strategy 3: reduction of the contact rate within each county (e.g., by wearing masks), but norestriction on intercounty movement: ρ ( t ) = 0 .
05 and d = 20 km;- Strategy 4: reduction of the contact rate within each county and restriction on intercountymovement: ρ ( t ) = 0 .
05 and d = 2 .
16 km.The daily number of deaths corresponding to each scenario is presented in Fig. 4. We estimate ineach case a value of the effective reproduction number R t over the whole country by fitting the globalnumber of infectious cases with an exponential function over the last 30 days. As expected, the morerestrictive the strategy, the less the number of deaths. After 30 days, the cumulative number ofdeaths with the first strategy is 17 271, and R t ≈
2. Restriction on intercounty movement (strategy2) leads to a 81% decrease in the cumulative number of deaths (3 281 deaths) and R t ≈ .
2; reducingthe contact rate within each county leads to a 88% decrease (strategy 3, 2 139 deaths) and R t ≈ . R t ≈ . Discussion
We show here that a parsimonious model can reproduce the local dynamics of the COVID-19epidemic in France with a relatively high goodness of fit. This is achieved despite the spatial het-erogeneity across French counties of some environmental factors potentially influencing the diseasepropagation. Indeed, our model only involves the initial spatial distribution of the infectious casesand spatially-homogeneous (i.e. countrywide) parameters. For instance, the mean temperatureduring the considered period ranges from 12 . ° C to 18 . ° C depending on the region. We do observea negative correlation between the mean temperature and the immunity rate (see Fig. S1), howeverit does not reflect a causality. Actually, our study shows that if there is an effect of local covariatessuch as the mean temperature on the spread of the disease once it emerged, its effect is of lowerimportance compared to the global processes at work at the country scale. Such covariates mightplay a major role in the emergence of the disease, but our work focuses on the disease dynamicsafter the emergence. 11igure 4: Daily number of deaths due to a new outbreak in logarithmic scale; comparison betweenfour management strategies. The number of deaths is computed over the whole country.12ence, we find that initial conditions and spatial diffusion are the main drivers of the spatialpattern of the COVID-19 epidemic. This result may rely on specific circumstances: e.g., mainlandFrance covers a relatively middle-sized area, with mixed urban and countryside populations acrossthe territory, a relatively homogeneous population age distribution, and a high level of centralismfor public decision (in particular regarding the disease-control strategy). Of course, these featuresare not universal. In other countries with more socio-environmental diversity within which envi-ronmental drivers and state decentralization could significantly induce spatial variations in diseasespread and, consequently, in which countrywide parameters would not be appropriate. Moreover, atthe global scale, the COVID-19 dynamics in different countries seem highly contrasted as illustratedby the data of Johns Hopkins University Center for Systems Science and Engineering [25] —seealso http://covid19-forecast.biosp.org/ [26]— and could probably not be explained with aunique time-varying contact rate parameter. Nonetheless, this model could be adapted to manyother situations at an appropriate geographical level, with a single well defined political decisionprocess.Herd immunity requires that a fraction 1 − / R ≈
70% of the population has been infected. Itis far from being reached at the country scale in France, but we observe that the fraction of immuneindividuals strongly varies across the territory, with possible local immunity effects. For instance,in the most impacted county the immunity rate is 16%, whereas it is less than 1% in less affectedcounties. At a thinner grain scale, even higher rates may be observed, for instance, by April 4 theproportion of people with confirmed SARS-CoV-2 infection based on antibody detection was 41%in a high-school located in Northern France [27].Real-time monitoring of the immunity level will be crucial to define efficient management poli-cies, if a new outbreak occurs. We propose such a tool which is based on the modeling approach M of this paper (see Immunity tab in http://covid19-forecast.biosp.org/ ). Remarkably,the estimated levels of immunity are comparable to those observed in Spain by population-basedserosurveys [28], with values ranging from 0 .
5% to 13 .
0% at the beginning of May at the provincialresolution. Such a large-scale serological testing campaign has not yet been carried out in France.However, if such data become available in France, our predictions could be evaluated and our modelupdated accordingly by including this new dataset in our estimation procedure. The mechanistic-statistical approach that we followed here can indeed be easily adapted to multiple observationprocesses.Our results indicate that —at the country scale— travel restriction alone, although they mayhave a significant effect on the cumulative number of deaths and the reproduction number over adefinite period, are less efficient than social distancing and other sanitary measures. Obviously,these results may strongly depend on the parameter values, which have been chosen here on thebasis of values estimated during the lockdown period. This is consistent with results for China [29],where travel restrictions to and from Wuhan have been shown to have a modest effect unless pairedwith other public health interventions and behavioral changes.The model selection criteria led to a strong evidence in favor of the selection of model M withnon-local transmission and spatially-constant contact rate. It is much more parsimonious thanthe fully heterogeneous model M and is therefore better suited to isolating key features of theepidemiological dynamics [4]. Despite important restrictions on movement during the consideredperiod (mandatory home confinement except for essential journeys until 11 May and a 100 kmtravel restriction until 2 June), the model M was also selected against the model M which doesnot take into account non-local transmission. This shows that intercounty transmission is one ofthese key-features that the non-local model manages to take into account.13ore generally, in France just as in Italy [8], the spatial pattern of COVID-19 incidence indicatesthat spatial processes play a key role. At this stage, only a few models can address this aspect.Some have adopted a detailed spatio-temporal modeling approach and use mobility data (see [8]for an SEIR-like model with 9 compartments). The framework we develop here, including thenon-local model and the associated estimation procedure, should be of broad interest in studyingthe spatial dynamics of epidemics, due to its theoretical and numerical simplicity and its abilityto accurately track the epidemics. This approach applies when geographic distance matters, whichmay not be the case at the scale of countries like the US. However, it does at more regional scales.Furthermore, we can envision natural extensions of the approach that would take into account longrange dispersal events in the interaction term. Data availability
Acknowledgments
This work was funded by INRAE (MEDIA network) and EHESS. We thank Jean-Fran¸cois Rey forassistance in developing the web app.
References [1] A Deslandes, et al., SARS-COV-2 was already spreading in France in late December 2019.
International Journal of Antimicrobial Agents , 106006 (2020).[2] J Demongeot, Y Flet-Berliac, H Seligmann, Temperature decreases spread parameters of thenew Covid-19 case dynamics.
Biology , 94 (2020).[3] L Roques, E Klein, J Papaix, A Sar, S Soubeyrand, Using early data to estimate the actualinfection fatality ratio from COVID-19 in France. Biology , 97 (2020).[4] AL Bertozzi, E Franco, G Mohler, MB Short, D Sledge, The challenges of modeling andforecasting the spread of COVID-19. Proceedings of the National Academy of Sciences https://doi.org/10.1073/pnas.2006520117 (2 July 2020).[5] BF Maier, D Brockmann, Effective containment explains subexponential growth in recentconfirmed COVID-19 cases in China.
Science , Issue 6492, 742–746 (2020).[6] K Prem, et al., The effect of control strategies to reduce social mixing on outcomes of theCOVID-19 epidemic in Wuhan, China: a modelling study.
The Lancet Public Health , Issue5, 261–270 (2020).[7] H Salje, et al., Estimating the burden of SARS-CoV-2 in France. Science , Issue 6500,208–211 (2020). 148] M Gatto, et al., Spread and dynamics of the COVID-19 epidemic in Italy: Effects of emergencycontainment measures.
Proceedings of the National Academy of Sciences , 10484–10491(2020).[9] AL Garc´ıa-Basteiro, et al., Monitoring the COVID-19 epidemic in the context of widespreadlocal transmission.
The Lancet Respiratory Medicine , 440–442 (2020).[10] L Roques, EK Klein, J Papaix, A Sar, S Soubeyrand, Impact of lockdownon the epidemic dynamics of COVID-19 in France. Frontiers in Medicine https://doi.org/10.3389/fmed.2020.00274 (5 June 2020).[11] J Zhang, et al., Age profile of susceptibility, mixing, and social distancing shapethe dynamics of the novel coronavirus disease 2019 outbreak in China. medRxiv https://dx.doi.org/10.1101%2F2020.03.19.20039107 (20 March 2020).[12] LM Berliner, Physical-statistical modeling in geophysics.
J Geophys Res , 8776 (2003).[13] S Soubeyrand, AL Laine, I Hanski, A Penttinen, Spatiotemporal structure of host-pathogeninteractions in a metapopulation.
American Naturalist , 308–320 (2009).[14] F Zhou, et al., Clinical course and risk factors for mortality of adult inpatients with COVID-19 in Wuhan, China: a retrospective cohort study.
The Lancet , Issue 10229, 1054–1062(2020).[15] X He, et al., Temporal dynamics in viral shedding and transmissibility of COVID-19.
NatureMedicine , , Issue 5, 672–675 (2020).[16] D Kendall, Discussion of ’Measles periodicity and community size’ by MS Bartlett. J. Roy.Stat. Soc. A , 64–76 (1957).[17] L Bonnasse-Gahot, H Berestycki, M-A Depuiset, M B Gordon, S Roch´e, N Rodriguez, J-PNadal, Epidemiological modelling of the 2005 French riots: a spreading wave and the role ofcontagion.
Scientific Reports , 1–20 (2018).[18] D Brockmann, L Hufnagel, T Geisel, The scaling laws of human travel. Nature , 462–465(2006).[19] S Meyer, L Held, , et al., Power-law models for infectious disease spread.
The Annals of AppliedStatistics , 1612–1639 (2014).[20] V Pareto, Cours d’Economie Politique, (Lausanne: F. Rouge, 1896).[21] J Gaudart, et al., Demography and diffusion in epidemics: malaria and black death spread. Acta Biotheoretica , 277–305 (2010).[22] H Akaike, A new look at the statistical model identification. IEEE transactions on automaticcontrol , 716–723 (1974).[23] G Schwarz, Estimating the dimension of a model. Ann Stat , 461–464 (1978).[24] H Nishiura, G Chowell, The effective reproduction number as a prelude to statistical estimationof time-dependent epidemic trends in Mathematical and statistical estimation approaches inepidemiology . (Springer), pp. 103–121 (2009).1525] E Dong, H Du, L Gardner, An interactive web-based dashboard to track COVID-19 in realtime.
The Lancet Infectious Diseases , Issue 5, 533–534 (2020).[26] S Soubeyrand, et al., The current COVID-19 wave will likely be mitigated in the second-lineEuropean countries. medRxiv https://doi.org/10.1101/2020.04.17.20069179 (22 April 2020).[27] A Fontanet, et al., Cluster of COVID-19 in northern France: A retrospective closed cohortstudy. medRxiv https://doi.org/10.1101/2020.04.18.20071134 (23 April 2020).[28] M Poll´an, et al., Prevalence of SARS-CoV-2 in Spain (ENE-COVID): a nationwide, population-based seroepidemiological study. The Lancet https://doi.org/10.1016/S0140-6736(20)31483-5(6 July 2020).[29] M Chinazzi, et al., The effect of travel restrictions on the spread of the 2019 novel coronavirus(COVID-19) outbreak.
Science , 395–400 (2020).[30] X Liu, S Zhang, Covid-19: Face masks and human-to-human transmission.
Influenza andOther Respiratory Viruses (2020). 16 upplementary MaterialSI 1. Supplementary figures • • In Fig. S2, we describe the timeline of the spatio-temporal dynamics of the immunity rateduring the observation period.Figure S1: (a) Mean temperature in each county over the period ranging from 30 March 2020 to11 June 2020. (b) Immunity rate vs mean temperature: we observe a negative correlation betweenthe mean temperature and the immunity rate (Pearson correlation coefficient: − . M , from 6 April 6 to 8 June.18 I 2. Auxiliary models
At each date t , the time-dependent parameters α ( t ) in M and M and α k ( t ) in M maximumlikelihood estimators (MLEs) in the window ( t − τ / , t + τ /
2) of the following auxiliary models: ˜ S (cid:48) ( s ) = − ˜ αN ˜ S ˜ I, ˜ I (cid:48) ( s ) = ˜ αN ˜ S ˜ I − ( β + γ ) ˜ I, ˜ R (cid:48) ( s ) = β ˜ I, ˜ D (cid:48) ( s ) = γ ˜ I, for s ∈ ( t − τ / , t + τ / , ( ˜ M ,t )and ˜ S (cid:48) k ( s ) = − ˜ αN k ˜ S k ˜ I k , ˜ I (cid:48) k ( s ) = ˜ αN k ˜ S k ˜ I k − ( β + γ ) ˜ I k , ˜ R (cid:48) k ( s ) = β ˜ I k , ˜ D (cid:48) k ( s ) = γ ˜ I k , for s ∈ ( t − τ / , t + τ / , ( ˜ M ,t )and ˜ S (cid:48) k ( s ) = − ˜ α k N k ˜ S k ˜ I k , ˜ I (cid:48) k ( s ) = ˜ α k N k ˜ S k ˜ I k − ( β + γ ) ˜ I k , ˜ R (cid:48) k ( s ) = β ˜ I k , ˜ D (cid:48) k ( s ) = γ ˜ I k , for s ∈ ( t − τ / , t + τ / . ( ˜ M ,t )The initial condition in these models is computed iteratively from the solutions of M , M and M , respectively, over the period [ t i , t − τ /τ /