Christophe Besse

Grégory Faye

Dynamics of epidemic spreading on connected graphs

Keywords: AMS classification: 34D05, 35Q92, 35B40, 92-10, 92D30 SIR model, graph, diffusion equation. 2

come

Introduction

Classical SIR compartment models are cornerstone models of epidemiology which allow one to study the evolution of an infected population at a given spatial scale (e.g. whole countries, regions, counties or cities). Such models date back to the pioneer work of Kermack and McKendrick [START_REF] Kermack | A contribution to the mathematical theory of epidemics[END_REF] and describe the evolution of susceptible (S) and infected (I) populations which eventually become removed (R) via systems of ordinary differential equations which typically take the form        S (t) = -τ S(t)I(t), I (t) = τ S(t)I(t) -ηI(t), R (t) = ηI(t), (1.1) where τ > 0 is a contact rate between susceptible and infected populations, and 1/η > 0 is the average infectious period; see [START_REF] Hethcote | The mathematics of infectious diseases[END_REF] for a review on SIR models. These models have been used in the past to reproduce data of epidemic outbreaks such as the bubonic plague [START_REF] Kermack | A contribution to the mathematical theory of epidemics[END_REF], malaria [START_REF] Mandal | Mathematical models of malaria -a review[END_REF], SARS influenza [START_REF] Centers For | Severe acute respiratory syndrome -Singapore, 2003[END_REF][START_REF] Magal | Final size of an epidemic for a two group SIR model[END_REF] and most recently COVID-19 [START_REF] Liu | Predicting the cumulative number of cases for the COVID-19 epidemic in China from early data[END_REF][START_REF] Magal | The parameter identification problem for SIR epidemic models: Identifying Unreported Cases[END_REF]21]; see also [START_REF] Magal | The parameter identification problem for SIR epidemic models: Identifying Unreported Cases[END_REF] for other applications.

• • • • • • • • • • • • • • • • • • • • • • • • • • • • •
In classical SIR models such as (1.1), interactions among the infected population are oversimplified, and when taken into account they typically involve transfer matrices of populations of infected between various uniform patches [START_REF] Magal | Final size of a multi-group SIR epidemic model: Irreducible and non-irreducible modes of transmission[END_REF][START_REF] Magal | Final size of an epidemic for a two group SIR model[END_REF][START_REF] Van Den Driessche | Reproduction numbers and sub-threshold endemic equilibria for compartmental models of disease transmission[END_REF]. Our interest lies in the understanding of the intricate interplay between spatial effects and heterogeneous interactions among infected populations. Schematically, we propose a model composed of cities linked by a given transportation network (roads, railroads or rivers), see Figure 1 for an illustration in the case of France. It will turn out that the appropriate theoretical framework will be graph theory where each vertices of the graph will be thought as the cities and the edges the lines of transportation. In a first approximation, we will assume that infected populations are only subject to spatial diffusion along the lines, as it is traditionally assumed in classical spatial SIR models [START_REF] Aronson | The asymptotic speed of propagation of a simple epidemic[END_REF][START_REF] Berestycki | Propagation of epidemics along lines with fast diffusion[END_REF][START_REF] Diekmann | Thresholds and travelling waves for the geographical spread of infection[END_REF][START_REF] Reluga | A two-phase epidemic driven by diffusion[END_REF]. As a consequence, in our model, the dynamics of the epidemic only takes place in the cities. Interactions are then modeled by flux exchanges between cities and lines where we assume that some fraction of infected individuals can either leave a city to be on a line, or leave a line and stop in a city, or pass from one line to another through a city. The typical question that we address here can easily be stated as follows. Given a connected graph of cities linked by roads and an initial configuration of infected individuals, how does the epidemic spread into the network and what is the eventual final configuration of the infected population? Our aim here is to gain insight into this spreading aspect at the fundamental mathematical level of a SIR type model that incorporates the possibility of infected individuals to travel along a specific given transportation network.

Our framework is at the crossroad of models that take into account lines of transportation such as recent reaction-diffusion models that study propagation of epidemics along lines with fast diffusion [START_REF] Berestycki | Propagation of epidemics along lines with fast diffusion[END_REF] and models that incorporate networks with more sophisticated interactions dynamics [START_REF] Ball | Epidemics on networks with preventive rewiring[END_REF]4,[START_REF] Britton | Epidemics on random graphs with tunable clustering[END_REF][START_REF] Spricer | An epidemic model on a weighted network[END_REF]. On a formal level, our proposed model can be thought of as being a one-dimensional version of the planar reaction-diffusion system of [START_REF] Berestycki | Propagation of epidemics along lines with fast diffusion[END_REF] with a line of fast diffusion in the case of one city and one line of transportation. Actually, the graph structure of the transportation network provides a natural embedding into a planar spatial domain. From a mathematical point of view, our model shares also some similarities with the PDE-ODE model of [START_REF] David | A novel approach to modelling the spatial spread of airborne diseases: an epidemic model with indirect transmission[END_REF] which studies the spread of airborne diseases where the movement of pathogens in the air is assumed to follow a linear diffusion.

Model formulation and main results

Throughout, we denote by G = (V, E) a compact metric graph, i.e. a collection of vertices V and edges E and further assume that G is finite and connected. Each edge e ∈ E is identified with a segment Ω e = [0, e ] with e ∈ (0, ∞), where e is the finite length of the edge. A real valued function u : G -→ R is a collection of one dimensional maps defined for each edge e ∈ E:

A SIR model on compact connected graph

Given a graph G, we let X v (t) := (S v (t), I v (t), R v (t)) ∈ R 3 , for each v ∈ V, where S v (t) represents the population of susceptible individuals, I v (t) the population of infected individuals and R v (t) the population of susceptible individuals at vertex v ∈ V and time t > 0. We assume that X v evolves according to a SIR model of the form

           S v (t) = -τ v S v (t)I v (t), I v (t) = τ v S v (t)I v (t) -η v I v (t) + e∼v α v e u e (t, v) -λ v I v (t), R v (t) = η v I v (t), (2.1)
where τ v , η v > 0 are the intrinsic parameters of the epidemic which may depend on the vertex v. The contribution -λ v I v (t) in the right-hand side of the equation for the infected population traduces the fact that infected individuals can leave the vertex v to incident edges whereas e∼v α v e u e (t, v) reflects the contribution of incoming infected individuals from incident edges. Here, e ∼ v denotes the edges incident to the vertex v and

λ v := e∼v λ v e ,
such that λ v e I v (t) infected individuals leave vertex v to edge e. We have assumed that only the infected population is subject to movement, and we think of S v being an ambiant population whose movement does not affect its distribution. We recover the standard SIR model (1.1) by considering the trivial graph G = ({v}, ∅). Throughout the manuscript, we will assume the following standing hypothesis on the coefficients α v e and λ v e in (2.1).

Hypothesis 2.1. For each (e, v) ∈ E × V we assume that α v e ∈ (0, 1) and λ v e ∈ (0, 1), together with e∼v λ v e ∈ (0, 1) and e∼v α v e ∈ (0, 1).

Next, for each e ∈ E, we let d e > 0 and we assume that u e evolves according to

∂ t u e (t, x) = d e ∂ 2 x u e (t, x), t > 0, x ∈ • Ω e . (2.2)
Assuming that infected individuals have local diffusion along the edges of the graph is a first approximation, and this can be viewed as a limiting Brownian movement of individuals. We shall come back to this modeling hypothesis later in the manuscript, but possible extensions could be to incorporate nonlocal diffusion or transport terms.

It now remains to model the exchanges of infected individuals at the vertices. Fo each v ∈ V, we associate an integer δ v ≥ 1 which we refer to as its degree (i.e. number of edges incident to the vertex v). We define u v (t) ∈ R δv as the column vector function

u v (t) := (u e (t, v)) e∼v ,
where we recall that e ∼ v denotes the edges incident to the vertex v, and thus u e (t, v) is the corresponding limit value of u e at x = v. Define also ∂ n u v (t) ∈ R δv as the column vector function

∂ n u v (t) := (∂ n u e (t, v)) e∼v ,
where ∂ n u e (t, v) is the outwardly normal derivative of u e at the vertex v. Our boundary conditions at the vertex v that link (2.1) and (2.2) are described by

D v ∂ n u v (t) + K v u v (t) = Λ v I v (t), (2.3)
where for each e ∼ v.

D v ∈ M δv (R) is the diagonal matrix D v = diag[(d e ) e∼v ] and K v ∈ M δv (R)
Remark 2.3. If the exchanges between the edges are symmetric, that is for each v ∈ V the matrices N v are symmetric, that is ν v e,e = ν v e ,e , ∀(e, e ) ∈ E × E, then Hypothesis 2.2 is automatically satisfied.

Initial configuration

Finally, we complement our coupled PDE-ODE (2.1)-(2.2)-(2.3) with some initial conditions. We assume that at t = 0, we have

u(t = 0, •) = u 0 ∈ BC(G, R),
such that for e ∈ E, u 0 e (x) ≥ 0, x ∈ Ω e .

On the other hand, for the ODE system (2.1), we suppose that

(S v (t = 0), I v (t = 0), R v (t = 0)) = (S 0 v , I 0 v , 0) ≥ 0, ∀v ∈ V.
We further assume that (2.3) is satisfied at t = 0

D v ∂ n u 0 v + K v u 0 v = Λ v I 0 v ,
with obvious notations u 0 v := (u 0 e (v)) e∼v and

I 0 v := (I 0 v , • • • , I 0 v ) t .
Last, we impose that the initial total population of infected individuals is strictly positive,

v∈V I 0 v > 0,
and that susceptible individuals are initially present at each vertex of the graph

S 0 v > 0, ∀v ∈ V.
This in turn implies that the total population is initially

M 0 := e∈E Ωe u 0 e (x)dx + v∈V S 0 v + I 0 v > 0.

Conservation of total population

Assuming that there is a solution to (u, (

X v ) v∈V ) to (2.1)-(2.2)-(2.
3), we have that the total mass of the system M (t) defined as

M (t) := e∈E Ωe u e (t, x)dx + v∈V (S v (t) + I v (t) + R v (t))
is a conserved quantity and thus independent of t.

To see that, we first remark that

S v (t) + I v (t) + R v (t) = e∼v α v e u e (t, v) -λ v I v (t) = A v u v (t), 1 δv -Λ v I v (t), 1 δv with 1 δv := (1, • • • , 1) t ∈ R δv ,
and •, • is the standard Euclidean inner product on R δv . On the other hand let us define

m(t) := e∈E Ωe u e (t, x)dx,
and assume that u is a classical solution of (2.2), which we will prove in the next section, and compute

m (t) = e∈E Ωe ∂ t u e (t, x)dx = e∈E d e [∂ x u e (t, x)] ∂Ωe = v∈V D v ∂ n u v (t), 1 δv = v∈V Λ v I v (t) -K v u v (t), 1 δv = v∈V Λ v I v (t) -A v u v (t), 1 δv - v∈V N v u v (t), 1 δv =0 = - v∈V S v (t) + I v (t) + R v (t) .
The fact that

v∈V N v u v (t), 1 δv = 0
is a direct consequence on the specific structure of each matrix N v and the fact that the sum of each column is zero. We therefore conclude that M (t) = 0 and

e∈E Ωe u e (t, x)dx + v∈V (S v (t) + I v (t) + R v (t)) = M 0 , ∀t ≥ 0.
Biological interpretation. Our model is thus consistent with the conservation of the total population as it is traditionally the case for SIR model in the case of zero natality/mortality rate. The exchanges between the vertices and the edges exactly compensate each other as is natural.

Main results and outline of the paper

We now present our main results regarding our model (2.1)-(2.2)- (2.3). At this stage of the presentation, we remain formal and refer to the following sections for precise statements and assumptions.

Main result 1: Existence and uniqueness of classical solutions. We prove in Theorem 1 below that for each well prepared initial condition our model (2.1)-(2.2)-(2.3) admits a unique positive classical solution which is global in time. We remark that the system (2.1)-(2.2)-(2.3) is not standard as it couples a system of PDEs to ODEs at each vertices through inhomogeneous Robin boundary conditions. As a consequence, the existence and uniqueness of classical solutions has to be proved. This analysis is conducted in Section 3.

Main result 2:

Long time behavior of the solutions. We fully characterize the long time behavior of the unique solution of our model. More precisely, we prove that the final total population of infected individuals at each vertex, denoted by I ∞ v , is a well defined quantity: 0

