^{1}

^{2}

^{3}

^{1}

^{4}

^{5}

^{6}

^{7}

^{8}

^{9}

^{1}

^{*}

^{1}

^{2}

^{3}

^{4}

^{5}

^{6}

^{7}

^{8}

^{9}

Edited by: Alemka Markotic, University Hospital for Infectious Diseases “Dr. Fran Mihaljevic”, Croatia

Reviewed by: Jan Clement, KU Leuven, Belgium; Vitaly V. Ganusov, The University of Tennessee, Knoxville, United States

This article was submitted to Virus and Host, a section of the journal Frontiers in Cellular and Infection Microbiology

This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

For wildlife diseases, one often relies on host density to predict host infection prevalence and the subsequent force of infection to humans in the case of zoonoses. Indeed, if transmission is mainly indirect, i.e., by way of the environment, the force of infection is expected to increase with host density, yet the laborious field data supporting this theoretical claim are often absent. Hantaviruses are among those zoonoses that have been studied extensively over the past decades, as they pose a significant threat to humans. In Europe, the most widespread hantavirus is the Puumala virus (PUUV), which is carried by the bank vole and causes nephropathia epidemica (NE) in humans. Extensive field campaigns have been carried out in Central Finland to shed light on this supposed relationship between bank vole density and PUUV prevalence and to identify other drivers for the infection dynamics. This resulted in the surprising observation that the relationship between bank vole density and PUUV prevalence is not purely monotonic on an annual basis, contrary to what previous models predicted: a higher vole density does not necessary result in a higher infection prevalence, nor in an increased number of humans reported having NE. Here, we advance a novel individual-based spatially-explicit model which takes into account the immunity provided by maternal antibodies and which simulates the spatial behavior of the host, both possible causes for this discrepancy that were not accounted for in previous models. We show that the reduced prevalence in peak years can be attributed to transient immunity, and that the density-dependent spatial vole behavior, i.e., the fact that home ranges are smaller in high density years, plays only a minor role. The applicability of the model is not limited to the study and prediction of PUUV (and NE) occurrence in Europe, as it could be easily adapted to model other rodent-borne diseases, either with indirect or direct transmission.

