Dynamics of epidemic spreading on connected graphs.

We propose a new model that describes the dynamics of epidemic spreading on connected graphs. Our model consists in a PDE-ODE system where at each vertex of the graph we have a standard SIR model and connections between vertices are given by heat equations on the edges supplemented with Robin like boundary conditions at the vertices modeling exchanges between incident edges and the associated vertex. We describe the main properties of the system, and also derive the final total population of infected individuals. We present a semi-implicit in time numerical scheme based on finite differences in space which preserves the main properties of the continuous model such as the uniqueness and positivity of solutions and the conservation of the total population. We also illustrate our results with a collection of numerical simulations for a selection of connected graphs.

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 (1927) 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 where is a contact rate between susceptible and infected populations, and is the average infectious period; see Hethcote (2000) 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 (Kermack and McKendrick 1927), malaria (Mandal et al. 2011), SARS influenza (Centers for Disease Control and Prevention 2003; Magal et al. 2016) and most recently COVID-19 (New England Journal of Medicine 2020; Magal and Webb 2018; Liu et al. 2020); see also Stolerman et al. (2015) and Magal and Webb (2018) 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 (Van den Driessche and Watmough 2002; Magal et al. 2016, 2018). 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 Fig. 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 (Aronson 1977; Diekmann 1978; Reluga 2004; Berestycki et al. 2020). 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.

Fig. 1

Map of France with an illustration of a connected graph connecting major cities

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 (Berestycki et al. 2020) and models that incorporate networks with more sophisticated interactions dynamics (Britton et al. 2008; Sahneh et al. 2013; Spricer and Britton 2019; Ball and Britton 2020; Bonnasse-Gahot et al. 2018). On a formal level, our proposed model can be thought of as being a one-dimensional version of the planar reaction-diffusion system of Berestycki et al. (2020) 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 David et al. (2020) which studies the spread of airborne diseases where the movement of pathogens in the air is assumed to follow a linear diffusion. Actually, one-dimensional PDE-ODE models have been successfully used in other biological contexts, and we refer to the line of works (Gou and Ward 2016; Gou et al. 2015; Paquin-Lefebvre et al. 2020) on cell-bulk models of cell signaling. Although these models share some similarities, the long time dynamics are different. In Gou and Ward (2016), Gou et al. (2015) and Paquin-Lefebvre et al. (2020), one typically observes oscillatory dynamics caused by the coupling at the boundary via the loss of stability of steady-states through Hopf bifurcations, whereas the long time behavior of the solutions of our model presents the characteristics of SIR compartment models with the convergence to a steady-state of the system. As it will be demonstrated, the selected steady-state is entirely determined by the initial configuration of individuals in our network.

Finally, we highlight that there is an intrinsic difference between the PDE-ODE model that we propose and the compartment models on fixed graphs that are prominent in the literature on the subject, see for example (Britton et al. 2008; Sahneh et al. 2013; Stolerman et al. 2015; Spricer and Britton 2019; Ball and Britton 2020; Bonnasse-Gahot et al. 2018) and references therein. In our setting, connections between vertices are not instantaneous as infected individuals have to diffuse along edges in order to be transported from one vertex to another. As a consequence, our PDE-ODE model intrinsically incorporates very subtle nonlocal in time behaviors that typical compartment models on fixed graph with local interactions do not take into account.

Model formulation and main results

Throughout, we denote by a compact metric graph, i.e. a collection of vertices and edges and further assume that is finite and connected. Each edge is identified with a segment with , where is the finite length of the edge. A real valued function is a collection of one dimensional maps defined for each edge :For future references, we define the space of bounded continuous functions on asand similarly with . We define the norm on for as

A SIR model on compact connected graph

Given a graph , we let , for each , where represents the population of susceptible individuals, the population of infected individuals and the population of removed individuals at vertex and time . We assume that evolves according to a SIR model of the formwhere are the intrinsic parameters of the epidemic which may depend on the vertex v. The contribution in the right-hand side of the equation for the infected population encodes the fact that infected individuals can leave the vertex v to incident edges whereas reflects the contribution of incoming infected individuals from incident edges. Here, denotes the edges incident to the vertex v andsuch that infected individuals leave vertex v to edge e. We have assumed that only the infected population is subject to movement, and we think of being an ambient population whose movement does not affect its distribution. We recover the standard SIR model (1.1) by considering the trivial graph . Throughout the manuscript, we will assume the following standing hypothesis on the coefficients and in (2.1).

