Propagation of Epidemics Along Lines with Fast Diffusion

It has long been known that epidemics can travel along communication lines, such as roads. In the current COVID-19 epidemic, it has been observed that major roads have enhanced its propagation in Italy. We propose a new simple model of propagation of epidemics which exhibits this effect and allows for a quantitative analysis. The model consists of a classical SIR model with diffusion, to which an additional compartment is added, formed by the infected individuals travelling on a line of fast diffusion. The line and the domain interact by constant exchanges of populations. A classical transformation allows us to reduce the proposed model to a system analogous to one we had previously introduced Berestycki et al. (J Math Biol 66:743–766, 2013) to describe the enhancement of biological invasions by lines of fast diffusion. We establish the existence of a minimal spreading speed, and we show that it may be quite large, even when the basic reproduction number \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$R_0$$\end{document}R0 is close to 1. We also prove here further qualitative features of the final state, showing the influence of the line.

Most of the models that are used rely on various extensions of the classical SIR cornerstone model of epidemiology. That is, they use population compartmental models that contain various additional compartments to the SIR ones to account for segments of the populations that are exposed, asymptomatic, presymptomatic, treated, etc. Such models yield the evolution of the infected population at a given level of territorial granularity (whole countries, regions, counties or cities). The spatial interplay aspect, mostly overlooked, is included by involving transfer matrices of populations and infected between various patches each of which being considered as uniform.
Yet, the propagation of COVID-19 exhibits remarkable spatial structure properties. Indeed, the spatial organization and spreading of epidemics in general reveal important features of the transmission process. This is especially true in relation with movements of individuals who carry with them infectious characteristics. It has long been known, since ancient times, that epidemics travel along lines of communications. In the black death epidemic of the 14th century, contagion advanced along roads connecting trade fair cities, and from there spread inwards leading to a front like invasion of Western Europe roughly from South to North, pulled by these roads. Such a mode of propagation was also at work for the propagation of rumours. (See for instance the propagation of the "big fear" in France, after the Revolution. Compare, for example, the presentation given by Siegfried (1960) and the analogies between these two phenomena.) In studying the early spread of HIV virus in Faria et al. (2014), Faria et al. pointed out that the virus mostly travelled along the main lines of communications (railways and waterways), see Fig. 1 below.
Some current studies are bringing to light a similar effect in the spread of the COVID-19 virus in Italy. Gatto et al. (2020) and Sebastiani (2020) have established that the coronavirus spread foremost along the main expressways. They argue indeed that cities located along the main North-South and East-West highways in Italy have faced an earlier and stronger contagion than cities with similar or higher population sizes but not located on these roads. The spatial signatures of the spreading of the epidemic are clearly revealed by Fig. 2 taken from the work by Gatto et al. (2020).
Then, Fig. 3 is a snapshot from a video in Gatto et al. (2020) which accounts for hospitalization. It highlights the radiation of the epidemic along highways and transportation infrastructures. Figure 4 below, taken from Sebastiani (2020), depicts the Italian cities with more than 1000 infected individuals reported as for April 5. For two of them, Piacenza and Cremona, which are located, respectively, at the crossroad of these two main highways and 40 Km away from it, infected people are above 1% of the population. We point out that the cities represented in this figure are by far not the most populated ones in Italy.
That this effect of roads still matters in such a global and "modern" epidemic as COVID-19 shows that spatial diffusion still plays an important role in the spreading of epidemics. At an early stage, long distance "jumps" through the air transport network played the key role in the global dissemination of the pandemic. However, the role of ground transportation became relevant in a second phase. This is even more true after the lockdown in many countries, that all but halted all flights. Even though ground travel was limited, it was never interrupted and, furthermore, ground transportation was essential in transporting all needed supplies. With respect to patch models that include transfer matrices, it thus transpires that the effect of lines is of a different nature. And as we have seen, these lines can be roads, railroads or rivers.
The aim of this paper is to propose a new model to account for such effects and then to study quantitatively how a line acts on the overall epidemics propagation. We thus introduce a model that we call the S I RT system, standing for Susceptibles,  Sebastiani (2020). The 33 Italian cities where, as of April 5, the total number of infected people reported is larger than 1,000; the streaked circles correspond to Piacenza and Cremona, where infected people are above 1% of the population. Reproduced with permission of G. Sebastiani (Color figure online) Infected, Recovered, and Travelling infected. This model takes explicitly into account the existence of a line along which infected individuals can travel with a specific diffusion coefficient. Our aim here is to gain insight into this spreading aspect at the fundamental mathematical level of a SIR type model that now incorporates the possibility of infected to travel along a specific line.
The SIR model dates back the fundamental paper of Kermack and McKendrick (1927), Kendall (1956) introduced the SIR model with spatial interaction in the dis-cussion of a statistical study of measles by Bartlett (1957). These models may be rewritten as nonlinear integro-differential equations in time and space variables, for which the study of spreading dates back to the 1970's with, in particular, the milestone papers of Aronson (1977) and Diekmann (1979). We choose in this paper to restrict our model to local (Brownian) diffusion and local interaction. Regarding the existence of travelling waves for the local (homogeneous) diffusion model, see the studies of Källén (1984), Ilyas (1994, 1995). For a survey, see Murray (2002).
The model we propose here aims at a fundamental mathematical understanding of this effect. We look at a stylized situation-a laboratory case as it were. The population of susceptibles diffuse in a homogeneous open territory, that we take to be a half plane. Then, the infected can travel along the line bounding this half plane. The populations of infected in the open territory and on the lines are in constant exchange. This is represented in our model by transmission parameters between these two populations (stable and travelling infected). We aim at understanding the effects of such lines on the speed of propagation of the epidemics, and how they affect the balance of the total number and locations of infected individuals.
As a benchmark, we use the classical SIR model with diffusion, for which we briefly recall the basic results and how to compute the speed of spreading. The latter coincides with that of Aronson and Weinberger (1978), we simply rewrite the formulas in terms of the basic reproduction number R 0 .
Let us mention that the SIR model falls into a more general class of systems in which one equation-the one for I here-exhibits a self-reinforcement mechanism when the other unknown function -S-is above some threshold level. Berestycki et al. (2020) develop a mathematical study of this class of systems, called Activity/Modulator. As described there, SIR models also arise in a variety of contexts to model contagion phenomena. They are particularly useful in the study of collective behaviours, such as social unrest. The work by Bonnasse-Gahot et al. (2018) developed such a system to model the spread of riots in France in 2005.
Of course, more realistic epidemiology models will incorporate networks rather than single lines and, likewise, the remaining propagation does not take place in a homogeneous territory. From this simplified model, one can nonetheless deduce a more realistic one involving a network of roads and a distribution of cities. Regarding the latter aspect, Bonnasse-Gahot et al. (2018) use explicitly such a network of cities.
This family of models lends itself to various extensions. But since we want to derive the mathematical properties of this system we only consider here the stylized model. We hope that this stylized model can shed some light on how lines influence the global unfolding of an epidemic.

