sPop: Age-structured discrete-time population dynamics model in C, Python, and R

R

Before we begin modelling, we load the albopictus package in R, and define the survival function as described in Equation 1.

R > library(albopictus)
R >
R > death <- function(pop) {
R +     if (nrow(pop)==0)
R +         return(data.frame(mean=480, sd=48.0))
R +     mn <- 480.0 - (48.0 * pop$devcycle)
R +     mn[mn < 240.0] <- 240.0
R +     return(data.frame(mean=mn, sd=0.1*mn))
R + }

The function returns a data.frame with a desired mean and standard deviation. Next, we initiate a population by calling the initiation routine of the spop class.

R > vec <- spop(stochastic=TRUE, prob="gamma")

With this line, we construct a stochastic population model with the gamma distribution as the basis of survival and development. Setting prob to nbinom selects the negative binomial distribution instead.

In order to introduce the first batch of individuals, we use the add method.

R > add(vec) <- data.frame(number=1000)

By default, age, development cycle, and the duration of development will be set to zero for all individuals. These can be customised by supplying additional fields to the data.frame: age to set age, devcycle to set the number of development cycles, and development to set the number of iterations the current development cycle has taken.

We can directly access the population structure of the spop class to inspect the number of individuals grouped with respect to age, development cycle, and the degree of development. Here, we will use these information to calculate the mean and standard deviation of expected lifetime for each age-development group.

R > tmp <- death(vec@pop)

The following step iterates the population for one time-unit by using the iterate method.

R > iterate(vec) <- data.frame(dev_mean = 50,
R +                            dev_sd = 10,
R +                            death_mean = tmp$mean,
R +                            death_sd = tmp$sd)

By defining dev_mean and dev_sd, we opt to use the gamma distribution to describe the probability of development. Setting dev_sd to zero results in the gamma distribution being discarded and a fixed number of iterations (indicated by dev_mean) being assigned for development. Instead, setting dev instead of dev_mean and dev_sd results in a daily constant development probability. The same principles apply for the survival process, where we provide the mean and the standard deviation of the gamma-distribution for each age- development group as calculated by the death function. After each iteration, the age and degree of development of the population are updated and the total number of individuals completing development is recorded together with the detailed account of the corresponding age-development groups. We access these data using the developed and devtable methods, respectively. In addition, the dead method returns the number of dead individuals following an iteration.

R > d <- developed(vec)
R > add(vec) <- devtable(vec)

In this example, we assume a periodic development process; therefore, we introduce all the individuals completing development back to the population using the add method.

Finally, we read the total size of the population using the size method.

R > s <- size(vec)

In the R implementation, we also provide an accessory method, perturb, to perform the same functions as the iterate method without updating age or development. By using perturb, the structure and size of a population can be modified to model the impact of an intervention or migration. The effect is immediate and in addition to the continuing survival and development processes. For instance, perturb can be used to model population control. If a certain agent kills a fraction of a population upon delivery, perturb can be used by specifying death to remove the affected individuals from the population. Otherwise, if the effect is age-dependent, death.mean and death.sd can be specified to calculate the fraction to be removed. If there is a need to migrate a subset of a population (population A) to a different population (population B), perturb can be used by specifying an appropriate development process (dev for age-independent, dev.mean and dev.sd for age-dependent migration). This results in keeping a detailed record of the age and stage of development of the individuals removed from population A. Essentially, this is the list of individuals selected for migration, and it can be used as desired. A simple example of migration is given below.

R > A <- spop(stochastic=FALSE, prob="gamma")
R > B <- spop(stochastic=FALSE, prob="gamma")
R > add(A) <- data.frame(number=1000)
R > perturb(A) <- data.frame(dev=0.5, death=0)
R > add(B) <- devtable(A)

The above code results in the migration of 500 individuals from population A to B while freezing the age, development cycle, and development counters.

Python

The Python implementation of the population dynamics model can be imported from the albopictus package.

Python >>> from albopictus.population import spop

We begin by declaring the survival function as in Equation 1.

Python >>> def death(pop):
Python ...     if pop.shape[1]==0:
Python ...         return [480.0, 48.0]
Python ...     mn = 480.0 - 48.0 * pop[:,1]
Python ...     mn[mn < 240.0] = 240.0
Python ...     return [mn, 0.1*mn]

