key: cord-0116585-6k1ky3hn
authors: Mammeri, Youcef
title: A reaction-diffusion system to better comprehend the unlockdown: Application of SEIR-type model with diffusion to the spatial spread of COVID-19 in France
date: 2020-05-06
journal: nan
DOI: nan
sha: 95b6ff2fe59c191deb73b0329b41c3532d2a744c
doc_id: 116585
cord_uid: 6k1ky3hn

A reaction-diffusion model was developed describing the spread of the COVID-19 virus considering the mean daily movement of susceptible, exposed and asymptomatic individuals. The model was calibrated using data on the confirmed infection and death from France as well as their initial spatial distribution. First, the system of partial differential equations is studied, then the basic reproduction number, R0 is derived. Second, numerical simulations, based on a combination of level-set and finite differences, shown the spatial spread of COVID-19 from March 16 to June 16. Finally, scenarios of unlockdown are compared according to variation of distancing, or partially spatial lockdown.

In late 2019, a disease outbreak emerged in the city of Wuhan, China. The culprit was a certain strain called Coronavirus Disease 2019 or COVID-19 in brief (Wor, 2020a) . This virus has been identified to cause fever, cough, shortness of breath, muscle ache, confusion, headache, sore throat, rhinorrhoea, chest pain, and nausea (Hui et al., 2020; Chen et al., 2020) . COVID-19 belongs to the Coronaviridae family. A family of coronaviruses that cause diseases in humans and animals, ranging from the common cold to severe diseases. Although only seven coronaviruses are known to cause disease in humans, three of these, COVID-19 included, can cause severe infection, and sometimes fatal to humans.

COVID-19 spreads fast. According to WHO (Wor, 2020b) , it only took 67 days from the beginning of the outbreak in China last December 2019 for the virus to infect the first 100,000 people worldwide. As of the 5th of May 2020, a cumulative total of 3,601,760 confirmed cases, while 251,910 deaths have been recorded for COVID-19 by World Health Organization (Wor, 2020c) .

Last 30th of January, WHO characterized COVID-19 as Public Health Emergency of International Concern (PHEI) and urge countries to put in place strong measures to detect disease early, isolate and treat cases, trace contacts, and promote social distancing measures commensurate with the risk (Wor, 2020d) . In response, the world implemented its actions to reduce the spread of the virus. Limitations on mobility, social distancing, and self-quarantine have been applied. More than 100 countries established full or partial lockdown. All these efforts have been made to reduce the transmission rate of the virus. For the time being, COVID-19 infection is still on the rise.

Many mathematical models have been proposed to help governments as an early warning device about the size of the outbreak, how quickly it will spread, and how effective control measures may be. Most of the model are (discrete or continuous) SIR-type and few are taken into account the spatial spread. Gardner (2020) implemented a metapopulation network model described by a discrete-time Susceptible-Exposed-Infected-Recovered (SEIR) compartmental model. The model gives an estimate of the expected number of cases in mainland China at the end of January 2020, as well as the global distribution of infected travelers. Wu et al. (2020) fused an SEIR metapopulation model to simulate epidemic. Danon et al. (2020) incorporated daily movements in an SEIR model, while Giuliani et al. (2020) proposed a statistical model to handle the diffusion of covid-19 in Italy.

Here, spatial propagation is translated by a diffusion and the reaction terms are deduced from an extension from the classical SEIR model by adding a compartment of asymptomatic infected (Arcede et al., 2020) . We aim at predicting the spread of COVID-19 by giving maps the basic reproductive numbers R 0 and its effective reproductive number R eff . Afterward, we also assess possible scenarios of unlockdown.

The rest of the paper is organized as follows. Section 2, outlines our methodology. Here the model was explained, where the data was taken, and its parameter estimates. Section 3 contained the qualitative analysis for the model. Here, we provide the reproductive number R 0 , then compare strategies to handle unlockdown. Finally, section 4 outlines our brief discussion on some measures to limit the outbreak.