Hypothesis 2.1

For each we assume thattogether with

Next, for each , we let and we assume that evolves according toAssuming 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. Possible extensions could be to incorporate nonlocal diffusion or transport terms.

It now remains to model the exchanges of infected individuals at the vertices. For each , we associate an integer which we refer to as its degree (i.e. number of edges incident to the vertex v). We define as the column vector functionwhere we recall that denotes the edges incident to the vertex v, and thus is the corresponding limit value of at . Define also as the column vector functionwhere is the outwardly normal derivative of at the vertex v. Our boundary conditions at the vertex v that link (2.1) and (2.2) are described bywhere is the diagonal matrix and whose structure will be specified below. Formally, (2.3) encodes the balance of fluxes of infected individuals at the vertex v, and we will demonstrate this heuristic rigorously by showing in the forthcoming Sect. 2.4 the conservation of total population.

Assumptions on the connectivity matrices

We now precise the form of the matrix entering in the boundary condition (2.3). Essentially, gathers two contributions. One contribution comes from the exchanges between infected individuals at the vertex with the incoming infected individuals for the incident edges. The second contribution models exchanges between edges. Indeed we allow infected individuals to pass from one edge to another one. More precisely, we have that splits into two partswhere the matrix is the diagonal matrix while the matrix is such that the sum of each column is zero. More precisely, if we label by the edges incident to the vertex v, we have that for all and for In the case , we getsee Fig. 2 for an illustration in that case.

Fig. 2

Schematic illustration of the exchanges at a given vertex v with

Furthermore, for the diagonal term we will use the shorthand notationThe fact that is such that the sum of each column is zero precisely encodes the fact that there is the conservation of infected individuals through exchanges between incident edges at each vertex. And, we remark that it implies that the matrix has a strict column diagonal dominance in the sense that for each because of this specific structure of . From now on we also assume that has a diagonal dominance with respect to its lines. This property will be crucial later on in the proof of existence of solutions. As a consequence, we impose the following running assumptions on the matrices .

Hypothesis 2.2

For each and , we assume thatFurthermore, we impose that for all together withfor each .

Remark 2.3

If the exchanges between the edges are symmetric, that is for each the matrices are symmetric, that isthen 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 , we havesuch that for ,On the other hand, for the ODE system (2.1), we suppose thatWe further assume that (2.3) is satisfied at with obvious notations and . Last, we impose that the initial total population of infected individuals is strictly positive,and that susceptible individuals are initially present at each vertex of the graphThis in turn implies that the total population is initially

Conservation of total population

Assuming that there is a solution to to (2.1)–(2.2)–(2.3), we have that the total mass of the system M(t) defined asis a conserved quantity and thus independent of t.

To see that, we first remark thatwithand is the standard Euclidean inner product on . On the other hand let us defineand assume that u is a classical solution of (2.2), which we will prove in the next section, and computeThe fact thatis a direct consequence on the specific structure of each matrix and the fact that the sum of each column is zero. We therefore conclude that andBiological 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 Sect. 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 , is a well defined quantity: for and are solutions of a system of implicit equations, where stands for the cardinal of , which belong to the parametrized submanifoldwhere 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 in the fully symmetric case where we obtain closed form formula (see Lemma 4.4) 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.5). 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 Sect. 4.

Main result 3: Amass preserving semi-implicit numerical schemeMain result 3: Amass 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 is symmetric for each , then our mass preserving semi-implicit numerical scheme is well-posed and preserves the positivity of the solutions. We refer to Sect. 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 Sect. 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 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 with , and that satisfy the boundary condition (2.3), there exists a unique positive global solution and .

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 with , and , that is for all , is bounded continuous on . 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 , , and are continuous functions of time satisfying , , and together with Hypotheses 2.1–2.2 verified at all times .

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 has the limit of a subsequence of solution of the following problemswithandstarting from and . Note that (3.1)–(3.2)–(3.3) is supplemented by the same initial condition at each step. We proceed along three main steps.