The SIR Model with Infected Diffusion
The classical SIR model divides the overall population in three compartments: susceptibles S, infected (and infectious) I and recovered R. Since the total population N is viewed as fixed, one derives the latter in a straightforward manner from the former two: R = N − S − I . Therefore, we do not mention explicitly this function henceforth. When taking into account spatial dependence, we view the unknowns S and I as depending on both time t and space location X ∈ R 2 . We consider that the individuals in the susceptible population S do not move around (or rather that movement does no affect its distribution). We can think of S as the ambient population. Thus, we assume that only the infected population I is subject to movement. We choose here to represent this movement as a pure local diffusion that can be viewed as a limiting Brownian movement of individuals. In the second part of this work [8], we consider the case of non-local diffusions. We denote by d the diffusion coefficient. We are thus led to the following spatial diffusion SIR model: The parameter β > 0 represents the transmission rate, while α > 0 is the recovery rate. The system is supplemented with initial conditions. We assume that the initial distribution of S, S(0, X ) ≡ S 0 is constant (and positive) and that I (0, X ) = I 0 (X ) is compactly supported (and nonnegative). The cumulative number of infected at location X , at time t, is given by t 0 I (t, X )dt. Several authors have considered this model or closely related ones with integral formulations. We mention some of the references in Sect. 3.1 below. For the sake of completeness, we derive here the main results concerning this system (1): existence and properties of a final state to which the solutions and the cumulative number of infections at each location converge, characterization of epidemic spreading, and asymptotic speed of propagation. These properties involve the classical basic reproduction number R 0 := S 0 β/α. As we will recall, the position of R 0 with respect to 1 determines here too a threshold for the epidemic to spread. These various properties will serve as benchmarks for our study of the effect of the presence of a road.

A "SIRT" Model in the Presence of a Road
Let us now introduce this system, which we call a road/field model. It represents a situation where there is a road on which infected individuals can travel with a specific diffusion coefficient D. We are especially interested in the case of a large D. The aim is to understand in what way such a road alters the epidemic spreading, in particular if it enhances its propagation and in what measure exactly. The starting point of our analysis consists in distinguishing the infected individuals which are moving on the line with fast diffusion from the ones present in the rest of the territory, by using two distinct density functions. This is the same modelling hypothesis we made in our previous paper Berestycki et al. (2013) to describe the dynamics of a single population in the presence of a road.
As before, we use the densities S(t, x, y) and I (t, x, y) of susceptible and infected individuals at time t ≥ 0 and position X = (x, y) ∈ R 2 . We still assume that we can ignore the movement of susceptibles. We introduce a new compartment of the population that we call T (t, x) standing for "travelling individuals". This is the density of infected individuals on the road, that is the line R. It is worth emphasizing that I (t, x, y) and T (t, x) are two different compartments, in particular T (t, x) is not the same as I (t, x, 0). Both populations are assumed to diffuse, but with different diffusion coefficients: D for T on the line and d for I in the plane. The two compartments interact by a constant exchange of individuals: at each time t > 0 and point x ∈ R, the line yields an amount μT (t, x) of individuals to the domain and receives an amount ν I (t, x, 0) of individuals from the domain.
To shed light on the effect of such a road, it is enough to consider the interplay of a half-plane with the road that bounds it. Indeed, results in this framework can be translated for results in the whole space by straightforward symmetry arguments. But it can also be seen that propagation properties in the plane in general can easily be deduced from the half-plane case. (For such a discussion, we refer to our earlier paper for the KPP equation with a road Berestycki et al. (2013)). We prefer to carry our analysis in the half-plane case for the sake of clarity and to simplify notations.
Thus, we consider the following system for the unknown S, I , T : By symmetry, we can reduce to study the problem in the upper half-space R × (0, +∞). We end up with the following system: We assume a uniform initial susceptible density: S(0, x, y) ≡ S 0 > 0. The initial density of infected individuals is assumed to be zero outside a bounded region: (T (0, x), I (0, x, y)) = (T 0 (x), I 0 (x, y)) are nonnegative and compactly supported in R and R × [0, +∞), respectively, with in addition I 0 ≡ 0, but possibly T 0 ≡ 0. Actually, to simplify matters, most of the time we will look at the case T 0 ≡ 0, not that it changes much in the proofs.

Organization of the Paper
In the next section, we will first analyse systems (1) to have sound benchmarks that will be useful to discuss the effects of the road. Then, the remaining of the paper will be devoted to the study of (2). We will start by discussing the stationary limiting state and see whether the epidemic spreads or not in Sect. 3.2. There we will also derive the asymptotic speed of propagation for the spatial spread of the epidemic. In Sect. 4, we discuss how the presence of the road affects the distribution of the total number of infected per location. Section 5 is devoted to the mathematical proofs of these various results. These depend on the various parameters, and we discuss their influence in Sect. 6.

Analysis of the Models
We are going to compare the behaviour of (2) to that of the classical SIR model in the whole plane, with no diffusion of the susceptible population, that is (1). The proofs of the results (in a slightly broader framework) are presented in Sect. 5.1.

