Epidemic model on a network: analysis and applications to COVID-19

TL;DR

This paper analyzes a network-based epidemic model with applications to COVID-19, introducing strategies for vaccination and containment, and providing insights into effective interventions and deconfinement scenarios.

Contribution

It presents a novel epidemic model on networks with a new vaccination strategy based on node degree and applies it to COVID-19 data.

Findings

Vaccinating high-degree nodes is most effective.

The model can evaluate deconfinement and second wave scenarios.

Few parameters allow fitting to COVID-19 data.

Abstract

We analyze an epidemic model on a network consisting of susceptible-infected-recovered equations at the nodes coupled by diffusion using a graph Laplacian. We introduce an epidemic criterion and examine different vaccination/containment strategies: we prove that it is most effective to vaccinate a node of highest degree. The model is also useful to evaluate deconfinement scenarios and prevent a so-called second wave. The model has few parameters enabling fitting to the data and the essential ingredient of importation of infected; these features are particularly important for the current COVID-19 epidemic.

Tables

Table 1. Table 1: Isolated node j 𝑗 j and associated eigenvalues of M 𝑀 M

$j$	degree	${\sigma }_3$	${\sigma }_2$	${\sigma }_1$
6	1	$-2.46$	$-1.30$	$1.22{10}^{-1}$
2	2	$-2.14$	$-7.85{10}^{-1}$	$3.83{10}^{-2}$
4	2	$-2.51$	$-4.79{10}^{-1}$	$-1.06{10}^{-2}$
3	3	$-1.44$	$-8.29{10}^{-1}$	$-3.71{10}^{-2}$
1	3	$-1.52$	$-6.41{10}^{-1}$	$-6.72{10}^{-2}$
5	3	$-1.52$	$-6.41{10}^{-1}$	$-6.72{10}^{-2}$

Table 2. Table 2: Isolated nodes i , j 𝑖 𝑗 i,j and associated eigenvalues of M 𝑀 M .

$i$	$j$	$d_i$	$d_j$	neighbors?	${\sigma }_3$	${\sigma }_2$	${\sigma }_1$
5	6	3	1	yes	$-2.97$	$-1.27$	$1.67{10}^{-1}$
2	3	2	3	yes	$-1.99$	$-9.49{10}^{-1}$	$1.56{10}^{-1}$
1	2	3	2	yes	$-2.57$	$-6.83{10}^{-1}$	$1.27{10}^{-1}$
4	6	2	1	no	$-3.00$	$-1.83$	$9.17{10}^{-2}$
1	3	3	3	yes	$-1.31$	$-9.66{10}^{-1}$	$7.53{10}^{-2}$
3	4	3	2	yes	$-2.26$	$-6.24{10}^{-1}$	$6.01{10}^{-2}$
2	4	2	2	no	$-2.68$	$-9.71{10}^{-1}$	$4.77{10}^{-2}$
4	5	2	3	yes	$-2.82$	$-3.94{10}^{-1}$	$-9.08{10}^{-3}$
1	4	3	2	no	$-2.61$	$-6.16{10}^{-1}$	$-3.48{10}^{-2}$
2	6	2	1	no	$-2.33$	$-1.89$	$-3.84{10}^{-2}$
1	6	3	1	no	$-2.60$	$-1.14$	$-1.09{10}^{-1}$
3	6	3	1	no	$-2.41$	$-1.06$	$-2.12{10}^{-1}$
1	5	3	3	yes	$-1.39$	$-4.24{10}^{-1}$	$-2.50{10}^{-1}$
2	5	2	3	no	$-1.89$	$-5.41{10}^{-1}$	$-2.91{10}^{-1}$
3	5	3	3	no	$-1.33$	$-5.97{10}^{-1}$	$-3.37{10}^{-1}$

Table 3. Table 3: Initial conditions and local R 0 subscript 𝑅 0 R_{0} for the plots of Fig. 4 .

node $j$	1	2	3	4	5	6
$S_j\left(t=0\right)$	0.26	0.14	0.55	0.16	0.5	0.18
$I_j\left(t=0\right)$	0	0	0.01	0	0	0
$R_0=\beta \frac{S_j^{}{}_0{}^{}}{\gamma }$	1.0296	0.5544	2.1780	0.6336	1.9800	0.7128
$I^{\ast }$	0.0001	0.0364	0.1109	0.0227	0.0750	0.0130

