Population.jl

This is the standalone Julia library of the dynamically-structured matrix population model sPop2. This version implements both age-dependent and accumulative processes.

Installation

Just type this in Julia:

    using Pkg
    Pkg.add("Population")

Alternatively, one could clone or download and extract the development version from this GitHub repository, and use the package as follows.

    using Pkg
    Pkg.add("Absolute path to the Population.jl directory")

    using Population

Using the library

The following creates a pseudo-structured population with 10 individuals and iterates it one step with 0 mortality and an Erlang-distributed development time of 20 ± 5 steps.

    pop = sPop2(PopDataSto(), AccErlang())

    AddPop(pop, 10)
    pr = (devmn=20.0, devsd=5.0)
    size, completed, poptabledone = StepPop(pop, pr)

See section Usage examples for further examples.

References

• Erguler, K., Mendel, J., Petrić, D.V. et al. A dynamically structured matrix population model for insect life histories observed under variable environmental conditions. Sci Rep 12, 11587 (2022). link
• Erguler K. sPop: Age-structured discrete-time population dynamics model in C, Python, and R [version 3; peer review: 2 approved]. F1000Research 2020, 7:1220. link

Usage examples

Let’s begin with a canonical example. Arguably, the straightforward way to model insect development is to use an ordinary differential equations system with exponentially distributed transition times.

A deterministic ODE model of larva development

Let’s assume that experimental observations of larva development yielded a round figure of 20 days as average development time. This is usually translated to the differential equations reals as an instantaneous rate of α = 1/20.

    using DifferentialEquations

    function develop!(du,u,p,t)
        du[1] = -(1.0/20.0)*u[1]
        du[2] = (1.0/20.0)*u[1]
    end

    u0 = [100.0; 0.0;]
    tspan = (0.0, 100.0)
    prob = ODEProblem(develop!,u0,tspan)
    sol = solve(prob)

Expectedly, this implies that larvae begin developing by the time they emerge from eggs, and keep producing pupae at the same constant rate until a negligible number of larvae is left.

A deterministic population with Erlang-distributed accumulative development

An alternative representation with a structured population can be constructed. For this, we need to know not the rate of pupa production, but the average duration of the larva stage and its variation. Let’s assume that stage duration follows a gamma distribution (Erlang to be precise), which implies that all individuals race to develop out of the larva stage, while some are faster and even more of them are slower.

    pop = sPop2(PopDataDet(), AccErlang())

    AddPop(pop, 10.0)
    out = [0 10.0 0.0]
    for n in 1:50
        pr1 = (devmn=20.0, devsd=5.0)
        ret = StepPop(pop, pr1)
        out = vcat(out, [n ret[1] ret[2]])
    end

Accounting for stochasticity

With a simple modification, we can simulate how the dynamics vary when each individual is given a daily chance of completing its development based on the Erlang-distributed development time assumption.

    pop = sPop2(PopDataSto(), AccErlang())

    AddPop(pop, 10)
    outst = [0 10 0.0]
    for n in 1:50
        pr1 = (devmn=20.0, devsd=5.0)
        ret = StepPop(pop, pr1)
        outst = vcat(outst, [n ret[1] ret[2]])
    end

Distributions and assumptions

The sPop2 framework allows for age-dependent and accumulative development times.

Age-dependent	Accumulative

Combining multiple processes

Multiple processes can be added to an sPop2 to represent more complex dynamics. For instance, we can represent survival with a daily constant rate and development with an Erlang-distributed accumulative process. The processes are executed in the order they are added to the sPop2.

a = sPop2(PopDataDet(), AgeConst(), AccErlang())
AddPop(a, 100.0)
ret = [0 GetPop(a) 0.0 0.0]
for i in 0:100
    pr1 = (prob=1.0/60.0,)
    pr2 = (devmn=30.0, devsd=5.0)
    out = StepPop(a, pr1, pr2)
    ret = vcat(ret, [i+1 out[1] out[2][1] out[2][2]])
end

The above instruction represent a life stage with 0.017 daily mortality and 30 ± 5 days of development time.

Gonotrophic cycle with a negative impact of ovipositioning

In this section, we will attempt to reproduce Figure 1 in the original age-dependent sPop model. We will represent an adult female mosquito with a lifetime of 20 ± 2 days. The female enters a cyclic process of obtaining bloodmeal and egg development, which takes about 2 days. However, at the end of each ovipositioning, her expected lifetime decreases by 2 days. I assure you that the actual physiology is far more complicated than this!

In order to represent the lifetime of this female mosquito, we will employ 3 processes:

• A custom age- and ovipositioning-dependent mortality process
• An age-dependent gamma-distributed egg development process
• A dummy process to count the number of ovipositioning events

    # Declare an sPop2 population with deterministic dynamics,
    # and define three processes in this order: Mortality, Gonotropic cycle, Ovipositioning
    a = sPop2(PopDataDet(), AgeCustom(custom, AgeStepper), AgeGamma(), AgeDummy())

The 3^rd process is a dummy (AgeDummy), which does not affect the population or even does not posess a time counter. The 2^nd process is the regular age-dependent gamma-distributed development process (AgeGamma).

The mortality process, on the other hand, is defined with a custom function and uses the AgeStepper. We included this stepper in process declaration, because we require that the status indicator be an age counter, i.e. the number of steps the sPop2 is iterated is kept in the status indicator (see qkey below). We declare the custom function as the following.

function custom(heval::Function, d::Number, q::Number, k::Number, theta::Number, qkey::Tuple)
    devmn = 480.0 - (qkey[3] > 4 ? 240.0 : 48.0 * qkey[3])
    devsd = 0.1 * devmn
    pr = (devmn=devmn, devsd=devsd)
    k, theta, stay = sPop2.age_gamma_pars(pr)
    return sPop2.age_hazard_calc(sPop2.age_gamma_haz, 0, qkey[1], k, theta, qkey)
end

With this function, we override an internal mechanism used by all other processes to calculate the probability of exit (need this be due to mortality, development, or something else) from the sPop2 population. This generic functional form takes the following parameters.

• heval is a function to calculate the cumulative probability density of an exit event. It may use some of the other parameters to do so.
• d refers to the number of days (more correctly time steps) elapsed in sPop2.
• q refers to the fraction of development (used in accumulative processes).
• k and theta are the parameters of the stage duration (development time) distribution.
• qkey is a tuple of status indicators (age for age-dependent and development fraction for accumulative processes)

Of all these, we use the status indicators to calculate the mean and standard deviation of mortality based on the third (dummy) process.

    devmn = 480.0 - (qkey[3] > 4 ? 240.0 : 48.0 * qkey[3]) # qkey[3]: Ovipositioning
    devsd = 0.1 * devmn
    pr = (devmn=devmn, devsd=devsd)

Then, we calculate the parameters of the corresponding gamma distribution using the internal function

    k, theta, stay = sPop2.age_gamma_pars(pr)

and, with these, we calculate the cumulative probability of daily mortality

    sPop2.age_hazard_calc(sPop2.age_gamma_haz, 0, qkey[1], k, theta, qkey) # qkey[1]: Mortality

Overall, the script to model the dynamics is given below.

    function custom(heval::Function, d::Number, q::Number, k::Number, theta::Number, qkey::Tuple)
        devmn = 480.0 - (qkey[3] > 4 ? 240.0 : 48.0 * qkey[3])
        devsd = 0.1 * devmn
        pr = (devmn=devmn, devsd=devsd)
        k, theta, stay = sPop2.age_gamma_pars(pr)
        return sPop2.age_hazard_calc(sPop2.age_gamma_haz, 0, qkey[1], k, theta, qkey)
    end

    # Stochastic dynamics, Mortality, Gonotropic cycle, Ovipositioning
    a = sPop2(PopDataDet(), AgeCustom(custom, AgeStepper), AgeGamma(), AgeDummy())
    AddPop(a,1000.0)
    ret = [0 0.0]
    for i in 0:480
        pr1 = NamedTuple{}()
        pr2 = (devmn=50.0, devsd=10.0)
        out = StepPop(a, pr1, pr2, pr1)
        for (q,n) in out[3][2]
            # Please note that this step is essential!
            AddPop(a, n, q.key[1], 0, q.key[3]+one(q.key[3]))
        end
        ret = vcat(ret, [i+1 out[2][2]])
    end
    return ret

At each step, StepPop returns the following tuple

• Current size of the population.
• Number of individuals completing each process in the given order.
• An array of (status indicators - number of individuals) pairs for each process.

The script adds all individuals completing a gonotrophic cycle (out[3][2]) back to the population one by one. While doing so, their status indicators are updated manually.

• n individuals are added
• Mortality indicator is left unchanged (q.key[1])
• Egg development indicator is reset
• Ovipositioning indicator is incremented (q.key[3] + one(q.key[3]))

Case studies

Nicholson’s blowflies

Contributed by Sean L. Wu @slwu89.

Nicholson’s blowflies are a classic example of a time-delayed stage-structured population model, and was presented as a case study in the sPop paper. We reproduce that example here using Population.jl.