< I ∞ v < ∞ for v ∈ V and (I ∞ v )
v∈V are solutions of a system of c V + 1 implicit equations, where c V stands for the cardinal of V, which belong to the parametrized submanifold

v∈V S 0 v e -τvI ∞ v + η v I ∞ v = M 0 ,
where M 0 is the initial total mass. We refer to Theorem 2 for a precise statement. We also present further qualitative results on the final total configuration (I ∞ v ) v∈V in the fully symmetric case where we obtain closed form formula (see Lemma 4.3) and in the case of two vertices where we manage to obtain sharp bounds on the final total populations of infected individuals (see Lemma 4.4). In each case, we manage to relate these quantities to standard basic and effective reproductive number for classical SIR model. The aforementioned results are proved in Section 4.

Main result 3: A mass preserving semi-implicit numerical scheme. We propose and analyze a semi-implicit in time numerical scheme based on finite differences in space which has the property to preserve a discrete total mass associated to the discretization. We prove that if the time discretization constant is smaller than a universal constant depending only on the parameters of the system (and not on the space discretization constant) and if N v is symmetric for each v ∈ V, then our mass preserving semi-implicit numerical scheme is well-posed and preserves the positivity of the solutions. We refer to Section 5 for a presentation of the numerical scheme and Theorem 3 for a precise statement of our main result.

Main result 4: Numerical results for various types of graphs. We illustrate our theoretical findings with selection of numerical simulations for various types of graphs in Section 6. We respectively study the case of 2 vertices and 1 edge, 3 vertices and 3 edges (closed graph), 4 vertices and 3 edges (star-shape graph) and N + 1 vertices and N edges with N being arbitrarily large (lattice graph). Most notably, in the last case, we show the propagation of the epidemics across the vertices of the graph in the form of a traveling wave.

The Cauchy problem: existence and uniqueness of classical solutions

This section is devoted to the proof of the following main theorem which guarantees that our model is well-posed.

Theorem 1. For each (S 0 v , I 0 v ) ≥ 0 with S 0 v > 0, v∈V I 0 v > 0 and u 0 ∈ BC(G, R + ) that satisfy the boundary condition (2.3), there exists a unique positive global solution

(S v , I v , R v ) ∈ C 1 (R + , R + × R + × R + ) and u ∈ C 1,2 (R + * × G, R + ).
The proof of Theorem 1 is divided into two parts. We first prove the existence of positive global classical solutions and then show that such constructed solutions are unique. We look for solutions that satisfy (2.1)-(2.2)-(2.3) in the classical sense, and we always assume that (S 0 v , I 0 v ) ≥ 0 with S 0 v > 0, v∈V I 0 v > 0 and u 0 ∈ BC(G, R + ), that is for all e ∈ E, u 0 e ≥ 0 is bounded continuous on Ω e . We further assume that the initial conditions satisfy the boundary condition (2.3). We remark that the system (2.1)-(2.2)-(2.3) is not standard as it couples a system PDEs to ODEs at each vertices through inhomogeneous Robin boundary conditions. As a consequence, the well-posedness of the Cauchy problem has to be proved. Remark 3.1. Our existence and uniqueness result extends trivially in the case that parameters τ v > 0, α v e ∈ (0, 1), λ v e ∈ (0, 1) and ν v e,e ∈ [0, 1) are continuous functions of time satisfying τ v (t) > 0, α v e (t) ∈ (0, 1), λ v e (t) ∈ (0, 1) and ν v e,e (t) ∈ [0, 1) together with Hypotheses 2.1-2.2 verified at all times t > 0.

Existence

In this section, we construct a classical solution to (2.1)-(2.2)-(2.3) through a limiting argument. We will obtain a solution (u, (X v ) v∈V ) has the limit of a subsequence of solution ((u n , (X n v ) v∈V )) n≥0 of the following problems

               dS n v (t) dt = -τ v S n v (t)I n v (t), dI n v (t) dt = τ v S n v (t)I n v (t) -(η v + λ v )I n v (t) + e∼v α v e u n-1 e (t, v), dR n v (t) dt = η v I n v (t), t > 0, ∀v ∈ V, (3.1)
with

D v ∂ n u n v (t) + K v u n v (t) = Λ v I n v (t), t > 0, ∀v ∈ V, (3.2)
and

∂ t u n e (t, x) = d e ∂ 2 x u n e (t, x), t > 0, x ∈ • Ω e , ∀e ∈ E. (3.3)
starting from u 0 ∈ BC(G, R + ) and (X 0 v ) v∈V . Note that (3.1)-(3.2)(3.3) is supplemented by the same initial condition (u 0 , (X 0 v ) v∈V ) at each step. We proceed along three main steps.

Step #1: solvability of (3. ) v∈V ) such that for each t → u n-1 e (t, v) is continuous, then we get the existence of a unique solution of (3.1) which is C 1 in time. Next we solve the system of PDEs (3.3)-(3.2) whose coupling comes from the boundary conditions and owing that now the right-hand side of (3.2) can be seen as given inhomogeneous term of class C 1 in time. As both D v and K v are invertible matrices, we get the existence of a classical solution u n ∈ C 1,2 which then ensures that t → u n e (t, v) is continuous.

Step #2: a priori estimates. Let 0 < T < 1 be fixed. We first show by a recursive argument that 0 < S n v , 0 ≤ I n v , 0 ≤ R n v for each v ∈ V and 0 ≤ u n e for each e ∈ E. It is trivial at n = 0. Let assume that is it true at n -1. We start from (3.1) and a direct integration gives Next, from the positivity of solutions, we obtain some uniform L ∞ bounds. More precisely, we claim that there exists a constant K > 0 depending only on (T,

S n v (t) = S 0 v e -τv t 0 I n v (s)ds > 0, I n v (t) = I 0 v e -(ηv+λv)t+ t 0 S n v (s)ds + e∼v α v e t 0 e -(ηv+λv)(t-s)+ t s S n v (τ )dτ u n-1 e (s, v)ds ≥ 0, R n v (t) = η v t 0 I n v (s)ds ≥ 0. Now owing that 0 ≤ I n v for each v ∈ V,
I 0 v , S 0 v , u 0 ∞ ) such that 0 ≤ S n v (t), I n v (t), R n v (t) ≤ K, and dS n v (t) dt , dI n v (t) dt , dR n v (t) dt ≤ K 0 < t ≤ T, and 0 ≤ u n e (t, x) ≤ K, 0 < t ≤ T, x ∈ Ω e .
First, using (3.1) we obtain that

dS n v (t) dt + dI n v (t) dt + dR n v (t) dt = -λ v I n v (t) + e∼v α v e u n-1 e (t, v), which gives that 0 ≤ S n v (t) + I n v (t) + R n v (t) ≤ S 0 v + I 0 v + T u n-1 ∞ , 0 < t ≤ T, together with dS n v (t) dt , dI n v (t) dt , dR n v (t) dt ≤ C u n-1 ∞ (1 + T + T 2 u n-1 ∞ ) 0 < t ≤ T,
which in turn implies that

u n ∞ ≤ CT u n-1 ∞ (1 + T + T 2 u n-1 ∞ ),
for C, C > 0 only depend on the initial condition (u 0 , (X 0 v ) v∈V ) and the parameters of the system. We now claim that by induction, we have for all 0 < t ≤ T ,

0 ≤ S n v (t) + I n v (t) + R n v (t) ≤ Ĉ an p=0 T p , dS n v (t) dt , dI n v (t) dt , dR n v (t) dt ≤ Ĉ 2an p=0 T p , u n ∞ ≤ Ĉ 2an p=0 T p+1 , with a n = 2 + 2a n-1 for n ≥ 2 and a 1 = 1 for some Ĉ > 0 only depending on (u 0 , (X 0 v ) v∈V ). As 0 < T < 1, we get that 0 ≤ S n v (t) + I n v (t) + R n v (t) ≤ ĈT , 0 < t ≤ T, together with u n ∞ , dS n v (t) dt , dI n v (t) dt , dR n v (t) dt ≤ ĈT ,
for some constant ĈT > 0 depending on (T, u 0 , (X 0 v ) v∈V ).

Step #3: existence of a solution. Parabolic Schauder estimates give that the time derivative and the space derivatives up to order 2 of u n are uniformly Hölder continuous in compact sets. As a consequence, we can use the Arzela-Ascoli theorem to show that (u n , (

X n v ) v∈V ) converges (up to sequences) toward (u, (X v ) v∈V ) in C 1,2 loc ((0, T ) × G) × C 1 loc ((0, T )) × C 1 loc ((0, T )) × C 1 loc ((0, T )). Passing to the limit n → +∞ in (3.1)-(3.2)-(3.3) we get that (u, (X v ) v∈V ) satisfies (2.1)-(2.2) subject to boundary conditions (2.3).
As a by product of the proof we get that for the just constructed solution (u, (X v ) v∈V ) we have the uniform bounds:

0 < S v (t) ≤ S 0 v , 0 < t ≤ T, v ∈ V, and 0 ≤ I v (t), R v (t) 0 < t ≤ T, together with 0 ≤ u e (t, x), 0 < t ≤ T, x ∈ Ω e , e ∈ E.
The fact that I v (t) ≥ 0 implies thanks to the strong maximum principle that in fact

0 < u e (t, x), 0 < t ≤ T, x ∈ Ω e , e ∈ E,
which in turn gives that I v (t) > 0 for each v ∈ V since

I n v (t) = I 0 v e -(ηv+λv)t+ t 0 S n v (s)ds + e∼v α v e t 0 e -(ηv+λv)(t-s)+ t s S n v (τ )dτ u n-1 e (s, v)ds > 0.
Finally, we use the conservation of mass which tells us that

e∈E Ωe u e (t, x)dx + v∈V (S v (t) + I v (t) + R v (t)) = M 0 > 0, 0 < t ≤ T,
such that both I v (t) and R v (t) are uniformly bounded in time, together with their derivatives. This also implies that there exists a constant M > 0, depending only (u 0 , (

X 0 v ) v∈V ) such that 0 < u(t, x) ≤ M, 0 < t ≤ M, x ∈ Ω e .
Using again parabolic regularity, we obtain the solution (u, (X v ) v∈V ) is global in time and satisfies

(2.1)-(2.2)-(2.
3) in the classical sense.

Uniqueness

Let assume that (u, (X v ) v∈V ) and (ũ, ( Xv ) v∈V ) are two classical solutions to (2.2)-( 2.3)-(2.1) starting from the same initial datum (u 0 , (X 0 v ) v∈V ). We denote (û, ( Xv ) v∈V ) where for each e ∈ E ûe = u eũe ,

and each v ∈ V Xv = ( Ŝv , Îv , Rv ) = (S v -Sv , I v -Ĩv , R v -Rv ).
By linearity, we get that for e ∈ E

∂ t ûe = d e ∂ 2 x ûe , t > 0, x ∈ • Ω e , together with D v ∂ n ûv (t) + K v ûv (t) = Λ v Îv (t), t > 0, v ∈ V.
On the other, one computes that Xv satisfies for each v ∈ V

           Ŝ v (t) = -τ v S v (t) Îv (t) + Ŝv (t) Ĩv (t) , Î v (t) = τ v S v (t) Îv (t) + Ŝv (t) Ĩv (t) -η v Îv (t) + e∼v α v e ûe (t, v) -λ v Îv (t), R v (t) = η v Îv (t),
We define the energy

E (t) := 1 2 e∈E Ωe (û e (t, x)) 2 dx + 1 2 v∈V Ŝv (t) 2 + Îv (t) 2 + Rv (t) 2 ,
and note that E (0) = 0 by definition. Next, differentiating E (t), we obtain

E (t) = e∈E Ωe ûe (t, x)∂ t ûe (t, x)dx + v∈V Ŝv (t) Ŝ v (t) + Îv (t) Î v (t) + Rv (t) R v (t) := E u (t) + E X (t).
On the one hand, we have

E u (t) = e∈E Ωe ûe (t, x)∂ t ûe (t, x)dx = e∈E d e Ωe ûe (t, x)∂ 2 x ûe (t, x)dx = - e∈E d e Ωe (∂ x ûe (t, x)) 2 dx + e∈E d e [û e (t, x)∂ t ûe (t, x)] ∂Ωe ≤ v∈V D v ∂ n ûv (t), ûv (t) = v∈V Λ v Îv (t) -K v ûv (t), ûv (t) = v∈V Λ v Îv (t) -A v ûv (t), ûv (t) - v∈V N v ûv (t), ûv (t) ≤ v∈V Λ v Îv (t) -A v ûv (t), ûv (t) ,
as N v is symmetric positive. On the other hand, we compute

