Understanding rabies persistence in low-density fox populations

ABSTRACT

Arctic fox (Vulpes lagopus) and its tundra habitat are a unique system for the study of rabies virus epidemics. Contrary to theoretical calculations reporting a critical density (K_T) of approximately 1 fox/km² for rabies endemicity, arctic rabies persists at densities below this. The calculation of K_T = 1 fox/km² assumes uniform fox density across the landscape and unrestricted mixing between susceptible and infected foxes. We hypothesize that spatial heterogeneity arising from resource distribution or social structure may result in regions where rabies is endemic, even though average fox densities at the regional or landscape-level are below K_T. To expand upon the limited body of research surrounding arctic rabies persistence, we examine arctic rabies via a two-patch structure. We find that arctic rabies can persist in heterogeneous landscapes where the mean carrying capacity is below the threshold carrying capacity required for endemicity in homogeneous landscapes. Rabies endemicity in low-carrying capacity regions within heterogeneous landscapes is further facilitated by high transmission rates, potentially due to non-breeding foxes (i.e. floaters), and when between-patch movement is restricted to latent and infected foxes. Our results suggest that rabies may persist in heterogeneous landscapes when the mean carrying capacity is as low as 0.25 foxes/km².

Résumé

Le renard arctique (Vulpes lagopus) et son habitat toundrique constituent un système unique pour l’étude des épidémies du virus de la rage. Contrairement aux calculs théoriques rapportant une densité critique (K_T) d’environ 1 renard/km² pour l’endémicité de la rage, la rage arctique persiste à des densités moindres. Le calcul de K_T = 1 renard/km² suppose une densité uniforme de renards dans le paysage et un libre mélange des renards susceptibles et infectés. Nous posons l’hypothèse que l’hétérogénéité spatiale due à la répartition des ressources ou à la structure sociale pourrait résulter en régions où la rage est endémique même si les densités moyennes de renards sont inférieures à K_T à l’échelle régionale ou du paysage. Pour ajouter aux connaissances limitées sur la persistance de la rage dans l’Arctique, nous examinons la rage arctique selon une structure bi-parcellaire. Nos résultats montrent que la rage arctique peut persister dans des paysages hétérogènes quand la capacité de charge moyenne est sous le seuil d’endémicité des habitats homogènes. L’endémicité de la rage dans les régions à faible capacité de charge dans des paysages hétérogènes est par ailleurs facilitée par des taux de transmission élevés, potentiellement attribuables aux renards non-reproducteurs et à la restriction des mouvements inter-parcellaires aux renards latents et infectés. Nos résultats suggèrent que la rage peut persister dans des paysages hétérogènes jusqu’à une capacité de charge moyenne aussi faible que 0,25 renard/km².

KEYWORDS:

• Rabies
• metapopulation
• endemicity
• source-sink
• arctic fox

Mots clés:

• Rage
• métapopulation
• endémicité
• source-puits
• renard arctique

Introduction

Thresholds for disease persistence derived from models assuming homogeneous mixing, and without spatial structure identify critical densities necessary for disease spread (Lloyd-Smith et al. 2005; Keeling and Rohani 2008). However, those commonly used models can predict densities not found in nature. For instance, Anderson et al. (1981) predicts rabies persistence in Europe when densities of its regional host, red fox (Vulpes vulpes), are as low as 1 fox/km², yet rabies persists in the tundra when densities of arctic fox (Vulpes lagopus) are well below this threshold (Simon et al. 2019). One way such a contradiction could exist is that environmental heterogeneity may lower threshold host densities for disease persistence. In this scenario, connectivity between populations can facilitate ‘rescue effects’, preventing disease extirpation that would otherwise occur in an equivalent homogenous environment (Hess 1996; Hagenaars et al. 2004). Here, we hypothesize that spatial structure reduces the threshold fox density for rabies persistence, to be consistent with reported arctic fox densities and observed rabies endemicity.

Previous results from metapopulation theory offer additional insights into how spatial structure will affect predicted rabies dynamics. For disease persistence in a metapopulation, at least one of the subpopulations must be able to support the disease independently (Hethcote 1976). For a two-patch model, where only one patch can support the disease independently, either the disease dies out or persists in both patches, depending on the between-patch connectivity (Wang and Zhao 2004). More generally, pairing two populations with different qualitative or quantitative behavior can lead to the emergence of new population dynamics and novel spatial patterns (Pedersen et al. 2016). Spatial structure and population connectivity affect persistence and threshold values (Bolker and Grenfell 1995; Wang and Mulone 2003; Wang and Zhao 2004), and once a disease becomes endemic, the intensity of outbreaks and fade-out dynamics can be influenced by the spatial arrangement (Post et al. 1983; Sattenspiel 1987). Considering arctic rabies within a metapopulation structure complements previous research that applies metapopulation theory to communicable diseases such as hepatitis A, gonorrhea, and HIV, which are all characterized by non-homogenous mixing and infections (Lajmanovich and Yorke 1976; Sattenspiel 1987; Jacquez et al. 1988). In these instances, the inclusion of spatial heterogeneity is justified in that humans, as a host, violate homogenous mixing assumptions, since most people use the same travel routes, work in the same areas, and reside in the same locations daily.

