Spreading speeds of rabies with territorial and diffusing rabid foxes

A mathematical model is formulated for the fox rabies epidemic that swept through large areas of Europe during parts of the last century. Differently from other models, both territorial and diffusing rabid foxes are included, which leads to a system of partial differential, functional differential and differential-integral equations. The system is reduced to a scalar Volterra-Hammerstein integral equation to which the theory of spreading speeds pioneered by Aronson and Weinberger is applied. The spreading speed is given by an implicit formula which involves the space-time Laplace transform of the integral kernel. This formula can be exploited to find the dependence of the spreading speed on the model ingredients, in particular on those describing the interplay between diffusing and territorial rabid foxes and on the distribution of the latent period.


1.
Introduction. In the last century, a rabies epizootic swept through large parts of Europe. Differently from the present rabies situation in African and Asian countries with the important role of domestic and stray dogs [42,69,70,71,73], it was mainly carried by foxes. It started in Poland in 1939 or 1940, south of Gdansk, and moved westward reaching the rivers Elbe in 1950 and Rhine in 1960, Denmark in 1964, Belgium, Luxembourg, and Austria in 1966, Switzerland in 1967, France in 1968, and Holland in 1974 [30,59]. It also moved to the east, to Hungary, the former Czechoslovakia and the former Soviet Union [59]. Its speed of spread ranged from 30 to 60 km per year according to [59,60] and from 20 to 60 km per year according to [29]. There were important variations where the epidemic front retreated at times in certain areas and moved up to 100 km in a given direction in other areas in a single year [59]. A study performed in three areas in the state of Baden-Württemberg (Germany) from January 1963 to March 31, 1971, found that the center of the frontwave moved at about 27 km per year [9] while the mean distance of new cases ahead of the frontline within a month was approximately 4.8 km [9,33].
Since the rabies virus appears in the saliva one or two days before symptoms are evident in the infected fox [29, p.247], the latent period is slightly shorter than the incubation period (transmission by biting). Fox rabies has a relatively long latent period with highly variable duration (from 12 to 110 days [4]) and a relatively short infection period (from 3 to 10 days [4]). While the literature seems to agree on the mean length of the infectious period (which is ended in most cases by the death of the fox), 5 days, there seems to be disagreement on the mean length of the latent period: 35 days ( [60] and the references therein) and 28 days [33] and 28 to 30 days [4], and 25 and 26.5 days and one month according to various sources cited in [59]. For simplicity, by slight abuse of language, we talk about foxes in the latent period as incubating foxes.
The above-mentioned fox epizootic has motivated quite a few mathematical studies [28,35,38,54,60]. Like this one, they all assume that susceptible and incubating foxes have home-ranges (unless they are migrating juveniles looking for new homeranges [38]).
Differently from previous studies, our model assumes that some of the rabid foxes essentially behave like susceptible and incubating foxes and keep their homeranges, while the other rabid foxes loose the attachment to their home-ranges and disperse by diffusion (Section 2). We call the first ones territorial rabid foxes and the second ones diffusing (wandering [59]) rabid foxes. Mathematical models have so far assumed that either all rabid foxes are territorial [31,60] or all rabid foxes diffuse [28,35,38] (see [34] for a comparison). The radio-tracking in [5] supports the existence of territorial rabid foxes, while the distances of new rabies cases in [9, Table 2] seem to support the existence of diffusing rabid foxes. Since the latent period is relative long and has a very variable duration, we allow it to have an arbitrary length distribution [28].
Our model does not take into account that juvenile foxes wander in search of home-ranges of their own [38] and may keep wandering if they do not find any [30,Chap.4].
Since our model assumes that some rabid foxes are territorial and that the other rabid foxes diffuse and includes a latent period of arbitrarily distributed length, it is a system of partial differential, functional differential, and differential integral equations. See Section 2. The movement of territorial foxes within their homeranges is modeled by integrals over space. This excludes using typical differential equations techniques (shooting methods, center unstable manifolds, e.g.) for finding traveling wave solutions and the associated speeds [14,21,23,35,38,42,70].
To address the spreading speeds of the rabies epidemic, we reduce this complex model to a scalar Volterra Hammerstein integral equation (Section 3 and Section 5) and apply the concept of asymptotic speeds of spread (Section 7 and [53,57]). This concept had been originally developed by Aronson and Weinberger [7] for partial differential equations in population dynamics and then used by Aronson [6] for certain epidemic models and extended to larger classes of population and epidemic models by Diekmann [13] and Thieme [53]. The spreading speed satisfies a precise, though implicit, formula involving a space-time Laplace transform of the integral kernel of the Volterra-Hammerstein equation. This formula can be used to study the dependence of the spreading speed on the model both analytically and numerically (Section 8, 10 and 11). The spreading speed coincides with the minimum speed associated with wave-form solutions [12,13,39,57].
The paper proceeds as follows. In Section 2, the rabies model is derived with both territorial and diffusing rabid foxes. It is reduced to a scalar Volterra Hammerstein integral equation in Section 3 [1].
In Section 4, the unique existence of solutions is addressed for the scalar integral equation and the rabies model. In Section 5, the model system is specialized to spatial homogeneity with the fox habitat being all of R n . This sets the stage for the concept of spreading speeds. The solutions of the system on R n as space dominate those on a bounded domain, with the diffusion equation being complemented by Dirichlet boundary conditions, as shown in Section 6. Numerical simulations of the model on a bounded domain can be found in [1,2] for one space dimension and in [3] for two space dimensions.
In Section 7, the main features of the spreading speed for Volterra Hammerstein integral equations are extracted from the literature [53,57]. Some results establishing the monotone dependence of the spreading speed on the space-time Laplace transform of the integral kernel that cannot be found in [53,57] are presented and proved in Section 8.
In Section 9, the results on spreading speeds are applied to the rabies model. The monotone and limiting dependence of the spreading speed on the model parameters is established in Section 10.
In Section 11, the dependence of the spreading speed on the initial density of susceptible foxes is investigated if the sizes of home-ranges depend on this density in a decreasing way as in [60,Sec.7.2].
In Section 12, the special case is considered that all rabid foxes diffuse and that susceptible foxes and incubating foxes always stay at the centers of their territories. The spreading speeds are numerically compared to the minimum wave speeds found for models with the same features that, in addition, incorporate the turnover of the fox population [35,36,37]. The comparison suggests that neglecting the fox population turnover does not change the qualitative behavior of the spreading speed and changes the quantitative behavior only moderately.
Numerical simulations of our model, which display the spread of rabies in the fox population in space and time for latent periods of fixed length and of exponential distribution [1,2,3], show expanding waves the speeds of which agree well with the spreading speeds determined from our analytic formulas. This paper is dedicated to the memory of Hans F. Weinberger (1928Weinberger ( -2017, who is one of the founders (with Donald G. Aronson [6,7]) and main promoters of the theory of spreading speeds ( [47,64,65], e.g., to mention a small sample).

2.
A rabies model with diffusing and territorial rabid foxes. To study the spread of rabies in a spatially distributed fox population, we consider an open subset Ω of R n which represents the habitat of the foxes.
We assume that susceptible and incubating foxes have home-ranges. Home ranges of settled foxes may partly or largely overlap those of other foxes though parts of them may be defended from encroachment by other foxes [29, p.149]. Therefore, it seems to be not only a convenient but also justifiable approximation of reality to model the locations of susceptible and incubating foxes by the centers of their home-ranges in a spatial continuum.
As for rabid foxes, we assume that some of them essentially behave like susceptible and incubating foxes and keep their home-ranges, while the others loose the attachment to their home-ranges and disperse by diffusion. We call the first territorial rabid foxes and the second diffusing rabid foxes.
Mathematical models have so far assumed that either all rabid foxes are territorial [60] or all rabid foxes diffuse [28,35,38]. The radio-tracking in [5] supports the existence of territorial rabid foxes, while the distances of new rabies cases in [9, Table 2] seems to support the existence of diffusing foxes.
In our model, the natural turnover of the fox population is ignored: No foxes are born, and the only deaths are those of rabid foxes dying from rabies. For the endemic dynamics of rabies, this would not be an appropriate assumption. For calculating the spreading speed, however, it may be admissible because the spreading speed will mainly be determined by the epidemic outbreak at the rabies front-line, and the outbreak should be short enough for the turnover to be irrelevant. This point will be revisited later (Section 12).

2.1.
Modeling the susceptible and rabid stages. Let S(x, t) denote the density of susceptible foxes (which are all territorial) at time t whose home-ranges center at location x ∈ Ω. Further R 1 (x, t) are the diffusing rabid foxes at location x and time t and R 2 (x, t) the territorial rabid foxes at time t whose home-ranges center at location x. Finally, L(x, t) are the incubating foxes with home-range center at x at time t (they are all territorial), and T (x, t) is the transition rate at which foxes in the latent period at location x ∈ Ω and at time t ≥ 0 become infectious.
Let κ 1 (x, z) denote the rate at which a fox with home-range center x visits the location z ∈ Ω, where κ 1 : Ω 2 → R + is continuous. The rate at which a susceptible fox with home-range center x meets a territorial rabid fox with home-range center z is given by which means that it is the rate at which they both visit some common point y ∈ Ω (Cf. [60, (5.1)]). The model takes the form, (2.2) If Ω = R n , the partial differential equation for R 1 is accompanied by a boundary condition. Before we continue with the derivation of the model in Sections 2.2 and 2.3, let us explain the equations and the parameters in (2.2). p 1 is the chance of a rabid fox to diffuse, and p 2 the chance to be territorial, p j ≥ 0 and p 1 + p 2 = 1. The radio-tracking in [5] supports the existence of territorial rabid foxes, and the distances of new rabies cases in [9, Table 2] supports the existence of diffusing foxes. As with many infectious diseases, different individuals react to the same disease in diverse ways. To model the two types of behavior, it seems the easiest to assign them at random when the foxes transit from the latent to the infectious state.
β > 0 is the transmission coefficient which incorporates fox activity and the chance that the meeting of a susceptible and a rabid fox leads to the infection of the susceptible fox. B(x, t) is the incidence of the disease, i.e., the number of new cases per unit of time.
The first differential equation in (2.2) describes how the new disease cases occur: A susceptible (S) fox with home-range center x is infected by meeting a diffusing rabid (R 1 ) fox that is momentarily located at z at rate κ 1 (x, z) or a territorial rabid (R 2 ) fox with home-range center at z at rate κ 2 (x, z) from (2.1).
The second differential equation in (2.2) describes the dynamics of diffusing rabid foxes (R 1 ). D is the diffusion constant and ∇ 2 x = n i=1 ∂ 2 xi the Laplace operator. Diffusing rabid foxes die at per capita mortality rate ν 1 > 0, and new diffusing rabid foxes are recruited from the latent stage at rate p 1 T (x, t).
The third differential equation in (2.2) describes the dynamics of territorial rabid foxes (R 2 ). Territorial rabid foxes die at per capita mortality rate ν 2 > 0, and new territorial rabid foxes are recruited from the latent stage at rate p 2 T (x, t).
The continuous functions S • , R • j : Ω → R + , j = 1, 2, denote the initial values of susceptible and of diffusing and territorial rabid foxes, respectively. Remark 2.1. We have modeled the duration of the rabid, infectious, stage in the simplest possible way, namely as exponentially distributed, i.e., by a constant per capita rate of dying from the disease. The duration of the latent period will be modeled by an arbitrary distribution because this will provide very interesting insights (Section 10.2). We could have done something similar for the rabid stage and let the disease-death rate and even the diffusion rate of diffusing rabid foxes depend on the time since becoming infectious [28]. We do not do so because this would add a layer of complexity that may obscure the interplay of diffusing and territorial rabid foxes. Instead of separating the latent and the rabid period, we could have looked only at the class of infected foxes and their infection age [15,28]. We did not do this in order to use the information about the lengths of the latent period and the rabid period that is available.
Alternatively to the diffusion equation, the wandering rabid foxes could be modeled by a differential-integral equation This would make it possible to include studying the effects of fat-tailed kernels κ 3 which should lead to accelerating solutions as in [10,20]. While in a few areas, rabies has spread with considerable higher speed than in most others [59], it is not clear whether this is caused by fat kernels (the main example for this seems to be rapid plant migration) or rather by spatial heterogeneities which are not considered in this paper. Yet another cause of anomalous high spreading speeds is the interaction of different species [65]. A mathematical reason for not venturing into the area of fattailed kernels is that they do not appear to lead to formulas for the spreading speed that are amenable to a qualitative analysis of how the spreading speed depends on rabies-relevant model parameters.
2.2. The cumulative infective force. We continue the derivation of the model by integrating the differential equation for S, are the time-cumulated densities of diffusing and territorial foxes, respectively. We call u(x, t) the cumulative infectious force at x up to time t. We also integrate the differential equations for the rabid foxes, The initial conditions are C j (x, 0) = 0. (2.9) 2.3. Modeling the passage through the latent period. Recall that, while T (x, t) is the rate of exiting the latent period, −∂ t S(x, t) = B(x, t) is the rate of entering the latent period. Since, by assumption, no foxes are born and susceptible and incubating foxes do not die, where L • (x) = L(x, 0) are the foxes at location x that undergo the latent stage of the disease at the beginning. By the fundamental theorem of calculus, U (x, t) is the time-cumulative rate of transition from latency to infectivity at x taken from 0 to time t. We stratify the foxes in the latency period along infection age, where E(x, t, a) are the foxes in the latent period that were infected (successfully exposed) a time units ago and so have the infection-age a. The infected foxes in the latent period with infection-age a are given by [55,Sec.13 are the foxes at location x ∈ Ω that are in the latent period at time 0 and have infection age a, is the conditional probability of still being in the latent period at infection age a under the condition of having been in the latent period already at age a − t. We substitute (2.12) into (2.11), with the convention that Υ(a) Υ(s) = 0 if a ≥ s and Υ(a) = 0. Since Υ is decreasing, it is of bounded variation; we integrate by parts and use a Stieltjes integral [55, Sec.B.1], L 0 (x, t) are the densities of foxes at location x at time t ≥ 0 that have been in the latent period at time 0 and still are in the latent period at time t. We substitute these formulas into (2.10), The equation for U can be rewritten as Finally, by (2.4), with u being given by (2.5) and 2.4. Summary of the model. In summary, we obtain the following system of two nonlinear functional partial differential equations, with the auxiliary equations As a reminder, C j (x, t) are the time-cumulated densities of diffusing (j = 1) and territorial (j = 2) rabid foxes at location x ∈ Ω up to time t ≥ 0 and u(x, t) is the time-cumulated density of the infective force. U (x, t) is the time-cumulative rate of transition from latency to infectivity at x taken from time 0 to time t and U 0 (x, t) is its part due to foxes that were in the latent period at time 0. The nonlinearity F is given by (2.21) and the spatial integral kernel κ 2 obtained from κ 1 by (2.1). Υ is the sojourn function of the latent rabies stage. For spreading speeds and/or traveling wave solutions for (not rabies specific) epidemic models that combine diffusion and spatial integrals in a different way, see [63,72]. Remark 2.3. Recall that U 0 is given by (2.18) and notice from that formula that U 0 (x, t) is an increasing function of t ≥ 0, U 0 (x, 0) = 0.
By (2.4) and (2.6), we can retrieve the densities of susceptible foxes and of diffusing and territorial rabid foxes via It is meaningful to assume that aΥ(a) → 0 as a → ∞ and Υ is integrable on R + . Then Extreme, but commonly used examples of the sojourn function Υ are those describing an exponential distribution of the length of the latent period, Υ(a) = e −a/τ , a ≥ 0, or a fixed length of the latent period, Υ(a) = 1, 0 ≤ a < τ , and Υ(a) = 0, a > τ , with τ as the mean length. Alternative derivations of our model for these two special cases that are closer to the numerical procedures used to simulate them can be found in [1,2]. If the latent period has fixed length τ , the model equations (2.2) can simply be completed by setting T (x, t) = B(x, t − τ ) for t ≥ τ and letting T (x, t) be prescribed for 0 ≤ t < τ , e.g., by T (x, t) = 0, which means that there are no incubating foxes at time 0 [1,2,3]. But, as we mentioned in the introduction, the length of the latent period is highly variable. In [1,2,3], also a derivation is given for differentiable Υ as for more realistic Gamma or lognormally distributed lengths of the latent stage. For this paper, we have chosen a derivation that works if Υ is an arbitrary sojourn function encompassing all these special cases. See Sections 10.2 and 11.1 for more about the influence of Υ on the spreading speeds of the epidemic.
3. Reduction to a scalar integral equation for the cumulative infectious force. If Ω = R n , there is a fundamental solution associated with the differential operator ∂ t − D∇ 2 x + ν 1 . If Ω is a bounded domain of R n , we assume that the boundary is smooth enough and the boundary conditions are of the kind that there is a Green's function G : Ω 2 ×(0, ∞) → (0, ∞) such that nonhomogeneous equations can be solved by the usual variation of constants formula [22]. Then, from the first equation in (2.22), G(x, y, t − r)U 0 (y, r)drdy.

(3.2)
Notice that the second integral expression can be written as t 0 Ω G(x, y, s) U 0 (y, t− s)dsdy after a change of variables and that C 10 (x, t) inherits from U 0 that it is an increasing function of t ≥ 0 (Remark 2.3). Further, C 10 (x, 0) = 0 for all x ∈ R n .
We substitute the fourth equation of (2.22) into (3.1), We change the order of integration, After a substitution, We change the order of integration another time, Similarly, and We substitute (3.6) and (3.8) into (2.5). After changing the order of integration, we have the following result. Recall (2.4) and (2.21).
Notice that u 0 (x, t) inherits from C 10 and C 20 that it is an increasing continuous function of t ≥ 0, u 0 (x, 0) = 0 for all x ∈ Ω. u 0 is the part of the cumulative infectious force that comes from the foxes that were either rabid or in the latent period at the beginning. We call u 0 the epidemic driving function.

Existence and uniqueness of solutions. The nonlinearity
is analytic, strictly increasing, strictly concave, Lipschitz continuous, and linearly dominated, F (u)/u is strictly decreasing in u > 0 and converges to 0 as u → ∞. Then there exists a unique continuous solution u : Ω × R + of (3.11) that is bounded on Ω × [0, c] for any c > 0.
The proof uses a standard application of the contraction mapping theorem (Banach's fixed point theorem) with the complete metric spaces C b The solution can be obtain by iteration. Define a sequence of functions u n : Ω × R + by the recursive definition x ∈ Ω, t ≥ 0, n ∈ N.
Since F is increasing, it follows by induction that After a change of variables, is an increasing function of t ≥ for all x ∈ Ω, the u n (x, t) inherit this property by induction and so does the limit u(x, t).
5. Spatial homogeneity. We now assume that Ω = R n and that the territorial foxes move in their home-ranges in a spatially homogeneous way, where | · | denotes the Euclidean norm in R n andκ 1 : R + → R + is assumed to be continuous. This means that the rate at which a fox with home-range center x visits the location y depends on the Euclidean distance between x and y, It is well-known from the theory of partial differential equations that the fundamental solution associated with the differential operator ∂ t − D∇ 2 x + ν 1 is given by x, y ∈ R n , t > 0.
We assume that the initial density of susceptible foxes is constant, S • (x) = S • . Then (3.11) takes the form Here, for x ∈ R n and s ≥ 0, Further, for x ∈ R n and t ≥ 0, (5.7) 6. Comparison of solutions on R n and on bounded Ω. If one wants to see how solutions of (2.2) look like, one cannot solve this system on R n but rather solves it on a bounded domain Ω of R n with sufficiently smooth boundary and Dirichlet boundary condition. A Dirichlet condition represents a hostile absorbing boundary, alternatively a Neumann boundary conditions would represent a nonhostile reflecting boundary. It is not quite clear to us whether a seashore is a hostile or reflecting boundary for foxes. We choose Dirichlet boundary conditions for mathematical reasons because we can prove that the solutions on a bounded domain with Dirichlet boundary condition are dominated by the solutions of the homogeneous solutions on R n while we have no good idea of how the solutions are related if the Dirichlet condition is replaced by a Neumann boundary condition. Let ξ : Ω 2 × (0, ∞) → R + be the integral kernel given by (3.12), (3.7), and (3.9) with G being the Green's function associated with ∂ t − D∇ 2 x − ν 1 and Dirichlet boundary conditions, G(x, y, t) = 0 for x ∈ ∂Ω, y ∈ Ω, t ∈ (0, ∞).
Letξ : R n × (0, ∞) → R + be the spatially homogeneous integral kernel given by (5.6) and the other formulas in Section 5 using the fundamental solution Γ n . Theorem 6.1. Let S • > 0 be constant and u : Ω × R + → R + be the solution to (3.11) andũ : R n × R + → R + be the solution to (5.5).
Proof. By the maximum principle, the Green's function G associated with Dirichlet boundary condition lies below the fundamental solution on its domain of definition, for the susceptibles S on the bounded domain Ω and the susceptiblesS on R n , we have S(x, t) ≥S(x, t) for x ∈ Ω, t ≥ 0. So the epidemic model on a bounded domain Ω with Dirichlet boundary conditions shows a less severe epidemic outbreak than the epidemic model on R n . This suggests that the spread of the rabies modeled on Ω is not so fast as the spread of the disease modeled on R n . There is no rigorous derivation of this suggestion, but it is confirmed by the numerical simulations in [2] (one space dimension) and [3] (two space dimensions). 7. Spreading speeds for spatially homogeneous Volterra Hammerstein integral equations. From now on, we only consider the spatially homogeneous situation described in Section 5. For notational simplicity, we drop superscript ∼, Here ξ : R n × (0, ∞) → R + is continuous and We first exclude the case in which there is no spatial spread.
We also derive an asymptotic bound in the case ξ * > 1.
We repeat the proof in [53] to be clear about that no more assumptions are needed than those we have made above. We prove the two results simultaneously. Notice u − u 0 is bounded by ξ * and the right hand side of this inequality is finite. By Fatou's lemma, applied to (7.2), where u has been extended to R n ×(−∞, 0) by 0. As F is increasing and lim sup n→∞ u(x n − z, t n − s) ≤ u for all n ∈ N, s ≥ 0 and z ∈ R n , lim sup We now impose additional assumptions on the integral kernel ξ.
Definition 7.4. We define the space-time Laplace transform of ξ by (7.3)), the candidate for the spreading speed is defined by Remark 7.5. Since ξ is nonnegative, Ξ can be defined without the assumption (7.3) and Assumption 7.3 (b) as long as the integral is taken in the Lebesgue sense and one accepts that Ξ(c, λ) = ∞ is possible. (7.3) can be rephrased in terms of Ξ as Ξ(c, 0) = Ξ(0, 0) < ∞, and Assumption 7.
We introduce a suitable restriction of the driving functions u 0 in (7.1).
Definition 7.6. The epidemic driving function u 0 : for all x ∈ R n and t ≥ 0.
Theorem 7.7. Let the Assumptions 7.3 be satisfied and ξ * > 1. Let u * > 0 be the unique solution of u * = ξ * F (u * ). Then c * ∈ (0, ∞) and c * is the spreading speed (aka asymptotic speed of spread) in the following sense: (a) Let u 0 be admissible and c > c * . Then, for any > 0, there exists some t > 0 such that u(x, t) ≤ whenever u is a solution of (7.1) and t ≥ t , |x| ≥ ct. (b) Let 0 < c < c * and let δ > 0 and t 2 > t 1 > 0 be such that u 0 (x, t) ≥ δ whenever t ∈ (t 1 , t 2 ) and x ∈ R n with |x| ≤ δ. Then, for any ∈ (0, u * ), there exists some t > 0 such that u(x, t) ≥ u * − whenever u is a solution of (7.1) and t ≥ t , |x| ≤ ct.
The next result suggests that c * is close to 0 if Ξ(0, 0) is slightly larger than 1.
Theorem 8.3. Let ξ j : R n × R + → R + be a family of functions satisfying the Assumptions 7.3 and let ξ : R n × R + → R + satisfy the Assumptions 7.3.
A series expansion of c * (α) in = √ ln α is derived in [60, Sec.6.1] using a so-called cumulant generating function. 9. Spreading speeds for the rabies model. How do the results of the previous section play out for the rabies model?
The integral kernel ξ in (7.1) is the sum of two integral kernels which, in turn, are temporal, spatial, or space-time convolutions of other functions. The space-time Laplace transform (7.4) converts these convolutions into products of appropriate Laplace transforms [57,Prop.4  Here E 1 (c, λ) is the space-time Laplace transform of η 1 andκ 1 is a spatial Laplace transform of κ 1 andΥ a temporal Laplace transform of Υ, By [57,Prop.4.2] or Lemma A.5 and by (5.6), if the denominator is positive; otherwise this integral is infinite. By (5.6), We substitute (9.3) and (9.4) into (9.1), νj is the average time a rabid fox has available for infecting others before diseaseinflicted death if it is territorial or diffusing, respectively. The weighted average of these, p1 ν1 + p2 ν2 is the time available for a typical rabid fox. S • is the density of susceptible foxes available to be infected, and β is the transmission rate. So ξ * can be interpreted as basic reproduction number of the disease and obtains the alternative notation R 0 .
We check Assumption 7.3 (a), (c). The isotropy of ξ i follows from Proposition A.4 and the fact that the ξ are convolutions of isotropic functions.

9.2.
Conditions for the admissibility of the epidemic driving function. Let R 0 > 1. The epidemic driving function u 0 in (5.7) is admissible (Definition 7.6) if, for any c, λ > 0 with Ξ(c, λ) < 1, there exists some γ ∈ (0, ∞) such that Here x · v is the Euclidean inner product.
So we make the following definition.
Definition 9.4. The initial data R • i , i = 1, 2, of rabid foxes and L • of foxes in the latent period are admissible if for any λ > 0 there exists some γ > 0 such that for all x ∈ R n . 9.3. Spreading speed and solution behavior. Theorem 7.1 and Theorem 7.2 do not need any translation to the rabies model except that now ξ * = R 0 is the basic rabies reproduction number.
Theorem 9.5. Let the Assumptions 9.2 be satisfied and R 0 > 1. Let u * > 0 be the unique solution of u * = R 0 F (u * ). Then c * ∈ (0, ∞) and c * is the spreading speed (aka asymptotic speed of spread) in the following sense: (a) Let the initial data R • 1 , R • 2 and L • be admissible and c > c * . Then, for any > 0, there exists some t > 0 such that u(x, t) ≤ whenever u is a solution of (7.1) and t ≥ t , |x| ≥ ct. (b) Let 0 < c < c * and let δ > 0 and be such that R Then, for any ∈ (0, u * ), there exists some t > 0 such that u(x, t) ≥ u * − whenever u is a solution of (7.1) and t ≥ t , |x| ≤ ct.

A final specialization for the numerical calculation of the spreading speed.
For the numerical computation of c * , we assume that the movements of territorial foxes about the center of their home-range are normally distributed, x ∈ R n , (9.12) where | · | is the Euclidean norm on R n , b > 0, and Γ n is the fundamental solutions associated with the differential operator ∂ t −∆ x for n space dimensions. By Lemma A.5,κ 1 (λ) = e bλ 2 = R n e −λxj κ 1 (|x|)dx, j = 1, . . . , n. (9.13) By (9.2), for the variance of movement in any direction, n R n |x| 2 κ 1 (|x|)dx, j = 1, . . . , n. (9.14) We have written the normal distribution in the somewhat unusual form (9.12) to have the formula for its Laplace transform as simple as possible. By (9.13) and (9.5), If R 0 > 1, then c * > 0 and the spreading speed satisfies the system [57,60] Ξ(c * , λ) = 1, ∂ ∂λ Ξ(c * , λ) = 0. (9.16) c * and λ are uniquely determined and can numerically computed using Newton's method as it is done to obtain some values for the spreading speed in Table 1, Figures 11.1, 11.2, 11.3, 13.1 and in [1,2,3]. Differently from [32,60], our qualitative analysis directly uses (7.6) rather than (9.16) in order to obtain lower and upper bounds and limit information for the spreading speed rather than approximations.
(9.12) is a special case of So σ 2 is some relative measure of how far the location of a fox varies from the center of its home-range in each direction. We compare (9.12) and (9.17) and find the relation between the two scaling parameters. By (9.2), and Ξ 0 is the space-time Laplace transform associated with a kernel resulting from a model where all territorial foxes stay a the center of their home-ranges and the diffusion rate is one squared space unit per unit of time. Letc 0 be the spreading speed associated with this model. By (7.6), c * ≥ c 0 =c 0 D 1/2 . Proposition 10.2. c * ≥c 0 D 1/2 wherec 0 is the spreading speed associated with the model where all territorial foxes always stay at the centers of their home-ranges and the diffusion rate of diffusing rabid foxes is one square kilometer per unit of time.
Notice that, for all k ∈ N, ThenΞ Since c * k are the spreading speeds associated with Ξ k ,c k = D −1/2 k c * k are the spreading speeds associated withΞ k .
The relation becomes more complicated if a class of sojourn functions depends on two-parameters, e.g., if they are are Gamma distributed or lognormally distributed. A Gamma distribution can be reparameterized over its mean and its variance, and its Laplace transform has the form where τ is the mean and σ the standard deviation (in our case the mean length of the latent period and the standard deviation of its distribution). Υ and, by Corollary 8.2, the associated spreading speeds are decreasing functions of the mean if the standard deviation is kept fixed and increasing functions of the standard deviation if the mean is kept fixed [25, Sec.6].
We were not able to prove something similar for the lognormal distribution (Section B) which has been fitted to the latent periods of quite a few infectious diseases [45,46].
If the mean is kept fixed and the variance tends to zero, the Gamma distributed sojourn functions approach the one for fixed length of the latent period which exemplifies that the spreading speed associated with fixed length is the smallest among the spreading speeds associated with length distributions of same mean length. This is a general fact, which is proved in the upcoming Section 10.2.1.
10.2.1. The special role of latent periods with fixed length. Among all arbitrary length distributions of the latent period with mean length τ , the distribution with fixed length τ is associated with the smallest spreading speed. This is a consequence of Corollary 8.2 as follows: Since the exponential function is convex, by Jensen's inequality [43,Thm Let Ξ andΞ be the space-time Laplace transforms with Υ andΥ, respectively (see (9.5), and c * andc * the associated spreading speeds (7.6). Then, Ξ(c, λ) ≥Ξ(c, λ) for all c, λ > 0. By Corollary 8.2, c * ≥c * . This means that the spatial spread is sped up if the length of the latent period has some variation rather than no variation at all. Whether the spreading speed is an increasing function of the variation of the length of the latent period seems an open problem. As we have seen in the previous section, this is the case for Gamma distributed lengths, but we could neither prove it in general nor find a counterexample.
10.2.2. Latent periods of positive minimum length. The latent period of fox rabies appears to have a minimum length of about 12 days [4]. This means that the sojourn function Υ has the form The mean length of the latent period is given by From (10.7) and Corollary 8.2 we obtain the following result.
Theorem 10.4. The spreading speed depends in a decreasing way on the minimum length τ of the latent period and the function Υ in (10.6).
In summary, we have the following result. 2 , the spreading speed decreases as a function of the proportion p of diffusing rabid foxes provided that R 0 > 1 is sufficiently close to 1.
2 , the spreading speed increases as a function of the proportion p of diffusing rabid foxes provided that R 0 > 1 is sufficiently close to 1.
Suppose the statement is false. Then there exists a sequence (R n ) in (1, ∞) with R n → 1 such that, for each n ∈ N, the spreading speed associated with (10.8) and R n replacing R 0 is not decreasing as a function of p. More precisely, this means the following.
(b) The case D ν > σ 2 1 2 is proved similarly using Lemma 10.5 (b). We consider the special case that the movements of territorial foxes about the center of their home-range are normally distributed, i.e., x ∈ R n . (10.14) Corollary 10.7. If bν > D, the spreading speed decreases as a function of the proportion p 1 of diffusing rabid foxes provided that R 0 > 1 is sufficiently close to 1. If bν < D, the spreading speed increases as a function of the proportion p 1 of diffusing rabid foxes provided that R 0 > 1 is sufficiently close to 1.
This result is remarkable because the right hand side of the inequality is positive for p = 0, while, if all rabid foxes are territorial, the spreading speed tends to 0 if the mean distance of foxes from their home-ranges tends to 0. See Section 10.5. Such a (discontinuity) phenomenon has already been observed in [32] for a population model that corresponds to our epidemic model without the latent period. Because of the complications raised by the inclusion of the latent period, we prefer to give a lower bound for the spreading speed rather than an approximate expression as in [32, (23.21)]. Notice that the lower bound in Theorem 10.8 holds for all p ∈ (0, 1], not just for small values of p.
11. Dependence of the spreading speed on the initial density of susceptible foxes if it influences home-range size. As mentioned in Theorem 10.1, the spreading speed is an increasing function of the initial density of susceptible foxes, S • , provided that all other parameters are kept fixed.
A study in three areas in the state of Baden-Württemberg of Germany finds that the mean distance of new rabies cases ahead of the monthly determined rabies frontline very slightly decreases if the hunting indicator of the fox density (foxes shot per km 2 per year) increases [9, Table 3]. The likely explanation that reconciles our mathematical result and this observation is that the home-range size is not independent of fox density but depends on it in a decreasing fashion as observed in [44] and assumed in [60,Sec.7.2].
In the following we show that, if that is the case, the spreading speed is still an increasing function of S • if S • is slightly larger than the critical density S at which the basic rabies reproduction number is one. So the decreasing dependence can only occur if S • is sufficiently larger than S .
Proof. Suppose not. Then there exists a sequence (S • j ) in (S , ∞) such that S • j → S as j → ∞, and there is some c > 0 such that c * j > c for all j ∈ N where c * j are the spreading speeds associated with S • j . Let ξ j be the associated integral kernels and Ξ j their space-time Laplace transforms given by (11.4) with S • j replacing S • . We apply Theorem 8.3. Notice that Ξ j (0, 0) = βS • j p1 ν1 + p2 ν2 → 1 as j → ∞. The other assumptions of Theorem 8.3 follow from (9.15) which implies that c * j → 0 as j → ∞, a contradiction.
Theorem 11.2. Under the assumptions of this section and 0 ≤ p 1 ≤ 1, there exists some S > S such that c * is an increasing function of S • ∈ (S , S ).
The next result gives a condition under which the behavior of the spreading speed as a function of S • is not compatible with the observations in Baden-Württemberg. Proof. We substitute λ = r cφ(S • ) , r ≥ 0, in (7.6) with Ξ given by (11.4), r) is convex in r > 0. Sinceκ is increasing (by the analog of Lemma 9.1),Ξ(c, r) is decreasing in c > 0. This is sufficient to make the proofs of Proposition 8.1 and Corollary 8.2 work withΞ replacing Ξ. It follows from our assumptions thatΞ is an increasing function of S • . So c * is an increasing function of S • . Example 11.4. Let φ(S • ) = √ S • and let the length of the latent period be exponentially distributed, Υ(a) = e −γa , a ≥ 0, with a constant γ > 0. Then c * is an increasing function of S • .
Proof. In this case, is an increasing function of S • . Apply Theorem 11.3.
All these results hold if there are only territorial rabid foxes or only diffusing rabid foxes or both territorial and diffusing rabid foxes. In case that there are territorial rabid foxes only, Example 11.4 shows that formula (7.8) and Figure 7 in [60] and the figure in Box 23.5 of [32] do not hold for exponentially distributed latent and infectious periods because they rely on the assumption that the lengths of these periods are highly concentrated about their mean lengths. Fig. 11.1 illustrates that c * is an increasing function of S • when the lengths of the latent and infectious periods are exponentially distributed and the home-range size is normally distributed with variance 2b = ω/S • (recall (9.14) and see Section 11.2 for more details). This figure is obtained by numerically solving system (9.16) with Ξ(c, λ) = e (ω/S • )λ 2 ν + λc βS • γ γ + cλ , c, λ > 0. (11.8) See the caption of Figure 11.1 for the numerical values of the parameters and Section 11.2 for how they were determined for ω and β. 1/γ and 1/ν are the mean lengths of the exponentially distributed latent and infectious periods.
In the next subsection, we will see how this picture changes if the latent period is not exponentially distributed and all rabid foxes are territorial. Compare Fig. 11.1 to Fig. 11.2 (a), where all rabid foxes are territorial and the length of the latent period is fixed. Notice the similarity of Fig. 11.2 (a)  11.1. Territorial rabid foxes only. Let us, for a moment, consider the case that all rabid foxes are territorial, p 2 = 1 and p 1 = 0. We look for scenarios in which c * → 0 as S • → ∞.
The space-time Laplace transform, (11.4), takes the form Suppose that c * does not converge to 0 as S • → ∞. Then there exists a sequence (S • j ) in (0, ∞) and some c > 0 such that S • j → ∞ as j → ∞ and 0 < c < c * j for all j ∈ N, where c * j are the spreading speeds associated with S • j . By (7.6), with Ξ j the Laplace transforms where S • j replaces S • , , j ∈ N. (11.11) We obtain a contradiction if we make assumption (11.12) in the next theorem, and we obtain the following result.
Theorem 11.5. Assume Then, under the assumptions of this section, c * → 0 as S • → ∞.
Example 11.6. In [60, (7.6)], it is assumed that the variance of fox movements about the home-range center in every direction is reciprocally proportional to the fox density S • , i.e., there is some constant (11.13) with constants , α > 0. By Theorem 11.5, c * → 0 as S • → ∞ if andκ(s) < ∞ for all s > 0. (11.14) is satisfied if α ≥ 1 or if the latent period has a positive minimum length (see (10.7)) or is lognormally distributed (Lemma B.1).
If the length of the latent period is Gamma distributed with mean length τ and standard deviation σ, then (11.14) is satisfied if (τ /σ) 2 > 1 α − 1. Recall (10.3). If α = 1/2 as in [60, (7.6)], (11.12) holds if the length is Gamma distributed with the mean being strictly larger than the standard deviation; it does not hold if the length is exponentially distributed, in which case c * is a strictly increasing function of S • (Example 11.4).
See [19,62] and the references therein for other epidemic and endemic models in which the systems' behavior changes when exponential stage-length distributions are replaced by more realistic ones.
We are not aware of any fits of standard probability distributions to fox rabies data for the length of the latent period. In [18], Weibull, Gamma and lognormal distributions are fitted to times from infection to disease death for tiger salamanders infected by ranaviruses. The lognormal distributions give a better fit than the Gamma distributions, which in turn yield a better fit than the Weibull distributions. For the best fitting distribution in each class, the mean is substantially larger than the standard deviation.
11.2. Territorial and diffusing rabid foxes. If some (and possibly all) rabid foxes diffuse, we have the following result for large initial densities of susceptible foxes. Proof. Suppose not. Then there exists a sequence (S • j ) in (0, ∞) such that S • j → ∞ as j → ∞, and there exists some c > 0 such that c * j < c for all j ∈ N, where c * j are the spreading speeds associated with S • j . Let Ξ j be the space-time Laplace transforms associated with S • j as given by (11.4). By (7.6), for each j ∈ N there exists some λ j > 0 such that Ξ j (c, λ j ) < 1. Since p 1 > 0 andκ 1 (λ) ≥ 1 (Lemma 9.1), by (11.4) the sequence (z j ) given by j is bounded and ν 2 + λ j c − Dλ 2 j > 0 for all j ∈ N. This implies that (λ j ) is bounded and (z j ) is unbounded because S • j → ∞ as j → ∞, a contradiction. If ν 1 = ν 2 , the uniformity of the divergence follows from Theorem 10.8.
For the rest of this section we assume that there are both diffusing and territorial rabid foxes. By combining Theorem 11.7 with the results in Section 11.1 and with Theorem 11.2, if p 1 > 0 is small enough, one expects that c * increases as long as S • is slightly larger than its critical value, then decreases for intermediate values of S • and finally tends to infinity as S • → ∞. This is confirmed by the numerical calculations shown in Fig. 11.2 (b), (c) and Fig. 11.3 (b), (c).
For these calculations, we assume that the movements of a territorial fox about the center of its territory are normally distributed and that the per capita rabiesdeath rates for territorial and diffusion rabid foxes are the same, ν 1 = ν 2 = ν.
We follow [60] and, recalling (9.14), we assume that 2b = ω/S • , (11.15) where ω = 5.3 is a dimensionless parameter taken from [27]. By (9.15), A direct estimate of the transmission coefficient β is difficult to come by. Let S be the threshold density of susceptible foxes at which the basic rabies reproduction number R 0 equals 1. By (9.6), For S , a range from 0.25 to 1.0 foxes km −2 is recorded from the literature by [60], while a range from 0.2 to 1.0 foxes km −2 is recorded by [36]. We follow [60] in choosing 18) an estimate taken from [51]. We also choose the average lengths of the latent and infectious periods to be very similar to those in [60, Sec.7.2], 33.44 and 5 days, respectively.
The behaviors of c * as a function of S • in Figure 11.2, (a)-(d), and in Figure 11 density (foxes shot per km 2 per year) increases [9, Table 3] and has an intermediate range. Notice that Figure 11.2(a) is independent of the choice of the diffusion constant D because all rabid foxes are territorial. However, in Figure 11.2, where the diffusion constant is D = 200 km 2 /year, as in the speed calculations in some pure diffusion rabies models [35,36,37] In order to have an educated guess about the impact of population turnover, we look at the special case of our model with the same assumptions, including that susceptible and incubating foxes stay at the center of their home-ranges all the time. (9.5) then takes the form In Table 1 we compare spreading speeds that have been determined by numerically solving system (9.16) with the minimum wave speeds calculated in [35,Sec.20.4] [36,37]. The results agree qualitatively and are not too different quantitatively. A non-spatial endemic model with non-seasonal births is compared to one with birth pulses in [40], and though there are differences in the solution behavior they are not too pronounced. This encourages us to believe that the qualitative behavior of the spreading speed is not affected by the omission of population turnover and that the quantitative results contain useful information as long as they are seen as approximations.
13. Discussion. Differently from other fox rabies models, the model presented here includes both territorial and diffusing rabid foxes. Territorial rabid foxes, like susceptible and incubating foxes, have fixed home-ranges through which they roam. Diffusing rabid foxes have lost the attachment to their home-ranges and randomly wander around. Mathematically, diffusing rabid foxes are modeled by a diffusion equation while an integral kernel describes the rate at which territorial foxes visit the various locations in their home-ranges.
The integral kernels preclude that certain standard differential equation techniques to study waveform solutions and their minimum speed (shooting methods, e.g.) can be used. Instead, the concept of spreading speed (aka asymptotic speed of spread) is employed which uses space-time Laplace transforms of the integral kernel. The spreading speed concept works best in this context if the rabies system can be reduced to a single integral equation (of Volterra Hammerstein type). To make this possible, the natural turnover of the fox population (births and non-rabies deaths) is ignored and density-dependent incidence is assumed. Comparisons in the special case that all rabid foxes are diffusing and susceptible and incubating foxes stay at the centers of their home-ranges show that the omission of the natural turnover hardly affects the qualitative features of the spreading speed and only moderately changes the quantitative features (Table 1). Frequency-dependent and other, more general, incidences and multiple groups have been handled in epidemic (not rabies specific) models by looking for traveling wave solutions using fixed point arguments [16,48,58,67]. It is not clear to us whether this approach would allow including the natural turnover of fox populations. In general (see [49,68] for exceptions), it is difficult to relate the existence of traveling waves to the actual spread of the epidemic while the concept of a spreading speed gives some mathematical sense to it (Theorem 9.5).
Our model has been tailored for a specific infectious agent, the rabies virus, and a specific host, the red fox, that was the carrier of the European rabies epizootic in the second half of the 20th century. It therefore does not exactly fit another infectious disease. But the modeling principles and techniques and the concept of spreading speed are very flexible and can be and have been used for the spatial spread of other infections like Salmonella in an industrial hen house [8] and of viral infections in bacteria [24,25,26].
Since the inclusion of both territorial and diffusing rabid foxes is the special feature of this model, it is of particular interest how the ratio between diffusing and territorial rabid foxes affects the spreading speed. If the mean length of the infectious period is the same for both types, L = 1/ν 1 = 1/ν 2 , and the movements of territorial foxes about the center of their home-ranges are normally distributed with variance 2b in every direction, this depends on the relation between b/L and D where D is the diffusion rate of diffusing rabid foxes. If the basic reproduction number R 0 is slightly larger than 1, then the spreading speed is an increasing function of the proportion of diffusing foxes if D − b/L > 0 and a decreasing function if b/L − D > 0 (Theorem 10.7). If R 0 is substantially larger than 1, numerical calculations ( Fig. 13.1) show that these monotonicities still hold if the differences above are sufficiently large.
If diffusing rabid foxes are added to a model with territorial rabid foxes, the following remarkable discontinuity phenomenon occurs, a population dynamics analog of which has been observed in [32,Sec.23.6]. As long as all rabid foxes are territorial, the spreading speed tends to 0 as the mean distance of foxes from their home-range centers tends to 0 (see Section 10.5). If there are diffusing rabid foxes, however small their proportion p > 0 among all rabid foxes may be, the spreading speed remains bounded away from zero as the mean distance of territorial foxes from their home-range centers tends to 0 and p → 0 (Theorem 10.8).
A study in three areas in the state of Baden-Württemberg of Germany finds that the mean distance of new rabies cases ahead of the monthly determined rabies frontline very slightly decreases if the hunting indicator of the fox density (foxes shot per km 2 per year) increases [9, Table 3]. According to our mathematical model, c * is an increasing function of the initial density S • of susceptible foxes when all other parameters are kept fixed (Theorem 10.1). The likely explanation that reconciles observation and model is that the home-range size is not independent of fox density but depends on it in a decreasing fashion as observed in [44] and assumed in [60,Sec.7.2]. If this is assumed and there are territorial rabid foxes only and if the latent period has a realistic (in particular not an exponential) length distribution, it can be shown that c * tends to 0 as the initial density of susceptible foxes tends to infinity S • ([60, Sec.7.2] and Theorem 11.5). If diffusing rabid foxes are added, numerical calculations suggest that c * depends in a decreasing fashion on intermediate values of S • unless the proportion of rabid foxes that diffuse or their diffusion rate are too large (Figures 11.2 and 11.3).
Numerical simulations of our model, which display the spread of rabies in the fox population in space and time for latent periods of fixed length and of exponential distribution, show expanding waves the speeds of which agree well with the spreading speeds determined from our analytic formulas [1,2,3] though the first are somewhat less than the second because of the zero Dirichlet boundary conditions that are imposed for the numerical calculations as explained in Section 6.
Lemma A.1. Let x, y ∈ R n such that |x| = |y| > 0. Then there exists an orthogonal real matrix A such that x = Ay.
Let z be the unit vector (1, 0, . . . , 0) * , i.e., written as a column vector. By the Gram-Schmid procedure, x can be expanded to an orthonormal basis of R n . Let A be the matrix with the basis vectors forming the columns and x the first column. Then Az = x.
By the same argument, there exists an orthogonal real matrix y = Bz with z = (1, 0, . . . , 0) * . Then z = B −1 y = B * y, with B * denoting the transposed real matrix. So x = AB * y with the orthogonal real matrix AB * .
Definition A.2. A function ψ : R n → R is called isotropic if ψ(x) = ψ(y) for all x, y ∈ R n with |x| = |y|.
A function ξ : R n × R + → R is called isotropic, if ξ(·, t) is isotropic for all t ≥ 0.
Lemma A.3. A function ψ : R n → R is isotropic if ψ(Ax) = ψ(x) for all x ∈ R n and all orthogonal real matrices A.
Proof. ⇒: Recall that |Ax| = |x| for all x ∈ R n and all real orthogonal matrices A. ⇐: Use Lemma A.1.