E u (t) = v∈V Ŝv (t) Ŝ v (t) + Îv (t) Î v (t) + Rv (t) R v (t) = v∈V τ v S v (t) Îv (t) Ŝv (t) + Ŝv (t) 2 Ĩv (t) + S v (t) Îv (t) 2 + Ŝv (t) Ĩv (t) Îv (t) - v∈V η v Îv (t) 2 - v∈V λ v Îv (t) 2 + v∈V Îv (t) e∼v α v e ûe (t, v) + v∈V η v Rv (t) Îv (t) ≤ CE (t) + v∈V Îv (t) e∼v α v e ûe (t, v),
where C > 0 is some large positive constant. Next, we see that

v∈V Îv (t) e∼v α v e ûe (t, v) = v∈V Îv (t), A v ûv (t) ,
such that we obtain

E (t) ≤ CE (t) + v∈V Λ v Îv (t) -A v ûv (t), ûv (t) + v∈V Îv (t), A v ûv (t) .
Next, if we denote ŵv (t

) := 1 2 A -1/2 v (Λ v + A v ) Îv (t) -A 1/2 v ûv (t), we compute 0 ≤ v∈V ŵv (t), ŵv (t) = 1 4 v∈V A -1 v (Λ v + A v ) 2 Îv (t), Îv (t) - v∈V (Λ v + A v ) Îv (t), ûv (t) + v∈V A v ûv (t), ûv (t) = 1 4 v∈V A -1 v (Λ v + A v ) 2 Îv (t), Îv (t) - v∈V Λ v Îv (t) -A v ûv (t), ûv (t) - v∈V Îv (t), A v ûv (t) .
As a consequence, we get

E (t) ≤ CE (t) + 1 4 v∈V A -1 v (Λ v + A v ) 2 Îv (t), Îv (t) - v∈V ŵv (t), ŵv (t) ≤ CE (t),
for some C > 0 and we conclude that E (t) = 0 for all time which then implies that û = 0 and Xv = 0.

4 Long-time behavior of the solutions

+ v∈V (S v (t) + I v (t) + R v (t)) = M 0 > 0, ∀t > 0.

Final total populations: general results

As 0 < S v (t) < S 0 v and S v (t) is strictly decreasing, it asymptotically converges towards a limit that we denote

S ∞ v := lim t→+∞ S v (t), v ∈ V.
Furthermore, as R v (t) is strictly increasing and uniformly bounded, it asymptotically converges towards a limit that is denoted

0 < R ∞ v := lim t→+∞ R v (t) < ∞, v ∈ V.
But as for each t > 0

R v (t) = η v t 0 I v (s)ds,
this implies that

I v (t) := t 0 I v (s)ds -→ I ∞ v = R ∞ v η v < ∞ as t → +∞,
which in turn proves that

I ∞ v = lim t→+∞ I v (t) = 0.
If one recall the notation m(t) for the total population on the edges then we have

m(t) = e∈E Ωe u e (t, x)dx,
and it verifies

v∈V (S v (t) + I v (t) + R v (t)) + m(t) = v∈V S 0 v + I 0 v + e∈E Ωe u 0 e (x)dx.
The above computations shows that m(t) has a limit as t -→ +∞, that we denote m ∞ and which satisfies v∈V

(S ∞ v + R ∞ v ) + m ∞ = v∈V S 0 v + I 0 v + e∈E Ωe u 0 e (x)dx. (4.1)
We shall also keep in mind that

S ∞ v = S 0 v e -τvI ∞ v = S 0 v e -τv ηv R ∞ v , or R ∞ v = - η v τ v ln S ∞ v S 0 v , v ∈ V
And so if we introduce the function Ψ v (x) := x -ηv τv ln x, then the above conservation of mass can be written as

v∈V Ψ v (S ∞ v ) + m ∞ = v∈V I 0 v + Ψ v (S 0 v ) + e∈E Ωe u 0 e (x)dx.
On the other, one can compute that

dm(t) dt = v∈V λ v I v (t) - v∈V A v u v (t), 1 δv , such that m(t) + v∈V t 0 A v u v (s), 1 δv ds = e∈E Ωe u 0 e (x)dx + v∈V λ v I v (t).
Now, as m(t) and each I v (t) are convergent we deduce that all

t 0 u v (s)ds are also convergent so that m ∞ + v∈V ∞ 0 A v u v (s), 1 δv ds = e∈E Ωe u 0 e (x)dx + v∈V λ v I ∞ v , (4.2)
and

∞ 0 u v (s)ds < ∞, v ∈ V, which proves that u v (t) -→ 0 as t → +∞, v ∈ V.
And the boundary conditions imply that

∂ n u v (t) -→ 0 as t → +∞, v ∈ V.
Next, we define the sequence of functions u n e (t, x) = u e (t + n, x) for each e ∈ E and X n v (t) = X v (t + n) for each v ∈ V which are uniformly bounded such that one can extract a convergent subsequence. On the one hand we have that lim

n→∞ X n v (t) = X ∞ v = (S ∞ v , 0, R ∞ v ) and on the other if u ∞ e (t, x) = lim n→∞ u n e (t, x) it is solution of ∂ t u ∞ e (t, x) = d e ∂ 2 x u ∞ e (t, x),
with the boundary conditions

∂ n u ∞ v (t) = u ∞ v (t) = 0 δv , v ∈ V.
This then shows that u ∞ e (t, x) = 0, t > 0 and x ∈

• Ω e for each e ∈ E. As there is unicity of the limit, we deduce that lim t→+∞ u e (t, x) = 0, e ∈ E.

From which we also get that m ∞ = 0 and that

v∈V Ψ v (S ∞ v ) = v∈V I 0 v + Ψ v (S 0 v ) + e∈E Ωe u 0 e (x)dx. This implies that each Ψ v (S ∞ v ) is bounded, we get that S ∞ v > 0 for all v ∈ V. We also get from (4.2), that v∈V ∞ 0 A v u v (s), 1 δv ds = e∈E Ωe u 0 e (x)dx + v∈V λ v I ∞ v .
Finally, we use the fact that

dI v (t) dt + dS v (t) dt - η v + λ v τ v d ln S v (t) dt = e∼v α v e u e (t, v), v ∈ V,
to obtain that

I v (t) + S v (t) - η v + λ v τ v ln S v (t) -I 0 v -S 0 v + η v + λ v τ v ln S 0 v = e∼v α v e t 0 u e (s, v)ds, v ∈ V.
As a consequence, the final total populations of infected individuals at each vertices satisfy the following scalar differential equation

dI v (t) dt = S 0 v 1 -e -τvIv(t) -η v I v (t) + I 0 v + e∼v α v e t 0 u e (s, v)ds -λ v I v (t), v ∈ V. (4.3)
Passing to the limit as t → +∞, we get

0 = S 0 v 1 -e -τvI ∞ v -η v I ∞ v + I 0 v + e∼v α v e ∞ 0 u e (s, v)ds -λ v I ∞ v , v ∈ V.
To summarize, we have proved the following result.

Theorem 2. For each (S 0 v , I 0 v ) ≥ 0 with S 0 v > 0, v∈V I 0 v > 0 and u 0 ∈ BC(G, R + ) that satisfy the boundary condition (2.3), the long time behavior of the unique corresponding solution (u, (X v ) v∈V ) is given by

lim t→+∞ u e (t, x) = 0, x ∈ Ω e , e ∈ E, with ∞ 0 u e (s, v)ds < +∞, (v, e) ∈ G,
and

lim t→+∞ (S v (t), I v (t), R v (t)) = S 0 v e -τvI ∞ v , 0, η v I ∞ v , v ∈ V, (4.4)
where the final total populations of infected individuals 0 < I ∞ v < ∞ at each vertices v ∈ V are solutions of the system

         S 0 v e -τvI ∞ v + η v I ∞ v = I 0 v + S 0 v + e∼v α v e ∞ 0 u e (s, v)ds -λ v I ∞ v , v ∈ V, v∈V ∞ 0 A v u v (s), 1 δv ds = e∈E Ωe u 0 e (x)dx + v∈V λ v I ∞ v .
(4.5)

As a consequence, (I ∞ v ) v∈V belongs to the parametrized submanifold given by

v∈V S 0 v e -τvI ∞ v + η v I ∞ v = M 0 . (4.6) Remark 4.1. Equivalently, (S ∞ v ) v∈V belongs to the parametrized submanifold given by v∈V S ∞ v - η v τ v log (S ∞ v ) + η v τ v log S 0 v = M 0 , (4.7)
and (R ∞ v ) v∈V belongs to the parametrized submanifold given by

v∈V S 0 v exp(-τ v /η v R ∞ v ) + R ∞ v = M 0 . (4.8)
The equations (4.6), (4.7), and (4.8)

also read v∈V (S ∞ v + R ∞ v ) = M 0 , which is nothing but (4.1) since we have proved that m ∞ = 0. Remark 4.2. If we assume that τ = τ v > 0 and η = η v > 0 are independent of v ∈ V and let S v = τ /η S v , R ∞ v = exp(-τ /η R ∞ v ) and I ∞ v = exp(-τ I ∞ v ).
Then, equations (4.6), (4.7), and (4.8) are respectively equivalent to

v∈V exp S ∞ v S 0 v S ∞ v = exp τ η M 0 , v∈V exp S 0 v R ∞ v R ∞ v = exp τ η M 0 ,
and

v∈V exp S 0 v I ∞ v I ∞ v = exp τ η M 0 .
The common right hand side features τ η M 0 that is nothing but the traditional basic reproductive number R 0 .

Final total populations of infected individuals: further properties

The aim of this section is to present further qualitative results on the final total configuration (I ∞ v ) v∈V in the fully symmetric case where one can obtain closed form formula and in the case of two vertices where we manage to obtain sharp bounds on the final total populations of infected individuals. In each case, we manage to relate these quantities to standard basic and effective reproductive number for classical SIR model [START_REF] Diekmann | On the definition and the computation of the basic reproduction ratio R0 in models for infectious diseases in heterogeneous populations[END_REF].

Fully symmetric case. We assume that the length e of every edge e ∈ E is equal to a reference length . For every e ∈ E, the diffusion coefficient d e is equal to d. We moreover suppose that for every vertex v ∈ V,

S 0 v = S 0 , I 0 v = I 0 and R 0 v = R 0 .
We also assume that τ = τ v > 0 and η = η v > 0 are independent of v ∈ V. In the same spirit, λ v e = λ and α v e = α for every e ∈ E and v ∈ V. We also assume ν v e i ,e j = ν for every edges incident to the vertex v. Finally, the components u 0 e of initial condition on each edges e ∈ E are supposed to be even with respect to the center of the interval Ω e = [0, ]. Thanks to all these assumptions, I ∞ v does not depend on the vertex v ∈ V and we set

I ∞ v = I ∞ .
Let us recall the notation c V for the cardinal of the set V. The parametrized submanifold given by (4.6) becomes

S ∞ + R ∞ = S 0 e -τ I ∞ + ηI ∞ = M 0 ,
where M 0 = M 0 /c V . We can transform this relation as

S 0 e -τ I ∞ + η τ τ I ∞ -M 0 = 0. (4.9)
Let I = -τ I ∞ . We have to solve

S 0 e I - η τ I -M 0 = 0.
The solutions are given in terms of Lambert W function that is the multivalued inverse relation of the function f (w) = we w for w ∈ C [START_REF] Corless | On the LambertW function[END_REF]. Let us recall how to compute the real solutions of the equation αe

x + βx + γ = 0 for (α, β, γ) ∈ R * × R * × R. Let ∆ = α/β exp(-γ/β) be the discriminant. If ∆ ≥ 0 or ∆ = -exp(-1
), the solution is unique and x = -W 0 (∆)γ/β where W 0 is the principal branch. If ∆ ∈ (-exp(-1), 0), there are two solutions x 0 = -W 0 (∆)γ/β and x -1 = -W -1 (∆)γ/β, where W -1 is another branch. When ∆ <exp(-1), there is no solution.

In our symmetric case, the discriminant writes

∆ = - S 0 τ η exp - M 0 τ η .
Since ∆ < 0, there exist solutions to (4.9) if ∆ ≥exp(-1), which is equivalent to

exp M 0 τ η -1 ≥ S 0 τ η . (4.10)
We recall that when we consider the standard SIR model (meaning in the context of this paper that we consider an isolated vertex), we can define the effective reproductive number R e and the basic reproductive number R 0 respectively given by R e := S 0 τ η , and

R 0 := M 0 τ η , (4.11)
see [START_REF] Diekmann | On the definition and the computation of the basic reproduction ratio R0 in models for infectious diseases in heterogeneous populations[END_REF][START_REF] Hethcote | The mathematics of infectious diseases[END_REF][START_REF] Van Den Driessche | Reproduction numbers and sub-threshold endemic equilibria for compartmental models of disease transmission[END_REF] for further properties of effective and basic reproductive numbers. If we denote R 0 = M 0 τ /η, the equation (4.10)

reads exp R 0 -1 ≥ R e .
This inequality is satisfied as long as S 0 ≤ M 0 , which is always true since M 0 = e∈E Ωe u 0 e (x)dx+ c V S 0 + I 0 ≥ c V S 0 . Since ∆ = -R e exp(-R 0 ), the solutions are

I 0,-1 = -W 0,-1 -R e exp(-R 0 ) -R 0 ,
and so

I ∞ 0,-1 = W 0,-1 -R e exp(-R 0 ) /τ + R 0 /τ.
Both W 0,-1 -R e exp(-R 0 ) < 0. However, we can show that I ∞ 0 > 0 and I ∞ -1 < 0. Thus, the only possibility is

I ∞ = W 0 -R e exp(-R 0 ) /τ + R 0 /τ.
We also have access to S ∞ and R ∞ thanks to (4.4). Since exp(-W 0 (x)) = W 0 (x)/x, we obtain

S ∞ = - η τ W 0 -R e exp(-R 0 ) , and
R ∞ = η τ W 0 -R e exp(-R 0 ) + M 0 .
We can summarize these results in the following lemma.

Lemma 4.3 (Fully symmetric case.). Assume that our model is fully symmetric, then the final total population of infected individuals as given by Theorem 2 is independent on the vertex that is

I ∞ v = I ∞ for each v ∈ V
, and I ∞ has the following closed form formula

I ∞ = W 0 (-R e exp(-R 0 /c V )) τ + R 0 c V τ ,
where R e and R 0 are respectively the effective and basic reproductive number defined in (4.11) and c V the cardinal of V. See Figure 3 for an illustration.

0 4 1 0 I v1 0 4 2 0 I v2 [I 0 v1 ; I 0 v2 ] [I 1 v1 ; I 1 v2 ] I v 2 = I v 1 ' !1 ' 0 S v1 ' !1 ' 0 S v2 [S 0 v1 ; S 0 v2 ] [S 1 v1 ; S 1 v2 ] S v 2 = S v 1 Figure 3: Schematic visualisation (red star) of I ∞ = I ∞ v1 = I ∞ v2 , resp. S ∞ = S ∞ v1 = S ∞ v2
, in the (I v1 , I v2 ) plane, resp. in the (S v1 , S v2 )-plane, in the fully symmetric case. The asymptotic value I ∞ , resp. S ∞ , lies at the intersection of the diagonal I v1 = I v2 , resp. S v1 = S v2 , and the implicit curve given by (4.6), resp. (4.7).

Case of two vertices. In this simple case, it is possible to build explicit formulas to deal with the implicit submanifold equations (4.6), (4.7), and (4.8

). Let R 0,v k := M 0 τ v k /η v k and R e,v k := S 0 v k τ v k /η v k , k = 1
, 2 be respectively the local to vertex v k basic and effective reproductive number. Then,

S ∞ v 1 = - η v 1 τ v 1 W   -exp (-R 0,v 1 ) R e,v 1 (R e,v 2 ) τv 1 ηv 2 τv 2 ηv 1 exp S ∞ v 2 τ v 1 /η v 1 S ∞ v 2 τ v 2 /η v 2 τv 1 ηv 2 τv 2 ηv 1   , (4.12)
where the Lambert W function W can be either W 0 or W -1 . Indeed, the argument of W being negative, two solutions have to be considered. We obviously also have

S ∞ v 2 = - η v 2 τ v 2 W   -exp (-R 0,v 2 ) R e,v 2 (R e,v 1 ) τv 2 ηv 1 τv 1 ηv 2 exp S ∞ v 1 τ v 2 /η v 2 S ∞ v 1 τ v 1 /η v 1 τv 2 ηv 1 τv 1 ηv 2   , (4.13)
Due to the definition of the domain of the Lambert W function, the argument has to be greater thanexp(-1). So, the following inequality must be satisfied for

S ∞ v 2 (respectively of S ∞ v 1 ) -exp (-R 0,v 1 ) R e,v 1 (R e,v 2 ) τv 1 ηv 2 τv 2 ηv 1 S ∞ v 2 τ v 2 /η v 2 τv 1 ηv 2 τv 2 ηv 1 exp -S ∞ v 2 τ v 1 /η v 1 ≥ -exp(-1).
Solving the equality part of this inequality, we find that

S ∞ v 2 = - η v 2 τ v 2 W -(R e,v 1 ) τv 2 ηv 1 τv 1 ηv 2 R e,v 2 exp τ v 2 η v 1 τ v 1 η v 2 (1 -R 0,v 1 )
.

This equation has to be verified both for W 0 and W -1 . Let Σ v 2 0,-1 be defined by

Σ v 2 0,-1 := - η v 2 τ v 2 W 0,-1 (A v 2 ) ,
where

A v 2 = (R e,v 1 ) τv 2 ηv 1 τv 1 ηv 2 R e,v 2 exp τ v 2 η v 1 τ v 1 η v 2 (1 -R 0,v 1 ) (4.14)
Then, the domain of

S ∞ v 1 as a function of S ∞ v 2 is S ∞ v 2 ∈ min Σ v 2 -1 , Σ v 2 0 , max Σ v 2 -1 , Σ v 2 0 . Concerning S ∞ v 2 as a function of S ∞ v 1 , we have S ∞ v 1 ∈ min Σ v 1 -1 , Σ v 1 0 , max Σ v 1 -1 , Σ v 1 0 , with Σ v 1 0,-1 := - η v 1 τ v 1 W 0,-1 (A v 1 ) ,
where

A v 1 = R e,v 1 (R e,v 2 ) τv 1 ηv 2 τv 2 ηv 1 exp τ v 1 η v 2 τ v 2 η v 1 (1 -R 0,v 2 ) . (4.15)
Thus,

S ∞ v 1 , S ∞ v 2 ∈ Ω S := min Σ v 1 -1 , Σ v 1 0 , max Σ v 1 -1 , Σ v 1 0 × min Σ v 2 -1 , Σ v 2 0 , max Σ v 2 -1 , Σ v 2 0 .
We present on Figure 4 (left) the functions W 0 and W -1 defining S ∞ v 2 as a function of S ∞ v 1 and the domain Ω S for a given set of the parameters and initial conditions. We refer to Section 5 for details regarding the numerical integration of the model and Section 6 for further numerical results on the case of two vertices.

Actually, we can reduce the domain of validity of (4.12)-(4.13) for S ∞ v 1 and S ∞ v 2 . Indeed, we know that S v k , k = 1, 2, decay with respect to time, so S ∞ v k < S v k . Moreover, the sum S ∞ v 1 + S ∞ v 2 < M 0 . Thus, we have

S ∞ v 1 , S ∞ v 2 ∈ ω S := [min Σ v 1 -1 , Σ v 1 0 , S 0 v 1 ] × [min Σ v 2 -1 , Σ v 2 0 , S 0 v 2 ] ∩ S ∞ v 1 + S ∞ v 2 < M 0 .
The domain ω S is drawn on Figure 4 (right).

Concerning I ∞ v 1 and I ∞ v 2 , we can perform the same analysis. Let

J ∞ v 1 = -R 0,v 1 + τ v 1 η v 2 τ v 2 η v 1 R e,v 2 exp -τv 2 I ∞ v 2 +τ v 2 I ∞ v 2 ,
and

J ∞ v 2 = -R 0,v 2 + τ v 2 η v 1 τ v 1 η v 2 R e,v 1 exp -τv 1 I ∞ v 1 +τ v 1 I ∞ v 1 .
We obtain for k = 1, 2,

I ∞ v k = 1 τ v k W -R e,v k exp J ∞ v k - J ∞ v k τ v k ,
still with W equal to W -1 and W 0 . Let ι 1 -1,0 and ι 2 -1,0 be defined by and

ι v 1 -1,0 = W -1,0 (-A v 1 ) τ v 1 + τv 1 ηv 2 τv 2 ηv 1 (R 0,v 2 -1 -log (R e,v 2 )) τ v 1 , (4.16) 0 ' v1 !1 ' v1 0 1 S v1 0 ' v2 !1 ' v2 0 1 S v2 + S W 0 W !1 0 ' v1 !1 ' v1 0 1 S v1 0 ' v2 !1 ' v2 0 1 S v2 [S 0 v1 ; S 0 v2 ] [S 1 v1 ; S 1 v2 ] ! S W 0 W !1 S 1 v1 + S 1 v2 6 1
ι v 2 -1,0 = W -1,0 (-A v 2 ) τ v 2 + τv 2 ηv 1 τv 1 ηv 2 (R 0,v 1 -1 -log (R e,v 1 )) τ v 2 , (4.17)
with A v 1 and A v 2 given by (4.14) and (4.15). Then,

(I ∞ v 1 , I ∞ v 2 ) ∈ [min ι v 1 -1 , ι v 1 0 , max ι v 1 -1 , ι v 1 0 ] × [min ι v 2 -1 , ι v 2 0 , max ι v 2 -1 , ι v 2 0 ].
We can show that min ι v k -1 , ι v k 0 < 0 for k = 1, 2. So, we can reduce this domain since

I ∞ v k > 0. So, we define the domain ω I (I ∞ v 1 , I ∞ v 2 ) ∈ ω I := [0, max ι v 1 -1 , ι v 1 0 ] × [0, max ι v 2 -1 , ι v 2 0 ].
As a consequence, we have proved the following lemma.

I ∞ v k = 1 τ v k W -R e,v k exp J ∞ v k - J ∞ v k τ v k , with J ∞ v k = -R 0,v k + τ v k η v j τ v j η v k R e,v j exp -τv j I ∞ v j +τ v j I ∞ v j , k = j ∈ {1, 2} , 4 v1 !1 4 v1 0 I v 1 4 v2 !1 4 v2 0 I v 2 ! I [I 1 v1 ; I 1 v2 ] W 0 W !1 ; v1 !1 ; v1 0 R v 1 ; v2 !1 ; v2 0 R v 2 ! R [R 1 v1 ; R 1 v2 ] W 0 W !1 Figure 5: Location of I ∞ k and R ∞ v k , k = 1, 2
, and visualisation of the domains ω I (left), and domain ω R (right). In both cases, where [START_REF] Liu | Predicting the cumulative number of cases for the COVID-19 epidemic in China from early data[END_REF])-(4.17). See Figure 5 for an illustration.

(I ∞ 1 , I ∞ 2 ) ∈ ω I and (R ∞ 1 , R ∞ 2 ) ∈ ω R
R 0,v k := M 0 τ v k /η v k and R e,v k := S 0 v k τ v k /η v k , k = 1, 2. Furthermore, we have the sharp bound (I ∞ v 1 , I ∞ v 2 ) ∈ ω I := [0, max ι v 1 -1 , ι v 1 0 ] × [0, max ι v 2 -1 , ι v 2 0 ], with ι v k -1,0 , k = 1, 2 defined in (4.
Remark 4.5. As the solutions R ∞ v k , k = 1, 2, are simply given by R

∞ v k = η v k I ∞ v k , if we let ρ v k -1,0 := ηι v k -1,0 then we have (R ∞ v 1 , R ∞ v 2 ) ∈ ω R := [0, max ρ v 1 -1 , ρ v 1 0 ] × [0, max ρ v 2 -1 , ρ v 2 0 ].
We represent on Figure 5 the domain ω R .

A semi-implicit numerical scheme which preserves total mass

In this section, we propose a semi-implicit in time numerical scheme based on finite differences in space which has the property to preserve the discrete total mass.

Notations

For each e ∈ E, we denote δx e > 0 the space discretization of each edge, and J e ∈ N the number of points of the corresponding discretization. For each i = 1, • • • , J e , the space grid on each edge is given by x i = (i -1)δx e with e = (J e -1)δx e . And we let J := e∈E J e ∈ N. Let δt > 0 be the time discretization and denote t m = mδt for m ≥ 0.

For a given function u ∈ C 1,2 (R + × G, R + ), its space-time discretization is given by some sequence of vectors

u ∼ (U m ) m≥0 , with U m = (U m 1 , • • • , U m J ) t ∈ R J .
For each e ∈ E, there exists an integer j e ∈ N such that

u e (t m , x i ) ∼ U m je+i , i = 1, • • • , J e , m ≥ 0.
We approximate the laplacian on each edge via finite differences. That is, for each e ∈ E,

∂ 2 x u e (t m , x i ) ∼ U m je+i-1 -2U m je+i + U m je+i+1 δx 2 e , i = 2, • • • , J e -1, m ≥ 1,
where we have only considered the interior points of the discretized domain. Let us now precise how we approximate the laplacian at a given vertex v ∈ V of the graph. So let v ∈ V such that there are δ v edges incident to the vertex. We locally label e ∼ v = (e 1 , • • • , e δv ) all these incident edges. For each v ∈ V, we introduce the map

σ v : {e 1 , • • • , e δv } → {1, • • • , J} such that σ v (e k )
corresponds to the global index of the grid discretization associated to the vertex v on edge e k . Finally, we denote by n(σ v (e k )) the global index of the nearest neighbor on edge e k to the vertex v.

Note that either n(σ v (e k )) = σ v (e k ) -1 or n(σ v (e k )) = σ v (e k ) + 1.
To approximate the laplacian at a given vertex v ∈ V on edge e k , we use the following formula

∂ 2 x u e k (t m , v) ∼ U * ,m σv(e k ) -2U m σv(e k ) + U m n(σv(e k ))
δx 2

e k := Z m v,k , k = 1, • • • , δ v , m ≥ 1.
The unknown U * ,m σv(e k ) can be expressed by discretization of the boundary condition as follows. For each v ∈ V with e ∼ v = (e 1 , • • • , e δv ), we approximate the normal derivative ∂ n u e k (t m , v) as

∂ n u e k (t m , v) ∼ U * ,m σv(e k ) -U m n(σv(e k )) 2δx e k , k = 1, • • • , δ v , m ≥ 1.
Using (2.3), and denoting I m v the time approximation of I v (t m ), we obtain the following expression for U * ,m

σv(e k ) U * ,m σv(e k ) = U m n(σv(e k )) - 2δx e k d e k α v e k U m σv(e k ) + δv l=1 (N v ) kl U m σv(e l ) -λ v e k I m v , k = 1, • • • , δ v , m ≥ 1.
As a consequence, we obtain that for each

k = 1, • • • , δ v and m ≥ 1 Z m v,k = 2U m n(σv(e k )) -2U m σv(e k ) δx 2 e k - 2 d e k δx e k α v e k U m σv(e k ) + δv l=1 (N v ) kl U m σv(e l ) -λ v e k I m v .

The semi-implicit numerical scheme

We introduce the following scheme for each m ≥ 0

                               U m+1 je+i = U m je+i + d e δt U m+1 je+i-1 -2U m+1 je+i + U m+1 je+i+1 δx 2 e , i = 2, • • • , J e -1, e ∈ E, U m+1 σv(e k ) = U m σv(e k ) + d e k δtZ m+1 v,k , k = 1, • • • , δ v , v ∈ V, S m+1 v = S m v -δtτ v S m+1 v I m v , I m+1 v = I m v + δt τ v S m+1 v I m v -η v I m+1 v + δt δv k=1 α v e k U m+1 σv(e k ) -λ v I m+1 v , R m+1 v = R m v + δtη v I m+1 v , (5.1)
initialized with U 0 ∈ R J and some (S 0 v , I 0 v , R 0 v ) v∈V . One can find similar semi-implicit discretization for the SIR part of the model in [START_REF] Sekiguchi | Dynamics of a discretized SIR epidemic model with pulse vaccination and time delay[END_REF].

Well-posedness and positivity. We prove that the numerical scheme defined through (5.1) is well defined and preserves positivity under some condition on δt. Indeed, we first remark that the equation for S m+1 v and I m+1 v in (5.1) can be used to obtain that

S m+1 v = S m v 1 + δtτ v I m v , I m+1 v = I m v + δtτ v I m v (S m v + I m v ) 1 + δt(η v + λ v ) (1 + δtτ v I m v ) + δt 1 + δt(η v + λ v ) δv l=1 α v e l U m+1 σv(e l ) ,
such that Z m+1 v,k can be expressed only in terms of elements of U m+1 as

Z m+1 v,k = 2U m+1 n(σv(e k )) -2U m+1 σv(e k ) δx 2 e k - 2 d e k δx e k α v e k U m+1 σv(e k ) + δv l=1 (N v ) kl U m+1 σv(e l ) - δtλ v e k 1 + δt(η v + λ v ) δv l=1 α v e l U m+1 σv(e l ) + 2λ v e k d e k δx e k I m v + δtτ v I m v (S m v + I m v ) 1 + δt(η v + λ v ) (1 + δtτ v I m v )
.

As a consequence, there exists a matrix A ∈ M J (R) such that

(I J + A) U m+1 = U m + Y m , where Y m ∈ R J is such that Y m j =      2δtλ v e k δx e k I m v + δtτ v I m v (S m v + I m v ) 1 + δt(η v + λ v ) (1 + δtτ v I m v ) , if j = σ v (e k ), 0
, otherwise.

Lemma 5.1. There exists a constant C 0 > 0, which only depends on the parameters of the system, such that if 0 < δt < C 0 then we have

• I J + A is invertible; • if N v is symmetric for each v ∈ V, then given V ∈ R J with V ≥ 0, the unique solution U ∈ R J of (I J + A)U = V also satisfies U ≥ 0.
Proof. Let U ∈ R J = 0 be such that (I J + A)U = 0. Without loss of generality, assume that

U j 0 = max j=1,••• ,J U j > 0.
If there exists e ∈ E such that j 0 = j e + i 0 for some i 0 ∈ {2, • • • , J e -1}, then we have

U j 0 + d e δt δx 2 e (2U j 0 -U j 0 -1 -U j 0 +1 ) = 0,
which is a contradiction by definition of

U j 0 . Next if j 0 is such that there is v ∈ V and k ∈ {1, • • • , δ v } such that j 0 = σ v (e k
), then we have

U σv(e k ) + 2d e k δt δx 2 e k (U σv(e k ) -U n(σv(e k )) ) = - 2δt δx e k α v e k U σv(e k ) + δv l=1 (N v ) kl U σv(e l ) - δtλ v e k 1 + δt(η v + λ v ) δv l=1 α v e l U σv(e l ) .
The left-hand side of the above equality is strictly positive and we claim that the right-hand side is negative. We use the fact that (N v ) kl = -ν v e l ,e k when k = l and (N v ) kk = j =k ν v e k ,e j δv l=1

(N v ) kl U σv(e l ) =   j =k ν v e k ,e j   U σv(e k ) - l =k ν v e l ,e k U σv(e l ) =   j =k ν v e k ,e j - l =k ν v e l ,e k   U σv(e k ) + l =k ν v e l ,e k U σv(e k ) -U σv(e l ) .
As a consequence, we deduce that

U v k := α v e k U σv(e k ) + δv l=1 (N v ) kl U σv(e l ) - δtλ v e k 1 + δt(η v + λ v ) δv l=1 α v e l U σv(e l ) =   α v e k + j =k ν v e k ,e j - l =k ν v e l ,e k - δtλ v e k 1 + δt(η v + λ v ) δv l=1 α v e l   U σv(e k ) + l =k ν v e l ,e k U σv(e k ) -U σv(e l ) + δtλ v e k 1 + δt(η v + λ v ) δv l=1 α v e l U σv(e k ) -U σv(e l ) .
The last two terms are positive by definition of U j 0 = U σv(e k ) = max j=1,••• ,J U j > 0. Now using Hypothesis 2.2, we have that

α v e k + j =k ν v e k ,e j - l =k ν v e l ,e k > 0,
such that the term in bracket is positive provided that

δtλ v e k 1 + δt(η v + λ v ) δv l=1 α v e l < α v e k + j =k ν v e k ,e j - l =k ν v e l ,e k , or equivalently δt   λ v e k δv l=1 α v e l -(η v + λ v )   α v e k + j =k ν v e k ,e j - l =k ν v e l ,e k     < α v e k + j =k ν v e k ,e j - l =k ν v e l ,e k .
As a consequence, we impose that

0 < δt < min v∈V min k=1,••• ,δv α v e k + j =k ν v e k ,e j -l =k ν v e l ,e k λ v e k δv l=1 α v e l -(η v + λ v ) α v e k + j =k ν v e k ,e j -l =k ν v e l ,e k + ,
where it is understood that when the positive part is zero there is no condition on δt. And we have reached a contradiction since

0 < U σv(e k ) + 2d e k δt δx 2 e k (U σv(e k ) -U n(σv(e k )) ) = - 2δt δx e k U v k < 0.
This shows that I J + A is invertible.

Next let U ∈ R J be the unique solution of (I J + A)U = V with V ≥ 0. We denote by U -∈ R J the vector with components given by

U - j = min(0, U j ), j = 1, • • • , J.
Our aim is to evaluate (I J + A)U, U - J where •, • J is the following scalar product on R J :

U, V J := e∈E Je-1 i=2 U je+i V je+i + 1 2 v∈V δv k=1 U σv(e k ) V σv(e k ) .
We divide (I J + A)U, U - J into three parts:

(I J + A)U, U - J = Q 1 + Q 2 + Q 3 ,
where

Q 1 := e∈E Je-1 i=2 U je+i + d e δt δx 2 e (2U je+i -U je+i-1 -U je+i+1 ) U - je+i , Q 2 := 1 2 v∈V δv k=1 U σv(e k ) + 2d e k δt δx 2 e k U σv(e k ) -U n(σv(e k )) U - σv(e k ) , Q 3 := δt v∈V δv k=1 1 δx e k α v e k U σv(e k ) + δv l=1 (N v ) kl U σv(e l ) - δtλ v e k 1 + δt(η v + λ v ) δv l=1 α v e l U σv(e l ) U - σv(e k ) .
The first and second terms are handled as follows

Q 1 + Q 2 = U, U - J + e∈E d e δt δx 2 e Je-1 i=1 (U je+i+1 -U je+i ) U - je+i+1 -U - je+i ≥ 0.
For the third term Q 3 , if we further assume that N v is symmetric, then the matrix

K v = A v + N v
is symmetric positive definite, and thus for each v ∈ V there exists some

β v > 0 such that δv k=1 1 δx e k α v e k U σv(e k ) + δv l=1 (N v ) kl U σv(e l ) U - σv(e k ) ≥ β v δv k=1 1 δx e k U σv(e k ) U - σv(e k ) ,
while there exists ω v > 0 such that

δv k=1 λ e k δx e k U - σv(e k ) δv l=1 α v e l U σv(e l ) ≤ ω v δv k=1 1 δx e k U σv(e k ) U - σv(e k ) .
And thus, we get an estimate for Q 3 of the form

Q 3 ≥ δt v∈V β v - δtω v 1 + δt(η v + λ v ) δv k=1 1 δx e k U σv(e k ) U - σv(e k ) ,
which is positive provided that δt is small enough. As a consequence, we have proved that

0 ≤ (I J + A)U, U - J = V, U - J ≤ 0,
which implies that U -= 0 and thus U ≥ 0.

The previous lemma demonstrates the well-posedness of our numerical scheme (5.1). It also ensures that if we start with positive initial conditions U 0 ≥ 0 and S 0 v > 0, I 0 v ≥ 0 with v∈V I 0 v > 0 and R 0 v ≥ 0, then for all m ≥ 1 we also have that U m ≥ 0, S m v > 0, I m v ≥ 0 and R m v ≥ 0, provided δt > 0 is small enough and N v is symmetric for each v ∈ V.

Preservation of total discrete mass. For any U ∈ R J , we define the following quantity

trap J (U ) := e∈E δx e Je-1 i=2 U je+i + 1 2 v∈V δx e k δv k=1 U σv(e k ) .
The expression trap J (U ) is simply the trapezoidal rule applied to the elements of U adapted to our graph G. From (5.1), we get that

trap J (U m+1 ) = trap J (U m ) + δt e∈E d e δx e Je-1 i=2 U m+1 je+i-1 -2U m+1 je+i + U m+1 je+i+1 + δt 2 v∈V δx e k d e k δv k=1 Z m+1 v,k
.

Upon denoting Z 1,m v,k the following quantity

Z 1,m v,k := 2U m n(σv(e k )) -2U m σv(e k )
δx 2 e k we get that 1 2

v∈V

δx e k d e k δv k=1 Z 1,m+1 v,k = v∈V d e k δx e k δv k=1 U m+1 n(σv(e k )) -U m+1 σv(e k ) .
Next, we observe that

e∈E d e δx e Je-1 i=2 U m+1 je+i-1 -2U m+1 je+i + U m+1 je+i+1 + v∈V d e k δx e k δv k=1 U m+1 n(σv(e k )) -U m+1 σv(e k ) = 0,
where the cancellation comes from the specific structure of the discretized laplacian through finite differences. As a consequence, we have that

trap J (U m+1 ) = trap J (U m ) -δt v∈V δv k=1 α v e k U m+1 σv(e k ) + δv l=1 (N v ) kl U m+1 σv(e l ) -λ v e k I m+1 v .
We also have that

δv k=1 δv l=1 (N v ) kl U m+1 σv(e l ) = δv l=1 δv k=1 (N v ) kl U m+1 σv(e l ) = 0,
as the sum over the lines of N v vanishes. And thus we get

trap J (U m+1 ) = trap J (U m ) -δt v∈V δv k=1 α v e k U m+1 σv(e k ) -λ v I m+1 v .
On the other hand, from (5.1) we also have

v∈V S m+1 v + I m+1 v + R m+1 v = v∈V (S m v + I m v + R m v ) + δt v∈V δv k=1 α v e k U m+1 σv(e k ) -λ v I m+1 v .
As a conclusion, we have proved the following result.

Lemma 5.2. Let (U m , S m v , I m v , R m v
) a solution of (5.1), then we have for each m ≥ 0

trap J (U m+1 ) + v∈V S m+1 v + I m+1 v + R m+1 v = trap J (U m ) + v∈V (S m v + I m v + R m v ) .
This is the discrete conter part of conservation of mass for the continuous model. Now, combining Lemma 5.1-5.2, we have proved the following theorem.

Theorem 3. There exists a constant C 0 > 0, which only depends on the parameters of the system, such that if 0 < δt < C 0 , then the numerical scheme (5.1) defines a unique sequence

(U m , S m v , I m v , R m v ) m≥0 .
If we further assume that N v is symmetric for each v ∈ V, then the numerical scheme (5.1) preserves the positivity of the initial condition. Finally, for each solution of (5.1), the total discrete mass is preserved, namely for each m ≥ 0, we have

trap J (U m+1 ) + v∈V S m+1 v + I m+1 v + R m+1 v = trap J (U m ) + v∈V (S m v + I m v + R m v ) .
6 Numerical results for a selection of graphs

In the present section, we illustrate our theoretical results with a collection of numerical simulations for various types of graphs. Throughout this section the time discretization is set to δt = 0.01 while the space discretization to δ xe = 0.01 for each e ∈ E.

6.1 Case of 2 vertices and 1 edge

v 1 v 2 • •
We first consider the case where c V = 2 and c E = 1, where c E denotes the cardinal of E. In this setting, we recall that our model reads as follows

∂ t u(t, x) = d∂ 2 x u(t, x), t > 0, x ∈ (0, ),
with boundary conditions

-d∂ x u(t, 0) + α 1 u(t, 0) = λ 1 I 1 (t), d∂ x u(t, ) + α 2 u(t, ) = λ 2 I 2 (t), t > 0,
where (S i (t),

I i (t), R i (t)), for i = 1, 2, solution of        S i (t) = -τ i S i (t)I i (t), I i (t) = τ i S i (t)I i (t) -η i I v (t) + α i u(t, v i ) -λ i I i (t), R i (t) = η i I i (t), t > 0,
where v 1 = 0 and v 2 = . This system is complemented by some initial condition (u 0 , S 0 i , I 0 i , R 0 i ) with S i > 0, I 0 1 + I 0 2 > 0, R 0 i = 0 and u 0 ≥ 0 such that the boundary condition is satisfied initially. Finally, we normalize the total mass as follows

M 0 = 0 u 0 (x)dx + 2 i=1 S 0 i + I 0 i = 1.
For the numerical simulations, we have fixed initial conditions to be of the form

u 0 (x) = λ 1 I 0 α 1 exp - α 2 x 2 2d , x ∈ [0, ],
with

(S 0 1 , I 0 1 , S 0 2 , I 0 2 ) = S 0 -I 0 - 0 u 0 (x)dx, I 0 , 1 -S 0 , 0 ,
where S 0 and I 0 may vary. In Figures 6-7-8, S 0 and I 0 are fixed to (S 0 , I 0 ) = (1/2, 10 -6 ), while in Figure 10, S 0 is allowed to vary and I 0 is fixed to I 0 = 10 -6 . ). We note that as d becomes smaller ∆T rapidly increases as increases.