The model comprises five distinct life-stages and exhibits stable quasi-cyclic oscillations. Although the model was originally constructed using delay differential equations, Population.jl adheres well to the observed dynamics despite the discrete time step. We demonstrate both deterministic and stochastic versions of the basic dynamical model.

    using Population
    using Plots
    using Distributions
    using Statistics

First we define parameters of the model.

    # Define daily propensity of death
    death = [0.07, 0.004, 0.003, 0.0025, 0.27]

    # Define the expected duration of development (in days)
    develop = [0.6, 5.0, 5.9, 4.1, Inf]

    # Define the rest of the parameters
    A0 = 600.0
    q = 8.5
    tau = 24
    # tau is a scaling factor to extend from daily to hourly iterations for more accurracy

Next we define a function to run the simulation for 300 days, in either stochastic or deterministic mode. The 5 life stages are stored in the vector vec. Each life stage has an AgeConst process for simulating mortality, which is independent of age in the blowflies model. It also has an AgeFixed process to simulate maturation, which has a deterministic (fixed) distribution in the blowflies model. Therefore variation in the stochastic version arises from the mortality process and egg laying, which is assumed to be drawn from a Poisson distribution centered on the deterministic mean.

On each day we calculate the terms required for the mortality and development processes (pr1 and pr2). Then we use StepPop to iterate the population for one day, and we store the output containing the total population size, dictionary of individuals completing development, and dictionary of all developing individuals for that stage in the vector out. Then we pass individuals who have completed development to the next stage life stage. Finally we compute the number of eggs produced by adults and add them to the egg stage.

    # this function is a wrapper routine to simulate a deterministic or a stochastic time series
    function sim_blowfly(stoch)
        output = zeros(300*tau+1, 5)
        # initiate 5 stages of the population and store them in a vector
        vec = [stoch ? sPop2(PopDataSto(), AgeConst(), AgeFixed()) : sPop2(PopDataDet(), AgeConst(), AgeFixed()) for n in 1:5]
        # add individuals to the initial stage (the egg stage)
        AddPop(vec[1], stoch ? 500 : 500.0)
        # record the initial conditions
        output[1, :] = [GetPop(v) for v in vec]
        # vector to store the output each day
        out = Vector{Any}(nothing, length(vec))
        # iterate for 300 days
        for n in 1:300*tau        
            # for each stage, repeat the following
            for m in 1:length(vec)
                # iterate the stage for one time unit
                pr1 = (prob = death[m]/tau,)
                pr2 = (devmn = develop[m]*tau,)
                out[m] = StepPop(vec[m], pr1, pr2)
                # and pass on to the next stage (if not the adult stage)
                if m > one(m)
                    for (q,s) in out[m-1][3][2]
                        AddPop(vec[m], s, 0.0, 0.0)
                    end
                end
            end
            # calculate the number of eggs to produce
            eggs = q * GetPop(vec[5]) * exp(-GetPop(vec[5])/A0)/tau
            if stoch
                eggs = rand(Poisson(eggs))
            end
            AddPop(vec[1], eggs)
            # record the current state
            output[n+1, :] = [GetPop(v) for v in vec]
        end
        # return the results
        return output
    end

We run the deterministic model, and sample 100 trajectories from the stochastic model. We plot results from the last 200 days of simulation.

    out = sim_blowfly(false)

    xtick = 100:(1/tau):300
    plot(xtick, out[100*tau+1:end,5], label = "Adults", color = :blue, ylim = [0,6000])
    plot!(xtick, out[100*tau+1:end,1], label = "Eggs", color = :red)

    outs = [sim_blowfly(true) for n in 1:100]
    outs = cat(outs..., dims=3)

    outs_adult_qt = [quantile(outs[i,5,:], [0.025, 0.5, 0.975]) for i in axes(outs,1)]
    outs_adult_qt = transpose(hcat(outs_adult_qt...))

    outs_eggs_qt = [quantile(outs[i,1,:], [0.025, 0.5, 0.975]) for i in axes(outs,1)]
    outs_eggs_qt = transpose(hcat(outs_eggs_qt...))

    xtick = (0:(size(outs_adult_qt)[1]-1))/24.0

    plot(xtick, outs_adult_qt[:,3], color="#0571b0", fillrange=outs_adult_qt[:,1], lw=1, alpha=0.5, label=:none)
    plot!(xtick, outs_adult_qt[:,2], color="#0571b0", lw=2, label=:none)
    plot!(xtick, outs_eggs_qt[:,3], color="#d6604d", fillrange=outs_eggs_qt[:,1], lw=1, alpha=0.5, label=:none)
    plot!(xtick, outs_eggs_qt[:,2], color="#d6604d", lw=2, label=:none)
    xlims!(100,300)