Could a rabies incursion spread in the northern Australian dingo population? Development of a spatial stochastic simulation model

Australia, home to the iconic dingo, is currently free from canine rabies. However northern Australia, including Indigenous communities with large free-roaming domestic dog populations, is at increased risk of rabies incursion from nearby Indonesia. We developed a novel agent-based stochastic spatial rabies spread model to evaluate the potential spread of rabies within the dingo population of the Northern Peninsula Area (NPA) region of northern Australia. The model incorporated spatio-temporal features specific to this host-environment system, including landscape heterogeneity, demographic fluctuations, dispersal movements and dingo ecological parameters—such as home range size and density—derived from NPA field studies. Rabies spread between dingo packs in nearly 60% of simulations. In such situations rabies would affect a median of 22 dingoes (approximately 14% of the population; 2.5–97.5 percentiles: 2–101 dingoes) within the study area which covered 1,131 km2, and spread 0.52 km/week for 191 days. Larger outbreaks occurred in scenarios in which an incursion was introduced during the dry season (vs. wet season), and close to communities (vs. areas with high risk of interaction between dingoes and hunting community dogs). Sensitivity analyses revealed that home range size and duration of infectious clinical period contributed most to the variance of outputs. Although conditions in the NPA would most likely not support a sustained propagation of the disease in the dingo population, due to the predicted number of infected dingoes following a rabies incursion and the proximity of Indigenous communities to dingo habitat, we conclude that the risk for human transmission could be substantial.


Probability of dying and probability of introduction of newly born dingoes
The study area, in northern Australia, is characterised by an equatorial climate, with a dry season (approximately May-October) and a wet season (approximately November-April). According to the camera-trap study results, the dingo density varies between the dry season (0.135 dingoes/ km 2 (CI: 0.127 -0.144)) and the wet season (0.147 dingoes/ km 2 (CI: 0.135 -0.159)). These values are interpreted as the mean density throughout each season (~6 months each). This temporal fluctuation of the population is incorporated into the rabies spread model by natural death and birth: 1. The daily probability of dying for each dingo is described by where δ is the death rate. It is assumed that the death rate is constant across all seasons (but varies between simulations).
2. The probability of introduction of new independent dingoes into the population occurs during a specific time period of the year, that is during the first half of the wet season. During this period of dingo introduction, the probability of introduction is set to be equal to a simple constant, as illustrated below: It is assumed that the population size is stable from one year to another. The number of days during the dry season is k, during the wet season is s and during the period of introduction is b.
On day 0 of the dry season, the population size is PD. At the end of the dry season, on day k, the population size is PD', which is also equal to the population size at the beginning of the wet season PW (day 0 of wet season). The population size at the midpoint of the dry season is PD_0.5. At the end of the wet season, on day s, the population size is PW', which is also equal to the population size at the beginning of the dry season PD. The population size at the end of the birth period is Pb, which corresponds approximately to the midpoint of the wet season.

DRY SEASON
On day 1, the population size is: On day 2, the population size is: By recursion, on day n, the population size is: The average daily density during the dry season must equal: The geometric series sums to: This yields the population size at day 0, at the beginning of the dry season.
Using equation (2) and (4), we can calculate the population size at the midpoint of the dry season, on day k/2 Similarly, the approximate population size on the last day of the dry season, on day k, becomes:

WET SEASON
Within the wet season, there are 2 different periods: 1) First part of the wet season (during b days): Period in which the probability of death AND the probability of introduction are occurring simultaneously. 2) Second part of the wet season (during s -b days): Period in which only the probability of death occurs (same as dry season).

Second part of the wet season (during s -b days)
Only death occurs during this period. Therefore, similar to the dry season, the population size on day n is: = − (7) And the population size at the end of the wet season is: The population at the last day of the wet season, PW', must equal to PD: ′ = Using equations (4) and (8), Using equation (7), the sum of daily densities during the second part of the wet season: On day 1, the population size is: Here, the first term on the right hand sider is the initial population of the dingoes subtracting the dingoes that died on day 1. The second term, CPb, is the number of dingoes born in day 1 of the wet season. On day 2, the population size is: By recursion, on day n, the population size is: The geometric series with index starting at zero converges according to: Since PW = PD' and using equation (2) The sum of daily densities during the first part of the wet season (from day 1 to day b) becomes: Using equation (3) for the geometric series, and equation (12) for C, yields The average daily density during the wet season must equal: Substituting the expressions for the sum densities, equations (10) and (13) yields The MeanD_wet may be expressed in terms of MeanD_dry by replacing the term Pb by equation (9).
For each simulation, the model randomly selects a value of the mean density during the dry season. This value is then multiplied by a factor of 1.0888889 to calculate the mean density during the wet season, as this corresponds to the ratio found between both mean density results from the camera-trap study. Using equation (14), with the selected and calculated values of mean densities during the dry and wet seasons, N = 365 days, s =180 days and b = 90 days, the model numerically calculates the value for the death rate, δ, in order to satisfy all above conditions and ensure that the population size remains stable from one year to another.

Home range distribution
For each dingo, based on its own selected 95% bivariate home range size estimate, a circular bivariate (2 dimensional) normal distribution will be fitted in view of obtaining a probability density function (referred to as home range distribution, in the manuscript). The probability density function for dingo i, ( , ), can be described by the following equation: We will consider the function to be centered at the origin xi = yi = 0. Since we want to describe a circular function, the correlation i = 0 and the standard deviations xi = yi =i. With these conditions, the probability density distribution simplifies to: The integral of ( , ) over the entire plane gives a probability of 1. The value of σ will be calculated by assuming that the probability of finding the dingo inside the home range estimate is 95%. Therefore, the probability of finding dingo i within a circle of radius R95, corresponding to the area of the 95% home range estimate is The standard error will be calculated for each dingo, based on its selected seasonal 95% home range value.

Probability of relocation in a new vacant home range area
We assume that dingo i will relocate to a new home range centroid if it found itself within an arc-shaped area of dimensions u around the location of the new home range centroid, and if the distance d between the dingo's home range centroid and the new location is larger than 95 . The probability of dingo i relocating is therefore calculated as the integral of the probability density function over the area formed by this arc, as followed: Here d equals the distance between the dingo i centroid and the new vacant centroid and equals the standard deviation of the home range distribution for dingo i. The probability of dingo i relocating ( i ) is based on the probability density function of the dingo. This function describes the probability of a dingo occurring at one location within a time frame which corresponds to the number of days, ni, needed for the dingo to explore most of its home range area (and therefore the time needed to build its probability density function). We will assume that relocation is an independent event from one day to another. Consequently, the daily probability of relocation therefore becomes: N.B. As the distance d increases to infinity, the arc-shaped area will slowly take the shape of a square of 4u 2 dimensions. u X u d Y Figure E. Representation of a very small subset of grid points and the network resulting from the connections between each point and its 32 nearest neighboring points in all directions. The connections represent the edges in which dingoes were allowed to travel from one point to another within the network.