at the second vertex decreases as λ 2 increases while at the first vertex it varies less significantly. As α 1 increases, the final total population of infected individuals at the first vertex increases while it decreases at the second vertex. This time the relative distance ∆T := T 2 -T 1 between time of maximal infection is monotonically increasing with α 1 . We get the opposite monotonicity properties as α 2 is varied.

In Figure 8, we investigate the joint effect of the diffusion coefficient d and the length of the edge on the dynamics of the epidemic at the vertices. Here, we focus on the delay between time of maximal infection T j in each infected population I j (t). As expected, when the diffusion coefficient is really small while the length is being kept at order one, ∆T takes large value: ∆T ∼ 10 4 when d = 10 -3 and = 1. Biologically, this means that when the diffusion coefficient is really small it takes more time for infected individuals from vertex one to reach the second vertex and start an epidemic. We also note that at fixed , ∆T monotonically decreases as d increases, while at fixed d, ∆T monotonically increases as increases.

In Figures 910, we vary respectively the initial population of susceptible individuals S 0 and infected individuals I 0 . We visualize the final total populations of infected individuals and corresponding final population of susceptible individuals on the parameterized surfaces (

I ∞ 1 , I ∞ 2 , S 0 ) and (S ∞ 1 , S ∞ 2 , S 0 ), respectively (I ∞ 1 , I ∞ 2 , I 0 ) and (S ∞ 1 , S ∞ 2 , I 0 )
, where the level sets of the parameterized surface are given by the conservation of total mass (4.6). We note that (I ∞ 1 , I ∞ 2 ) and (S ∞ 1 , S ∞ 2 ) are almost independent of I 0 when I 0 ≤ 10 -3 with sensible variations only occurring for larger values of I 0 . On the other hand, we observe that as S 0 is increased the final total population of infected individuals increases at the first vertex while it decreases at the second one. The dependence of (S ∞ 1 , S ∞ 2 ) as a function of S 0 is more subtile and is presented in Figure 11. In the same figure, we also show the location of max I j (t) and its amplitude. We observe a strong nonlinear dependence with respect to S 0 . As S 0 increases, we first see that the time at which I 1 (t) is maximal increases and then decreases, while max I 1 (t) is monotonically increasing. The converse is observed 6.2 Case of 3 vertices and 3 edges

v 1 v 2 v 3 • • • A B C
Next, we consider the case of 3 vertices and 3 edges arranged in a triangular configuration. For the numerical simulations presented in Figure 12, we have assumed full symmetry in the parameters that is

( e , d e ) = ( , d), e ∈ E, (τ v , η v ) = (τ, η), v ∈ V, (α v e , λ v e ) = (α, λ), (e, v) ∈ E × V, ν v e,e = ν, (e, e , v) ∈ E × E × V.
Regarding the initial condition, we have chosen

(S 0 1 , I 0 1 , S 0 2 , I 0 2 , S 0 3 , I 0 3 ) = S 0 -I 0 , I 0 , S 0 , 0, S 0 , 0 v ∈ V,
for a given (S 0 , I 0 ), while for each e ∈ E we have set u 0 e (x) = 0 on Ω e . Note that, we have initially a boundary layer as our initial condition does not satisfy (2.3) for small times. We remark that the final total populations of infected individuals and corresponding final population of susceptible individuals belong to a surface as provided by (4.6)-(4.7) from Theorem 2.

In Figure 13, we tested a different configuration. Upon labeling by A the edge between vertices v 1 and v 2 , B the edge between vertices v 2 and v 3 and C the edge between vertices v 1 and v 3 , we have set the parameters to

α 1 A = α 2 B = α 3 C = 0, and α 2 A = α 3 B = α 2 C = 1/10, while λ 1 A = λ 2 B = λ 3 C = 1/20, and λ 2 A = λ 3 B = λ 2 C = 0, and ν 1 A,C = ν 2 B,A = ν 3 C,B = 0, and ν 1 C,A = ν 2 A,B = ν 3 B,C = 1/30
. The length of each edge is fixed e = = 1 and (τ v , η v ) = (1, 1/7) at each vertex v ∈ V. Finally, we have set different coefficients on each edge, namely d A = 1, d B = 10 -2 and d C = 10 -3 . Initially, we assume that infected individuals are only present at vertex v 1 and each vertex has the same number of susceptible individuals fixed to 1/3. Finally, for each e ∈ E we have set u 0 e (x) = 0 on Ω e . We see in Figure 13 that such a configuration can generate a second wave of infection at the first and second vertices showing that transient dynamics can be complex with multiple bumps of infection.

Figure 12: Final total populations of infected individuals (left) and corresponding final population of susceptible individuals (right) as ν is varied from 0.05 to 0.95. The dark blue surfaces are the location of Figure 13: Time plot of infected populations in the case of 3 vertices and 3 edges in a triangular configuration between times [0, 500] (left) and a zoom for times between [150, 400] (right). We observe a second wave of infection at the first vertex resulting from incoming infected individuals that have successively passed through the two other vertices. This second wave is also present at the second vertex with a slight increase of I 2 (t) after the second wave has reached the first vertex. Parameters values are set in the text. Next, we consider a star-shape graph with 4 vertices and 3 edges where one vertex is connected to the three others. In this configuration, we assume that our parameters may vary with respect to time, modeling locked down strategies for example [START_REF] Griette | Unreported cases for Age Dependent COVID-19 Outbreak in Japan[END_REF][START_REF] Liu | Predicting the cumulative number of cases for the COVID-19 epidemic in China from early data[END_REF]. More precisely, we will assume that there exists T lock and µ lock such that the transmission rates can be written as

(I ∞ 1 , I ∞ 2 , I ∞ 3 ) respectively (S ∞ 1 , S ∞ 2 , S ∞

Case of 4 vertices and 3 edges

τ v (t) = τ , t ∈ [0, T lock ], τ exp(-µ lock (t-T lock ))+τ lock 1+exp(-µ lock (t-T lock ))
, t > T lock , for each v ∈ V and for a given 0 < τ lock < τ . We will assume that the four vertices are at equal distance such that e = for each e ∈ E and that the coefficient diffusion are equal on each edge, d e = d, e ∈ E. We further assume that at the central vertex v 2 exchanges are no longer allowed. That is, we impose that Finally, we set η v = η for all v ∈ V. Regarding the initial condition, we work with (S 0 1 , I 0 1 , S 0 2 , I 0 2 , S 0 3 , I 0 3 , S 0 4 , I 0 4 ) = (S 0 -I 0 , I 0 , S 0 , 0, S 0 -, 0, S 0 + , 0) , for given (S 0 , I 0 , ), while for each e ∈ E we have set u 0 e (x) = 0 on Ω e . In Figure 14, we report the location of the time of maximal infection T max for each vertex together with the corresponding amplitude I j max = max t≥0 I j (t) as a function of τ lock . We observe that below a critical value of τ lock , the time of maximal infection always occurs at t = T max = T lock traducing the fact that the locked down strategy has no effect on the dynamics of the epidemic. At each vertex, we observe the same pattern: as τ lock is decreased the corresponding I j max is decreasing while T max is increasing up to some value of τ lock where we observe a sudden turning point (see the right panel of Figure 14). We observe that τ tp lock,v k , the value of the turning point, is well approximated (actually always bounded by below) by the value at which the effective reproduction number of each vertex is equal to 1. Indeed we have R We also point out that when τ lock is below the turning point τ tp lock,v k , the corresponding value of I j max is below 10 -3 . On the other hand, in Figure 15, we present similar results but this time τ lock is fixed and T lock varies. Above some critical value of T lock , I j max saturates to a fixed value independent of T lock traducing the fact that the locked down strategy has no effect on the dynamics of the epidemic if it occurs to late in time. Depending on the initial configuration of susceptible populations at each vertex, we observe intricate nonlinear relationships on the location of the time of maximal infection T max .

α 2 e = α , t ∈ [0, T lock ], α exp(-µ lock (t -T lock )) , t > T lock , e ∈
6.4 Case of N + 1 vertices and N edges

v 1 v 2 v 3 v 4 v N v N +1 • • • • • •
In our final example, we have considered a network of N + 1 vertices and N edges arranged in a lattice, in the sense that vertex v j is only connected to vertices v j-1 and v j+1 via two different edges.

Figure 16 shows the time evolution of the infected population I v j (t) and susceptible populations S v j (t) at each vertex for several different initial conditions when the length and diffusion coefficient of each edge are equal. In the first case (top panel), we assume that I 0 v 1 > 0 while I 0 v j = 0 for all other vertices, and observe a propagation of burst of activity among infected and susceptible populations. In the second case (middle panel), we assume that I 0 v N/2 > 0 while I 0 v j = 0 for all other vertices, and we see the propagation of two bursts of activity among infected and susceptible populations going leftwards and rightwards. In the last case (bottom panel), we assume that I 0 v 1 = I 0 v N +1 > 0 while I 0 v j = 0 for all other vertices, and we note the propagation of two waves activity which collide at the middle vertex v N/2 . For very small values of the diffusion coefficient d, this burst of epidemic activity seems to travel coherently and forms a coherent traveling wave, as can be seen in Figure 17 where we represent the location of max t>0 I v j (t) at each vertex. Such a traveling wave of epidemic activity share similarities with traveling waves in excitable media such as the propagation of electrical activity along a nerve cell [START_REF] Hodgkin | A quantitative description of membrane current and its application to conduction and excitation in nerve[END_REF][START_REF] Hupkes | Traveling pulse solutions for the discrete FitzHugh-Nagumo system[END_REF] or calcium waves [START_REF] Sneyd | Tutorials in Mathematical Biosciences II[END_REF]. When d = 10 -3 , they are all aligned on the same line, where for smaller values d ∈ 10 -1 , 10 -2 the location is a nonlinear curve. We also demonstrate that larger diffusion coefficient leads to a faster propagation of epidemic burst across vertices. Finally, we also remark that if I d max,1 denotes the maximum as a function of d at the first vertex, we have I Regarding the initial condition on the edge, we have set u 0 e (x) = 0 on Ω e for each e ∈ E.

Figure 1 :

 1 Figure 1: Map of France with an illustration of a connected graph connecting major cities.

  u e : Ω e -→ R. For future references, we define BC(G, R) the space of bounded continuous functions on G as BC(G, R) := e∈E BC(Ω e , R), and similarly BC k (G, R) with k ≥ 1. We define the L ∞ norm on G for u ∈ BC(G, R) as u ∞ := max e∈E sup x∈Ωe |u e (x)|.

  whose structure will be specified below. Formally, (2.3) traduces the balance of fluxes of infected individuals at the vertex v, and we will demonstrate this heuristic rigorously by showing in the forthcoming Subsection 2.4 the conservation of total population. together with (K v ) e,e := α v e + e =e ν v e,e ∈ (0, 1),

Figure 4 :

 4 Figure 4: Location of S ∞ v1 and S ∞ v2 together with the visualisation of the domains Ω S (left) and ω S (right). The final configuration of susceptible individuals (S ∞ v1 , S ∞ v2 ) lies on the closed curve parametrized by the two branches of the Lambert W function (blue and red curve). We note that S ∞ v1 , S ∞ v2 ∈ ω S as indicated by the red star on the right figure. Values of the parameters are d = 10 -3 , λ 1 = λ 2 = 6/10, α 1 = α 2 = 1/8, τ v1 = 1, τ v2 = 9/10, η v1 = 2/5, η v2 = 2/6, and initial conditions are set to: I 0 1 = I 0 2 = 10 -6 , S 1 0 = 3/4 -I 0 1 , S 2 0 = 1/4 -I 0 1 and u 0 (x) = 0. The mass M 0 is therefore equal to 1.

Lemma 4. 4 .

 4 Case of two vertices. Assume that |V| = 2 and |E| = 1. The final total population of infected individuals at each vertex I ∞ v k , k = 1, 2 can be expressed as

  Figure 5: Location of I ∞ k and R ∞ v k , k = 1, 2, and visualisation of the domains ω I (left), and domain ω R (right). In both cases, (I ∞ 1 , I ∞ 2 ) ∈ ω I and (R ∞ 1 , R ∞ 2 ) ∈ ω R are represented by a red star. Values of the parameters and initial conditions are similar to Figure 4.

Figure 7 :

 7 Figure7: Final total populations of infected individuals (left) and corresponding final population of susceptible individuals (middle) as one parameter is varied from 0.05 to 0.95 while all parameters are fixed. The blue curve is the location of(I ∞ 1 , I ∞ 2 ) respectively (S ∞ 1 , S ∞ 2 )while the dark red circles indicate the numerically computed values. Right: relative distance ∆T := T 2 -T 1 between time of maximal infection T j in each population, indicated by dark red circles, as the parameter is varied from 0.05 to 0.95. Varying parameters: λ 1 (top panel), λ 2 (second panel), α 1 (third panel) and α 2 (bottom panel).

Figure 8 :

 8 Figure 8: Log-plot of the relative distance ∆T = T 2 -T 1 between time of maximal infection T j in each population I j (t) as the diffusion coefficient d and the length of the edge are varied while all other parameters are fixed to λ 1 = λ 2 = 1/10, α 1 = α 2 = 1/4, and τ 1 = τ 2 = 1 with η 1 = η 2 = 1/3. For the initial condition we have (S 0 , I 0 ) = (1/2, 10 -6). We note that as d becomes smaller ∆T rapidly increases as increases.

Figure 9 :

 9 Figure 9: Final total populations of infected individuals (left) and corresponding final population of susceptible individuals (right) as the initial population of susceptible individuals I 0 is varied from 10 -7 to 10 -1 in logscale while S 0 = 1/2 is fixed. The dark blue curves are the location of (I ∞ 1 , I ∞ 2 ) respectively (S ∞ 1 , S ∞ 2 ) for each value of S 0 , while the dark red circles indicate the numerically computed values. Each dark blue curve is a level set of the parameterized surface given by the conservation of total mass (4.6). All other parameters are fixed to d = = 1, λ 1 = λ 2 = 1/10, α 1 = α 2 = 1/4, and τ 1 = τ 2 = 1 with η 1 = η 2 = 1/3.

Figure 10 :

 10 Figure 10: Final total populations of infected individuals (left) and corresponding final population of susceptible individuals (right) as the initial population of susceptible individuals S 0 is varied from 0.05 to 0.95 while I 0 = 10 -6 is fixed. The dark blue curves are the location of (I ∞ 1 , I ∞ 2 ) respectively (S ∞ 1 , S ∞ 2 ) for each value of S 0 , while the dark red circles indicate the numerically computed values. Each dark blue curve is a level set of the parameterized surface given by the conservation of total mass (4.6). All other parameters are fixed to d = = 1, λ 1 = λ 2 = 1/10, α 1 = α 2 = 1/4, and τ 1 = τ 2 = 1 with η 1 = η 2 = 1/3.

Figure 11 :

 11 Figure 11: Locations of S ∞ 1,2 (top) and max I 1,2 (t) (bottom left and right) as functions of S 0 1,2 . Values of all other parameters are similar to Figure 10. The initial condition is of the form (S 0 , 1 -S 0 ) with S 0 ∈ [0, 1] and to each initial configuration is associated a color code from blue to red. The curve in the top right panel is a projection on the (S ∞ 1 , S ∞ 2 )-plane of the parametrized a curve from Figure 10, right panel.

3 )

 3 Figure12: Final total populations of infected individuals (left) and corresponding final population of susceptible individuals (right) as ν is varied from 0.05 to 0.95. The dark blue surfaces are the location of(I ∞ 1 , I ∞ 2 , I ∞ 3 ) respectively (S ∞ 1 , S ∞ 2 , S∞3 ) for each value of S 0 , while the dark red circles indicate the numerically computed values. Parameters were set to = d = 1, (τ, η) = (1, 1/6), and (α, λ) = (1/8, 1/10), while the initial condition is (S 0 , I 0 ) = (1, 10 -6 ).

Figure 14 :

 14 Figure 14: Location of the time of maximal infection T max for each vertex together with the corresponding amplitude I j max as a function of τ lock (left) with its projection in the (T max , τ lock )-plane (middle) and a zoom near the turning points (right). Other parameters are set to = 1, d = 0.1, η = 1/8, (α, λ, ν) = (1/8, 1/20, 1/20), T lock = 50 and µ lock = 100 with (S 0 , I 0 , ) = (1/4, 10 -6 , 10 -2 ).

  E, while α j e = α for j = 2 and e ∈ E, together withλ 2 e = λ , t ∈ [0, T lock ], λ exp(-µ lock (t -T lock )) , t > T lock , e ∈ E,while λ j e = λ for j = 2 and e ∈ E, and alsoν 2 e,e = ν , t ∈ [0, T lock ], ν exp(-µ lock (t -T lock )) , t > T lock ,(e, e ) ∈ E × E.

Figure 15 :

 15 Figure 15: Location of the time of maximal infection T max for each vertex together with the corresponding amplitude I j max as a function of T lock for two configurations of initial susceptible populations at vertices v 3 and v 4 , with = 10 -1 (left) and = 10 -2 (right). Other parameters are set to = 1, d = 0.1, η = 1/8, (α, λ, ν) = (1/8, 1/20, 1/20), T lock = 50 and η lock = 0.6 with (S 0 , I 0 ) = (1/4, 10 -6 ).

  e,v k = 1 if and only if τ c v k =

  Figure16shows the time evolution of the infected population I v j (t) and susceptible populations S v j (t) at each vertex for several different initial conditions when the length and diffusion coefficient of each edge are equal. In the first case (top panel), we assume that I 0 v 1 > 0 while I 0 v j = 0 for all other vertices, and observe a propagation of burst of activity among infected and susceptible populations. In the second case (middle panel), we assume that I 0 v N/2 > 0 while I 0 v j = 0 for all other vertices, and we see the propagation of two bursts of activity among infected and susceptible populations going leftwards and rightwards. In the last case (bottom panel), we assume that I 0 v 1 = I 0 v N +1 > 0 while I 0 v j = 0 for all other vertices, and we note the propagation of two waves activity which collide at the middle vertex v N/2 . For very small values of the diffusion coefficient d, this burst of epidemic activity seems to travel coherently and forms a coherent traveling wave, as can be seen in Figure17where we represent the location of max t>0 I v j (t) at each vertex. Such a traveling wave of epidemic activity share similarities with traveling waves in excitable media such as the propagation of electrical activity along a nerve cell[START_REF] Hodgkin | A quantitative description of membrane current and its application to conduction and excitation in nerve[END_REF][START_REF] Hupkes | Traveling pulse solutions for the discrete FitzHugh-Nagumo system[END_REF] or calcium waves[START_REF] Sneyd | Tutorials in Mathematical Biosciences II[END_REF]. When d = 10 -3 , they are all aligned on the same line, where for smaller values d ∈ 10 -1 , 10 -2 the location is a nonlinear curve. We also demonstrate that larger diffusion coefficient leads to a faster propagation of epidemic burst across vertices. Finally, we also remark that if I d max,1 denotes the maximum as a function of d at the first vertex, we have I d 1 max,1 ≤ I d 2 max,1 for d 1 ≤ d 2 while for larger vertices j ≥ 6 we have the reverse ordering I d 1 max,j ≥ I d 2 max,j for d 1 ≤ d 2 . For the numerical simulations presented in Figures 16-17, we have assumed full symmetry in the parameters that is

Figure 16 :

 16 Figure 16: Time evolution of the infected (left) population I j (t) and susceptible (right) populations S j (t) at each vertex for several different initial conditions and N = 24. Top: infected individuals are initially present only at vertex. Middle: infected individuals are initially present only at the middle vertex. Bottom: infected individuals are initially present only at the first and last vertices. We observe a traveling wave of infectious activity propagating though the vertices. Parameters were set to = 1, d = 10 -3 , (τ, η) = (1, 1/75), and (α, λ, ν) = (1/8, 1/10, 1/20), while the initial condition is (S 0 , I 0 ) = (1/25, 10 -6 ).

  the maximum principle implies that u n e ≥ 0 for each e ∈ E. Assume by contradiction that e * ∈ E is the component which reaches a negative minimum, namely u n e * (t * , x * ) = -δ < 0 with u n e * (t, x) > -δ for t < t * and x ∈ Ω e * and for each e = e * we have u n e (t, x) > -δ for t ≤ t * and x ∈ Ω e . We know that x * ∈ ∂Ω e * and let denote v * = x * ∈ V the vertex where this occurs. The Hopf lemma implies that ∂ n u n e * (t * , v * ) < 0. Inspecting the boundary condition (3.2) at v * , we obtain that

	d e * ∂ n u n e * (t * , v * ) + α v * e * u n e * (t * , v * ) +	e∼e *	ν v * e * ,e u n e * (t e∼e *	ν v * e,e * u n e (t * , v * ) = λ v * e * I n v * (t),
	which writes
	0 > d e * ∂ n u n e * (t e∼e *	ν v * e,e * -α v * e * -	e∼e *	ν v * e * ,e -	e∼e *	ν v * e,e * (δ + u n e (t

* , v * ) -* , v * ) + δ * , v * )) = λ v * e * I n v * (t) ≥ 0,