Unlike the R implementation, the population structure is stored in a numpy.ndarray with the following order of columns: age, development cycle, degree of development, and number. Although the initiation step is similar to the R implementation, a two-dimensional list or a numpy.ndarray should be supplied to intriduce batches of individuals to a population.

Python >>> vec = spop(stochastic=True,prob="gamma")
Python >>> vec.add( [ [0, 0, 0, 1000] ] )

By using the add method above, we introduce 1000 individuals with zero age, development cycle, and degree of development. The population structure is directly accessible, which enables us to calculate a different mean and standard deviation for the gamma-distributed development of each age-development group.

Python >>> tmp = death(vec.pop)

Next, we iterate the population for one time-unit using the iterate method.

Python >>> vec.iterate(dev_mean = 50,
Python ...             dev_sd = 10,
Python ...             death_mean = tmp[0],
Python ...             death_sd = tmp[1])

In order to read the total number of individuals completing development, the detailed account of the corresponding age-development groups, and the total size of the population, we access the developed, devtable, and size attributes of the spop class. Please note that these attributes are overwritten each time the iterate method is called. Here, we record the total number of developed individuals and the population size, and reintroduce the developed individuals to the population for the next round of development.

Python >>> d = vec.developed
Python >>> vec.add(vec.devtable)
Python >>> s = vec.size

The Python implementation of the iterate method accepts an additional logical indicator pause to prevent updating age and development. If this parameter is supplied and if it is false, the iterate method acts as the perturb method of the R implementation.

C

The C implementation of the sPop package is further optimised for speed. The source code resides in the albopictus package of Python, and it needs to be compiled with the GNU Scientific Library (version 2.1 or later). We begin by locating the package directory and compiling three source files into the object code. Assuming that the file name of our model is test_spop.c, we produce the executable with the following.

$ gcc -c -o ran_gen.o ${pkgdir}/ran_gen.c
$ gcc -c -o gamma.o ${pkgdir}/gamma.c
$ gcc -c -o spop.o ${pkgdir}/spop.c
$ gcc -I${pkgdir} -lgsl -o test_spop ran_gen.o
    gamma.o spop.o test_spop.c

where $pkgdir is a bash variable holding the package directory. In order to use the package, we need to include the following header files in test_spop.c.

C   #include "ran_gen.h"
C   #include "gamma.h"
C   #include "spop.h"

The first header file defines the routines required for random number generation, and the second one defines the routines for the gamma and negative binomial distributions. The last header file defines the spop population structure and the associated functions for initialisation, modification, and garbage collection.

Each age-development group is stored in the individual_st data structure,

C   typedef struct individual_st {
C     unsigned int age;
C     unsigned int devcycle;
C     unsigned int development;
C     sdnum number;
C   } individual_data;

where the age, development cycle, degree of development, and the number of individuals in each age-development group are stored in age, devcycle, development, and number variables in the same order. The sdnum is a union data structure holding an unsigned int for a stochastic population or a double for a deterministic population. The spop data structure holds an array of individuals (individuals), population size (size), the number of dead and developed individuals following an iteration (dead and developed, respectively), a detailed account of developed individuals (devtable), an indicator for the probability distribution of age dependence (gamma_mode), a logical indicator for a stochastic or a deterministic model (stochastic), and two counters to manage the dynamic size of individuals (ncat and cat).

C   typedef struct population_st {
C     individual_data *individuals;
C     sdnum size;
C     sdnum dead;
C     sdnum developed;
C     void *devtable;
C     unsigned char gamma_mode;
C     unsigned char stochastic;
C     unsigned int ncat;
C     unsigned int cat;
C   } *spop;

Following the procedure in previous sections, we begin implementing the model in test_spop.c by declaring the survival function in Equation 1.

C   void death(const individual_data *ind,
C              double *death_prob,
C              double *death_mean,
C              double *death_sd) {
C     (*death_prob) = 0;
C     (*death_mean) = 480.0 - (ind->devcycle > 4 ?
     240.0 : 48.0 * ind->devcycle);
C     (*death_sd) = 0.1 * (*death_mean);
C   }

Please note that the C implementation handles a single age-development group at a time; therefore, the survival function is redesigned accordingly.

Next, we initiate a stochastic model with the gamma distribution as the basis of survival and development using the spop_init function with the first parameter set to a logical true.

C   vec = spop_init(1,MODE_GAMMA_HASH);