For wildlife diseases, to make an inference on infection prevalence and the force of infection to reservoir and spillover hosts (e.g., humans), host population density is often used as a proxy (e.g., Davis et al.,

Because of this complexity, it can be difficult to understand the relationship between host density, infection prevalence and force of infection to humans. Moreover, in wildlife disease systems there is often a lack of longitudinal infection data with sufficient (temporal) resolution, as this requires long-term labor-intensive sampling of a large number of host animals in the field (Allen et al.,

Hantaviruses are among the best studied (re-)emerging zoonotic pathogens, whose transmission is considered to be density-dependent (Mills et al.,

While the relationship between prevalence and density may be complex on a short timescale (within the year) (Adler et al.,

PUUV is a hantavirus carried by bank voles (

The density of bank voles (black) and the density of the different infection status categories: susceptible (and exposed) for PUUV infection (blue), infectious (red), and carrying maternal antibodies (green) on the core 5.8 ha study grid (data from Voutilainen et al.,

This observation is supported by another independent record: the number of human NE cases reported during the same time period. Although the number of cases can be expected to be proportional to bank vole density, it is not higher in a bank vole peak year but often even slightly lower than in increase years (Kallio et al.,

In order to explain this discrepancy, Kallio et al. (

Another explanation may be that the spatial behavior of bank voles varies with density. In Finland, it has been observed that in peak years, when the abundance is high already at the start of the breeding season, breeding female voles have smaller territories compared to low-abundance years (Koskela et al.,

In this paper we tested two hypotheses as possible causes for the reduced density of infectious voles (and the reduced number of human NE-cases) in a peak year: the transient immunity of maternal antibodies and the spatial behavior of bank voles during a peak year. To this end, we advanced a new spatially-explicit model for PUUV dynamics that includes maternally derived waning immunity and reduced home ranges in high density years. We also considered different shedding patterns (constant/decreasing) and different transmission routes (direct/indirect) in the simulations, since these system parameters are difficult to assess and not well-agreed upon. This allowed us to test to what extent the conclusions drawn are sensitive to these parameters. Understanding the relationship between density and prevalence is not only of practical use as it would help to predict disease outbreaks and the possible infection risk for humans, but, as we will see further, this also allows us to deduce important model parameters, hence giving insight into the virus-host interactions driving the infection dynamics.

Earlier hantavirus models have dealt mainly with hantaviruses in the New world, e.g., the Sin Nombre virus carried by deer mice, as they pose a significant threat to humans. A basic population dynamics model was introduced by Abramson and Kenkre (

For models of the Old world PUUV (Sauvage et al.,

In this paper, we developed a novel spatially-explicit individual-based model and since we hypothesize that the varying spatial behavior of the bank vole is important for the virus transmission dynamics, we put great effort into modeling the bank vole movements, which depend on sex, maturation stage and the season. Therefore, each of the bank voles has to be simulated individually and the home range of each individual needs to be able to change over time and space. Because of the fact that the virus transmission is likely to be indirect (via feces, urine, saliva, Kallio et al.,

Bank vole movements are modeled on an

We consider four maturation stages: juveniles, subadults, dispersing adults and settled adults, of which only the settled adults can reproduce (Bondrup-Nielsen and Karlsson,

In the model, female settled adult bank voles have a time-dependent probability to produce a five offspring litter given by

see

Bank vole demography as it is implemented as a 3-years cycle in the model simulations. These same conditions shown in are repeated every 3 years.

After birth, juveniles stay close to their maternal site for 1 month, after which they are ready for maturation (Viitala et al.,

Dispersing bank voles roam around for about 2 weeks, before settling in a new place. The transition from dispersing adult to settled adult is governed by the settling rate μ = 26/year, which is the reciprocal duration of the dispersing stage. After settling, adults start reproducing (Bujalska,

In Central-Finland, bank voles live for about 7.5 months in the wild, and hence the adult mortality rate was set to ^{*}(

We assumed the population density fluctuates with a 3 years cycle, as is typically the case for Central-Finland (see ^{*}(^{2}(^{2}(

and have chosen _{1} = 655, _{2} = 0.185, and _{3} = 1.3 such that the resulting simulated population dynamics curve, shown in

The simulated bank vole density during a 3 years cycle for different life stages (averaged over 100 cycles; the width equals twice the standard deviation). Note that adult density is limited to ≈10 bank voles/ha. In winter, all juveniles have turned either to adults or subadults; the latter mature (and turn into new adults) in spring.

The bank vole density and adult density vary at a local scale within the larger grid. Therefore, every bank vole “experiences” its proper bank vole (adult) density, which we calculated within a circular area with radius

If we define respectively the distribution of bank voles on the grid

where ^{2}:

Note that this also allows for the carrying capacity (and consequently the mortality) to be spatially-dependent. Indeed, if the environment is heterogeneous, e.g., the vegetation varies throughout the patch, this will be reflected in the local carrying capacity, which will change across the grid. This is beyond the scope of this paper, though, and we have assumed the carrying capacity to be constant over the grid.

Bank voles have a home burrow and their movements are limited to their home range. During winter, all bank voles display more or less the same spatial behavior: they all have rather small home ranges of ~0.2 ha (Eccard et al.,

Adults and juveniles move around diffusively, but are bound to their territory at position _{A}. The probability density function (pdf) of this adult (or juvenile) can then be written as the stationary solution of the harmonic model,

which is essentially a Gaussian pdf, as derived by Abramson et al. (

Short-term movement is the movement during a time step and is modeled by a Gaussian distribution with standard deviation σ_{Δt} (see

Schematic illustration of how bank vole movement is implemented in the model.

In order to model the spatial behavior on the longer term and to arrive at a long term pdf with standard deviation σ, we simulate the spatial behavior as follows: for every time step Δ

Apart from the fact that this approach gives the flexibility necessary to vary the home range size (depending on the season, sex, maturation status), it simplifies (and speeds up) the simulations significantly, while capturing the essential spatial dynamics. Moreover, it also makes sense in a biological way. Indeed, for large home ranges, the bank vole cannot travel its entire home range during a single time step, as its movement is spatially limited due to the short time interval and its limited mobility. Therefore, as it can not sample its total home range, the long-term movement has to be modeled stochastically.

Using the common definition of home range as the area in which the vole spends 95% of its time, we can then easily deduce σ which corresponds to this area from Equation (4):

such that a home range (area) of 0.2 ha corresponds to σ ≈ 10.3 m.

During winter we assume that the home range size is equal for all voles, irrespective of their sex and maturation status (σ = 10 m). At the start of the breeding season, the home range size increases for adults, and the home range is different for females (0.5 ha, σ = 16 m) and males (2 ha, σ = 32 m). At the end of the breeding season, the home range of the remaining adults again shrink to the smaller home range sizes (0.2 ha, σ = 10 m). For simplicity, we chose the variance of the short term kernel equal to the variance corresponding to the smallest home ranges, i.e., σ_{Δt} = 10 m. Hence, we only need to model the long term movements in the case of larger home ranges.

The above simulation applies to all settled voles, adults and juveniles, whose home range is fixed. But after maturation, new adults disperse and start moving around in search for a place to settle. We model this non-localized emigration as a random walk with turning angle distribution: at each discrete time step, the dispersing new adult travels a straight fixed distance Δ

such that the vole has a higher probability to keep on moving in the same direction. We know the voles travel about 1 km during their 2 weeks period of dispersal (Viitala et al.,

In our model, these new adults roam around for on average 2 weeks and those that have the least home range overlap with other adults mature first.

The infectiousness of a vole (the amount of virus shed per time step) over time is modeled by the infection curve

i.e., we assume that voles undergo an acute phase of about a month, followed by a chronic infectious stage with lower shedding (10%) which lasts for the rest of their life.

In recent years, an extensive (capture-mark-recapture) dataset was collected to assess the virus shedding pattern of voles in the field (Voutilainen et al.,

to test the sensitivity of the simulation results to the actual shedding pattern.

Because at a very young age, juveniles spend most of their time inside the burrow, under protection of a territorial mother, their virus contamination in the environment is likely to be very limited, hence we assume that infected juveniles only start contaminating the environment at the age of 2 weeks.

A newborn is born free of infection, as there is no vertical transmission of PUUV. But if the mother is infected, the newborn will be temporarily protected by MatAb (Gavrilovskaya et al.,

In a longitudinal study in which offspring from infected mothers were exposed to virus infected bedding at different ages, Kallio et al. (

The virus is transmitted through infected saliva, feces and urine. An animal can get infected through a direct (hostile or friendly) encounter with an infected animal or through contact with infected ground. Since bank voles are non-territorial during the main transmission season (winter), when transmission via fighting related wounds is improbable (Voutilainen et al.,

The first term accounts for the viral memory of the contaminated ground, whose viral load decays at a rate of ξ. The second term describes the accumulation of virus during the current time step, due to summed virus shedding of voles _{j}, _{j}) visit the grid point (_{j}, _{j}), which is the value of the short term pdf of vole _{j}(_{1} in the equation (fitting method described below).

In the model, virus survival in the environment is governed by decay rate ξ. Kallio et al. (

Susceptible voles move within the grid of contaminated ground _{i}, _{i}), will be infected during a single time step, one has to calculate to what extent the vole gets in contact with virus contaminated ground. Again making use of the short-term pdf describing the bank vole movement during a single time step, this can be calculated as

where β_{2} is again an unknown factor, describing the transmission efficiency from infected ground to vole (fitting method described below). The kernel _{i}, _{i}) describes to what extent the susceptible vole

In order to simplify and speed up calculations, we define the infectiousness matrix as

which is essentially the distribution of the infectious voles _{j}(

so that the calculation of the second term is reduced to the 2D convolution of the infectiousness matrix

where ^{*}(^{2}). Hence, including the movement of the susceptible vole itself does not add an extra computational burden, since we simply have to convolve the infectiousness matrix ^{*}. On the contrary: since the width of the kernel is slightly larger (a factor _{2} can now be merged into one single unknown parameter β. This fitting parameter will be adapted to fit the simulation results to the field data, as is discussed below. The grid element size Δ^{*} is sampled sufficiently.

Using this model on a standard pc, we are able to model patches as large as 10 × 10 km, keeping track of and modeling the interaction of more than a million voles.

In the above equations, we have assumed that indirect transmission,via virus accumulated in the environment is the most dominant route. As mentioned before, it is also possible that voles get infected through direct contact, e.g., by mating, biting, exploratory sniffing, etc. If one wants to study only direct transmission, one should use:

instead of Equation (11).

In the following we show simulation results for different model scenarios. For every model scenario, there is only one unknown parameter left which cannot be estimated from the literature: the parameter β, which is related to the efficiency of virus transmission from an infectious bank vole to the ground and from the ground to a susceptible bank vole. In our analyses, we chose to fit β for every model scenario such that the mean peak density of infectious voles in an increase year matched that measured in the field (with resolution of 0.1 voles/ha). The mean peak density of infectious voles in the field was obtained by averaging the two peak densities of infectious voles in the increase years 2004 and 2007 (see

In order to compare the simulated infection dynamics with the field data and to assess the “quality” of the model, we focused on a selected set of features of the infection dynamics. Our primary feature was the mean peak density of infectious bank voles in the peak abundance year following an increase year. From

In the analysis, we focus on these three features and discuss to what extent these are in agreement with the field data (using the quality measures). We deliberately refrained from an elaborate statistical evaluation of the quality of the different models, because we have only two full 3-years cycles of field data to compare the model output with, and consequently any statistical test would lack power. One could argue that one cycle consists of multiple data points, but these data points are highly correlated in time. Moreover, the simulated bank vole density shown in

All simulations were done on an

The first column shows the different values for the model parameter β such that the maximal density of infectious bank voles in the increase year equals 8.4 bank voles/ha (marked by “+” in

^{−1}vole^{−1}) |
|||||
---|---|---|---|---|---|

16,400 | 16.6 | 1.66 | 1.9 | 1.0 | |

16,400 | 16.2 | 1.62 | 1.9 | 1.0 | |

59,000 | 11.3 | 0.99 | 2.0 | 1.4 | |

59,000 | 11.0 | 0.96 | 2.0 | 1.5 | |

59,000 | 11.0 | 0.96 | 2.0 | 1.5 | |

300,000 | 10.2 | 1.33 | 1.2 | 1.0 | |

29,000 | 10.1 | 0.64 | 3.4 | 3.7 |

First, we ran simulations without including immunity or density-dependent spatial behavior. The immunity period was set to zero, such that antibodies were cleared after the first day (i.e., effectively removing the effect of maternal antibodies). The spatial behavior only varied within the year (σ = 10 m for all voles, except during the breeding season, when σ for males and females temporarily increased to σ = 16 and 32 m, respectively) but did not differ between years. We assumed indirect virus transmission (via the environment) and an acute and a chronic infectious phase where virus shedding of infected individuals is lower (10%) in the latter. The simulation results are shown in

Simulation results for a 3-years cycle (averaged over 100 cycles; the width equals twice the standard deviation) for the different scenarios, see text:

Next, we considered the effect of inter-annual variation in density-dependent spatial behavior of bank voles on the infection dynamics. Again, we assumed that there is no maternally-derived immunity, but now, in a peak year, we reduced the home ranges of female and male voles, respectively to σ = 10 and 16 m during the breeding season. The results in

To study the effect of immunity on the infection dynamics, we assumed that the voles are protected for about 80 days, as found by Kallio et al. (

If, in addition to immunity, we also include the density-dependent spatial behavior, the results in

Comparing the simulation results to the field data shown in

In order to test what the effect of this would be, we have plotted in

Simulation results (averaged over 100 cycles; the width equals twice the standard deviation) for different scenarios, see text:

In the above simulations, we have considered only the indirect route (Kallio et al.,

In all previous simulations, we have assumed that the voles go through an acute phase of about a month, followed by a chronic infectious stage which lasts for the rest of their life. In case the amount of virus shed does not decrease over time but remains constant, the results are shown in

Incorporating the effect of maternal antibodies on zoonotic infection dynamics has long been neglected (Boulinier and Staszewski,

Our model simulations also help explaining the occurrence of NE in the human population in Central-Finland during the same time period. In line with the pattern observed for the abundance of infectious voles, Kallio et al. (

Although the density of infectious voles in a peak year matched that measured in the field, the simulated density of immune voles was much higher than measured in the field. This difference may be due to the limited sensitivity of the antibody assay to detect MatAb immune voles (Gilbert et al.,

Bank voles can transmit an infection through direct or indirect contact (Kallio et al.,

Related to this, there is still no conclusive evidence on the exact shedding pattern that has to be taken into account (Hardestam et al.,

It is plausible that in reality both indirect and direct routes are to some extent used for virus transmission and the “effective” pattern of virus shedding is somewhere in between the two considered shedding patterns. In fact, there are likely seasonal differences in the role of direct (breeding in summer) and indirect (non-breeding in winter) transmission and since the temporal change in virus presence is different in blood, feces, urine and saliva, it is to be expected that the shedding patterns too may be different for both routes.

The model used here improved on previous models by allowing complex spatial dynamics of the bank vole/virus interaction: it models the spatial behavior (dispersal and varying home ranges) of the bank voles explicitly, in an environment that is being contaminated with (and gradually cleared from) virus. Apart from the particular situation in Central-Finland, the model can be used for studying the distribution of PUUV throughout Europe. The model, attuned to (the best documented situation in) Central-Finland, can be used as a baseline to investigate to what extent the variation of different system (model) parameters alters the infection dynamics, and to identify processes responsible for the heterogeneous distribution of PUUV across Europe. Finally, the applicability of the model is not limited to PUUV, as it could easily be adapted to model the spatial transmission dynamics of other pathogens. It is especially suitable to study infections in which the variable spatial dynamics of the host are thought to play an important role in the transmission dynamics.

In this paper, we have advanced a novel spatially-explicit model to simulate PUUV infection dynamics in a population of bank voles. We have applied the model to the specific situation in Central-Finland, where PUUV infection occurrence is generally high and poses a significant threat to humans (Vaheri et al.,

The datasets analyzed in this article are not publicly available. Requests to access the datasets should be directed to

JR, KT, BB, and NH designed the model, based on input on bank vole and Puumala virus ecology provided by LV, HH, and HL. JR implemented the model, performed the simulations, and wrote the first draft of the manuscript. LV and HH provided field data. All authors contributed to, read, and approved the submitted version.

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

We thank Frank Sauvage for valuable discussions.

The Supplementary Material for this article can be found online at: