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. 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.

In addition to iterate, the perturbate method performs the same procedure to iterate the population for a day, however, leaves age and development unchanged.

This method can be used for modelling external influences such as harvesting and population control, and also emmigration.

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

R > s <- size(vec)

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 perturbate 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 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^11,16. 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.

Figure 2. 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.

As a result, the output of the model is almost identical to the output of the original model (compare Figure 2(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 (Figure2(b)). However, the amplitude of the oscillations and the relative heights of the minor cycles in each oscillation change drastically at successive generations. Similar fluctuations were also reported by Nicholson (1954)¹⁷ in a laboratory culture of the Australian sheep blowfly (also presented in Gurney et al. 1983¹⁶ Figure 1).

Figure 3. 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 green: life expectancy of precisely 10 generations; solid red, teal, and purple: gamma-distributed life expectancy with mean 10 and standard deviation 1.25, 2.5, and 5 generations, respectively.

Age-structured host-parasite interactions

Another classic example of age dependency in population dynamics was proposed by 18 as a variation of the host-parasite interaction model of 19. 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 introduced age-structure to the host population and assumed that only juvenile hosts are targeted by parasites.

Both the original Nicholson-Bailey model and its age-structured version are implemented in Supplementary File 4 (case_studies.py). As shown in Figure 3(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 3(a) (solidlines).

Hastings (1984)¹⁸ also discusses the disruptive effect of age-dependent host survival in stability. This is evident in Figure 3(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 resembles the Hastings’ model with no age limit (solid lines in Figure 3(a) and the dashed line in Figure 3(b)). When mortality is considerably higher in older individuals (red line in Figure 3(b)) or a strict age limit is introduced (green line in Figure 3(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 1965–1966²⁰. 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.

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 4, the number of infectious cases with respect to the number of susceptible individuals, the S-I plane, (Figure 4(a)) and the time trajectory of the outbreak (Figure 4(b)) closely follow the data presented in Brauer et al. 2001²¹ Table 9.1.

Figure 4. The age-structured SIR model of the Great Plague in Elam.

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 5 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 red lines). In the opposite scenario, where the length of the infectious case is precise, the entire outbreak resolves rapidly as seen in the green 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 5. The effect of age-structure in modelling epidemics.

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