The macro MODE_GAMMA_HASH refers to the optimised implementation of the gamma distribution. Alternatively, MODE_NBINOM_RAW and MODE_GAMMA_RAW refer to the unoptimised implementations of the negative binomial and the gamma distributions. Optimisation involves recording previously-used values in a hash table for reuse, however, is memory intensive and should be used with caution. Faster more efficient implementations of the probability distributions are the main concern for future releases.

Having initiated vec, we introduce 1000 individuals of zero age with the spop_add function.

C   spop_add(vec,0,0,0,1000);

spop_add accepts parameters in the following order:

In order to iterate the population for one time interval, we use the spop_iterate function.

C   spop_iterate(vec,
C                0,
C                50.0, 10.0,
C                0,
C                0,
C                0, 0,
C                death,
C                0);

spop_iterate accepts the following parameters in the given order:

Following each iteration, the list of age-development groups that completed their development is stored in the devtable variable of the spop data structure. In order to reintroduce these individuals back to the population, we use the spop_popadd function.

C   spop_popadd(vec,vec->devtable);

It is possible to obtain a summary output of the population structure by using the spop_print function, which takes the spop data structure as the only parameter. spop can be recycled by emptying its contents with the spop_empty function.

C   spop_empty(vec);

Finally, in order to clear the memory used by vec, we supply its address to the spop_destroy function.

C   spop_destroy(&vec);

Model output