Step #1: solvability of (3.1)–(3.2)–(3.3)

We first show that (3.1)–(3.2)–(3.3) admits a unique solution. It can be done by induction. Assume that at step , we have constructed a solution such that for each is continuous, then we get the existence of a unique solution of (3.1) which is 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 in time. As both and are invertible matrices, we get the existence of a classical solution which then ensures that is continuous.

Step #2: a priori estimates. Let be fixed. We first show by a recursive argument that , , for each and for each . It is trivial at . Let assume that is it true at . We start from (3.1) and a direct integration givesNow owing that for each , the maximum principle implies that for each . Assume by contradiction that is the component which reaches a negative minimum, namely with for and and for each we have for and . We know that and let denote the vertex where this occurs. The Hopf lemma implies that . Inspecting the boundary condition (3.2) at , we obtain thatwhich writesand leads to a contradiction. Here we have used the fact thatfrom Hypothesis 2.2 on the matrices .

Next, from the positivity of solutions, we obtain some uniform bounds. More precisely, we claim that there exists a constant depending only on such thatandFirst, using (3.1) we obtain thatwhich gives thattogether withwhich in turn implies thatfor only depend on the initial condition and the parameters of the system. We now claim that by induction, we have for all ,with for and for some only depending on . As , we get thattogether withfor some constant depending on .

Step #3: existence of a solution. Parabolic Schauder estimates give that the time derivative and the space derivatives up to order 2 of are uniformly Hölder continuous in compact sets. As a consequence, we can use the Arzela–Ascoli theorem to show that converges (up to sequences) toward in . Passing to the limit in (3.1)–(3.2)–(3.3) we get that 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 we have the uniform bounds:andtogether withThe fact that implies thanks to the strong maximum principle that in factwhich in turn gives that for each sinceFinally, we use the conservation of mass which tells us thatsuch that both and are uniformly bounded in time, together with their derivatives. This also implies that there exists a constant , depending only such thatUsing again parabolic regularity, we obtain the solution is global in time and satisfies (2.1)–(2.2)–(2.3) in the classical sense.

Uniqueness

Let assume that and are two classical solutions to (2.2)–(2.3)–(2.1) starting from the same initial datum . We denote where for each and each By linearity, we get that for together withOn the other, one computes that satisfies for each We define the energyand note that by definition. Next, differentiating , we obtainOn the one hand, we haveas is symmetric positive. On the other hand, we computewhere is some large positive constant. Next, we see thatsuch that we obtainNext, if we denote , we computeAs a consequence, we getfor some and we conclude that for all time which then implies that and .

Long-time behavior of the solutions

Throughout this section, we denote by the unique positive bounded classical solution of the Cauchy problem (2.1)–(2.2)–(2.3) as given by Theorem 1 and which further satisfies the conservation of total population, namely

Final total populations: general results

As and is strictly decreasing, it asymptotically converges towards a limit that we denoteFurthermore, as is strictly increasing and uniformly bounded, it asymptotically converges towards a limit that is denotedBut as for each this implies thatwhich in turn proves thatpending that we prove that the limit of as exists. To see that, we use the fact thatto obtain thatWe now argue that the left-hand side of the previous equation is uniformly bounded as , which ensures that has a limit as by positivity of . This in turn implies that the limit of as exists.

If one recall the notation m(t) for the total population on the edges then we haveand it verifiesThe above computations shows that m(t) has a limit as , that we denote and which satisfiesWe shall also keep in mind thatAnd so if we introduce the function , then the above conservation of mass can be written asOn the other, one can compute thatsuch thatNow, as m(t) and each are convergent we deduce that all are also convergent so thatandwhich proves thatAnd the boundary conditions imply thatNext, we define the sequence of functions for each and for each which are uniformly bounded such that one can extract a convergent subsequence. On the one hand we have that and on the other if it is solution ofwith the boundary conditionsThis then shows that , and for each . As there is unicity of the limit, we deduce thatFrom which we also get that and thatWe also get from (4.3), thatFinally, as a consequence of (4.1), the final total populations of infected individuals at each vertices satisfy the following scalar differential equationPassing to the limit as , we getTo summarize, we have proved the following result.

