Example description

  • Interacting age-structured population with genetically variable traits (from Ferrière and Tran (2009)).
  • Individuals characterized by their age $a \in [0,2]$ and size at birth $x_0 \in [0,4]$.
  • Two types of events :
    • Births
    • Deaths

Birth events

  • Each individual can give birth to an offspring, with an intensity


$$b(x_0) = \alpha (4 - x_0) \leq 4\alpha$$
  • New individual have the same size than his parent with probability $1-p$.
  • A mutation can occur with probability $p$ and then the birth size of the new individual is
$$x_0' = \min(\max(0, x_0 + G), 4),$$

where $G$ is a Gaussian random variable $\mathcal{N}(0,\sigma^2)$.

Death events

  • Size of an individual of age $a$ : $x_0 + ga$.

  • Competition : If a individual of size $x$ encounters an individual of size $y$, then it can die with the intensity $ U (x, y)$.

$$U(x,y) = \beta \left( 1- \frac{1}{1+ c\exp(-4(x-y))}\right) \leq \beta$$
  • The death intensity of an individual of size $x_0 + ga$ at time $t$ :


$$d(x_0,a,t,pop) = \sum_{(x_0', a') \in pop} U (x_0 + g a, x_0' + g a').$$

Population creation

  • Initial population is a dataframe with three columns:

    • birth (mandatory) : Date of birth of individuals in initial population (here the initial time is $t_0=0$).
    • death (mandatory) : Date of death (NA if alive).
    • birth_size corresponding to the size at birth.
# Generate population
N <- 900  # Number of individuals in the initial population
x0 <- 1.06
agemin <- 0.
agemax <- 2.

pop_init <- data.frame(
  "birth" = -runif(N, agemin, agemax), # Age of each individual chosen uniformly in [0,2]
  "death" = as.double(NA),
  "birth_size" = x0 # All individuals have initially the same birth size x0.
)
head(pop_init)
A data.frame: 6 × 3
birthdeathbirth_size
<dbl><dbl><dbl>
1-1.1390462NA1.06
2-0.1650035NA1.06
3-0.7701448NA1.06
4-1.1031752NA1.06
5-1.7765681NA1.06
6-0.2058391NA1.06

Model parameters

params_birth <- list( # parameters for birth event.
  "p" = 0.03,
  "sigma" = sqrt(0.01),
  "alpha" = 1)

params_death <- list(
  "g" = 1,
  "beta" = 2./300.,
  "c" = 1.2
)

Events creation

Birth event creation

  • Birth intensities : no interactions $\rightarrow$ created using function mk_event_individual.
  • Intensities functions have to be implemented using small C++ code chunks.
  • Facilitated by available parameters.
  • Individual is called I : age(I,t), I.birth_size....
  • Characteristics of individuals entering the population have to be specified in the kernel_code. New individual is called ... newI.
birth_event <- mk_event_individual( type = "birth",
  intensity_code = 'result = alpha * (4 - I.birth_size);',
  kernel_code = 'if (CUnif() < p)
                     newI.birth_size = min(max(0., CNorm(I.birth_size, sigma)), 4.);
                 else
                     newI.birth_size = I.birth_size;')

Death event creation

  • Death intensities : interactions → created using function mk_event_interaction.
  • Drawn individual I interacts with individual named J
death_event <- mk_event_interaction( # Event with intensity of type interaction
  type = "death",
  interaction_code = "double x_I = I.birth_size + g * age(I,t);
                      double x_J = J.birth_size + g * age(J,t);
                      result = beta * ( 1.- 1./(1. + c * exp(-4. * (x_I-x_J))));" #  U
)

Model creation

model <- mk_model(
  characteristics = get_characteristics(pop_init),
  events = list(birth_event, death_event),
  parameters = c(params_birth, params_death) # Model parameters
)
summary(model)

Events:
#1: individual event of type birth
#2: interaction event of type death
--------------------------------------- 
Individual description:
names:  birth death birth_size 
R types:  double double double 
C types:  double double double
--------------------------------------- 
R parameters available in C++ code:
names:  p sigma alpha g beta c 
R types:  double double double double double double 
C types:  double double double double double double

Simulation

Computation of event intensity bounds


$$b(x_0) = \alpha (4 - x_0) \leq 4\alpha$$


$$U(x,y) \leq \beta$$ 
birth_intensity_max <- 4*params_birth$alpha
U_max <- params_death$beta

Simulation

sim_out <- popsim(model = model,
  population = pop_init,
  events_bounds = c('birth'=birth_intensity_max, 'death'=U_max),
  parameters = c(params_birth, params_death),
  age_max = 2,
  time = 500)

Simulation on  [0, 500] 
sim_out$logs[["effective_events"]] / sim_out$logs[["proposed_events"]]
sim_out$logs[["duration_main_algorithm"]]
0.448581523713725
0.219941

Outputs

Population dataframe

  • Dataframe sim_out$population: all individuals present in the population over the period $[0, 500]$. 
pop_out  <- sim_out$population
tail(pop_out)
A data.frame: 6 × 3
birthdeathbirth_size
<dbl><dbl><dbl>
304712497.9870499.98702.754715
304713498.6642499.99192.564276
304714499.1359499.98852.578342
304715499.7084499.99612.578342
304716499.8724499.98932.675707
304717499.9892500.00072.948514

Age pyramids

pyr1 <- age_pyramid(pop_out, ages = seq(0,2,by=0.2), time = 500)

Change of parameters

params_death$g <- 0.3

sim_out <- popsim(model = model,
  population = pop_init,
  events_bounds = c('birth'=birth_intensity_max, 'death'=U_max),
  parameters = c(params_birth, params_death),
  age_max = 2,
  time = 500)

Simulation on  [0, 500]

Increase in initial population size

N <- 2000
pop_init_big <- data.frame(
  "birth" = -runif(N, agemin, agemax), # Age of each individual chosen uniformly in [0,2]
  "death" = as.double(NA),
  "birth_size" = x0 # All individuals have initially the same birth size x0.
)

params_birth$alpha <- 6
params_birth$p <- 0.01 # Mutation probability
params_death$beta <- 1/100
params_death$g <- 1

Simulation

birth_intensity_max <- 4*params_birth$alpha
U_max <-  params_death$beta

sim_out <- popsim(
  model = model,
  population = pop_init_big,
  events_bounds = c('birth'=birth_intensity_max, 'death'=U_max),
  parameters = c(params_birth, params_death),
  age_max = 2,
  time = 500)

Simulation on  [0, 500] 
sim_out$logs[["duration_main_algorithm"]]
4.904684