Study site and demographic data
The study was carried out in Volcanoes National Park, the Rwandan part of the Virunga massif, which is further shared with Uganda and the Democratic Republic of the Congo. We focused on habituated mountain gorilla groups monitored by the Dian Fossey Gorilla Fund’s Karisoke Research Center, often referred to as the Karisoke subpopulation. Since 1967, groups in this subpopulation have been followed on a near daily basis. Through the mid-2000s, the Karisoke groups generally numbered three but over the last decade, group fission events and new group formations resulted in an average of ten groups in the region (see42,43). During daily observations, detailed demographic data are recorded, such as group composition, birthdate and death date, group transfers (for further details see Strier et al.50). The data used for this study covers demographic data from 1967 to 2018 and includes 396 recognized individuals.
Epidemiological data
We obtained published data on four variables that control the disease dynamics of COVID-19 in humans, namely (a) the basic reproductive number (R0)34,35, (b) the infection fatality rate (IFR) based on estimates from China and Italy24,25,36,37, (c) the probability of developing immunity and (d) the duration of immunity37,38,39,41.
Stochastic projection model
We used the stochastic projection model proposed by Colchero et al.51, that models population dynamics for both sexes on fully age-dependent demographic rates. The model incorporates the yearly variance–covariance between demographic rates, while it accounts for infanticide as a function of the number of silverbacks (mature males > 12 years old) in the population51. Because of this relationship between infanticide and number of silverbacks, this source of mortality changes in time and cannot be assumed to be part of the infant mortality rate. To explore the extinction probability for the Karisoke subpopulation as a function of different diseases, we gathered information from the model on the proportion of individuals that died for each disease and the frequency of outbreaks (i.e., how often outbreaks occurred).
Demographic-epidemiological projection model for COVID-19
We constructed a predictive population model that combines the species’ baseline demographic rates with a model based on the susceptible-infected-recovered-susceptible (SIRS) framework. As the baseline demographic rates, we used the age-specific mortality and fecundity estimated by Colchero et al.51 for mountain gorillas (Karisoke subpopulation). We defined four epidemiological stages, namely (a) susceptible, (b) infected, (c) immune and (d) dead, each of which we further divided into a fully age-specific structure (Fig. 1). Based on recent research on COVID-19 on humans, we assumed that the dynamics of the model allowed for the recovered individuals to be divided into either susceptible or immune37,38,39,41. Furthermore, we incorporated the potential age-specific infection fatality rate (IFR) based on current estimates from medical and epidemiological research24,25,36,37, adjusted to the lifespan of the gorillas by means of the logistic function
$$qleft(xright)=frac{{q}_{M}}{1+{text{exp}}left[-0.2left(x-25right)right]},$$
(1)
where qM is the maximum infected mortality probability. Similarly, we modeled the probability of developing immunity as a function of the strength of the disease, which, based on recent research, we measured as mirroring Eq. (1) as
$$mleft(xright)=frac{{M}_{I}}{1+{text{exp}}left[-0.2left(x-25right)right]},$$
(2)
where MI is the maximum immunity probability (Fig. 2B).
To explore the potential impact of COVID-19 on the growth rate of the Karisoke mountain gorilla subpopulation, we varied four of the critical epidemiological variables, namely (a) the basic reproductive number, R0, from 0.5 to 6 (which helps to simulate factors such as increased group density, which may increase the likelihood of transmission), (b) the maximum infected mortality probability, qM = (0.3, 0.6) (Fig. 2A), (c) the immunity duration, TI to 1, 3, 6, and 12 months, and (d) the maximum immunity probability, MI, from 0.2 to 0.8 (Fig. 2B). As time units we used year fractions in half months (i.e., t1 − t0 = 0.5/12), which allowed us to simplify the model, based on current information on the average time of serial interval and incubation period in humans21. This implementation assumes that susceptible individuals could become infected at the beginning of the time interval, while infected individuals in time interval t would either recover (immune or susceptible) or die in t + 1.
The deterministic structure of the model implies that the number of individuals in each sex, age and epidemiological stage was given by the possible contribution from the other stages 1/2 month before. This is, the number of susceptible individuals of age x at time t is given by the difference equation
$$begin{aligned} n_{s,x,t} & = p_{x – 1} left{ {n_{s,x – 1,t – 1} + n_{i,x – 1,t – 1} left[ {1 – qleft( {x – 1} right)} right]left[ {1 – mleft( {x – 1} right)} right]} right} & quad + n_{{m,x – T_{i} ,t – T_{i} }} prodlimits_{{j = x – T_{i} :j > 0}}^{x – 1} {p_{j} – n_{i,x,t} } , end{aligned}$$
where the ns,x,t is the number of susceptible individuals of age x at time t, and subscripts i and m refer to infected and immune individuals, respectively. For simplicity of notation, we do not include a subscript for sex, although the model does distinguish between sexes. The probability px is the age-specific survival probability. Functions q(x) and m(x) are as in Eqs. (1) and (2). Similarly, the number of immune individuals at time t and age x are
$${n}_{m,x,t}={n}_{i,x-1,t-1}left[1-qleft(x-1right)right]mleft(xright)+sum_{{j:0le jle {T}_{i}wedge x-j>0}}{p}_{x-j}{n}_{i,x-j,t-j}.$$
We incorporated this mechanistic structure into a stochastic model, where all contributions from time t to t + 1 were drawn from binomial or Poisson distributions. For instance, the total new number of infected individuals, Ni,t, was obtained as a random draw from a Poisson distribution with expected value
$$Eleft[{N}_{i,t}right]={text{min}}left[{{R}_{0}N}_{i,t-1},{N}_{t}right],$$
where Nt is the total number of individuals in the study subpopulation. We then distributed randomly these individuals into different available ages and sex corresponding to the term ni,x,t, in the susceptible equation above. The number of newborns, Bx,t, at each age for which there were available females at time t was drawn from a binomial distribution with expected value
$$Eleft[{B}_{x,t}right]=left({n}_{s,x,t}+{n}_{m,x,t}right){f}_{x}$$
where fx is the age-specific average female fecundity rate and ns,x,t and nm,x,t refers to the number of susceptible and immune females, respectively, of age x at time t. The sex of each newborn was then determined by means of a Bernoulli draw with probability given by the proportion of males in the population. Thus, if the draw produced 1 for that individual, it became a male, and if 0 a female.
For each scenario, we ran stochastic simulations for 2000 iterations for 10 years and recorded the average number of individuals at each age–sex and epidemiological state at every month. We then ran long-term stochastic simulations for four scenarios with R0 = 3 and maximum immunity probability MI = 0.2. For these, we recorded also the number of subpopulations that went extinct at each month.
Source: Ecology - nature.com