Theorem 2

For each with , and that satisfy the boundary condition (2.3), the long time behavior of the unique corresponding solution is given byandwhere the final total populations of infected individuals at each vertices are solutions of the systemAs a consequence, belongs to the parametrized submanifold given by

Remark 4.1

Equivalently, belongs to the parametrized submanifold given byand belongs to the parametrized submanifold given byThe equations (4.7), (4.8), and (4.9) also read , which is nothing but (4.2) since we have proved that .

Remark 4.2

If we assume that and are independent of and let , and . Then, equations (4.7), (4.8), and (4.9) are respectively equivalent toandThe common right hand side features that is nothing but the traditional basic reproductive number .

Remark 4.3

The above equation (4.4) can be interpreted as the counter-part on graphs to equations that were already derived in a spatially continuous setting (Diekmann 1978; Berestycki et al. 2020). In our case, we can rewrite (4.4) aswhere and . As a consequence, Eq. (4.4) can be thought as a discrete reaction-diffusion equation set on the graph where is a heterogenous diffusion process which takes into account the connectivity of the graph, while encodes the nonlinear reaction terms. As in Diekmann (1978) and Berestycki et al. (2020), the nonlinearity has some concavity properties asWe do not pursue this direction as it would require a thorough study of the heterogeneous diffusion operator which is beyond the scope of the present manuscript.

Final total populations of infected individuals: further properties

The aim of this section is to present further qualitative results on the final total configuration 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 (Diekmann et al. 1990). We also refer to Stolerman et al. (2015) in the case of networks without diffusion on the edges where properties of the basic reproduction number are linked to the geometry and heterogeneity of the network.

Fully symmetric case

We assume that the length of every edge is equal to a reference length . For every , the diffusion coefficient is equal to d. We moreover suppose that for every vertex , , and . We also assume that and are independent of . In the same spirit, and for every and . We also assume for every edges incident to the vertex v. Finally, the components of initial condition on each edges are supposed to be even with respect to the center of the interval . Thanks to all these assumptions, does not depend on the vertex and we set . Let us recall the notation for the cardinal of the set . The parametrized submanifold given by (4.7) becomeswhere . We can transform this relation asLet . We have to solveThe solutions are given in terms of Lambert W function that is the multivalued inverse relation of the function for (Corless et al. 1996). Let us recall how to compute the real solutions of the equation for . Let be the discriminant. If or , the solution is unique and where is the principal branch. If , there are two solutions and , where is another branch. When , there is no solution.

In our symmetric case, the discriminant writesSince , there exist solutions to (4.10) if , which is equivalent toWe 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 and the basic reproductive number respectively given bysee Diekmann et al. (1990), Van den Driessche and Watmough (2002) and Hethcote (2000) for further properties of effective and basic reproductive numbers. If we denote , the equation (4.11) readsThis inequality is satisfied as long as , which is always true since . Since , the solutions areand soBoth . However, we can show that and . Thus, the only possibility isWe also have access to and thanks to (4.5). Since , we obtainandWe can summarize these results in the following lemma.

Lemma 4.4