Of the various host species and their habitable environments, the arctic fox and its tundra habitat have proved to be a unique system for the rabies virus. Contrary to the well-documented disease density threshold (K_T) of approximately 1 fox/km² (Anderson et al. 1981), arctic rabies persists endemically at densities below this. Indeed, landscape-level arctic fox densities rarely exceed 0.3 breeders/km² (Angerbjörn et al. 1999; Eide et al. 2004; Simon et al. 2019). When these numbers are scaled into total number of foxes per square kilometer, by including non-breeding foxes and juvenile survival, Arctic densities are still unlikely to exceed an average of 1 fox/km². While most Arctic regions have average densities below K_T, some local areas have carrying capacities that exceed the K_T threshold. These areas are, for instance, town dump sites, migratory bird colonies, and carrion along caribou migratory pathways, all of which are potentially regions within a metapopulation that can independently support rabies (Harris 1981; Trewhella et al. 1988; Savory et al. 2014), thus meeting Hethcote’s (1976) requirement for disease persistence in a metapopulation. In addition, ‘floater’ foxes move between social groups, having high mobility rates, and are an overlooked aspect of arctic rabies persistence (Lindström 1989).

Previous rabies modeling research provides meaningful insights and shows some agreement with empirical observations. Specifically, previous models predict spread rates (Källén et al. 1985; Murray et al. 1986; Smith et al. 2002) and define thresholds for vaccination regimes and efficient implementation (Broadfoot et al. 2001; Russell et al. 2006; Asano et al. 2008; Clayton et al. 2010; Neilan and Lenhart 2011). When the carrying capacity of all foxes in an area is greater than the threshold density (K>K_T), either dampened oscillations or limit cycles result, and these characteristics have been well documented in the epidemiological data of red fox populations across Europe and most of North America (Anderson et al. 1981). Still, many critical questions, such as the spread and maintenance of rabies in the Arctic, are largely unresolved (Mørk and Prestrud 2004). Simon et al. (2019) extended Anderson et al. (1981) to consider high transmission rates, short incubation periods, prolonged infectious periods, periodicity in the birth rate, and interaction with red foxes. With these modifications to Anderson’s model, rabies can persist in the Arctic with fox densities lower than 0.15 fox/km², yet at low densities, immigration will cause sporadic outbreaks of rabies, and this leaves open the question of dispersal’s role in arctic rabies endemicity (Simon et al. 2019).

Without spatially structured source-sink disease dynamics, host interactions are not partitioned into a biologically representative scheme but are represented as an average; without spatial structure the potential for rescue effects and disease reintroduction is eliminated. Here, we formulate and parameterize a two-patch disease model to explore rabies persistence in the Arctic and disease dynamics in low-density regions. We identify the necessary conditions for rabies endemicity in a metapopulation and consider how rabies can persist at landscape-level densities below K_T, given the assumptions our model makes about space, connectivity, and individuals.

Methods

Two-patch model

We use a two-patch deterministic model (Figure 1). One of the two patches is a low-carrying capacity patch (K₁<K_T) that is characteristic of many Arctic areas (hereafter referred to as the ‘sink patch’, and where tundra-specific variables and parameters are indexed with a subscript 1). The sink patch is coupled with a higher-carrying capacity patch (K₂ > K₁) to represent the area surrounding a town, migratory pathway, or a migratory bird colony (hereafter referred to as the ‘source patch’, where specific variables and parameters are indexed with a subscript 2).

Figure 1. A two-patch model describing rabies dynamics. The epidemiological status of foxes are susceptible, S₁ and S₂; latent, E₁ and E₂; and infected, I₁ and I₂, where the subscript 1 indicates residence on a sink (low-carrying capacity) patch and the subscript 2 indicates residence on a source patch (high-carrying capacity). The transmission rate is β, the rate of disease progression from exposed to infected is p, and the disease-induced mortality rate is v. The inter-patch movement rate, from a patch i to j is m_ij and the epidemiological status of foxes does not change while travelling between patches. The model assumes that each of these parameters (β, p, and v) are the same for all individuals. The figure does not show reproduction, mortality, or density-dependent constraints for visual clarity. See Equations 1(1) $S_1={\mathrm{r}\mathrm{S}}_1-{\mu }_1{\mathrm{S}}_1{\mathrm{N}}_1-\beta {\mathrm{S}}_1{\mathrm{I}}_1\hspace{0.17em}-\hspace{0.17em}{\mathrm{S}}_1{m}_{12}+{\mathrm{S}}_2{m}_{21}$(1) –6(6) $I_2=\hspace{0.17em}p{\mathrm{E}}_2\hspace{0.17em}-\hspace{0.17em}{\mathrm{I}}_2\left(v+d\right)-{\mu }_2{\mathrm{I}}_2{\mathrm{N}}_2-\hspace{0.17em}{\mathrm{I}}_2{m}_{21}+{\mathrm{I}}_2{m}_{12}$(6) for the complete model