Table 4. Table 4: The twelve regions of mainland France, the main cities and the number of patients in hospital on April 5 2020.

region	main cities	number of patients
Ile de France	Paris	11474
Grand est	Strasbourg, Mulhouse	4739
Auvergne-Rhone-Alpes	Lyon, Grenoble, Clermont-Ferrand	2983
Provence-Alpes-Cote d’Azur	Marseille, Nice	1722
Hauts de France	Lille, Valenciennes	2119
Bourgogne-Franche Comté	Dijon, Besancon	1120
Occitanie	Toulouse, Montpellier	1024
Nouvelle Acquitaine	Bordeaux	762
Centre Val de Loire	Orleans, Tours	751
Normandie	Rouen, Caen, Le Havre	710
Pays de la Loire	Nantes	667
Bretagne	Rennes, Brest	423

Equations

$$\left\{\begin{array}[]{ll}{\dot{S}}=-\beta SI,\\ {\dot{I}}=\beta SI-\gamma I\\ {\dot{R}}=\gamma I\\ \end{array}\right.$$

$$R_0\equiv \frac{\beta }{\gamma },$$

$$S^{\ast }=\frac{1}{R_0},$$

$$I^{\ast }=I_0+S_0-\frac{1}{R_0}(1+\log (R_0{S}_0))$$

$${\dot{S}}_j=-\beta \frac{S_j}{N_j}{I}_j+ϵ\sum _{k\sim j}(S_k-S_j),$$

$${\dot{I}}_j=\beta \frac{S_j}{N_j}{I}_j-\gamma I_j+ϵ\sum _{k\sim j}(I_k-I_j),$$

$$\left\{\begin{array}[]{ll}{\dot{S}}=\epsilon\Delta S-\beta S~{}I,\\ {\dot{I}}=\epsilon\Delta I+\beta SI-\gamma I.\\ \end{array}\right.$$

$${\Delta }_{kl}=1\mathrm{if}kl\mathrm{connected},0\mathrm{otherwise};{\Delta }_{kk}=-\sum _{l\ne k}{w}_{kl}.$$

$$0=-{\omega }_1^2\ge -{\omega }_2^2\ge \cdots \ge -{\omega }_n^2.$$

$$\Delta v^j=-{\omega }_j^2{v}^j.$$

$$\Omega =\{(S,I):S\ge 0;I\ge 0\}.$$

$$\dot{S}=ϵ\Delta S$$

$$\dot{I}=ϵ\Delta I-\gamma I$$

$$\dot{I}=(ϵ\Delta -\gamma )I+\beta SI.$$

$$\dot{I}=MI$$

$$M=ϵ\Delta -\gamma {Id}_n+\mathrm{diag}(\beta S_1,\beta S_2,\dots ,\beta S_n).$$

$$I=\sum _{k=1}^n{\gamma }_k{v}^k,$$

$${\dot{\gamma }}_k=(-{\omega }_k^2-\gamma ){\gamma }_k+<\beta SI|v^k>.$$

$${\dot{\gamma }}_k=(-{\omega }_k^2-\gamma +\beta S){\gamma }_k.$$

$${\dot{\gamma }}_k=(-{\omega }_k^2-\gamma ){\gamma }_k+\beta \sum _{l=1}^n{\gamma }_l\left(\sum _{j=1}^n{S}_j{v}_j^l{v}_j^k\right).$$

$$ϵ=1,\beta =6.5,\gamma =\mathrm{6.25.}$$

$$ϵ=1,\beta =6.75,\gamma =6.25,$$

$${\sigma }_1=\underset{\parallel X\parallel =1}{\sup }{X}^{T}MX.$$

$$X=(1,0\dots 0{)}^T,X=(0,1,0\dots 0{)}^T,\dots$$

$${\sigma }_1$$

$$\dots$$

$${\sigma }_1$$

$${\sigma }_1\ge -{\min }_k{d}_k+\beta S-\gamma .$$

$$M=\left(\begin{array}{ccccc}-1& 1& 0& \dots & 0\\ 1& -1& 0& \dots & 0\\ 0& \dots & \dots & 0& \end{array}\right)$$

$${\sigma }_1^{\prime }\le {\sigma }_1.$$

$$\beta =1.98,\gamma =0.5,ϵ=10^{-3}$$

$$\beta =0.33,\gamma =0.13,$$

Full text

1 Introduction

Many models of the propagation of an epidemic such as the current COVID-19 [1] involve a network. This can be a contact network between individuals. Then, the network is oriented and is used to understand how a given individual can infect others at the very early stages. The models are typically probabilistic, see [2] for example. Once the epidemic is established, the geographical network becomes important. There, nodes represent locations and edges the means of communication; for COVID-19 these are the airline routes [3]. Such a network is non oriented and the important nodes are the ones that are most connected.

One of simplest models of a disease is the Kermack-McKendrick system of equations [4] involving three populations of susceptible, infected and recovered individuals $(S,I,R)$. Using this model together with a probability transition matrix [5] for the geographic coupling, Brockman and Helbling [6] performed a remarkable study of the propagation of well-known epidemics like SARS or H1N1 due to airline travel. They emphasized that the fluxes between the nodes govern the propagation of the epidemic.

In this article, we consider $(S,I,R)$ Kermack-McKendrick equations coupled to a network through a graph Laplacian matrix [7]. The combination of the simple SIR dynamics with the diffusion yields the essential ingredients to model and understand an epidemic, such as the COVID-19. In particular,

• •

there are few parameters so that fitting to data can be successful,

• •

it contains the essential ingredient of importation of infected subjects from country to country.

The epidemic front is controlled by the availability of susceptibles. If susceptibles are large enough, the front cannot be stopped. The number of susceptibles varies from node to node. Reducing this number at a given location can be done through isolation. This is expensive and cannot be done for the whole network. It is therefore important to address the question: what nodes are more useful to isolate to mitigate the epidemic?

Using this model together with the detailed data available [10] [11], we predicted the onset of the COVID-19 epidemic in Mexico [3]. The present article is devoted to the detailed analysis of the model. We first prove that it is well-posed and that solutions remain positive. We introduce an epidemic criterion that generalizes the well-known $R_{0}$ of the scalar case. For small diffusion, nodes are almost decoupled and an outbreak occurs at a node if the local $R_{0}$ is larger than one. When the diffusion is moderate, the epidemic criterion depends on the network and when there is an outbreak, it starts synchronously on the network. Using this criterion, we define an isolation policy. We find that it is most useful to isolate the high connectivity nodes and not efficient to isolate neighbors. For the particular case of the COVID-19 we discuss the effect of deconfining; the model shows that allowing circulation between heavily and weakly infected areas will prolong the outbreak in the latter.

The article is organized as follows. In section 2, we introduce the model, discuss its main features and present the epidemic criterion. Section 3 shows a simple six node network based on the country of Mexico; there the effect of isolation is discussed. The COVID-19 disease is studied in section 4 and we show the estimation of the time of outbreak in Mexico. The important issue of deconfinement is studied in section 5. We conclude in section 6.

2 The model and epidemic criterion

One of the main models to describe the time evolution of the outbreak of an epidemic is the Kermack-McKendrick model [4]

[TABLE]

where the dynamics of transmission depends of the frequency and intensity of the interactions between (healthy) susceptible $S$ and infected individuals $I$ and produce recovered individuals $R$. The parameters $\beta$ and $\gamma$ are the infection rate and the recovery rate. The model conserves $N=S+I+R$ the total number of individuals. Note that $R$ is essentially the integral of $I$ and therefore plays no role in the dynamics. We will omit it below and only discuss $S$ and $I$.

An epidemic occurs if $\beta S-\gamma>0$ [4]. At $t=0$, $S=1$ so that an infection occurs if the infection factor defined as

[TABLE]

is greater than one. An important moment in the time evolution of $S$ and $I$ is when the number of infected is maximum. The corresponding values $(S^{*},I^{*})$ can be calculated easily; we give the derivation in the Appendix. The expressions are

[TABLE]

Note that $I^{*}$ and $S^{*}$ depend strongly on $R_{0}$. Take for example $\gamma=0.625$ and different values of $\beta$.

[TABLE]

The value of $I^{*}$ depends also on the initial number of susceptibles $S_{0}$ which is smaller than the total number $N$. Generally, a large $N$ gives a large $S_{0}$ and $I^{*}$.

2.1 SIR on a network

We consider a geographic network of cities connected by roads or airline routes. This introduces a spatial component so that $(S,I)$ become vectors; we also drop $R$. This is similar to Murray’s model where he introduces spatial dispersion in an SI model using a continuous Laplacian term [9]. The evolution at a node $j$ in a network of $n$ nodes reads

[TABLE]

where $N_{j}$ is the population at node $j$, the $\sum_{k\sim j}$ is the exchange with the neighboring nodes $k$ of $j$ and where $\epsilon$ is a constant. The main difference with the model of [6] is that we assume symmetry in the exchanges.

The equations (5) can be written concisely as

[TABLE]

where $S=(S_{1},S_{2},\dots,S_{n})^{T},~{}~{}I=(I_{1},I_{2},\dots,I_{n})^{T},~{}~{}\beta\equiv(\beta/N_{1},\beta/N_{2},\dots,\beta/N_{n})^{T}$, $\Delta$ is the graph Laplacian matrix [7] and we denote by $SI$ the vector

$(S_{1}I_{1},S_{2}I_{2},\dots,S_{n}I_{n})^{T}$. The infection rate $\beta$ can vary from one geographical site to another while the recovery rate $\gamma$ depends only on the disease. The diffusion $\epsilon$ should be small so that the populations involved in that process remain much smaller than the node populations $N_{j}$. Another point is that the diffusion could act only on the infected population. We chose to put the diffusion on both $S$ and $I$ for symmetry reasons.

The graph Laplacian ${\Delta}$ is the real symmetric negative semi-definite matrix, defined as

[TABLE]

The graph Laplacian has important properties, see ref. [7], in particular it is a finite difference approximation of the continuous Laplacian [8]. The eigenvalues of $\Delta$ are the $n$ non positive real numbers ordered and denoted as follows:

[TABLE]

The eigenvectors $\{v^{1},\dots,v^{n}\}$ satisfy

[TABLE]

and can be chosen to be orthonormal with respect to the scalar product in $\mathbb{R}^{n}$, i.e. ${v}^{i}\cdot{v}^{j}=\delta_{i,j}$ where $\delta_{i,j}$ is the Kronecker symbol.

2.2 Well posedness and positivity

The model (7) is well posed in the sense that the solution remains bounded. We show this in the Appendix using standard techniques.

The biological domain of the system is

[TABLE]

Let us show that $\Omega$ is an invariant set for (7) so that the model makes sense in biology. Consider the different axes $S_{j}=0$ and $I_{j}=0,~{}j=1,\dots n$. First assume $I_{j}=0,~{}j=1,\dots n$, then equation (7) reduces to

[TABLE]

which conserves the positivity of $S$. Similarly when $S=0$, we get

[TABLE]

and again the positivity of $I$ is preserved.

2.3 Epidemic criterion

Here we extend the 1D epidemic criterion of Kermack-McKendrick [4] to our graph model. Initially, the vector $I$ will follow the second equation of (7)

[TABLE]

Equation (11) describes the onset of the epidemic on the network. It can be written

[TABLE]

where $M$ is the symmetric matrix

[TABLE]

The eigenvalues of $M$ $\sigma_{1},\dots,\sigma_{n}$ are real. If one of them is positive, then the solution $I(t)$ increases exponentially and the epidemic occurs. We can then write

**Epidemic criterion : ** there is an onset of the epidemic if one eigenvalue $\sigma_{i}$ of $M$ is positive.

Two situations occur, depending whether the diffusion is small or moderate. For small diffusion, the contribution of $\Delta$ to $M$ can be neglected. Then each node will develop independently from the others. We will have outbreaks in some and not in others.

When the diffusion is moderate, the Laplacian contributes to $M$. Since $M$ is symmetric the eigenvalues of $M$ remain in the same order as the ones of $\Delta$. This is the interlacing property [7]. Then $\sigma_{1}$ will tend to 0 for $\gamma,\beta\to 0$. Note also that since $\parallel S\parallel$ decreases with time, the estimate given by the eigenvalues of $M$ indicates the size of the epidemic i.e. ${\rm max}~{}\parallel I\parallel$. Then, the eigenvector of $M$ for the eigenvalue $\sigma_{1}$ will be almost constant and the epidemic will start synchronously on the network.

The analysis of the moderate diffusion case can be extended when $\beta$ is constant. Expanding $I$ on an orthonormal basis of eigenvectors ($v^{k}$) of $\Delta$

[TABLE]

we get

[TABLE]

Assume that the susceptible population is constant on the network. Then ${\rm diag}(S_{1},S_{2},\dots,S_{n})=S{\rm Id}_{n}$ so that equation (16) reduces to

[TABLE]

The epidemic starts if $-\gamma+\beta S>0$ which is a simple generalization of the criterion in the scalar case.

When the population of susceptibles is inhomogeneous and $\beta$ is homogeneous, equation (14) becomes

[TABLE]

Then the eigenvectors and the geometry of the network play a role.

3 A simple example

We illustrate the results given above on a 6 node network inspired by the geographical map of Mexico, with six main cities surrounding Mexico city, see Fig. 1. A node represents a city and an edge is a road link between two cities. For simplicity, here we assume that $N_{j}$ is independant of $j$ so that $I$ and $S$ are given in percentages. This will be the case throughout the article unless specified.

The graph Laplacian is

\Delta=\left[\begin{array}[]{cccccc}-3&1&1&0&1&0\\ 1&-2&1&0&0&0\\ 1&1&-3&1&0&0\\ 0&0&1&-2&1&0\\ 1&0&0&1&-3&1\\ 0&0&0&0&1&-1\end{array}\right].

The eigenvalues of this graph laplacian are

[TABLE]

The corresponding eigenvectors are

[TABLE]

3.1 Influence of the diffusion

The variable $\epsilon$ measures the intensity of the diffusion of $S$ and $I$ on the network. When $\epsilon<<1$ the diffusion is very weak and the evolution at each node can be decoupled from the one of its neighbors. For larger $\epsilon$, the diffusion and reaction occur on similar time periods and need to be analyzed together. To see the influence of the diffusion, we plot in Fig. 2 the evolution of $I_{k}(t),~{}k=1,\dots,6$ for $\epsilon=0.1$ (left panel) and $\epsilon=10^{-7}$ (right panel).

Note the times of arrival of the infection, first in node 5 the neighbor of node 6, then nodes 1 and 4 and finally nodes 3 and 2. For the large diffusion (left panel of Fig. 2) the peaks are very close and there is a strong influence between the nodes. On the other hand for a small diffusion, the peaks are well separated and the nodes are decoupled. The maximum of $I_{k}$ is given by the estimate (4).

Infecting node 2 changes the time of arrival of the outbreak as shown in Fig. 3. It reaches first nodes 1 and 3, then nodes 4 and 5 and finally node 6.

3.2 Isolation policies for large diffusion

We first consider that the diffusion and the nonlinear terms have similar orders of magnitude. We will address the case of weak diffusion in the next section.

When the diffusion is large, one should consider the epidemic on the network as a whole and use the topology of the network to reduce the strength of the outbreak. From the amplitude equations (16), one can devise a strategy of isolation. By this we mean reducing $\beta$ so that the maximal eigenvalue of $M$ is minimum.

We choose

[TABLE]

Table 1 shows the eigenvalues $\sigma_{1},\dots,\sigma_{n}$ of $M$ from (12) when isolating a node of the network, i.e. setting $S_{j}=0$ at a specific node $j$ and keeping the other nodes the same. We chose $S=(1,1,1,1,1,1)^{T}$.

In the absence of isolation, the largest eigenvalue of the matrix $M$ is $2.5~{}10^{-1}$. The table shows that it is most effective to isolate nodes 1,3 and 5. These nodes have the highest degree of the network.

We now isolate two cities in the network. The results are presented in table 3. We chose

[TABLE]

and $S=(1,1,1,1,1,1)^{T}$.

Again the high degree nodes 1,3 and 5 are the ones that reduce $\sigma_{1}$ the most and are therefore the most effective when applying isolation. It is also not effective to isolate neighboring nodes.

The results shown in tables 1 and 3 can be explained in part from the properties of the matrix $M$ and the graph Laplacian $\Delta$. The maximal eigenvalue $\sigma_{1}$ of $M$ verifies [7]

[TABLE]

We can find inequalities for $\sigma_{1}$ by choosing

[TABLE]

Denoting $d_{i}$ the degree of node $i$, we get

[TABLE]

so that

[TABLE]

This relation shows that isolating a node that has not smallest degree does not change the estimate. Conversely, if there is a unique node of minimal degree and we isolate it, then the bound changes.

Using similar arguments, it can be shown that isolating two neighboring nodes, say 1 and 2 will be less effective than isolating two non neighboring nodes.

Now we look at what happens if we cut a link, which corresponds to condemning a road for example. Let $\Delta^{\prime}$ be the Laplacian of the new graph obtained by deleting a link. Without loss of generality we can assume this link to be between vertices 1 and 2. Then $\Delta^{\prime}=\Delta-M$ where

[TABLE]

$M$ has all eigenvalues equal to [math] except one which has value $-2$. Applying the Courant-Weyl inequalities, see for example [7], we get the following result for the maximum eigenvalue of $\Delta^{\prime}$

[TABLE]

Note that equality is possible: when $S$ is homogeneous, the maximum eigenvalue of $M$ will always be $-\gamma+\beta S$. In such a case, cutting a link is ineffective.

3.3 Small diffusion and isolation policies

When diffusion is small, the nodes evolve almost independently so that the simple dynamics of the scalar SIR model apply. Then some nodes can be isolated and the epidemic is not seen there. Fig. 4 shows such a situation. We choose

[TABLE]

and the initial conditions are given in the table 3 below. The local $R_{0}$ is also computed and one sees that an outbreak will occur at nodes 1,3 and 5 and not at nodes 2,4 and 6. Fig. 4 shows the peaks for $I_{3}$ and $I_{5}$ and the maxima of $I_{3}$ and $I_{5}$ are close to the ones predicted by the SIR formulas (4).

4 Propagation of COVID-19

We now consider the propagation of the Corona virus COVID-19 on a network consisting of a complete graph of 7 nodes with an additional link to an 8th node, see top left panel of Fig. 5. The 7 nodes correspond to the following cities or regions, Hubei, Beijing, Shanghai, Japan, western Europe, eastern USA and western USA and the 8th node is Mexico. The links are the main airline routes. We assume a complete graph because airline routes connect any two of these regions.

The parameters chosen are $\beta=0.5,\gamma=0.2$ and $\epsilon=10^{-6}{\rm days}^{-2}$. These were suggested by very early estimations of the outbreak in Wuhan. Using the data from the John Hopkins website [10] and our model, we estimated the starting time of the outbreak in Mexico to be from March 20 to March 30 2020 [3]. These results were shared with the Ministry of health of Mexico at the end of February 2020 so that preparations could be made. In the present article, we analyze the model; we therefore comment briefly the evolutions obtained in [3] as an interesting case.

The top right panel of Fig. 5 shows the time evolution in days of the infected in the different nodes. The simulation is started at node 1 (Hubei in red in Fig. 5) with $I_{1}=0.1$ and the other nodes are set at 0. The susceptibles are set to 1 everywhere. As expected the maximum $I_{1}^{*}=0.32$ and the subsequent maxima $I_{j}^{*}=0.21,~{}~{}j=2,\dots 8$. Notice how the nodes 2-7 start simultaneously while node 8 is delayed. Communications with the province of Hubei were restricted at the end of December 2019. To model this, we now consider that node 1 is only connected to node 2 which forms a complete graph with nodes 3-7, see bottom left panel of Fig. 5. The infected are shown in the bottom right panel; as expected the epidemic first arrives in node 2 then synchronously in nodes 3-7 and then in node 8.

To test the effect of having different $\beta$ s at each node we increased $\beta_{3}$, reduced $\beta_{4}$ and kept the other $\beta$ s the same. Then the epidemic arrives sooner at node 3 and later at node 4. As expected from the formulas (4) $I_{3}^{*}>I_{j}^{*}>I_{4}^{*}$ for $j=5,6,7$, see [3].

5 Confinement and deconfinement

Since there is no vaccine for COVID-19 disease and the mortality is relatively high, many countries put in place a confinement or measures to reduce the movement of the population.

China confined the Hubei region around January 22, Italy confined its population on March 9, France on March 17 and so on. In the middle of the epidemic, Spain and France reached a situation like the one shown in Fig. 6. This picture shows the number of patients in hospitals on April 5 2020 for the 12 different regions of mainland France (see Table 4), the data was obtained from the website [12]. Note how some regions are highly infected while others have many fewer cases. To use this regional data of France, one could establish a graph of the main roads and railways connecting the main cities, very much like the one for Mexico in Fig. 1.

Let us now consider the confinement. It can be implemented by

• (i)

reducing the contact ratio $\beta_{j}$ of each node $j$

• (ii)

reducing the diffusion $\epsilon$, ie the travel between nodes

When deconfining the population once the peak of the epidemic has passed the two options (i) and (ii) need to be relaxed, so as to avoid a so-called second wave. This happens in particular when the epidemic affected a small fraction of the total number $N$. Then, relaxing $\beta$ or equivalently increasing $N$ causes a number of new susceptibles to enter the reaction and therefore produce a second peak of infection. The model (5) allows to analyze the effect of the two options. We do this separately but note that in reality both act together and effects can cancel.

Consider option (i). For this, we study the situation at a single node. We choose the parameters

[TABLE]

and the computation is started at $t=0$ with $S_{0}=1,~{}I_{0}=0.01$. Here, as before, the units of $I$ are percentages. In Fig. 7, we show the evolution of the infected for a sudden increase of $\beta$ from 0.33 to 0.5 at $t=20$ before the peak (left panel, line b ) and $t=30$ after the peak (right panel, line b). Clearly, the deconfinement before the peak causes many more infections and could saturate the hospitals. Deconfining after the peak as shown on the right panel of Fig. 7 is harmless. Only one node is involved.

Another way of deconfining is to relax option (ii), that is allowing travel from one node to another. This corresponds to increasing the diffusion which can bring infected from large centers to small centers. To illustrate this, consider the graph of two nodes shown in Fig. 8. It corresponds for the example of France Fig. 6, to allowing travel between the large urban area of Paris and the much less populated Normandy. The capacity of node 1 (Paris) is $N_{1}=20~{}10^{6}$ while $N_{2}=10^{6}$. Now, $I_{1}$ and $I_{2}$ are actual numbers and not percentages. We consider equations (5) where $S_{j},I_{j}$ have been normalized by $10^{6}$. The parameters are

[TABLE]

The initial conditions are

[TABLE]

in millions. To model option (ii) we suddenly increase the diffusion parameter $\epsilon$ to $\epsilon=10^{-2}$.

Fig. 9 shows $I_{1}(t)$ (top panels) and $I_{2}(t)$ (bottom panels) when deconfining before the peak (left panels) and after the peak (right panels). The curves $I_{1}(t),I_{2}(t)$ before and after deconfining are indicated by $a$ and $b$ respectively. The $y$ scale are in millions. Deconfining before the peak is not so harmful for node 1 while it is catastrophic for node 2, $I_{2}$ is multiplied by 3. After the peak, a sudden deconfinement prolongs the epidemic at node 2.

6 Discusion and conclusion

For the COVID-19 epidemic, few articles have addressed the coupling of an epidemic model to the geographical landscape. The subject is difficult as mobility is not well understood at microscopic level. It is also difficult to extract a model from the available data.

In the reference [13], the authors estimated how mobility and transmissibility affect the onset of the epidemic in the cities adjacent to Wuhan in the territory of China. They found that reducing the mobility between cities by some factor did not change the time of onset of the epidemic peaks. This is unexpected as seen in the present article where reducing the mobility by orders of magnitude delays the onset of the epidemic. Also the world data [10] and [11] shows clearly the epidemic peaking first in China, then in Iran, Italy and so on. Therefore, the result of [13] might be due to the fact that small changes are effected on the mobility. It might also be due to the model of mobility chosen.

For Brazil, the study [14] considers a six equation compartimental model coupled to a complex mobility scheme. The authors show that the epidemic curves vary ”enormously” over different geographic scales. Outbreaks can start in big cities and propagate to the countryside or there might be multiple foci of infection. It is difficult to get a view of the epidemic at the scale of the network with such a complex model.

In our study, we used the simplest susceptible-infected equations at the nodes coupled by a geographic diffusion term. This contains the essential ingredient of importation of infected subjects from country to country. We have kept the number of parameters to a minimum so that fitting the data can be successful. This is particularly important for the present epidemic of COVID-19 where one wants to get a global picture at the level of the network.

The Laplacian models well the symmetric flow from one node to another. This is a first order approximation of the spreading of infected in terms of diffusion. One can refine the approximation using a kernel of an anomalous or stochastic diffusion but then more careful measures of the ways the infection travels are needed.

On the analysis side, using this model, we generalized the well-known epidemic criterion of Kermack-McKendrick. For small diffusion, outbreaks occur at different times as the disease advances through the network. A larger diffusion will cause the outbreak to occur synchronously on the network. Using this criterion, we designed an isolation policy: we find it best to isolate high degree nodes and not efficient to isolate neighbors.

We also discussed the important aspect of deconfining a region after the outbreak. Circulation between highly infected regions and less impacted areas should be reduced to prevent the spread of infected to the latter. Large clusters should be carefully controlled.

Finally the study points out the usefulness of having accurate data at country and local levels (cities, neighborhoods and hospitals).

Acknowledgements

This work is part of the XTerM project, co-financed by the European Union with the European regional development fund (ERDF) and by the Normandie Regional Council.

Appendix A Analysis of the SIR model

The expressions (4) may easily obtained in terms of $R_{0}$ as it is shown below.

From $\dot{I}=0$ we obtain

[TABLE]

Dividing the second equation of (1) by the first, we get

[TABLE]

which can be integrated to yield

[TABLE]

where we assumed $S(t=0)=S_{0}$ and $I(t=0)=I_{0}$. Then one can compute $I^{*}$

[TABLE]

Assuming $S_{0}=1$, equations (23,25) can written in terms of $R_{0}$ as

[TABLE]

The time $t^{*}$ corresponding to $S^{*},I^{*}$ can be calculated in the following way.

From the second equation of (1) one can write

[TABLE]

where we have substituted $I(S)$ from (24). Integrating this expression from $S^{*}$ to $S_{0}$ yields the value $t^{*}$

[TABLE]

This expression can be used to predict the time $t^{*}$ from data.

Appendix B Well-posedness of the model

To prove the well-posedness, we rewrite the system (7) as the following abstract differential equation:

[TABLE]

where x:=\left(\begin{array}[]{c}s\\ i\end{array}\right), $A$ is the matrix given by

[TABLE]

and $f:R^{n}\times R^{n}\longrightarrow R^{2n}$ defined by

[TABLE]

and x_{0}:=\left(\begin{array}[]{c}S_{0}\\ I_{0}\end{array}\right).

It is clear that, the function $f$ is $L_{f}$-lipschitzian with $L_{f}$ depends only on $\beta$ and $\gamma$. Now, we formulate the well-posedness theorem, which is the main theorem of this section:

Theorem B.1

Given $x_{0}\in R^{n}$. Then, the equation (28) has a unique solution satisfying the following formula:

[TABLE]

**proof

**

Let $x_{0}\in R^{n}$ and $T>0$. Consider the mapping $\Gamma:C\longrightarrow C$ given by

[TABLE]

where $C:=C([0,T],R^{n})$. Let us prove that $\Gamma$ is a contraction. Indeed, let $u,v\in C$, then

[TABLE]

On the other hand

[TABLE]

Hence, by iterating for $n\geq 1$, we conclude that

[TABLE]

Now, for $n$ large enough,

[TABLE]

The mapping $\Gamma^{n}$ is a contraction. Therefore, by using the iterating fixed point theorem $\Gamma$ is also a contraction. Consequently, the system (13) has a unique solution which is given by (14). end  proof