(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 for each , and has the following closed form formulawhere and are respectively the effective and basic reproductive number defined in (4.12) and the cardinal of . See Fig. 3 for an illustration.

Fig. 3

Schematic visualisation (red star) of , resp. , in the plane, resp. in the -plane, in the fully symmetric case. The asymptotic value , resp. , lies at the intersection of the diagonal , resp. , and the implicit curve given by (4.7), resp. (4.8) (color figure online)

Case of two vertices

In this simple case, it is possible to build explicit formulas to deal with the implicit submanifold equations (4.7), (4.8), and (4.9). Let and , be respectively the local to vertex basic and effective reproductive number. Then,where the Lambert W function W can be either or . Indeed, the argument of W being negative, two solutions have to be considered. We obviously also haveDue to the definition of the domain of the Lambert W function, the argument has to be greater than . So, the following inequality must be satisfied for (respectively of )Solving the equality part of this inequality, we find thatThis equation has to be verified both for and . Let be defined bywhereThen, the domain of as a function of isConcerning as a function of , we havewithwhereThus,We present on Fig. 4 (left) the functions and defining as a function of and the domain for a given set of the parameters and initial conditions. We refer to Sect. 5 for details regarding the numerical integration of the model and Sect. 6 for further numerical results on the case of two vertices.

Fig. 4

Location of and together with the visualisation of the domains (left) and (right). The final configuration of susceptible individuals lies on the closed curve parametrized by the two branches of the Lambert W function (blue and red curve). We note that as indicated by the red star on the right figure. Values of the parameters are , , , , , , , and initial conditions are set to: , , and . The mass is therefore equal to 1 (color figure online)

Actually, we can reduce the domain of validity of (4.13)–(4.14) for and . Indeed, we know that , , decay with respect to time, so . Moreover, the sum . Thus, we haveThe domain is drawn on Fig. 4 (right).

Concerning and , we can perform the same analysis. LetandWe obtain for ,still with W equal to and . Let and be defined byandwith and given by (4.15) and (4.16). Then,We can show that for . So, we can reduce this domain since . So, we define the domain As a consequence, we have proved the following lemma.

Lemma 4.5

Case of two vertices. Assume that and , where is the cardinal of . The final total population of infected individuals at each vertex , can be expressed aswithwhere and , . Furthermore, we have the sharp boundwith , defined in (4.17)–(4.18). See Fig. 5 for an illustration.

Fig. 5

Location of and , , and visualisation of the domains (left), and domain (right). In both cases, and are represented by a red star. Values of the parameters and initial conditions are similar to Fig. 4 (color figure online)

Remark 4.6

As the solutions , , are simply given by , if we let then we haveWe represent on Fig. 5 the domain .

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 , we denote the space discretization of each edge, and the number of points of the corresponding discretization. For each , the space grid on each edge is given by with . And we let . Let be the time discretization and denote for .

For a given function , its space-time discretization is given by some sequence of vectorsFor each , there exists an integer such thatWe approximate the laplacian on each edge via finite differences. That is, for each ,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 of the graph. So let such that there are edges incident to the vertex. We locally label all these incident edges. For each , we introduce the map such that corresponds to the global index of the grid discretization associated to the vertex v on edge . Finally, we denote by the global index of the nearest neighbor on edge to the vertex v. Note that either or . To approximate the laplacian at a given vertex on edge , we use the following formulaThe unknown can be expressed by discretization of the boundary condition as follows. For each with , we approximate the normal derivative asUsing (2.3), and denoting the time approximation of , we obtain the following expression for As a consequence, we obtain that for each and

The semi-implicit numerical scheme

We introduce the following scheme for each initialized with and some . One can find similar semi-implicit discretization for the SIR part of the model in Sekiguchi and Emiko (2011).

Well-posedness and positivity

We prove that the numerical scheme defined through (5.1) is well defined and preserves positivity under some condition on . Indeed, we first remark that the equation for and in (5.1) can be used to obtain thatsuch that can be expressed only in terms of elements of asAs a consequence, there exists a matrix such thatwhere is such that

Lemma 5.1

There exists a constant , which only depends on the parameters of the system, such that if then we have

is invertible;

if is symmetric for each , then given with , the unique solution of also satisfies .

Proof

Let be such that . Without loss of generality, assume that . If there exists such that for some , then we havewhich is a contradiction by definition of . Next if is such that there is and such that , then we haveThe 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 when and As a consequence, we deduce thatThe last two terms are positive by definition of . Now using Hypothesis 2.2, we have thatsuch that the term in bracket is positive provided thator equivalentlyAs a consequence, we impose thatwhere it is understood that when the positive part is zero there is no condition on . And we have reached a contradiction sinceThis shows that is invertible.

Next let be the unique solution of with . We denote by the vector with components given byOur aim is to evaluate where is the following scalar product on :We divide into three parts:whereThe first and second terms are handled as followsFor the third term , if we further assume that is symmetric, then the matrix is symmetric positive definite, and thus for each there exists some such thatwhile there exists such thatAnd thus, we get an estimate for of the formwhich is positive provided that is small enough. As a consequence, we have proved thatwhich implies that and thus .

The previous lemma demonstrates the well-posedness of our numerical scheme (5.1). It also ensures that if we start with positive initial conditions and , with and , then for all we also have that , , and , provided is small enough and is symmetric for each .

Preservation of total discrete mass

For any , we define the following quantityThe expression is simply the trapezoidal rule applied to the elements of U adapted to our graph . From (5.1), we get thatUpon denoting the following quantitywe get thatNext, we observe thatwhere the cancellation comes from the specific structure of the discretized laplacian through finite differences. As a consequence, we have thatWe also have thatas the sum over the lines of vanishes. And thus we getOn the other hand, from (5.1) we also haveAs a conclusion, we have proved the following result.

Lemma 5.2

Let a solution of (5.1), then we have for each

This is the discrete conter part of conservation of mass for the continuous model.

Now, combining Lemmas 5.1–5.2, we have proved the following theorem.

Theorem 3

There exists a constant , which only depends on the parameters of the system, such that if , then the numerical scheme (5.1) defines a unique sequence . If we further assume that is symmetric for each , 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 , we have

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 while the space discretization to for each .

Case of 2 vertices and 1 edge

We first consider the case where and , where denotes the cardinal of . In this setting, we recall that our model reads as followswith boundary conditionswhere , for , solution ofwhere and . This system is complemented by some initial condition with , , and such that the boundary condition is satisfied initially. Finally, we normalize the total mass as followsFor the numerical simulations, we have fixed initial conditions to be of the formwithwhere and may vary. In Figs. 6, 7 and 8, and are fixed to , while in Fig. 10, is allowed to vary and is fixed to .

Fig. 6

Profiles of the solutions together with the total population on the edge and the total mass of the system M(t) as the parameter is varied from 0.05 to 0.95. All other parameters are fixed and set to , , , and with . For the initial condition we have

Fig. 7

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 respectively while the dark red circles indicate the numerically computed values. Right: relative distance between time of maximal infection in each population, indicated by dark red circles, as the parameter is varied from 0.05 to 0.95. Varying parameters: (top panel), (second panel), (third panel) and (bottom panel) (color figure online)

Fig. 8

Log-plot of the relative distance between time of maximal infection in each population as the diffusion coefficient d and the length of the edge are varied while all other parameters are fixed to , , and with . For the initial condition we have . We note that as d becomes smaller rapidly increases as increases

Fig. 10

Final total populations of infected individuals (left) and corresponding final population of susceptible individuals (right) as the initial population of susceptible individuals is varied from 0.05 to 0.95 while is fixed. The dark blue curves are the location of respectively for each value of , 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.7). All other parameters are fixed to , , , and with (color figure online)

In Fig. 6, we report the profiles of the solutions together with the total population on the edge and the total mass of the system M(t) as the parameter 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 while it has a significant impact on the dynamics at the first vertex. Indeed, as is increased, the maximum of infected individuals is decreased. In the last panel of the figure, we also illustrate the conservation of total population where the fluctuations around is of order . In the top panel of Fig. 7, we present the final total populations of infected individuals and corresponding final population of susceptible individuals as is varied. The blue curve is the location of respectively while the dark red circles indicate the numerically computed values. We recover the fact that has a more significant impact on the final total populations at the first vertex than it has at the second vertex. To get a better understanding of the intricate dynamics between the epidemic at the two vertices, we also present the relative distance between time of maximal infection in each population as is varied. We observe that is not monotone in , as it first decreases and then increases. But we also note that for traducing the fact that the pick of the epidemic occurs at the second vertex before it does at the first vertex, although initially . This illustrates the effect of the diffusion of infected individuals along the edge.

Similarly, in Fig. 7, we report the final total populations of infected individuals and corresponding final population of susceptible individuals as (second panel), (third panel) and (bottom panel) are varied from 0.05 to 0.95. As expected, the final total population of infected individuals at the second vertex decreases as increases while at the first vertex it varies less significantly. As 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 between time of maximal infection is monotonically increasing with . We get the opposite monotonicity properties as is varied.

In Fig. 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 in each infected population . As expected, when the diffusion coefficient is really small while the length is being kept at order one, takes large value: when and . 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 , monotonically decreases as d increases, while at fixed d, monotonically increases as increases.

In Figs. 9 and 10, we vary respectively the initial population of susceptible individuals and infected individuals . We visualize the final total populations of infected individuals and corresponding final population of susceptible individuals on the parameterized surfaces and , respectively and , where the level sets of the parameterized surface are given by the conservation of total mass (4.7). We note that and are almost independent of when with sensible variations only occurring for larger values of . On the other hand, we observe that as is increased the final total population of infected individuals increases at the first vertex while it decreases at the second one. The dependence of as a function of is more subtile and is presented in Fig. 11. In the same figure, we also show the location of and its amplitude. We observe a strong nonlinear dependence with respect to . As increases, we first see that the time at which is maximal increases and then decreases, while is monotonically increasing. The converse is observed at the second vertex.

Fig. 9

Final total populations of infected individuals (left) and corresponding final population of susceptible individuals (right) as the initial population of susceptible individuals is varied from to in log-scale while is fixed. The dark blue curves are the location of respectively for each value of , 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.7). All other parameters are fixed to , , , and with (color figure online)

Fig. 11

Locations of (top) and (bottom left and right) as functions of . Values of all other parameters are similar to Fig. 10. The initial condition is of the form with 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 -plane of the parametrized a curve from Fig. 10, right panel (color figure online)

Case of 3 vertices and 3 edges

Next, we consider the case of 3 vertices and 3 edges arranged in a triangular configuration. For the numerical simulations presented in Fig. 12, we have assumed full symmetry in the parameters that isRegarding the initial condition, we have chosenfor a given , while for each we have set on . 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.7)–(4.8) from Theorem 2.

Fig. 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 respectively for each value of , while the dark red circles indicate the numerically computed values. Parameters were set to , , and , while the initial condition is (color figure online)

In Fig. 13, we tested a different configuration. Upon labeling by A the edge between vertices and , B the edge between vertices and and C the edge between vertices and , we have set the parameters towhileandThe length of each edge is fixed and at each vertex . Finally, we have set different coefficients on each edge, namely , and . Initially, we assume that infected individuals are only present at vertex and each vertex has the same number of susceptible individuals fixed to 1/3. Finally, for each we have set on . We see in Fig. 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.

Fig. 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 after the second wave has reached the first vertex. Parameters values are set in the text

Case of 4 vertices and 3 edges

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 (Griette et al. 2020; Liu et al. 2020). More precisely, we will assume that there exists and such that the transmission rates can be written asfor each and for a given . We will assume that the four vertices are at equal distance such that for each and that the coefficient diffusion are equal on each edge, , . We further assume that at the central vertex exchanges are no longer allowed after locked down. That is, we impose thatwhile for and , together withwhile for and , and alsoFinally, we set for all . Regarding the initial condition, we work withfor given , while for each we have set on .

In Fig. 14, we report the location of the time of maximal infection for each vertex together with the corresponding amplitude as a function of . We observe that below a critical value of , the time of maximal infection always occurs at 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 is decreased the corresponding is decreasing while is increasing up to some value of where we observe a sudden turning point (see the right panel of Fig. 14). We observe that , 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 if and only if , and we findwith our specific values of the initial condition, while we have computedWe also point out that when is below the turning point , the corresponding value of is below . On the other hand, in Fig. 15, we present similar results but this time is fixed and varies. Above some critical value of , saturates to a fixed value independent of 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 .

Fig. 14

Location of the time of maximal infection for each vertex together with the corresponding amplitude as a function of (left) with its projection in the -plane (middle) and a zoom near the turning points (right). Other parameters are set to , , , , and with

Fig. 15

Location of the time of maximal infection for each vertex together with the corresponding amplitude as a function of for two configurations of initial susceptible populations at vertices and , with (left) and (right). Other parameters are set to , , , , and with

Case of vertices and N edges

In our final example, we have considered a network of vertices and N edges arranged in a lattice, in the sense that vertex is only connected to vertices and via two different edges. Figure 16 shows the time evolution of the infected population and susceptible populations 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 while 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 while 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 while for all other vertices, and we note the propagation of two waves activity which collide at the middle vertex . 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 Fig. 17 where we represent the location of 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 (Hodgkin and Huxley 1952; Hupkes and Sandstede 2010) or calcium waves Sneyd (2005). When , they are all aligned on the same line, where for smaller values 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 denotes the maximum as a function of d at the first vertex, we have for while for larger vertices we have the reverse ordering for .