The within-patch epidemiological dynamics, including density-dependent population growth, are based on Anderson et al. (1981). Within a patch, i, the model variables describe the density of susceptible, S_i, latent or exposed, E_i, and infected, I_i, foxes (fox/km²), and the model does not consider recovery as rabies is almost always fatal. The density on patch i is $N_i=S_i+E_i+I_i$, and rabies percent prevalence is defined as $\left(\left(E_i+I_i\right)/N_i\right)\cdot 100$ which is a percent ranging between 0 and 100. The rate that each susceptible fox is exposed to rabies per infected fox is β (foxes/km²)⁻¹ yr⁻¹ (transmission dynamics assume mass action). Foxes in the latent compartment (E_i) have the virus but are not yet contagious to others. The rate of disease progression from E_i to I_i is p, and the rate of disease-induced mortality for infected individuals is v (both with units yr⁻¹). Our formulation assumes that only foxes susceptible to rabies (S₁ and S₂) can reproduce, as the degenerative effects of rabies make it unlikely that foxes of any other epidemiological status would be able to reproduce, and assumes that pups are born susceptible. Foxes disperse between patches at the rate m_ji from patch j to i. Parameter values are the same on both patches except interpatch dispersal, m_ji, and the carrying capacity, K_i, which appears in Equations 1(1) $S_1={\mathrm{r}\mathrm{S}}_1-{\mu }_1{\mathrm{S}}_1{\mathrm{N}}_1-\beta {\mathrm{S}}_1{\mathrm{I}}_1\hspace{0.17em}-\hspace{0.17em}{\mathrm{S}}_1{m}_{12}+{\mathrm{S}}_2{m}_{21}$(1) –6(6) $I_2=\hspace{0.17em}p{\mathrm{E}}_2\hspace{0.17em}-\hspace{0.17em}{\mathrm{I}}_2\left(v+d\right)-{\mu }_2{\mathrm{I}}_2{\mathrm{N}}_2-\hspace{0.17em}{\mathrm{I}}_2{m}_{21}+{\mathrm{I}}_2{m}_{12}$(6) via µ_i, since µ_i = r/K_i (i.e. see Table 1). The model dynamics are,

(1) $S_1={\mathrm{r}\mathrm{S}}_1-{\mu }_1{\mathrm{S}}_1{\mathrm{N}}_1-\beta {\mathrm{S}}_1{\mathrm{I}}_1\hspace{0.17em}-\hspace{0.17em}{\mathrm{S}}_1{m}_{12}+{\mathrm{S}}_2{m}_{21}$(1)

(2) $E_1=\hspace{0.17em}\beta {\mathrm{S}}_1{\mathrm{I}}_1\hspace{0.17em}-\hspace{0.17em}{\mathrm{E}}_1\left(p+d\right)-{\mu }_1{\mathrm{E}}_1{\mathrm{N}}_1-\hspace{0.17em}{\mathrm{E}}_1{m}_{12}+{\mathrm{E}}_2{m}_{21}$(2)

(3) $I_1=\hspace{0.17em}p{\mathrm{E}}_1\hspace{0.17em}-\hspace{0.17em}{\mathrm{I}}_1\left(v+d\right)-{\mu }_1{\mathrm{I}}_1{\mathrm{N}}_1-\hspace{0.17em}{\mathrm{I}}_1{m}_{12}+{\mathrm{I}}_2{m}_{21}$(3)

(4) $S_2=r{\mathrm{S}}_2-{\mu }_2{\mathrm{S}}_2{\mathrm{N}}_2-\beta {\mathrm{S}}_2{\mathrm{I}}_2\hspace{0.17em}-\hspace{0.17em}{\mathrm{S}}_2{m}_{21}+{\mathrm{S}}_2{m}_{12}$(4)

(5) $E_2=\hspace{0.17em}\beta {\mathrm{S}}_2{\mathrm{I}}_2\hspace{0.17em}-\hspace{0.17em}{\mathrm{E}}_2\left(p+d\right)-{\mu }_2{\mathrm{E}}_2{\mathrm{N}}_2-\hspace{0.17em}{\mathrm{E}}_2{m}_{21}+{\mathrm{E}}_2{m}_{12}$(5)

(6) $I_2=\hspace{0.17em}p{\mathrm{E}}_2\hspace{0.17em}-\hspace{0.17em}{\mathrm{I}}_2\left(v+d\right)-{\mu }_2{\mathrm{I}}_2{\mathrm{N}}_2-\hspace{0.17em}{\mathrm{I}}_2{m}_{21}+{\mathrm{I}}_2{m}_{12}$(6)

Table 1. Parameter descriptions for the two-patch rabies model (Equations 1(1) $S_1={\mathrm{r}\mathrm{S}}_1-{\mu }_1{\mathrm{S}}_1{\mathrm{N}}_1-\beta {\mathrm{S}}_1{\mathrm{I}}_1\hspace{0.17em}-\hspace{0.17em}{\mathrm{S}}_1{m}_{12}+{\mathrm{S}}_2{m}_{21}$(1) –6(6) $I_2=\hspace{0.17em}p{\mathrm{E}}_2\hspace{0.17em}-\hspace{0.17em}{\mathrm{I}}_2\left(v+d\right)-{\mu }_2{\mathrm{I}}_2{\mathrm{N}}_2-\hspace{0.17em}{\mathrm{I}}_2{m}_{21}+{\mathrm{I}}_2{m}_{12}$(6) ). Parameters are the same as Anderson et al. (1981) and Simon et al. (2019), but where subscripts denote patch-specific values. The parameter values are the same on each patch, except for carrying capacity and dispersal

Definition	Parameter	Value
Birth rate	a	1 yr^−1
Mortality	d	0.5 yr^−1
Net population growth rate at low densities	r =a-d	0.5 yr^−1
Patch 1 dispersal	$m_{12}$	0.1–0.25 yr^−1
Patch 2 dispersal	$m_{21}$	0.1–80 yr^−1
Patch 1 carrying capacity	K1	0.0–1.0 fox/km^2
Patch 2 carrying capacity	K2	0.1–5.0 fox/km^2
Latency	p	13 yr^−1
Disease-induced mortality	v	73 yr^−1
Transmission coefficient	$\beta$	80 km^2/fox$\cdot$yr
Density-dependent constraints	µi =r/Ki	Varied km^2/fox$\cdot$yr

Parameter values provided in Table 1 are based on Anderson et al. (1981) and Simon et al. (2019), and parameter values were also cross-referenced with Mørk and Prestrud (2004) to ensure biological relevance for arctic fox populations. The parameter set, derived from demographic studies of foxes, yields a theoretically simulated disease density threshold of 1 fox/km². Anderson et al. (1981) notes that observed epidemiological evidence suggests K_T is between 0.25 and 1.0 fox/km²; data which reflect realized thresholds rather than a fundamental threshold. We move forward with the above parameter set to provide comparable results to equivalent spatially homogenous models, and changes in K_T via our different patch and demographic configurations are relative to the initial threshold, K_T = 1 fox/km², and those baseline parameters.

To estimate carrying capacity, we assumed that the estimated densities of foxes are near carrying capacity, and we included breeding pairs, juveniles that remain at the den, and adult non-breeding foxes (floaters). The density of breeders is estimated as 0.02–0.3 breeders/km² (Angerbjörn et al. 1999). The regression given by Strand et al. (1995) estimates litter size from placental scars, giving an average litter of 9 pups for a breeding pair. Pup to juvenile survival is about 10% (Meijer et al. 2008), so there are approximately 0.45 juvenile foxes/breeder/year. Next, we consider floaters, which have been documented as up to 25% of the population (Lindström 1989), and for this study we assume floaters are on average 20% of the population. The carrying capacity estimate, per breeder, is (1$+0.45)/0.8=1.81$. As the density of breeding pairs spans a range that is 15x greater than its lower bound, we will estimate low-, mid- and high-carrying capacity values for the sink patch to acknowledge this uncertainty. Our estimates for the carrying capacity on the sink patch are:

$K_1=0.02\cdot 1.81=0.036\mathrm{f}\mathrm{o}\mathrm{x}/k{m}^2$ (Low),

$K_1=0.16\cdot 1.81=0.29\mathrm{f}\mathrm{o}\mathrm{x}/k{m}^2$ (Mid),

$K_1=0.3\cdot 1.81=0.54\mathrm{f}\mathrm{o}\mathrm{x}/k{m}^2$ (High).

The transmission rate for rabies is estimated to be $\beta =80k{m}^2/\mathrm{f}\mathrm{o}\mathrm{x}\cdot \mathrm{y}\mathrm{r}$ (Anderson et al. 1981; Llyod 1980). However, it is unclear whether this estimate is calculated for resident foxes or also considers highly mobile ‘floaters’, which may be 7–10 times more mobile (Lindström 1989). When floaters are excluded from the $\beta$ estimate, and given our previous assumption that floaters are 20% of the population, an estimate of the transmission rate that considers floaters would range from $\beta =176to224k{m}^2/\mathrm{f}\mathrm{o}\mathrm{x}\cdot \mathrm{y}\mathrm{r}$. Further explanation of this calculation and its assumptions can be seen in Appendix 2.