and leads to a contradiction. Here we have used the fact that e∼e * ν v * e,e * ≤ α v * e * + e∼e * ν v * e * ,e , from Hypothesis 2.2 on the matrices (K v ) v∈V .

Acknowledgment

This works was partially supported by Labex CIMI under grant agreement ANR-11-LABX-0040.

Figure 6: Profiles of the solutions (S i (t), I i (t)) together with the total population on the edge 0 u(t, x)dx and the total mass of the system M (t) as the parameter λ 1 is varied from 0.05 to 0.95. All other parameters are fixed and set to d = = 1, λ 2 = 1/10, α 1 = α 2 = 1/4, and τ 1 = τ 2 = 1 with η 1 = η 2 = 1/3. For the initial condition we have (S 0 , I 0 ) = (1/2, 10 -6 ).

In Figure 6, we report the profiles of the solutions (S i (t), I i (t)) together with the total population on the edge 0 u(t, x)dx and the total mass of the system M (t) as the parameter λ 1 is varied from 0.05 to 0.95, while all other parameters are being kept fixed. We observe that the dynamics of the epidemic at the second vertex is almost independent of the parameter λ 1 while it has a significant impact on the dynamics at the first vertex. Indeed, as λ 1 is increased, the maximum of infected individuals max t≥0 I 1 (t) is decreased. In the last panel of the figure, we also illustrate the conservation of total population where the fluctuations around M 0 = 1 is of order 10 -12 . In the top panel of Figure 7, we present the final total populations of infected individuals and corresponding final population of susceptible individuals as λ 1 is varied. The blue curve is the location of (I ∞ 1 , I ∞ 2 ) respectively (S ∞ 1 , S ∞ 2 ) while the dark red circles indicate the numerically computed values. We recover the fact that λ 1 has a more significant impact on the final total populations at the first vertex than it has at the second vertex. The get a better understanding of the intricate dynamics between the epidemic at the two vertices, we also present the relative distance ∆T := T 2 -T 1 between time of maximal infection T j in each population as λ 1 is varied. We observe that ∆T is not monotone in λ 1 , as it first decreases and then increases. But we also note that ∆T < 0 for λ 1 ≥ 0.1 traducing the fact that the pick of the epidemic occurs at the second vertex before it does at the first vertex, although initially I 0 2 = 0. This illustrates the effect of the diffusion of infected individuals along the edge.

Similarly, in Figure 7, we report the final total populations of infected individuals and corresponding final population of susceptible individuals as λ 2 (second panel), α 1 (third panel) and α 2 (bottom panel) are varied from 0.05 to 0.95. As expected, the final total population of infected individuals


MD Isabelle Mosnier
email: isabelle.mosnier@aphp.fr

MD Olivier Sterkers

PhD Yann Nguyen

PhD Ghizlene Lahlou

Benefits in noise from sound processor upgrade in thirty-three cochlear implant users for more than 20 years

Keywords: long-term benefit, hearing in noise, upgrading, speech processor

Purpose: Some oldest patients rehabilitated with a cochlear implant more than 20 years ago could still be upgraded with new generations of speech processor (SP). The aim of this study was to show the benefit of a recent generation of SP in this population.

Methods: A monocentric prospective study was designed to evaluate the performance of 33 ancient CI22M users implanted between 1989 and 1997 and upgraded with the late compatible sound processor CP900. Performance were evaluated in quiet and noise with Framatix, an automated adaptative test.

Results: Performance using Framatix significantly improved with the CP900, with a decrease of the median speech perception threshold of 6 dB in quiet (p<0.05) and 5,3 dB in noise (p<0.0005). No subjective benefit using the APHAB questionnaire was observed.

Conclusion:

Upgrading of cochlear implant recipients who were implanted more than 20 years ago with recent compatible and new technological SP provide benefit in speech recognition in noise.

Introduction

Deafness has been established to play a key role on cognitive functions all along the life, as early deafness restoration permits development of language [START_REF] Kral | Neurocognitive factors in sensory restoration of early deafness: a connectome model[END_REF], deafness at mid-age is the main modifiable factor for dementia [START_REF] Livingston | Dementia prevention, intervention, and care[END_REF], and profoundly deaf rehabilitation by cochlear implant (CI) in elderly reduces dramatically cognitive impairment and its evolution to dementia [START_REF] Mosnier | Long-term cognitive prognosis of profoundly deaf older adults after hearing rehabilitation using cochlear implants: cognitive prognosis after hearing rehabilitation[END_REF]. Multichannel CI available since the late 1980s has changed the life of severe to profoundly deaf people [4, who would further gain benefit from new technologies by upgrading the external sound processor (SP). New SP generation has improved dramatically CI performance in noise intelligibility for most patients [START_REF] Mosnier | Benefits from upgrade to the CP810 Sound Processor for Nucleus 24 Cochlear implant recipients[END_REF][START_REF] Mosnier | Benefit of the UltraZoom beamforming technology in noise in cochlear implant users[END_REF][START_REF] Franco-Vidal | The Oticon Medical Neuro Zti cochlear implant and the Neuro 2 sound processor: multicentric evaluation of outcomes in adults and children[END_REF]. However, some oldest CI implanted more than 20 years ago could not be upgraded with new SPs due to the lack of compatibility for some brands of CI. It raised the question of surgical reimplantation of the internal part in this population which would have concerned at most more than 600 recipients implanted before 2000 in France. The aim of this study was to show the benefit of a recent generation of SP in the oldest CI population.

Methods

A monocentric study was prospectively conducted on CI recipients, implanted between 1989 and 1997 with the first generation of CI (CI22M, Cochlear, Sydney, Australia), to test the benefit of the more recent compatible SP (CP 900). Among 63 implanted patients, 17 were either lost of follow-up or deceased, 11 others were not eligible for a reimbursement of a new SP, and 2 had incomplete data. Finally 33 daily CI22M users for 21±2.2 years (range: 19-26) were included. Among them, 16 pre-lingual, 2 peri-lingual, and 15 post-lingual subjects with a profound to total hearing loss were implanted at 5 years (range: 3-16), at 37 and 38 years, and at 36 years (range: 16-60), respectively. Demographic data of the population are described in Table 1. The current SP at the time of the upgrade was a Freedom and an ESPrit 3G SP for 31 and 2 subjects respectively, and all of them used the SPEAK coding strategy.

Characteristics of the three generations of SP were indicated in Table 2. The CP 900 SP was fitted with the same SPEAK strategy and the same MAP characteristics as was used in the original processor. Three programs were available: 2 new automatic signal processing algorithm programs (SCAN) and a non-SCAN program with the same options used as in the original SP (Table 1). CI users were tested with their current revised SP, and at 2 months post-upgrade with the CP900 with their preferred program, in a soundproof room with two lists of monosyllabic words presented at 60 dB SPL in quiet (Lafon list) and for the best performers, with an automated adaptive test in quiet and noise (Framatrix 8), with the current revised SP and the CP900 at 2 months post-upgrade. Subjective perceptions were evaluated using the Abbreviated Profile of Hearing Aid Benefit (APHAB) questionnaire 9.

The results for each test session were compared independently. Scores for words in quiet, the Framatrix test and the APHAB questionnaire were not normally distributed, so a nonparametric Wilcoxon paired test were used. A p value of less than 0.05 was considered to be significant.

Results

Two months after upgrading first generations of internal receiver (CI22M) with the CP900 SP, 31 out of the 33 patients used the SCAN programs all day long (Table 1). No change was observed in quiet for median scores of words (41% vs 43%) and phonemes (72% vs 67%) recognition of the monosyllabic words list. Using adaptive test,speech reception threshold measurement in quiet was feasible before upgrading, in 18 out of 33 CI 22M recipients in quiet (54%) and 13 out of these 18 patients in noise (39%). In this population, the median speech reception threshold was decreased by 6dB (range: -27, +5; p<0.05) with CP900 as compared to previous SP (Fig. 1A). In noise, it decreased with CP900 by 5.3dB (range: -16, -0.4; p<0.0005) as compared to older SP, and even became measurable at +12.4dB in 6 other subjects (Fig. 1B). The benefit in noise was observed in pre-lingual and in post-lingual CI22M recipients. APHAB scores were available before and after upgrading in 30 out of 33 CI 22M recipients and was unchanged with the CP900 SP (Fig. 2).

Discussion

Over the last 15-20 years, the series of SPs from Cochlear Limited has introduced a number of refinements designed to enhance CI performance mainly in noise. The processor CP810, introduced with the 5 system implant was a technological breakthrough for hearing in noisy conditions by an electrical stimulation of the deaf cochlea. Being compatible with previous CI, upgrading Nucleus 24 with CP810 noise program improved by more than 20% performance in 77% of users, thanks to the two adaptive omni-directional microphones and new front-end processing options5. In the Nucleus 6 (CP 900 series), an automatic signal processing algorithm (SCAN) has been added allowing automatic transition between six scenes based on the analysis of environmental signal. Because CI is a highly reliable implanted device (cumulative survival percentage of 92.1% over 29 years for CI22M reported by the company)10, reimplantation would not be necessary to upgrade ancient CI22M devices which have been implanted more than 20 years ago with new technologies if also compatible. Once made possible, upgrading the old CI with the compatible new generation CP900 SP, yielded to significant improvement in noise intelligibility for more than half of CI22M users in a short period of time. The number of CI22M implantees who reached the level of intelligibility in noise using the Framatix test increased from 39% to 58% with previous SPs and CP900 SP, respectively and most patients prefer new automatic algorithms.

No change was observed in quiet regardless of the test used. These results confirm a study including two groups of 15 and 24 CI22M recipients recruited from six North American clinics, showing no change in quiet for CNC words but speech intelligibility improvement in noise with the CP900 compared with the Freedom SP, for AzBio sentences at fixed signal-tonoise ratio [START_REF] Biever | Upgrade to Nucleus 6 in previous generation Cochlear sound processor recipients[END_REF]. This suggests that after such a long duration of cochlear implantation, the frequency mapping should be located in the auditory cortex which may allow a rapid adaptation to new coding strategy and audiological technology. However, the patients implanted more than 20 years ago were different from those who have been implanted more recently and have already benefit from new technologies. Early CI candidates with a prelingual deafnesswere implanted later than nowadays, a critical period for hearing maturation [START_REF] Bruijnzeel | Evaluation of pediatric cochlear implant care throughout Europe: Is European pediatric cochlear implant care performed according to guidelines?[END_REF]. Further post-lingual CI users had a more profound deafness and for longer time than during the last 10 years [START_REF] Holder | Current profile of adults presenting for preoperative cochlear implant evaluation[END_REF]. This may explain that, before upgrade, only 54% of CI recipients in quiet and 39% in noise were able to perform adaptive tests which are more sensitive to show a benefit with the new SP but also more difficult. Moreover, the benefit in noise of upgrading CI22M with CP900 appears to be less striking than upgrading Nucleus 24 with CP810 5 and APHAB questionnaire was probably not sensitive enough to capture a weaker subjective benefit.

Conclusion

More than half of the upgraded old CI recipients with a compatible recent SP performed better the difficult Framatix test in quiet and noise in a short period of time. Although it was thought that reimplantation of the internal receiver should be necessary several times during life-spanning of the majority of profound to total deaf patients, this study evidences that cochlear implantation would be performed only once, especially in children operated on nowadays before 1 year old.

Figure legend

Figure 1 :

 1 Figure 1: Performance using Framatrix in quiet (A, n=18) and in noise (B, n=19) of CI 22 M

Figure 2 :

 2 Figure 2: APHAB scores of CI 22 M recipients upgraded with CP900. The box plots show the

No funding. Written informed consent was obtained from each subject before enrolment into the study (CNIL n° 2040853v0 ). help of Solange Lator for assistance in data collection.