Permafrost dynamics and the risk of anthrax transmission: a modelling study

In this section, we first present the general formulation of our anthrax transmission model following both a deterministic and a stochastic approach. The latter seems particularly suitable in this case because of its ability to capture the dynamics of discrete events furthering pathogen spread, especially in the case of a small host population or episodic disease transmission. Then, we illustrate the methods used to derive conditions for the establishment of sustained disease transmission, in particular considering seasonal variations of environmental forcings and herding practice.

The anthrax transmission model

Our formulation builds on a compartmental model describing the epidemiological dynamics that affect a target population (composed of susceptible and infected individuals) being exposed to environmental contamination. We focus on domestic herbivores because they are both the most at risk and the most socio-economically valuable for Arctic communities. We thus neglect any direct interaction between infected carcasses and carnivores or scavengers, and consider environmental spores as the only source of infection (i.e. by ingestion while herds graze). As little is known yet on how age and sex of the animals may influence anthrax transmission, we suppose that all animals are equally vulnerable to the infection.

Let S(t), I(t) and R(t) be the total abundances of susceptible, infected and temporarily immune animals at time t, respectively, and let H be the total size of the animal population. Differently from previous formulations, we introduce two different reservoirs of spores. The first (with abundance (B_1 (t))) accounts for fresh spores that are released after the death of infected hosts and that are immediately available on the soil surface. As these spores infiltrate, are washed away or get buried, they enter the second reservoir (with abundance (B_2 (t))), which describes long-term storage in the soil. Anthrax transmission can thus be described by the following system of ordinary differential equations:

As for susceptible animals (Eq. 1), we assume that pastoralist practices keep herd size under controlled demographic growth, with (mu H) being the constant recruitment rate compensating the natural (non disease-induced) mortality occurring at rate (mu). Susceptible animals may become infected at a rate expressed by the total force of infection, F(t), which will be described in details later. A fraction (sigma) of animals that have been exposed to anthrax spores develops symptoms and enters the infected compartment I (Eq. 2). Once infected, symptomatic animals may die as a result of the anthrax infection at rate (alpha), or die for other causes not related to the disease at rate (mu). The remaining fraction of exposed animals ((1-sigma), e.g. animals that have been exposed to lower doses of spores) may not exhibit symptom and develop a temporary immunity³⁸. As these individuals do not shed spores, we assume that they enter directly the immune compartment R (Eq. 3). These animals lose their immunity and return to the susceptible class at rate (rho). When infected animals die of anthrax disease, spores proliferate in the host carcass. We assume that each death produces a constant number (theta) of spores, which are then released from the carcass and contaminate the surrounding environment, whose areal extent is A (Eq. 4). Both (B_1 (t)) and (B_2 (t)) are environmental densities of spores per unit area (# spores m(^{-2})). The freshly released spores decay at a rate (delta _1), or may be removed from the surface soil layer and stored in the active layer reservoir at a rate (chi) (Eq. 5). The latter rate thus conceptually encapsulates the combined effect of infiltration, runoff, and the burying of infected carcasses without appropriate sanitary precautions. Spores stored in the active layer decay at a rate (delta _2). The main processes involved in anthrax transmission dynamics have been conceptualized in Fig. 1.

The force of infection F(t) , which controls the rate at which susceptible animals get infected, depends on the concentrations of spores ((B_1), (B_2)) and the rate of exposure ((beta (t))) according to

where the fraction (B_i/(K+B_i)) (for (i=1,2)) is the probability of becoming infected after being exposed to a certain density (B_i) of spores, K being the half-saturation constant (i.e. the dose of spores for which infection risk is half of its maximum value). As mentioned before, because of the processes involving active layer thaw, including e.g. cryoturbation, soil cracking, and solifluction, we assume that spores (B_2) may become available to grazing animals. The exposure to spores (B_2) is therefore influenced by active layer thawing, which has a significant seasonal component. The interaction between thawing and the release of spores is expressed through the parameter (eta (t)), which quantifies the probability of being exposed to spores (B_2) relatively to that of being exposed to freshly released spores ((B_1)). Because all the processes mentioned above are more likely to occur with thawing, we assume the probability (eta (t)) to be proportional to the depth of the active layer. We will later relax this simple assumption and investigate more complex relationships. We initially mimic the annual cycle of active layer thawing with a simple sinusoidal function (which also simplifies stability analysis via Floquet theory), so that (eta (t)) can be expressed as:

with (epsilon _{eta }) indicating the maximum amplitude of seasonal fluctuations, i.e. the maximum probability for susceptibles to be exposed to spores (B_2) relative to spores (B_1). Note that the soil thaws only during the warmer months, during which susceptibles are potentially exposed to spores (B_2) ((eta (t)>0)). Later on, we model more realistically the annual cycle of active layer depth, relating it to a real record of surface soil temperatures via the Stefan equation^{39, 40}, so as to better analyze anthrax risk in the Arctic environment.

Herding practices and grazing activity might vary seasonally as well, favouring increased exposure during warmer months. Therefore, we set

where (beta _0) is the average value of (beta (t)), (epsilon _{beta }) is the maximum amplitude of seasonal grazing fluctuations, while (0le phi le 1) is the temporal lag between the phases of (beta (t)) and (eta (t)). Note that t is expressed in days.

Finally, to reduce the number of model parameters, we introduce the dimensionless spore concentrations (B^{*}_1={B_1}/{K}) and (B^{*}_2={B_2}/{K}) (see equations S1–S5 in the Supplementary Information). This substitution allows the aggregation of parameters (theta), A, and K into a single one, namely (theta ^{*}=theta /(A K)).

Stochastic formulation

To build a stochastic version of our anthrax transmission model, we rely on an extension of the classic exact stochastic simulator algorithm (SSA)⁴¹ that has recently been proposed to describe the Haiti cholera epidemic⁴². In the SSA, each animal is considered individually, i.e. the abundances of susceptible and infected animals are treated as discrete variables, ({mathscr {S}}(t)) and ({mathscr {I}}(t)). Accordingly, each individual experiences stochastic events (i.e. birth, death, infection, anthrax-related death, etc.; see Table 1) that occur at different rates, (e_k), where k indicates a generic event, depending on the state of the system. The overall occurrence of events is modeled as a Poisson point process whose rate e is defined as the sum of the rates of occurrence of all possible events, i.e.

The inter-arrival time between two subsequent events is thus an exponentially distributed random variable with mean 1/e, and the next event to occur is selected according to the probability (e_k/e)⁴¹.

The concentrations of anthrax spores, ({mathscr {B}}^{*}_1(t)) and ({mathscr {B}}^{*}_2(t)), are instead treated as continuous stochastic variables, because they are typically large enough to allow a continuous description. At each anthrax-related death event, ({mathscr {B}}^{*}_1(t)) undergoes a step increase of size (theta ^{*}), whereas between events spore concentrations are updated using the analytical solution of equations S4–S5 (see the Supplementary Information), with (theta ^{*}=0). In analogy with the deterministic formulation (Eq. 6), the force of infection reads

Finally, a Monte Carlo approach, in which many different trajectories (realizations) of the SSA are evaluated, is used to study the long-term behaviour of the stochastic formulation of the anthrax transmission model.

Derivation of disease transmission conditions

Linear stability analysis of time-invariant systems

Conditions for long-term pathogen invasion and sustained transmission (endemicity) are first derived in the absence of seasonal fluctuations. To that end, we consider the exposure rate and the probability to be infected by spores (B_2) to be constant over time (i.e. (beta (t)=const=beta _0) and (eta (t)=const=eta _0), respectively).

Endemic anthrax transmission is possible if the disease-free equilibrium (DFE), a state of system 1–5 where ((S,I,R,B_1,B_2)=(H,0,0,0,0)), is asymptotically unstable. Linear stability analysis is used to determine a threshold condition based on the basic reproduction number⁴³

Specifically, the DFE is asymptotically stable when (R_0<1) (hence, disease cannot spread in the long run), unstable otherwise (hence, endemic anthrax transmission is possible). The conditions for positivity and asymptotic stability of the endemic equilibrium (EE) of model 1–5 can also be evaluated based on Eq. 9. Specifically, the EE is endowed with strictly positive components and asymptotically stable if (R_0>1), unfeasible and unstable otherwise. Clearly, (R_0=1) represents a bifurcation point, where the two equilibria collide and exchange their stability (transcritical bifurcation). For further mathematical details, the reader may refer to the Supplementary Information.

In the absence of the long-term spore reservoir (B_2), the basic reproduction number ({tilde{R}}_{0}) reads:

This definition will become useful in the following section.

Periodic systems: Floquet analysis

Conditions for endemic anthrax transmission to occur in a seasonally forced environment can be studied by applying Floquet theory^36,37. The disease-free equilibrium of model 1–5 subject to periodic fluctuations is unstable when its maximum Floquet exponent, (xi _{max}), is positive (for further theoretical details see Supplementary Information). To compare transmission dynamics between periodic and time-independent conditions we calculated also ({overline{R}}_0) by assuming (eta (t)) and (beta (t)) to be constant and equal to their average value. For any parameter set, ({overline{R}}_0) provides information regarding the stability conditions of the system if temporal fluctuations of parameters were neglected.

Note that other parameters may vary seasonally: for instance, the spore transition rate (chi) and the decay rates of the spores stored in the two reservoirs may vary over time because of fluctuations in the environment, temperature, and freezing or thawing conditions. However, for the sake of simplicity, and also due to the lack of detailed information, in the following we limit our analyses on the coupled effect of (beta (t)) and (eta (t)) (Eqs. 8 and 7, respectively).

Model setting and data

Most of the model parameters have been estimated according to reference values proposed in the literature, as shown in Table 2. The average lifespan of domestic livestock varies widely (on average between 5 and 20 years, or even more⁴⁴), depending on the animal species and herding management. Since animals with shorter life expectancy are more likely to be infected (see Supplementary Fig. S1), we assumed an average mortality rate of 0.2 years(^{-1}) as a representative case, i.e. domestic cattle with an average lifespan of 5 years. Given high mortality rates among infected herbivores¹⁴, we assumed that about 70% of infected animals develop symptoms. The remaining 30% grow a temporary immune response, ensuring animal immunity for about 6 months³⁸. Typically, infection with anthrax bacterium leads symptomatic animals to death in about 14 days¹⁴. Then, as the spores are released from infected carcasses, we assumed they remain directly available for about 10 days before their removal from the soil surface^15,32. Spores can be viable for decades¹⁴, thus we assumed an average viability of 10 years. Finally, we assumed that the probability (eta) can vary between 0 and 1, thus implying that animals can be equally exposed to the two spore reservoirs during the periods of maximum thawing. The parameters (beta) and (theta ^{*}), which quantify overall exposure and contamination, respectively, are critical in determining transmission dynamics. However, the lack of suitable epidemiological records prevents a proper estimation. Therefore, we explored different combinations of these parameters to compare different scenarios for anthrax transmission dynamics and discuss the results. While the maximum exposure rate (beta) has an easily interpretable physical meaning, the parameter (theta ^{*}) has a less immediate interpretation. We therefore illustrate results in terms of the corresponding ({tilde{R}}_{0}) (Eq. 10), that is, the basic reproduction number of the simplified model that does not account for the long-term spore reservoir (B_2). All simulations have been run with a total population size of (H=10^4) animals.

Finally, we exploited real data on the current variability of climate and permafrost dynamics to investigate the relationship between warm years (and related deeper active layers) and the risk of anthrax outbreaks. To that end, we run model simulations using the stochastic formulation and realistic forcing. To obtain the latter, we exploited a 17-year-long (2002–2018) dataset of thawing depth available at the Samoylov monitoring site (Lena River delta, northern Siberia)⁴⁵ which we combined with records of soil surface temperature ((hbox {T}_{S})). We then modeled the yearly cycle of the active layer depth Z via the Stefan equation^39,40 according to which (Z=Esqrt{C_S}), where E is the edaphic factor taking into account soil properties and (C_S) is the cumulative soil surface temperature, calculated when the top soil temperature is above (0,^{circ }hbox {C}). The estimated value of parameter E is equal to (2.58,hbox {cm}^circ hbox {C}^{-0.5}), when Z is in cm and (C_S) in (^{circ }hbox {C}) (see Supplementary Fig. S2). To produce synthetic time series of active layer depth to be used in the simulations, we first calculated the mean annual soil surface temperature and fitted a normal probability distribution to the 17 records. We then produced 200-year-long time-series of daily soil surface temperature randomly sampling the mean annual temperature from the normal distribution and assigning an annual pattern obtained shifting the trajectory of the average year (i.e. the year whose daily values are the averages across the available record for that specific day). Soil surface temperature is then transformed into active layer depth using the calibrated Stefan equation. In each 200-year-long model simulation, we discarded the first 100 years, which were used as model spin-up period, and retained the last 100 years for analysis. We run 100 replicas of the process, thus obtaining a total of 10,000 years of simulated anthrax incidence. Note that a 100-year-long simulation should not be interpreted as a future projection for which the hypothesis of steady climate is hardly justifiable, but rather as a computational way to obtain a large sample of simulations exploring the current climate variability without the need to repeat the spin-up phase of the model.

In the analysis described so far, we assumed the probability of contact between animals and spores (B_2), i.e. the parameter (eta (t)), to be proportional to the active layer depth. This implicitly assumes that the underground concentration of spores is uniform. However, as the potential sources of spores are on the surface, a negative gradient of spore concentration for increasing depth could be expected. Mathematically, this can be mimicked assuming a saturating relationship between the probability (eta (t)) and the active layer depth Z(t) so that the marginal increase of risk associated with a unit increase of Z decreases with the depth of the active layer itself. We have therefore explored two scenarios: in the first one (case 1) we assumed a linear relationship, i.e. (eta (t)propto Z(t)); in the second (case 2) a saturating relationship of the type (eta (t)propto Z(t)/(Z(t)+Z_0)), where (Z_0) represents the semi-saturation depth, which has been set to 0.2 m. We then scaled (eta (t)) so that the maximum value for case 1 is equal to 0.2. Accordingly, (eta (t)) in case 2 has been scaled so that it has the same long-term mean of case 1.

Source: Ecology - nature.com