All three implementations of the model are given in Supplementary File 1–Supplementary File 3 (test_spop.R, test_spop.py, and test_spop.c. The resulting distribution of the number of individuals completing a development cycle during the first 20 days of simulation is given in Figure 1.

Figure 1. Number of individuals completing a development cycle in 20 days.

Thick solid line indicates the mean trajectory from the deterministic simulation, while the thin solid line and the dotted lines indicate the median and the 95% range of the stochastic simulation output. The timestep for each iteration is one hour.

Five cycles of development are clearly seen from the figure, while the population survives for less than 15 days with the survival function defined in Equation 1. Blending of development cycles is apparent and progressive due to the uncertainty in the duration of development (50 hours on average with a standard deviation of 10 hours).

Nicholson’s blowflies

We begin with Nicholson’s Blowflies, a classic example of time-delayed stage-structured population model^16,19. The model comprises five distinct life-stages and exhibits stable quasi-cyclic oscillations. Although, originally the model was constructed using continuous time-delay equations, we will demonstrate that the sPop model adheres well to the observed dynamics and the implementation presented in the StagePop package¹⁶. Furthermore, we will present a stochastic version of the model, which helps to improve our understanding of the observed variation.

Both the deterministic and stochastic versions of the model are presented in Supplementary File 4 (case_studies.py). The adaptation assumes fixed daily survival and strict development durations, values of which are the same as the original model (Figure 2 in Gurney et al. 1983¹⁹). A scaling factor is introduced to calculate hourly instead of daily propensities to improve accuracy on a par with the continuous-time simulations.

As a result, the output of the model is almost identical to the output of the original model (compare Figure 3(a) with the Figure 3a of Gurney et al. 1983¹⁹). The six peaks shown between generations 100 and 300 are matched by the stochastic version of the model (Figure3(b)). Evidently, the highest variability in stochastic simulations is in peak amplitude, whereas the frequency of the oscillations is well conserved. During each peak, two sub-peaks are observed, which are separated by a trough. The heights of each peak and trough largely vary suggesting that the two peaks may not be resolved in the observations of natural populations. Nevertheless, similar fluctuations were observed in the laboratory culture of the Australian sheep blowfly reported in Nicholson 1954²⁰ (also presented in Gurney et al. 1983¹⁹ Figure 1).

Figure 3. Age- and stage-structured model of Nicholson’s Blowflies in discrete time.

In (a), the deterministic trajectories of eggs (dashed line) and mature adults (solid line) are shown for the corresponding generations. In (b), five stochastic trajectories of mature adults are shown for the same duration.

Age-structured host-parasite interactions

Another classic example of age dependency in population dynamics was proposed by 21 as a variation of the host-parasite interaction model of 22. The Nicholson-Bailey model considers dynamics in discrete generations where parasites traverse a given area in search of a host. As a result, the number of parasites in the subsequent generation corresponds to the number of hosts parasitised. Hastings used this model for prey-predator interactions where he represented preys as hosts and predators as parasites. He introduced age-structure to the prey population and assumed that only juvenile preys are targeted by predators. Following Hastings, we use the terms host and parasite as analogous to prey and predator, respectively, regarding the emerging dynamics.

Both the original Nicholson-Bailey model and its age-structured version are implemented in Supplementary File 4 (case_studies.py). As shown in Figure 4(a), the original model without age dependency (dashed lines) exhibits oscillations with increasing amplitude around an unstable steady state. The model output matches Figure 10 in Nicholson (1935)²². Introducing age-structure with a fixed survival rate for host results in the stabilisation of dynamics as reported by Hastings (1984)²¹ and seen in Figure 4(a) (solid lines).

Figure 4. Age-structured host-parasite interactions of Nicholson-Bailey-Hastings.

In (a), the Nicholson-Bailey model (dashed lines) is compared with its age-structured version (solid lines) where parasites choose host in an age-dependent manner. The number of parasites is scaled down to 25% to aid visualisation. In (b), the effect of age-dependent host survival on stability is shown. Mortality rate is given in the inset with respect to age (the number of generations). Dashed line: age-independent mortality; solid orange: life expectancy of precisely 10 generations; solid green, red, and purple: gamma-distributed life expectancy with mean 10 and standard deviation 1.25, 2.5, and 5 generations, respectively.

Hastings (1984)²¹ also discusses the disruptive effect of age-dependent host survival in stability. This is evident in Figure 4(b), where the survival imbalance between young and old individuals drives the dynamics away from the steady state, eventually rendering it unstable. When the age-dependent mortality curve is close to being horizontal, the dynamics closely resemble the Hastings’ model with no age limit (solid lines in Figure 4(a) and the dashed line in Figure 4(b)). When mortality is considerably higher in older individuals (green line in Figure 4(b)) or a strict age limit is introduced (orange line in Figure 4(b)), the stability is lost and the amplitude of oscillations increase in time.

The Great Plague in Eyam

The final case we study is the severe outbreak of bubonic plague in Eyam (Sheffield, UK) in 1665–1666²³. The outbreak was initially modelled by Raggett et al. 1982²³, and later, a simple deterministic susceptible-infectious-removed (SIR) epidemic model was developed by Brauer et al. 2001²⁴ to study the outbreak.

Although canonical SIR models are not age-structured, the underlying transmission dynamics is intricately time-dependent. This becomes clear when, for instance, the (intrinsic) incubation period is observed to be strictly more than a couple of days but not less. While a system of ordinary differential (or difference) equations describes an uninterrupted flow between the incubation period and the infectious state, the sPop age-structured model is able to track the length of the incubation period for each individual and trigger state change only for the right individuals at the right times.

An age-structured version of the SIR epidemic model where the infectious stage duration is modelled with a gamma distribution is presented in Supplementary File 4 (case_studies.py). As shown in Figure 5, the number of infectious cases with respect to the number of susceptible individuals, the S-I plane, (Figure 5(a)) and the time trajectory of the outbreak (Figure 5(b)) closely follow the data presented in Brauer et al. 2001²⁴ Table 9.1.

Figure 5. The age-structured SIR model of the Great Plague in Eyam.

In (a), the simulated numbers of infectious versus susceptible individuals are shown together with the outbreak data. In (b), the number of susceptibles and infectious cases are plotted with respect to time for the duration of the outbreak. Model simulations (solid lines) are compared with the outbreak data (dots).

Figure 6 demonstrates the effect of age-structure on outbreak dynamics by comparing model output with different characteristics of the infectious period. The blue lines indicate the trajectory matching the Great Plague in Eyam. Please note that the corresponding infectious period is only slightly time-dependent where there is a minor difference between the mortalities of newly infected and long-time infectious individuals (plotted in the inset). If the rate of exit from the infectious stage is completely independent from the duration spent in the stage (as is the case with canonical SIR models), the outbreak duration increases (the green lines). In the opposite scenario, where the length of the infectious case is precise, the entire outbreak resolves rapidly as seen in the orange trajectory. Since time-dependence has a significant impact on outbreak trajectory, using realistic mathematical formulations to infer outbreak parameters, such as the incubation period and the rate of infection, becomes a critical step in developing predictive epidemic models.

Figure 6. The effect of age-structure in modelling epidemics.

The outbreak trajectory shown in Figure 5(b) (blue lines) is compared with alternative forms of mortality. Green lines indicate time-independence where the gamma distributed infectious period has σ = µ. Orange lines indicate that the infectious period has a precise length of µ = 11.71 days (σ = 0).