In this study, we used the publicly available dataset of COVID-19 provided by the Santé Publique France. This dataset includes daily count of confirmed infected cases, recovered cases, hospitalizations and deaths. Data can be downloaded from https://www.data.gouv.fr/fr/datasets/donnees-hospitalieresrelatives-a-lepidemie-de-covid-19/. These data are collected by the National Health Agency and are directly reported public and unidentified patient data, so ethical approval is not required. The map of population density are from Géoportail (https://www.geoportail.gouv.fr) established by the National Geographic Institute. Data concerning transport are extracted from National Institute of Statistics and Economic Studies (https://statistiques-locales.insee.fr/). We focus our study on six components of the epidemic flow (Figure 1 ), i.e. the densities of Susceptible individual (S), Exposed individual (E), symptomatic Infected individual (I s ), asymptomatic Infected individual (I a ), Under treatment individual (U ) and Removed individual (R). To simplify the readings, treatments are not distinguished between quarantine, hospitalization or medicine. To build the mathematical model, we followed the standard strategy developed in the literature concerning SIR model (Diekmann and Heesterbeek, 2000; Brauer and Castillo-Chavez, 2012; Arcede et al., 2020) . We assumed that susceptible can be infected by exposed and by infected individuals. We suppose that only susceptible, exposed and asymptotic individuals are moving. The dynamics is governed by a system of three partial differential equations (PDE) and three ordinary differential equations (ODE) as follows, for

(1)

Frontiers being now closed, homogeneous Neumann boundary conditions is imposed. The total living population is N = S + E + I a + I s + U + R and the death is D = µ(I s + U ) No new recruit is added. The parameters are described in Table A of Figure 3 .

Latency period and infection period have been estimated as 5 days and 7 days respectively (Lauer et al., 2020) , and thus δ = 1/5, γ = 1/7. To account for the lockdown and unlockdown, the average number of contacts is updated as follows )

while the diffusion coefficient is set up to

Here bol denotes for beginning of lockdown and eol for end of lockdown. Unlockdown is assumed to be faster than lockdown. The parameter 0 ≤ η ≤ 1 is a varying coefficient translating respect for distancing. From INSEE, the daily commuting in France is around 25km, the value of d 0 is fixed equal to 25 2 16 (Okubo, 1980; Shigesada and Kawasaki, 1997) . Six parameters θ = (ρ, β e , β s , β a , p, µ) remain to be determined. Given, for N days, the observations I s,obs (t i ) and D obs (t i ), the cost function consists of the nonlinear least square function

denote the output of the mathematical model at time t i computed with the parameters θ. The optimization problem is solved using Approximate Bayesian Computation combined with a quasi-Newton method (Csilléry et al., 2010) . 

From the map of the country, the level-set (Osher and Fedkiw, 2002; Sethian, 1999) is defined as the function φ such that the territory

and its boundary is the zero level of φ The exterior normal is n = ∇φ ||∇φ|| . The computation domain consists in a cartesian grid given by the image pixels (Figure 2) . Then the map of population density allows to build initial population density N 0 (x). From data of confirmed infecteds, the spatial distribution of initial infection I s,0 (x) with respect to department number is created ( Figure  2 ).

Finally, Runge-Kutta 4 is used for the time discretization and central finite difference for the space discretization.

Let us suppose that for all time, d 1 ≤ d(t) ≤ d 0 . We first prove that the model is globally well posed.

Then there exists a unique global in time weak solution (S, E, I a , I s , U, R) ∈ L ∞ (R + , L ∞ (Ω)) 6 , of the initial boundary value problem. Moreover, the solution is nonnegative and S + E + I a + I s + U + R ≤ N 0 .

Proof. Let Y = (S, E, I a , I s , U, R). Thanks to the comparison principle and according to the Duhamel formulation, we look for a time T > 0 such that the map

is a contraction mapping from the closed ball

Here G d , K e , and K s are the kernels of the respective operators ∂ t − d∆, ∂ t − d∆ + δ, and ∂ t − d∆ + γ for d = d 0 or d 1 . According to Ouhabaz (2005) , there exists a constant C Ω > 0 depending only on Ω such that the kernels satisfy

Combining with the fact that the integral terms of the right-hand-side are locally Lipschitz, choosing T 1 CΩ allows to apply Picard's fixed point theorem (Henry, 1981 ).

If f = (f 1 , . . . , f 6 ) denotes the right-hand-side of the system (1), and Y = (S, E, I a , I s , U, R), since f i (Y i = 0) ≥ 0, we deduce that the solution is nonnegative if the initial datum is nonnegative. Finally, maximum principle provides boundedness of solution.

We give a condition on parameters such that the disease has an exponential initial growth.

Theorem 3.2. Let (S 0 , E 0 , I a,0 , I s,0 , 0, 0) be a nonnegative initial datum. If the basic reproduction number

then (E, I a , I s ) exponentially grows.

This number has an epidemiological meaning. The term βe δ represents the transmission rate by exposed during the average latency period 1/δ. The term

(1−p)βa γ is the transmission rate by asymptomatic during the average infection period 1/γ, and the last one is the part of symptomatic. Proof. A linearization around (S 0 , E 0 , I a,0 , I s,0 , 0, 0) is written as the linear system of differential equations

Let (v k ) k≥1 be an orthonormal basis of eigenfunctions of the Laplace operator with the homogeneous Neumann boundary condition, i.e. −∆v k = k 2 v k . Therefore, the characteristic polynomial of the matrix

is P (x) = x 3 + a 2 x 2 + a 1 x + a 0 , with a 0 = (γ + µ)(d 0 k 2 + γ)(d 0 k 2 + δ) (1 − R k ) and

If R k > 1, there is at least one positive eigenvalue that coincides with an initial exponential growth rate of solutions.

To reflect the spatiotemporal dynamic of the disease, we consider the effective reproduction number

and its mean with respect to the domain Ω

The value of R 0 is computed in Table A of Figure 3 . We now establish the asymptotic behavior of solution.

Theorem 3.3. With the same assumptions as Theorem 3.1. Suppose moreover ω 0 β e ≤ δ. Then the solution converges almost everywhere to the Disease Free Equilibrium (S * , 0, 0, 0, R * ) with S * + R * = N * .

Proof. From the last differential equation in system (1), we deduce that R is an increasing function bounded by N (0). Thus R(t) converges to R * a.e. as t goes to +∞. Then integrating over time this equation provides

which is finite. Furthermore, I s , I a , U also go to 0 a.e. as t → +∞ thanks to the positivity of the solution. Multiplying the second equation by E and integrating over Ω give

Since 0 ≤ S N ≤ 1 and ω ≤ ω 0 , Young's inequality followed by Poincaré inequality provide

Since I s and I a go to 0 when t → +∞, it is enough to choose ε > 0 such that dC Ω + δ − ω 0 β e − ε 2 > 0, to conclude that E → 0 a.e.

Remark 3.4. The basic reproduction number R 0 can be computed thanks to the next generation matrix of the model without diffusion as in van den Driessche and Watmough (2000) . Since the infected individuals are in E, I a and I s , new infections (F) and transitions between compartments (V) can be rewritten as

Thus, R 0 = ρ(F V −1 ) of the next generation matrix

No treatment is applied, then ν = 0 Calibration of the model is done from January 24, 2020, the day of first confirmed infection, to April 30, i.e. 97 days. Since R 0 = 3.425257 and 1 − 1 R0 = 0.708051, 70% of the population is set as susceptible to the infection due to the herd immunity. The objective function J is computed to provide a relative error of order less than 10 −2 . In Figure 3 , Table A shows estimated parameters. Remark that ω 0 β e ≤ δ. The rest of the Figure compares the data and the fitted total symptomatic infected and death of the posterior distribution. 

Simulations are performed from January 24, 2020, the day of first confirmed infection, to June 16. Images are 957 × 984 pixels, time step is chosen to verify the stability condition preserving positivity, and η = 0.6. Figure 4 -A presents the day before lockdown, the disease is mainly located in regions North (Hauts-de-France,Île-de-France), , East (Alsace), Center-East (Rhône) and South-East (Bouches-du-Rhône). When lockdown ends on May 10, Figure 4 -B1 shows that the disease remains located into these regions, except for the South where it almost vanishes. Provided that distancing is predominantly respected (here η = 0.6), the density of symptomatic infecteds is decreasing until June 16 everywhere (Figure 4-B2) .

Figures 4-C represents what would happen if nothing has been done. On May 10, the disease is located in the same regions, but with a greater number of infected. As a consequence, the spread continues from East to West. While the disease vanishes in the Eastern part of the country on June 16, the West is peaking. It is important to see that, without intervention, the total number of infected is high and the whole country is deeply affected. 

A naive strategy is to unlock the regions where R eff (x, .) is the lowest. Looking at the maps in Figures 5, we see that the effective reproduction number before the lockdown is between 3.34 and 3.42 (Figures 5-A) . On May 10, it goes down and is between 0.38 and 0.42 (Figures 5-B) . After May 10, with a distancing mostly respected (η = 0.6), R eff (x, .) increases again. It exceeds 1 with lower value (1.04) in the regions that have been most affected upstream (Figures 5-C) . A less forced strategy is to believe in the respect of distancing after lockdown. With this set of parameters, R eff goes to 1 when η = 0.63. As shown in Figure  6 -A, provided that at least 63% of the individuals respect the distancing rules, the number of symptomatic infected individuals is controlled. Figure 6 -B shows that in this case, the effective reproduction number R eff remains less than 1. Below 63% regarding distancing, the number of symptomatic infected increases again and R eff becomes greater than 1.

A third strategy, more constraining for some, is to continue the lockdown in the most affected areas. Here, lockdown is continued on the eastern half, while the western half is free. Then the number of symptomatic infected and the effective reproduction number pursue their decay as shown in Figure 6 (dotted lines). It is important to note that in the simulation no individuals move from a confined to an unconfined region. A homogeneous Neumann condition, equal to 0, is imposed along this fictitious frontier. 

Lockdown has reduced both the number of infections and the spread. The propagation focused in the Eastern half and kept the West intact. Without lockdown, the whole country would be affected and more severely. Lockdown caused a significant reduction of the reproduction number, from 3.42 to 0.38. Since the number of susceptible individuals remain large, as soon as contacts increase, the effective reproduction number R eff grows and can exceed 1.

With a lack of treatment, social distancing remains the most effective means. We notice that it has to be highly respected (here at over 63%). As shown in Figure 6 -A, the number of symptomatic infected individuals, therefore potentially hospitalized, is restrained as soon as at least 63% of the individuals respect the distancing, Nevertheless, this constraint can be relaxed since it can be imposed only in the most infectious areas.

In summary, to obtain a possible unlockdown map, the local value of the effective reproduction number should be taken into account, as well as the number of infected individuals and the direction of the spread of the disease.

Of course, this is a simplified model based only the population density and mean daily commuting. For example, the model could be improved by considering a larger diffusion along major axes of travel (i.e. d = d(x, t)), by taking local effects of distancing (e.g η = η(x, t) where η(x, t) = 0 in closed schools), or by opening of the frontiers (by changing the boundary conditions and adding new recruits). Supplementary video Spatiotemporal propagation of COVID-19 from March 16 to June 16 with lockdown occurring from March 17 to May 11 (left) and without intervention (right).

Accounting for symptomatic and asymptomatic in a seir-type model of covid-19

Epidemiological and clinical characteristics of 99 cases of 2019 novel coronavirus pneumonia in wuhan, china: A descriptive study

Approximate bayesian computation (abc) in practice

A spatial model of covid-19 transmission in england and wales: early spread and peak timing. medRxiv

Mathematical Epidemiology of Infectious Diseases: Model Building, Analysis and Interpretation

Modeling the spreading risk of 2019-ncov

Modelling and predicting the spatiotemporal spread of coronavirus disease 2019 (covid-19) in italy. arXiv

Geometric theory of semilinear parabolic equations

The continuing 2019-ncov epidemic threat of novel coronaviruses to global health -the latest 2019 novel coronavirus outbreak in wuhan, china

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

Predicting the cumulative number of cases for the covid-19 epidemic in china from early data

Diffusion and Ecological Problems: Mathematical Models

Level Set Methods and Dynamic Implicit Surfaces

Analysis of Heat Equations on Domains

Level Set Methods and Fast Marching Methods: Evolving Interfaces in Computational Geometry, Fluid Mechanics, Computer Vision, and Materials Science

Biological Invasions: Theory and Practice

Reproduction numbers and sub-threshold endemic equilibria for compartmental models of disease transmission

Naming the coronavirus disease (COVID-19) and the virus that causes it. World Health Organization

World Health Organization

Emergency Committee Regarding the Outbreak of Novel Coronavirus (2019-nCoV)

World Health Organization

detail/30-01-2020-statementon-the-secondmeeting-of-the-international-health-regulations-(2005)-emergencycommittee-regarding-the-outbreak-of-novel-coronavirus

Nowcasting and forecasting the potential domestic and international spread of the 2019-ncov outbreak originating in wuhan, china: A modelling study