To perform our analyses, we numerically solved Equations (1(1) $S_1={\mathrm{r}\mathrm{S}}_1-{\mu }_1{\mathrm{S}}_1{\mathrm{N}}_1-\beta {\mathrm{S}}_1{\mathrm{I}}_1\hspace{0.17em}-\hspace{0.17em}{\mathrm{S}}_1{m}_{12}+{\mathrm{S}}_2{m}_{21}$(1) –6(6) $I_2=\hspace{0.17em}p{\mathrm{E}}_2\hspace{0.17em}-\hspace{0.17em}{\mathrm{I}}_2\left(v+d\right)-{\mu }_2{\mathrm{I}}_2{\mathrm{N}}_2-\hspace{0.17em}{\mathrm{I}}_2{m}_{21}+{\mathrm{I}}_2{m}_{12}$(6) ) using the ode45 function in MATLAB 2018a until the total population size reached an equilibrium state. The system of equations and the computer code are available via the following link: https://doi.org/10.6084/m9.figshare.c.4840710

Results

Low to intermediate levels of dispersal in heterogeneous environments can support rabies in low-carrying capacity ‘sink’ patches, where rabies would be absent if the patch was isolated, or if the landscape was homogeneous with the carrying capacity equal to this low value, K₁ (Figure 2). As the carrying capacity on the sink patch increases toward the threshold carrying capacity in a homogeneous environment, K_T, the sink patch requires less augmentation via dispersal from the source patch to maintain disease, which is reflected by the expanding parameter space for endemicity as shown for the low- (Figure 2a), mid- (Figure 2b), and high- (Figure 2c) estimates for K₁. Disease is supported in the sink patch until dispersal from the source patch removes too many susceptible and infected foxes, such that disease dynamics cannot be maintained on the source patch, at which point, the disease dies out in both patches (Figure 3).

Figure 2. Rabies is endemic to the sink patch (low-carrying capacity), where disease would otherwise be absent, when low to intermediate levels of movement couples the disease dynamics between the sink and source (high-carrying capacity) patches. Panels show rabies prevalence, (100(E₁+I₁)/N₁), on the sink patch (% prev (1)) for our three estimates of carrying capacity with K₁ equal to low: 0.036 (a); mid: 0.29 (b), and high: 0.54 (c) fox/km². Rabies prevalence (%) on the sink patch is highest for large values of the carrying capacity in the source (large K₂), large values of the carrying capacity in the source (c; high estimate of K₁), and for intermediate movement rates to the sink from the source patch (m₂₁), a pattern that is explained further in Figure 3. Parameter values are given in Table 1 and m₁₂ = 0.25

Figure 3. Rabies prevalence (%) on the sink patch (a) peaks at intermediate movement rates because high movement rates eradicate the infection in the source patch (b). For high levels of movement, the source patch cannot maintain its function as a disease source when the density of susceptible and infected foxes are depleted through movement and become too few to sustain the epidemic. Parameter values are given in Table 1, with m₁₂ = 0.25, K₁ = 0.54 fox/km², and K₂ = 5 fox/km²

In a heterogeneous two-patch landscape, rabies can persist when the landscape-level mean carrying capacity is below the threshold carrying capacity for rabies endemicity in a homogeneous landscape, K_T (Figure 4). When the landscape-level mean carrying capacity is fixed at 0.9 foxes/km², and different combinations of K₁ and K₂ are considered, we find that rabies prevalence is highest when the carrying capacities on the sink and source patches are the most different (Figure 4). However, for lower values of the landscape-level mean carrying capacity, for example, $\bar{\mathrm{K}}$ = 0.5 fox/km², rabies endemicity is not possible (see Appendix 3), as infection prevalence on the source patch is too low to both sufficiently subsidize the sink patch and maintain disease locally.

Figure 4. Rabies can be endemic in a heterogeneous landscape where the mean carry capacity is less than the threshold carrying capacity for rabies endemicity in a homogeneous landscape. For our parameter values (Table 1), rabies is endemic in a homogenous landscape when the carrying capacity is greater than 1 fox/km² (K>K_T = 1; (Andersen and May (1981)). We set the landscape-level mean carrying capacity for our two-patch model to $\bar{\mathrm{K}}$ = (K₁ +K₂)/2 = 0.9 foxes/km². When K₁ =K₂ =$\bar{\mathrm{K}}$ =0.9< K_T (far right on the x-axis), no disease occurs on either patch since the landscape is homogeneous; however, as the variance between the two carrying capacities on each patch increases (toward the left on the x-axis), rabies becomes established on both patches (red and blue lines). The between-patch movement rates are m₁₂ =m₂₁ =0.25, and all other parameter values are as given in Table 1

We estimate that rabies can persist in fox populations with heterogeneous environments where the mean carrying capacity is as low as ~0.25 foxes/km² and the carrying capacity in the sink patch is as low as 0.036 foxes/km² (Figure 5). To generate this lower bound, we considered the maximum feasible estimate of β = 224 km²/fox$\cdot$yr and our lowest estimate of fox carrying capacity.

Figure 5. Higher transmission rates, for example, due to highly mobile ‘floater’ foxes, allow for rabies endemicity for landscape-level carrying capacities as low as 0.25 foxes/km². The transmission rate is elevated by floaters or higher mobility foxes (we consider a range of values beginning from the baseline value of 80 and increase to 224 km² foxes⁻¹ yr⁻¹). With the carrying capacity on the sink patch set to its lowest estimate: K₁ =.036 fox/km², we find that rabies can persist when the carrying capacity on the source patch is ~0.46 foxes/km² (the value of K₂ for the blue contour when β = 224 km²/foxes⁻¹ yr⁻¹) corresponding to a landscape-level mean carrying capacity of ~0.25 foxes/km², as seen on the right y-axis. Parameter values are m₁₂ =m₂₁ =0.1 yr⁻¹, and all other parameters are as described in Table 1

We found that rabies persistence is sensitive to our assumptions restricting the epidemiological status of dispersing foxes (Figure 6). We set K₁ = 0.54 foxes/km² and K₂ = 2 foxes/km² and found that irrespective of the level of dispersal from the source to the sink patch, rabies could never persist if all epidemiological compartments disperse equally, or if only susceptible foxes disperse, for dispersal rates >1 yr⁻¹ (Figure 6, black lines). By contrast, rabies did persist when only latent or infected foxes disperse (Figure 6, blue and red lines). Rabies persistence occurred over the widest range of ecological conditions when only infected foxes dispersed; however, rabies prevalence reached its highest value when only latent foxes disperse.

Figure 6. Rabies can be endemic on the sink patch, when it might otherwise be eliminated (black), if between-patch movement occurs for only infected foxes (red) or only latent foxes (blue). Assuming that between-patch movement occurs for foxes with all epidemiological statuses or only for susceptible foxes, then rabies cannot persist for a wide range of movement rates from the source to the sink patch (black dashed and dotted lines), but when assuming only infected foxes move between patches rabies is present on the sink patch for a wide range of movement rates (red). Assuming that only latent foxes move between patches, rabies prevalence on the sink patch can be high, but results in the extinction of the fox population for high movement rates (blue). The model formulation is described in Appendix 1 and parameters are given in Table 1, with K₂ = 2 fox/km², K₁ = 0.54 fox/km², and m₁₂ = 0.25 yr⁻¹

Discussion

Low fox densities and spatial discontinuity heavily influence rabies disease dynamics in the Arctic. Spatially homogenous disease models, which assume uniform density and mixing, typically give a useful simplification of infection dynamics, as seen in Anderson et al. (1981). The increased resolution that heterogeneity gives a model is especially important near threshold values, where small changes can shift the qualitative outcome of a deterministic model. Using spatial heterogeneity as a lens to view arctic rabies, we extended the work of Simon et al. (2019) whereby infected animals subsidize the system via an unspecified source without any dynamics. Doing so, we determined the conditions in which rabies can persist in the Arctic under a two-patch structure, where infectious individuals move between source and sink patches.

We showed endemicity in the low-carrying capacity patch when paired with a higher carrying capacity disease source patch. This was observed across a range of low-densities and dispersal values (Figure 2), and when latent only or infected only individuals move between patches (Figure 6). We found that rabies can persist when the average carrying capacity, $\bar{\mathrm{K}}$ = (K₁+K₂)/2, across the two patches in the metapopulation is less than K_T = 1 fox/km² (Figure 4). Our results are consistent with previous studies. Like Hethcote (1976), we found that for a low-carrying capacity patch to support disease, the subsidization from the source patch must not be so much that it suppresses disease dynamics in that patch (Figure 3). We also found equilibrium rabies prevalence in the metapopulation peaks with low-level symmetric bi-directional dispersal (Figure 2) as noted in Gurarie and Seto (2008), and further extend those results by showing that high levels of dispersal can maximize disease prevalence when only latent, or infectious individuals disperse (Figure 6).

We examined the potential for higher mobility of arctic foxes via the transmission coefficient, β, as floater foxes may play a major role in disease transmission. We showed endemicity in all low-carrying capacity patches across a wide range of β values, and for larger β values, even when the source patch carrying capacity (K₂) was well below K_T (Figure 5). Figures 2 and 4 show the plausibility of rabies on low-carrying capacity patches; however, if we assume that the neurological effects of rabies make its host the most likely demographic to disperse (Figure 6), then there is a much larger parameter space that will allow for rabies persistence in the Arctic than seen in Figure 2. This is particularly relevant, given that specific strains of the virus can influence the dispersal and habitat permeability to infected individuals (Scott 1988), and rabies may alter movement patterns (Barton et al. 2010), which is in contrast to the well-studied territoriality and demographic dispersal events observed in fox populations.

We independently considered several mechanisms that facilitated rabies endemicity at low carrying capacity, yet it is likely that many of these mechanisms are working in combination to create an environment where rabies is endemic at extremely low densities. Furthermore, we used a two-patch system; however, additional patches would allow for a lower threshold for endemicity, as a single patch could potentially experience the additive effects of multiple patches contributing to or subsidizing that single patch population and disease dynamics. With the presence of a disease reservoir, i.e., around towns or goose colonies, floater foxes with higher transmission rates, and dispersal of foxes with latent or clinical rabies infection, we suggest it is feasible to have rabies endemicity in landscapes where the average landscape density is 0.25 fox/km², and where sink patches have fox densities of 0.036 fox/km² (see Figure 5 for details). Whereas we consider patch dynamics in the context of spatial heterogeneity and the related assumptions, it has also been noted that positive temporal autocorrelation and dispersal can enhance metapopulation persistence (Matthews and Gonzalez 2007; Roy et al. 2005), even when that metapopulation is composed entirely of sinks (Manojit et al. 2005).

A discrete-patch metapopulation ensemble resembles the structure of the Arctic: the discrete patches are disconnected by several land masses (e.g. mainland Canada), island systems (e.g. Greenland), and archipelagos (e.g. Svalbard) connected by sea ice. Genetic studies show that strains of rabies are spread between these patches, as consistent with the assumptions of metapopulation epidemic models for arctic foxes (Raundrup et al. 2015; Hanke et al. 2016). The geography of Svalbard, an endemically infected area, is consistent with the discrete space assumption of a disease metapopulation model. The most likely origin of rabies in Svalbard is via the migration of arctic foxes from Greenland or the Siberian islands (Mørk et al. 2011). On the island of Spitsbergen in the Svalbard archipelago, the density of breeding foxes is approximately 0.1–0.15 fox/km² with prevalence values of ~0.3% (Eide 2002; Mansfield et al. 2006; Mørk et al. 2011). These results closely mirror those seen in our low-carrying capacity patches when connected to a disease source.

Rabies is usually absent on the island of Newfoundland, but in the spring of 2002, Newfoundland saw its first outbreak of rabies in 14 years (Nadin-Davis et al. 2008). It is assumed that the disease was introduced to the island from an infected mainland fox that travelled on an ice sheet across the 9-mile-wide Strait of Belle Isle from Labrador. This outbreak captures the essence of arctic rabies; it is a disease that is largely governed by the spatial arrangement of the Arctic and its connectivity. Here, we examined arctic rabies through this lens using a two-patch structure and found that rabies can persist endemically in the Arctic via source-sink dynamics, partitioning of densities, selective dispersal, and increased mobility. These all provided realistic parameter space for rabies endemicity, further supported by the disease dynamics observed in Arctic regions characterized by a metapopulation structure. Our study expanded upon the limited body of research surrounding the persistence of arctic rabies and showed the feasibility of endemicity at low densities by considering spatial dynamics.

Acknowledgments

We would like to thank J. Marino, A. Fofana, A. McLeod, and M. Noonan for their helpful comments on previous drafts of this paper.

Disclosure statement

The authors report no financial interests or conflict of interests.

Additional information

Funding

A. Hurford was supported by an NSERC Discovery Grant (RGPIN 2014-05413), and N. Lecomte by the Canada Research Chairs program; Natural Sciences and Engineering Research Council of Canada [RGPIN 2014-05413];

Appendix 1

The following model is a modified version of the previously stated two-patch model. In addition to our prior formulation, each equation here has two additional terms. The general dispersal from patch to patch, m₁₂ and m₂₁ have been modified to $m_{12}^S$, $m_{12}^E$, $m_{12}^I$, and $m_{21}^S$, $m_{21}^E$,$m_{21}^I$. The addition of S, E, and I superscripts to each dispersal term represents the epidemiological compartment that is dispersing. This formulation allowed us to selectively disperse foxes of different disease statuses.

(1) $S_1=r{S}_1-{\mu }_1{S}_1{N}_1-\beta S_1{I}_1-S_1{m}_{12}-\hspace{0.17em}{S}_1{m}_{12}^S+S_2{m}_{21}+S_2{m}_{21}^S$(1)

(2) $E_1=\hspace{0.17em}\beta S_1{I}_1-E_1\left(p+d\right)-{\mu }_1{E}_1{N}_1-\hspace{0.17em}{E}_1{m}_{12}-\hspace{0.17em}{E}_1{m}_{12}^E+E_2{m}_{21}+E_2{m}_{21}^E$(2)

(3) $I_1=\hspace{0.17em}p{E}_1-I_1\left(v+d\right)-{\mu }_1{I}_1{N}_1-\hspace{0.17em}{I}_1{m}_{12}-\hspace{0.17em}{I}_1{m}_{12}^I+I_2{m}_{21}+I_2{m}_{21}^I$(3)

(4) $S_2=r{S}_2-{\mu }_2{S}_2{N}_2-\beta S_2{I}_2-S_2{m}_{21}-\hspace{0.17em}{S}_2{m}_{21}^S+S_2{m}_{12}+S_2{m}_{12}^S$(4)

(5) $E_2=\hspace{0.17em}\beta S_2{I}_2-E_2\left(p+d\right)-{\mu }_2{E}_2{N}_2-\hspace{0.17em}{E}_2{m}_{21}-\hspace{0.17em}{E}_2{m}_{21}^E+E_2{m}_{12}+E_2{m}_{12}^E$(5)

(6) $I_2=\hspace{0.17em}p{E}_2-I_2\left(v+d\right)-{\mu }_2{I}_2{N}_2-\hspace{0.17em}{I}_2{m}_{21}-\hspace{0.17em}{I}_2{m}_{21}^I+I_2{m}_{12}+I_2{m}_{12}^I$(6)

Appendix 2

This model uses the same formulation as previously used, with the exception of the transmission term, which is further discussed here.

Blackwood formulation: βS(I + φ) where the parameter φ represents the additive disease effect from floaters. This model assumes that when the infected animals on a patch go to zero, the floaters (ie. φ) continues to infect the population. Therefore, when disease dynamics lead to a fade out or burnout, the equilibrium prevalence of a population is proportional to φ. New infections are transmitted from floaters in every timestep, regardless of the population dynamics of residents, forcing an unstable disease-free equilibrium (Brauer and van den Driessche 2001).

Our formulation: $\left(\beta \cdot b_f\right)\left(I\cdot fp\right)S$, where b_f and fp are parameters that represent the efficacy of floater transmission relative to residents, and the proportion of infected floaters relative to residents. With this modification, we are able to elevate the population level transmission coefficient to account for floater-driven interactions. Although this model includes higher mobility through an additive infection term, we multiply the additive effect of floaters by the resident population; therefore, when S or I go to zero, so do the effects of floaters. This allows us to relax the assumptions of Blackwood et al. (2013) and Simon et al. (2019) that if the infected population on a patch goes extinct, then floaters can still infect the resident population by reemerging every timestep in a proportion equal to φ. Because of this assumption, the Blackwood formulation forces an unstable disease-free equilibrium because new infections are transmitted from floaters in every timestep, regardless of resident population dynamics, specifically K< K_T. This formulation is not seen in the paper, as we just report the beta values that result from this multiplication.

Appendix 2 References

Blackwood JC, Streicker DG, Altizer S, Rohani P. 2013. Resolving the roles of immunity, pathogenesis, and immigration for rabies persistence in vampire bats. P Natl Acad Sci USA. 110:20,837–20,842.

Brauer F, van den Driessche P. 2001. Models for transmission of disease with immigration of infectives. Math Biosci. 171: 143–154.

Simon A, Hurford A, Lecomte N, Belanger D, Leighton PA. 2019. Dynamics and persistence of rabies in the Arctic. Polar Res. 38: 33-66.

Appendix 3

This figure shows that in our two-patch system, with lower values of the landscape-level mean carrying capacity, rabies endemicity is not possible.

Figure A3. Very low landscape level averages $\bar{K}=0.7fox/k{m}^2$ (a) and $\bar{K}=0.5fox/k{m}^2$ (b) provide limited parameter space for endemicity (a) and can be exhausted of its disease by dispersal without effectively infecting the other patch (b). is taken across the two patches in the landscape. The landscape-level average prevalence (%) of rabies, when the mean carrying capacity, $\bar{\mathrm{K}}$= 0.7 fox/km² (a), and, $\bar{\mathrm{K}}$= 0.5 fox/km² (b) is for an ascending value of the low-carrying capacity patch 1, K₁, and the carrying capacity on patch 2 is 2$\bar{\mathrm{K}}$ – K₁. When $K_1=K_2=\bar{K}/2<K_T$, no disease occurs on either patch; however, as the variance in the K_i between the two patches increases (toward the left side), rabies becomes established on both patches in (a), although at extremely low prevalences. Dispersal is unidirectional with m₁₂ = 0 and m₂₁ = 0.1. Unless otherwise stated, all parameter values are as given in Table 1