Fig. 16

Time evolution of the infected (left) population and susceptible (right) populations at each vertex for several different initial conditions and . 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 , , , and , while the initial condition is

Fig. 17

Left: Location of the maxima of the infected population at each vertex for several different diffusion coefficients . Right: Zoom of the left figure for small times. We observe that for very small d the location of the maxima is a long a line while it is curved for larger values of d. We also remark that if denotes the maximum as a function of d at the first vertex, we have for while for larger vertices we have the reverse ordering for . Other parameters were set to , , and , while the initial condition is

For the numerical simulations presented in Figs. 16 and 17, we have assumed full symmetry in the parameters that isRegarding the initial condition on the edge, we have set on for each .

Discussion

Summary of main results

In this work, we have proposed a new model that describes the dynamics of epidemic spreading on connected graphs. Our model consists in a PDE-ODE system where at each vertex of the graph we have a standard SIR model and connections between vertices are given by heat equations on the edges supplemented with Robin like boundary conditions at the vertices modeling exchanges between incident edges and the associated vertex. Our first main result is the existence and uniqueness of classical, global in time, solutions of our PDE-ODE model. Our second main result is a complete characterization of the long time behavior of the unique solution of our model. We proved that the final total populations of infected individuals at each vertex are well defined quantities and solutions of a system implicit equations. We also managed to obtain further qualitative properties in the fully symmetric case by exhibiting closed form formula and relate these quantities to standard basic and effective reproductive number for classical SIR model. Next, we proposed and analyzed 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 have further proved 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 the exchanges of fluxes are symmetric at each vertex, then our mass preserving semi-implicit numerical scheme is well-posed and preserves the positivity of the solutions. And finally, we illustrated our theoretical findings with selection of numerical simulations for various types of graphs showing very interesting dynamics such traveling waves.

Biological limitations of the model and possible extensions

The first biological limitation of our model comes from the purely diffusive behavior imposed on each edge of the graph. It would be biologically relevant to incorporate some kind of directed or ballistic motion (drift term in (2.2)) given that many individuals set off on transportation lines travel from one specific city to another preferentially, see Bertaglia and Pareschi (2021) for a purely hyperbolic model. This would result in an equation on each edge of the formfor some speed whose sign will determine the preferred direction of transportation. Such a drift term modifies the boundary conditions at each vertex according towhere is given by where . We expect that the long time dynamics of the system with such a directed motion, which is still mass preserving, will be somehow similar to the purely diffusive case, and that only transient dynamics will be affected. However, we leave it as an open modeling problem at the moment.

A second biological limitation comes from our assumption that the movement of individuals in the susceptible population does not affect its distribution and thus only the infected population is subject to movement. Considering the susceptible population as an ambient population was a first step, and it would be natural to extend our model to the case that individuals in the susceptible population can also move along our transportation network. In the case of spatially extended systems of reaction-diffusion type, it is notorious that allowing the susceptible population to diffuse is more challenging from a theoretical point of view as monotonicity properties of the solutions are lost (Berestycki et al. 2020). We leave such an analysis for a future work.

Scaling limits

Several scaling limits could be considered and we present two directions which seem natural and very promising. First, it would be very interesting to investigate the limit of large diffusivity along the edges of the graph. One should be able to recover a class of network-based ODE models that have been studied in the literature. Indeed, we expect that in this limit, one should obtain an ODE for some (now independent of ) along the edges that depends only on the population levels at the two connecting nodes. This should lead to an ODE system at the vertices of the graph that is globally coupled to .

Another natural limit to investigate is the one where the number of vertices of the graph goes to infinity along with the length of the edges which shall converge to zero. One would expect to obtain in the limit a PDE system describing the evolution of spatially continuous populations of susceptible, infected and removed populations. Depending on the initial geometry of the graph the limiting PDE will either be posed on the real line or on some two-dimensional domain. The exchange term between infected individuals at the edges, given by , is likely to produce in the limit a spatial diffusion term and we expect to recover reaction-diffusion like models such as the ones in Aronson (1977), Berestycki et al. (2020) and Berestycki et al. (2020).