Benchmark: The SIR Model
System (1) can be reduced to the classical Fisher-KPP equation by noticing that S(t, X ) is easily computed from I (t, X ). Indeed, calling v the cumulative number of infected, integrating the second equation in (1) yields Then, taking the sum of the equations in (1) and integrating in time in (0, t) we find the equation for v: Transform (3) is a well-known particular case of a broader class (see for instance Aronson (1977), Diekmann (1978)) that reduces the SIR model with nonlocal interactions, and with no diffusion on the susceptible individuals, to nonlinear integral equations. In this particular case, we retrieve a parabolic equation. Pulsating waves and spreading speeds for periodic S 0 are studied for (1) by Ducrot and Giletti (2014). The integral equation resulting from the model with nonlocal interactions is studied by Ducasse (2020). The nonlinearity f is concave and vanishes at 0, it is then of the KPP type. Only the presence of the "source" term I 0 differs from the standard Fisher-KPP equation. However, being I 0 compactly supported, the dynamics of the equation is still governed by the sign of f (0). Using a notation commonly employed in the literature, we write The quantity R 0 can be viewed as the classical basic reproduction number, see for instance Källén (1984) and the discussion on such number and its interpretation in Diekmann et al. (1990). This parameter plays a crucial role and is widely mentioned by experts, decision makers and media for the analysis of the present COVID-19 pandemic. The sign of f (0) is determined by the position of R 0 with respect to the threshold value R 0 = 1. The presence of I 0 in (4) prevents the existence of constant steady states, making the dynamics of the equation more complex. Nevertheless, the following Liouville-type result holds.
Theorem 3.1 The equation in (4) admits a unique positive, bounded, stationary solu- where v * is the unique positive zero of f .
The proof of this result, as well as the others in this section, are presented in Sect. 5.1. They can also be found in Ducrot and Giletti (2014), at least for the case R 0 > 1, in the more general framework of periodic coefficients α(X ), β(X ) and distribution S 0 (X ). Some additional qualitative properties of v ∞ are contained in Theorem 5.1 below. Namely, v ∞ is radially decreasing outside the support of I 0 and it has the shape of a bump above the value 0 or v * .
It turns out that the steady state v ∞ is a global attractor for the dynamics of (4).
Since the above convergence holds true for the time derivatives (see the proof in Sect. 5.1), a first consequence we derive is: that is, the number of infected individuals drops to 0 asymptotically in time.
One can also read Theorem 3.2 in terms of number of remaining susceptibles at time t and position X , which is S(t, X ) = S 0 e −βv(t,X ) . Hence, the amount of people that will be infected by the virus at a given place X , throughout the whole course of the epidemic, is given by We point out that the fact that this function is not constant is a consequence of the hypothesis that susceptibles do not diffuse. Otherwise, the diffusion tends to "flatten" the density S, and this mechanism occurs for rather general systems, see, e.g. (Berestycki et al. 2020, Theorem 4.3). For this reason, the description of the function v ∞ is extremely important in the study of the epidemic. With this regard, Theorem 3.1 leads to the following crucial dichotomy for the total amount of infected people far from the epicentre: This corresponds to two opposite scenarios. If R 0 ≤ 1, the epidemic does not propagate and areas very far from the initial outburst (the support of I 0 ) will be essentially not infected. Conversely, if R 0 > 1, the epidemic propagates across the territory, and, even though its impact will be weaker on places far from the epicentre, the total infected people there will be a portion 1 − e −βv * of the overall population.
What is important to determine in the case R 0 > 1 is the speed at which the epidemic spreads. This is provided by the following.
The quantity c S I R is the asymptotic speed of spreading of the epidemic wave. It coincides from one hand with the speed 2 d f (0) of the standard Fisher-KPP equation, and from the other with the minimal speed of the waves for the SIR model (1) obtained by Källén (1984).
A word must be said about the proof, presented in Sect. 5.1. It is essentially a direct consequence of Theorem 3.2 above and Aronson and Weinberger (1978). Indeed, the former describes the behaviour of the solution in compact regions, while the latter provides the estimates far from the origin, where the influence of I 0 is negligible.

Dynamics of the Model in the Presence of the Line
It occurs that the very same transform works for the model with the line (2): with, as before, This is the same system studied in Berestycki et al. (2013), except for the "source" terms I 0 , T 0 . We recall that these are assumed here to be nonnegative and compactly supported, with I 0 ≡ 0 (and typically T 0 ≡ 0).
We start with the Liouville-type result, analogous to Theorem 3.1. (7) admits a unique positive, bounded, stationary solution

Theorem 3.4 System
where v * is the unique positive zero of f .
The long-time behaviour for (7) is described by the following result, which is the counterpart of Theorem 3.2 about the model without the line.
At this stage, the picture is identical to the standard SIR model described in the previous section. The total amount of infected individuals after the passage of the epidemic wave, at a given point (x, y), is with v r ∞ exhibiting two qualitatively distinct behaviours depending on whether R 0 ≤ 1 or R 0 > 1. Namely, the total number of infected people very far from the epicentre of the epidemic is lim |(x,y)|→∞ which is the same as in the case without the line.
Next, we investigate the speed at which the epidemic spreads across the territory, along the line.
locally uniformly with respect to y ≥ 0.

In addition, the spreading speed c T S I R satisfies
The asymptotic speed of spreading c T S I R coincides with the one for the homogeneous model, i.e. with I 0 ≡ 0 and T 0 ≡ 0, which is provided by (Berestycki et al. 2013, Theorem 1.1). From a mathematical point of view, this means that the presence of the compact perturbation does not affect the speed of propagation, as in the case of the SIR model (4). The fact that the speed cannot decrease if one adds the perturbation is a straightforward consequence of the comparison principle. Instead, to derive the opposite inequality, we need to go into the proof of (Berestycki et al. 2013, Theorem 1.1) and use the supersolutions constructed there. This is done in Sect. 5.2 below.
Theorem 3.6 shows that the presence of the line has a true impact on the speed at which the epidemic spreads. Indeed, if the diffusion coefficient on the line D is larger than twice the one in the field d, the asymptotic speed of spreading in the direction of the line is enhanced, compared with the standard one c S I R . How much the spreading is enhanced is discussed in Sect. 6. One can then wonder what is the effect of the line on the other directions. In the case I 0 , T 0 ≡ 0, we have shown in Berestycki et al. (2016a) that the enhancement of the speed occurs in a cone of directions around the line, with an associated critical angle. We suspect the same scenario to hold true for system (7).
We conclude this section by showing that our model accounts for a true epidemic wave, in the following sense. At every point of the domain under consideration, the number of infected individuals, that was initially close to 0, raises to a nontrivial level around a certain time, then decays back to 0.
Proposition 3.7 Assume that R 0 > 1. There is a constant T * > 0, a positive function I * : [0, +∞) → R and a function τ * : R → R satisfying for which the following is true for any solution to (2).
2. (decay far from τ * (x)). The following limits hold uniformly in x ∈ R and locally uniform in y ≥ 0: This proposition is proved in Sect. 5.2. We note that, for the true SIR model, one has a much stronger result, as one may relate the maximum of I to the maximum of the derivative of the one-dimensional Fisher-KPP wave. While we do have the existence of travelling wave here (see Berestycki et al. (2016b)), we do not have precise convergence results for the solution of the Cauchy Problem. This will be done in a future study.

Effects of the Line on the Total Number of Infected Per Location
We recall that the total number of infected people at a given location (x, y), for both models without and with the road, is given by where v ∞ is either v ∞ or v r ∞ , that is, (the v component of) the unique positive, bounded, stationary solution of the problem. A pointwise comparison between v ∞ and v r ∞ would then tell us in which way the presence of the road affects the total number of infected, depending on the location. In order to compare the two problems with and without the road, we consider for the SIR model (1) an initial datum I 0 (x, y) in R 2 which is symmetric with respect to the x-axis, and for the S I RT model (2) the same I 0 restricted to R × [0, +∞) and T 0 ≡ 0. Here is our result.
Theorem 4.1 Suppose that R 0 = 1, that I 0 (x, y) is an even function of y and that T 0 ≡ 0. Then, there exist two subsets This result shows that the presence of the road increases the number of infected people in some areas and decreases it in some others. The next theorem states a further property of the limiting state.

Theorem 4.2
The stationary solution to the problem with the road (7) satisfies Moreover, a * is independent of S 0 , I 0 , T 0 and satisfies This theorem will allow us to specify the region where the enhancement due to the road takes place. It will be seen to happen close to the road, far from the epicentre of the epidemic. Recall that v ∞ and v r ∞ have the same limit at infinity, namely 0 if R 0 ≤ 1 and v * , the positive zero of f , if R 0 > 1. The two solutions differ by their rate of exponential decay towards the limit state. Indeed, Theorem 5.1(ii)-(iii) below implies that the rate of exponential decay of v ∞ is whereas Theorem 4.2 shows that the decay of (u r ∞ , v r ∞ ) is strictly slower. Observe that the decay of v ∞ is the natural one provided by the linearized equation (outside supp I 0 ) and it is obtained through rather standard arguments. On the contrary, the decay of v r ∞ is slower than the natural one. The proof of Theorem 4.2, presented in Sect. 5.2, relies on some ideas from Berestycki et al. (2013), the keystone consisting in the construction of suitable sub-and supersolutions.
Let us sketch the argument. To start with, we look for a stationary solution to (7) as a perturbation of the limit at infinity: where, for notational simplicity, we have set v * := 0 in the cases R 0 ≤ 1. Dropping the o(h) terms, outside the supports of I 0 and T 0 we get the linearized system Observe that f (v * ) ≤ 0. We seek for solutions in the form with a, b, γ > 0. After some computation, the problem reduces to the algebraic system If R 0 = 1, i.e. f (v * ) < 0, this system admits a unique positive solution, that we denote (a * , b * , γ * ), whose first two components are represented in Fig. 5 (observe that a * < − f (v * )/d, which is the radius of the disk). If instead R 0 = 0, i.e. f (v * ) = 0, the unique solution of (14) is (a * , b * , γ * ) := (0, 0, μ/ν). In both cases, the pair ( u, v) will serve as a supersolution to (7) and will provide the upper bound in Theorem 4.2. For the lower bound, we need to modify the construction in order to obtain a subsolution supported in a strip. First of all, we penalize the reaction term in the field by considering a parameter ζ > − f (v * ) ≥ 0. Next, we modify the above function v, considering now the pair where 0 < ε < 1. The linear system (12) with f (v * ) replaced by −ζ reduces to the new algebraic system This system converges to (14) as (ζ, ε) → (− f (v * ), 0), and thus its unique positive solution, denoted by (a This leads to a family of bounded subsolutions with a decay in the x-variable arbitrarily close to a * .

The SIR Model
We give here the proofs of the results of Sect. 3.1, and we also derive some further qualitative properties of the solution. We consider a more general framework: we set problems (1), (4) in arbitrary space dimension N ≥ 1, i.e. X ∈ R N , and for initial data v(0, X ) not necessarily identically equal to zero, but rather nonnegative and compactly supported. This does not entail any additional difficulty.
Proof (Theorem 3.1) On one hand, the function identically equal to 0 is a subsolution to the stationary equation On the other hand, the constant functionv ≡ s is a supersolution of this equation provided s > 0 is sufficiently large, because I 0 is bounded and f (+∞) = −∞. The existence of a solution 0 ≤ v ∞ ≤ s then follows from the standard sub-/supersolution method. Actually, as 0 is not a solution, the elliptic strong maximum principle yields v ∞ > 0 in R N . Let us now derive the limit as |X | → ∞. By elliptic estimates, any sequence of translations v ∞ (X + X n ), with (X n ) n∈N diverging, converges (up to subsequences) towards a nonnegative, bounded solution v to If R 0 ≤ 1, i.e. f (0) ≤ 0, then f < 0 on (0, +∞), from which one readily derives that necessarily v ≡ 0 (for instance by comparison with the ODEu = f (u)). This shows that v ∞ (X ) → 0 as |X | → ∞ when R 0 ≤ 1.
In the case R 0 > 1, we remark that v ∞ is a supersolution to the classical Fisher-KPP equation for which the "hair-trigger" effect holds, see Aronson and Weinberger (1978): any solution with a positive, bounded initial datum converges as t → +∞, locally uniformly in space, to the positive zero v * of f . We infer by comparison that v ∞ ≥ v * in R N . This shows in particular that v is positive, and thus applying the "hair-trigger" To prove the uniqueness, we distinguish again the cases R 0 > 1 and R 0 ≤ 1. Case R 0 > 1. We need to show that, for any pair of positive, bounded solutions v 1 , v 2 to (18), there holds that v 1 ≤ v 2 in R N . Assume by way of contradiction that Since v 1 , v 2 → v * at infinity, the above supremum is a maximum, attained at some pointX . Subtracting the equations, we get which, evaluated at the pointX (where v 1 = kv 2 and Δ(v 1 − kv 2 ) ≤ 0) yields But this is impossible because, being concave and vanishing at 0, f satisfies Case R 0 ≤ 1. Now, any given positive, bounded solutions v 1 , v 2 to (18) tend to 0 at infinity. Take ε > 0 and call v ε Assuming by contradiction that v 1 > v ε 2 somewhere and repeating the same arguments as before with v 2 replaced by v ε 2 , we end up with the inequality at some pointX and for some k > 1. As seen before, this is impossible. This means that v 1 ≤ v ε 2 in R N , which, by the arbitrariness of ε, entails the desired inequality v 1 ≤ v 2 .
The following result contains some additional properties of v ∞ .

Theorem 5.1 The stationary solution v ∞ (X ) satisfies the following properties:
(i) v ∞ is radially decreasing outside the support of I 0 , i.e.
where δ is such that X · e ≤ δ for all X ∈ supp I 0 ;

Proof Statement (i).
We make use of a reflection argument due to Jones (1983). Fix a direction e ∈ S N −1 and take r > δ, where δ is such that X · e ≤ δ for all X ∈ supp I 0 . Let T be the reflection with respect to the hyperplane {X · e = r }, and define v := v ∞ • T . In the half-space {X · e < r } (⊃ supp I 0 ), the function v solves dΔ v + f ( v) = 0, hence it is a subsolution of the equation satisfied by v ∞ . In addition, v ≡ v ∞ on {X · e = r }.
Repeating the comparison argument in the proof of Theorem 3.1 with v 1 = v and v 2 = v ∞ , but in {X · e < r }, and observing that the pointX cannot belong to the boundary {X · e = r }, one infers that v ≤ v ∞ in {X · e < r }. Actually, the strong maximum principle and Hopf's lemma imply that the inequality is strict and moreover Since ∇ v(re) · e = −∇v ∞ (re) · e, this gives the desired monotonicity.

Statements (ii)-(iii). We have shown in the proof of Theorem 3.1 that, in the case
The strict inequality follows from the strong maximum principle.
We simultaneously derive properties (20)-(21). To do so, we set v * = 0 in the case R 0 ≤ 1. For given e ∈ S N −1 , k > 0 and λ = − f (v * )/d, consider the function w(X ) := ke −λX ·e . Using the concavity of f , we derive that is, v * + w is a supersolution to (18) outside supp I 0 . Let ρ > 0 be such that supp I 0 ⊂ B ρ . We choose k := e λρ max v ∞ , so that v * +w ≥ v ∞ in B ρ ⊃ supp I 0 . As a consequence, the function min v ∞ , v * +w is a generalized supersolution to (18). It then follows from the sub-/supersolution method, that there exists a solution 0 ≤ v ≤ v ∞ , v * + w , which is actually positive by the strong maximum principle. Then, v ≡ v ∞ thanks to the Liouville result of Theorem 3.1. We have thereby shown that This being true for all e ∈ S N −1 , we obtain the upper bound Let us derive the lower bound. Take and define z := e −λ|X | . This function satisfies, for |X | ≥ dλ(N −1) It follows that, for h > 0 sufficiently small, the function v * +hz is a subsolution to (18) for +hz is a generalized subsolution to (18). By the sub-/supersolution method, there exists a solution v ≤ v ≤ s, where s is such that f (s) < − max I 0 . Theorem 3.1 eventually yields v ≡ v ∞ . This shows the lower bound Call The estimates (22)-(23) yield, for |X | sufficiently large, from which the limits in (20)-(21) follow due to the arbitrariness of ε > 0.
We now turn to the results about the Cauchy problem (4), with an initial datum v(0, X ) nonnegative and compactly supported.
Proof (Theorem 3.2) Let v, v be the solutions to the Cauchy problem emerging from the initial data identically equal to 0 and s, respectively, with s > 0 large enough so that f (s) + max I 0 < 0. We further take s > max v(0, ·). The parabolic comparison principle yields Since the initial data of v, v are, respectively, a sub-and a supersolution of the problem, the comparison principle implies that v, v are, respectively, increasing and decreasing in t and then, by parabolic estimates, they converge to two stationary solutions 0 < V ≤ V < s. The convergences occur locally uniformly in space and hold true for the time derivatives. Theorem 3.1 implies that V ≡ V ≡ v ∞ . The proof is complete.
Proof (Theorem 3.3) Take ε ∈ (0, c S I R ) and consider a sequence (t n ) n∈N diverging to +∞ and a sequence (X n ) n∈N in R N such that |X n | ≤ (c S I R − ε)t n . If (X n ) n∈N is bounded, we already know from Theorem 3.2 that v(t n , X n ) − v ∞ (X n ) → 0 as n → ∞. Suppose now that (X n ) n∈N diverges (up to subsequences). Recall that v is a supersolution to the standard Fisher-KPP equation (19), for which spreading occurs with the asymptotic speed 2 d f (0), that is exactly c S I R . Because v(1, ·) > 0, we infer that To derive the upper bound, we consider the same function w(X ) := ke −λX ·e as in the proof of Theorem 5.1(ii), with λ = − f (v * )/d and e ∈ S N −1 , k > 0. We have seen that v * + w is a supersolution to (18) outside supp I 0 . We then take k large enough, independently of e, so that v * + w(X ) > v(t, X ) for all t > 0 and X ∈ supp I 0 ∪ supp v(0, ·). Hence, by comparison, v(t, X ) < v * + w(X ) for all t > 0, X ∈ R N . This being true for any e ∈ S N −1 , with k independent of e, yields v(t, The proof of the first limit stated in the theorem is achieved.
Let us deal with the second limit. For c = c S I R = 2 d f (0), λ = c SI R 2d and any given e ∈ S N −1 and k > 0, the function w(t, Hence, by the concavity of f , it is a supersolution to (4) outside supp I 0 . We choose k large enough, independently of e, in such a way that w(0, X ) > v(t, X ) for all t > 0 and X ∈ supp I 0 ∪ supp v(0, ·). Hence, by comparison, v(t, X ) < w(t, X ) for all t > 0, X ∈ R N , and therefore, letting e vary in S N −1 , we get v(t, x) ≤ ke −λ(|X |−c SI R t) , which gives the desired limit.
Using the same reflection argument as in the proof of Theorem 5.1, one can show that the solution v(t, X ) to (4) is radially decreasing with respect to X outside the support of I 0 .
We conclude this section with a monotonicity result with respect to the diffusion coefficient d. We assume now for simplicity that the initial datum v(0, ·) is identically equal to 0.

Proposition 5.2 Under the assumption that I 0 is radially decreasing, the values v ∞ (0)
and v(t, 0), for any t > 0, are decreasing with respect to the diffusion coefficient d.

Proof
Let v 1 , v 2 be the solutions to (4) associated with the coefficients d = d 1 and d = d 2 , respectively, with 0 < d 1 < d 2 . The function v j (t, Thus, owing to the monotonicity of I 0 , the comparison principle yields v 1 (t, X ) > v 2 (t, X ) for all t > 0, X ∈ R N , which, evaluated at X = 0, gives v 1 (t, 0) > v 2 (t, 0). Moreover, because v j converges towards the steady state v j ∞ as t → +∞, thanks to Theorem 3.2, one further derives This inequality is actually strict due to the strong maximum principle. We infer that v 1

The SIRT Model
We will make repeatedly use of the weak and strong comparison principles for the road-field system. They are provided by (Berestycki et al. 2013, Proposition 3.2) in the case I 0 ≡ 0, and one can check that the presence of the bounded source term I 0 does not affect their proofs. When applied to a stationary subsolution (u, v) and supersolution (u, v) satisfying (u, v) ≤ (u, v), the strong comparison principle implies that the inequality is strict unless (u, v) ≡ (u, v). Here and in the sequel, inequalities are understood component-wise.
Proof (Theorem 3.4) We start with constructing a stationary supersolution to (7). Take 0 < σ < √ α/d and define the functions where K is a positive constant that will be fixed later. The pair (u, v) satisfies the last equation of (7). Moreover, we compute We then choose K sufficiently large so that the above right-hand sides are larger than max T 0 and max I 0 , respectively. Then, (u, v) is a supersolution to (7). The existence of a positive, stationary solution follows from the sub-/supersolution method, applied with (0, 0) as a subsolution and (u, v) as a supersolution. Owing to the strong maximum principle, this provides us with a solution We derive the uniqueness result, as well as the limit at infinity, by distinguishing the cases R 0 ≤ 1 and R 0 > 1.
Case R 0 > 1. The positive steady state (u r ∞ , v r ∞ ) is a supersolution to the problem with I 0 , T 0 ≡ 0, which reduces to the system studied in Berestycki et al. (2013). For such system, we know from (Berestycki et al. 2013, Theorem 4.1) that positive solutions converge as t → +∞ to the steady state (ν/μ, 1)v * (v * = 1 for the f in Berestycki et al. (2013)), whence by comparison, Consider a sequence ((x n , y n )) n∈N in R×R + , with (y n ) n∈N diverging to +∞. Then, by elliptic estimates, as n → ∞, v r ∞ (x + x n , y + y n ) converges locally uniformly (up to subsequences) towards a solution v of the equation (25). Then, as seen in the Proof of Theorem 3.1, we necessarily have v ≡ v * . This proves the second limit of the theorem in the case R 0 > 1.
Take now a diverging sequence (x n ) n∈N in R, and let ( u(x), v(x, y)) be the limit of (a subsequence of) (u r ∞ (x + x n ), v r ∞ (x + x n , y)), whose existence is guaranteed by elliptic estimates up to the boundary. Thus, ( u, v) is a stationary solution of (7) with I 0 , T 0 ≡ 0, which is positive due to (25). It follows from (Berestycki et al. 2013, Theorem 4.1) that ( u, v) ≡ (ν/μ, 1)v * . This proves the first limit stated in the theorem (the uniformity in y following from the first limit).
It remains to prove the uniqueness. Let (u 1 , v 1 ) and (u 2 , v 2 ) be two pairs of positive, bounded, stationary solutions to (7). Assume by way of contradiction that Because of the limits we have just proved, one of the following situations necessarily occurs: Suppose we are in the latter case. Then, exactly as in the proof of ..., the concavity of f prevents the maximum from being achieved in the interior of R × R + . Then, it is achieved at some point (x, 0), and Hopf's lemma yields Using the third equation in (7), together with v 1 (x, 0) = kv 2 (x, 0), we find that which contradicts the definition of k. Consider the remaining case: Computing the difference of the equations satisfied by ku 2 and u 1 at a pointx where this maximum is achieved, we derive This is impossible because v 1 < kv 2 .
We have thereby shown that k ≤ 1, that is, (u 1 , v 1 ) ≤ (u 2 , v 2 ). Exchanging the roles of the solutions yields the uniqueness result.
Case R 0 ≤ 1. We start with the uniqueness result. We derive it in the more general framework of nonnegative solutions, with possibly I 0 ≡ 0. We need to show that, for any two pairs (u 1 , v 1 ), (u 2 , v 2 ) of nonnegative, bounded, stationary solutions to (7), there holds that (u 1 , v 1 ) ≤ (u 2 , v 2 ). Assume by contradiction that, on the contrary, Suppose first that sup R×R + (v 1 − v 2 ) = h, and let ((x n , y n )) n∈N be a maximizing sequence. If (y n ) n∈N is bounded from below away from 0, then the functions v j (x + x n , y + y n ) converge locally uniformly (up to subsequences) towards two solutions v j of the equation −dΔ v j = f ( v j ) in a neighbourhood of the origin. Moreover, This is impossible, because f (0) = α(R 0 − 1) ≤ 0 and hence f is decreasing on R + . If, instead, y n → 0 (up to subsequences), then the pairs (u j (x + x n ), v j (x + x n , y)) converge locally uniformly (up to subsequences) towards two solutions ( u j , v j ) of the same system, which is of the form (7) with I 0 either translated by some vector (ξ, 0), or replaced by 0. Moreover, If the maximum is also attained at some interior point, we get the same contradiction as before. Therefore, Hopf's lemma yields

This contradicts the definition of h.
Suppose now that Considering now a maximizing sequence (x n ) n∈N for u 1 − u 2 , then the limits (up to subsequences) ( u j , v j ) of the translations (u j (x + x n ), v j (x + x n , y)), which, once again, satisfy a system analogous to (7). The difference u 1 − u 2 attains its maximum ν μ h at the origin, whence

This contradicts sup R×R
The proof of the uniqueness is concluded. Let us pass to the limits at infinity. The limit v r ∞ (x, y) → 0 as y → +∞, uniformly with respect to x, follows from the negativity of f , exactly as in the above proof of Theorem 3.1. Consider now a diverging sequence (x n ) n∈N in R. The sequence of translations (u r ∞ (x +x n ), v r ∞ (x +x n , y)) converges locally uniformly (up to subsequences) towards a bounded, stationary solution ( u, v) to (7) with I 0 ≡ 0. We have seen above that the Liouville-type result holds for nonnegative stationary solutions of such system, hence, necessarily, ( u, v) ≡ (0, 0). This concludes the proof of the theorem.
We now turn to the set of results on the solution of the Cauchy problem (7)-(8).
Proof (Theorem 3.5) Since the initial datum (0, 0) is a subsolution to (7) which is not a solution, the comparison principle implies that the solution (u, v) is strictly increasing in t. It further implies that (u, v) is smaller than the supersolution constructed in the proof of Theorem 3.4, defined by (24). It follows that, as t → +∞, (u, v) converges locally uniformly to a positive, bounded, stationary solution, which necessarily coincides with (u r ∞ , v r ∞ ) due to Theorem 3.4. Proof (Theorem 3.6) In the case I 0 , T 0 ≡ 0, the result reduces to (Berestycki et al. 2013. Theorem 1.1), with the only difference that the initial datum there must not be identically equal to (0, 0) (otherwise the solution remains (0, 0) for all times). In that case, the limit state is simply (u r ∞ , v r ∞ ) ≡ (ν/μ, 1)v * . We call c T S I R the speed provided by (Berestycki et al. 2013, Theorem 1.1).
Take ε ∈ (0, c T S I R ) and consider a sequence (t n ) n∈N diverging to +∞ and a sequence (x n ) n∈N in R such that |x n | ≤ (c T S I R − ε)t n . If (x n ) n∈N is bounded, then the convergence of (u(x n ), v(x n , y)) towards the steady state follows from Theorem 3.5.
Suppose that (x n ) n∈N diverges (up to subsequences). By the strong maximum principle, the solution (u, v) is strictly larger than (0, 0) at, say, t = 1. Fitting a compactly supported datum below it, and applying the spreading result from (Berestycki et al. 2013, Theorem 1.1) to the solution of (7) with I 0 , T 0 ≡ 0, emerging from such datum, we infer by comparison that where we have also used the limit given by Theorem 3.4. On the other hand, by comparison, (u, v) ≤ (u r ∞ , v r ∞ ), for all t ≥ 0. Thus, the first limit in Theorem 3.6 is proved.
We turn to the second limit. We restrict to x > 0, the case x < 0 being obtained by a specular argument. Let us briefly remind how the asymptotic speed-named here c T S I R -is obtained in Berestycki et al. (2013). To start with, one linearizes (7) around (0, 0), i.e. replaces f (v) with f (0)v = α(R 0 − 1)v and set I 0 , T 0 to 0. One then looks plane wave solutions of the form ϕ(t, x), ψ(t, x, y) = e −a(x−ct) (1, γ e −by ), a, b, c, γ > 0. This leads to the following algebraic system for the unknowns (a, b, c) (the coefficient γ is easily computed from the exchange condition): The analysis of this system is somewhat involved, and we refer to Berestycki et al. (2013) for a detailed discussion. Following this analysis, one sees that the system admits a solution starting from a minimal value of c , called c T S I R . Hence, for c = c T S I R , the pair (ϕ, ψ) above is a solution to the linearized system and therefore, by the concavity of f , it is a supersolution to the original system (7) outside the supports of I 0 , T 0 . The second limit in Theorem 3.6 for x > 0 then follows by comparison with k(ϕ, ψ), with k sufficiently large so that k(ϕ, This is possible due to Theorem 3.6 and (25). These also entail that the function τ * (x) is locally bounded and satisfies (9). It is also clear that inf τ * > 0, because v identically vanishes at t = 0 and it is uniformly continuous by parabolic estimates. Assume (10) to be false. Namely, there exists a sequence (x n ) n∈N in R and a bounded sequence (y n ) n∈N in [0, +∞) such that one of the following situations occurs: Because τ * (x) is bounded from below away from zero and locally bounded from above, and I , T are positive for t > 0, we necessarily have that (τ * (x n )) n∈N diverges. Set from parabolic estimates, a subsequence of (S n , I n , T n ) n∈N (that we may assume without loss of generality to be the whole sequence) converges to an entire (i.e. for all t ∈ R) solution (S ∞ , I ∞ , T ∞ ) of the S I RT system (2). We may also assume (y n ) n∈N to converge to some y ∞ ≥ 0; by (26) we have that either I ∞ (0, 0, y ∞ ) = 0 or T ∞ (0, 0) = 0. In the latter case, we deduce from the last equation in (2) that T ∞ ≡ 0 for t ≤ 0, and then I ∞ ≡ 0 too; in the former case, the same conclusion follows from the strong maximum principle and Hopf's lemma, using the first and third equations in (2). Set u n (t, x) = u(τ * (x n ) + t, x n + x), v n (t, x, y) = v(τ * (x n ) + t, x n + x, y); as T n = ∂ t u n and I n = ∂ t v n , a subsequence (that we may once again to be the whole sequence) converges to a t-independent pair (u ∞ (x), v ∞ (x, y)), as their time derivatives converge to 0 as n → ∞. So, (u ∞ , v ∞ ) is a stationary solution of (7) with I 0 , T 0 ≡ 0. But we know from (Berestycki et al. 2013, Theorem 4.1) that the Liouvilleproperty holds for such system, that is, either by the definition of τ * (x).
We now prove the result about the rate of decay of the steady state.
Proof (Theorem 4.2) We derive these limits in the case x → +∞, the limits at −∞ being analogous. Upper bound. In order to simultaneously treat the cases R 0 ≤ 1 and R 0 > 1, we set in the former v * := 0. It follows that f (v * ) ≤ 0 for all R 0 > 0. Consider the pair ( u, v) given by (13) Namely, ( u, v) satisfies the linearized problem (12). By the concavity of f , we see that, for any h > 0, −ζ replaced by f (v * ), as correctly observed by the referee and therefore the pair ( ν μ v * + h u, v * + h v) is a supersolution to (7) outside the supports of T 0 , I 0 . We take h large enough so that, inside these supports, such pair is larger than the supersolution (u, v) used in the proof of Theorem 3.4, defined by (24). As a consequence, the pair is a generalized supersolution to (7). Hence, we can find a stationary solution between (0, 0) and such supersolution, which is necessarily positive and thus coincides with (u r ∞ , v r ∞ ) thanks to the Liouville result of Theorem 3.4. We have thereby shown that and therefore the desired upper bounds. Lower bound. Fix ζ > − f (v * ) ≥ 0. We consider now the pair ( u, v) from (15), with 0 < ε < 1 and (a, b, γ ) = (a ζ,ε * , b ζ,ε * , γ ζ,ε * ) unique positive solution of (16). The function v(x, y) vanishes at y = y ε := − log ε 2b ζ,ε * , it is bounded in the half-strip S ε := {x > 0, 0 < y < y ε }, and, by construction, satisfies there −dΔ v = −ζ v. Hence, because f (v * ) = 0 and ζ < − f (v * ), the function v * + h v is a subsolution to the second equation of (7) in S ε provided h > 0 is sufficiently small, depending on ζ . As a consequence, choosing h = h ζ small, we have that the pair is a subsolution to (7) in S ε . Up to replacing h ζ with a smaller quantity h ζ,ε if need be, we can also require and, in addition, that in the half-strip S ε , (u, v) is smaller than the supersolution (u, v) defined by (24). Let us consider the solution (u, v) of the evolution problem (7) having (u, v) as initial datum. We know from the Liouville-type result that (u, v) (u r ∞ , v r ∞ ) as t → +∞. It follows that, for all t > 0, and moreover v(t, x, y ε ) > 0 = v(x, y ε ) for x ≥ 0. We can therefore apply the comparison principle for the road-field system in the half-strip S ε and deduce that, there, (u, v) remains larger than (u, v) for all t > 0. Thus, we infer that Owing to the arbitrariness of ζ > − f (v * ) and 0 < ε < 1, one then gets the desired lower bound using (17) and noticing that y ε → +∞ as ε → 0.
We conclude with the comparison between the steady state for the models without and with the line.
Proof (Theorem 4.1) The existence of the set E + , of the form |x| ≥ ρ, y < h, directly follows from Theorems 5.1 and 4.2. Let us show the existence of E − .
Because I 0 (x, y) is an even function of y, the stationary solution of (4) v ∞ is even in y too, by uniqueness. Hence, ∂ y v ∞ (x, 0) = 0 for all x ∈ R. Thus, integrating the equation −dΔv ∞ = f (v ∞ ) + I 0 (x, y) on R × (0, +∞) yields On the other hand, integrating the equations (7) for (u r ∞ , v r ∞ ) shows 1

R×(0,+∞)
f (v r ∞ ) + I 0 (x, y) dx dy = d interpret this as saying that the epidemic does not propagate. On the contrary, when R 0 > 1, the steady state converges to some positive constant and the epidemic has propagated. Thus, the position of R 0 relative to 1 still governs the propagation or dying out of the epidemic. At this stage, the scenario is exactly the same as for the standard SIR model. 2. We compute an asymptotic speed of propagation for this model, that we call c T S I R , and compare it with c S I R , the speed of propagation for the classical SIR model with diffusion (1). We show that c T S I R > c S I R if D > 2d, where D and d are the diffusion coefficients on the road and in the rest of the territory, respectively. In the case D ≤ 2d, then, c T S I R = c S I R . Thus, the presence of a road with fast diffusion enhances the speed of propagation of the epidemic. 3. We show that the S I RT system is governed by four parameters: the basic reproduction number R 0 (from which the classical SIR speed of propagation is computed), the reduced transmission coefficientsμ andν, and the ratio D = D/d between the diffusion on the road and the diffusion in the field. We find that, even if R 0 is very close to 1, the diffusion on the road may trigger a wave of contamination spreading at high speed. Even though the growth of the infection at each location may be slow, the propagation along the road may be fast. This may lead to situations where an epidemic may seem dormant and thus innocuous while it spreads seeds far away bringing about outbreaks and clusters apparently unconnected to the regions with a significant prevalence. 4. We compare the total cumulative number of infected individuals per location, I tot (x), in the cases with and without the road. We show that, compared with the standard SIR model, the presence of the road increases I tot (x) in the range of x large, that is, far from the epicentre of the epidemic. This result is not intuitive a priori. Indeed, by the enhancement of the speed of spreading of the epidemic wave due to the road, it also follows that at any location the epidemic peak lasts less than without the road. Therefore, one might have thought that, moving faster, the total number of infected by the epidemic would go down. However, far from the epicentre, the contrary happens. We also prove that while the total number of infected is higher far away, there is also a region E − , presumably close to the epicentre, where the total number of infected I tot (x) for the model with the road is smaller than the corresponding one for the standard SIR model. It would be interesting to characterize such a set E − in some specific cases, establishing for instance whether it actually contains the epicentre of the epidemic. We leave this as an open question. The lower number of infected near the epicentre due to the presence of the road could be related to another phenomenon that we observe on the standard SIR model. In Proposition 5.2 above, we prove that the quantity I tot (x) evaluated at the epicentre of the epidemic x = 0 is a decreasing function of the diffusion coefficient d. This reflects the fact that a higher diffusion coefficient "scatters" more quickly the infected individuals far from the epicentre. The same mechanism could be at work for our model S I RT , where the road allows for more infected individuals to move away from the centre. 5. In a general manner, the S I RT -type of models we introduce here adds a compartment in the population dynamics and fills a gap at an intermediate scale. Indeed, for instance, most models presently used in monitoring the COVID-19 epidemic rely on the detailed analysis of two types of networks: the global network of air travel or the microscopic socio-economic or family and schools network. As our analysis shows, the intermediate network of roads and railways plays an important role in the spatial spreading of epidemics at levels of a region or a country.
This model and its generalizations open the way to many open problems. For instance, a similar study would have to be carried out when there is diffusion not only of the infected but also of the susceptibles. This would lead to a system: The model we have presented and analysed in this paper sheds light on the effect of a road within an environment of slow diffusion for the spreading of epidemics. It allows us to explain some observations and to uncover various effects. The simplified structure allows us to carry a fairly complete mathematical analysis. Yet, this model could lend itself to more practical developments. Indeed, involving a network of roads should yield more precise results. To take into account roads in a more realistic fashion in discrete models is an important perspective.