# Population with genetically variable traits

This document provides an example of usage of the package IBMPopSim, for simulating an interacting age-structured population with genetically variable traits, based on in [@MR2562651].

See vignette('IBMPopSim') for a detailed presentation of the package.

# Example description

We recall here the example 1 of [@MR2562651].

Individuals are characterized by their body size at birth $x_0 \in [0,4]$, which is a heritable trait subject to mutation, and by their physical age $a \in [0,2]$. The body size is an increasing function of age, and the size of an individual of age $a$ is

$$x=x_0 + ga,$$

where $g$ is the growth rate, which is assumed to be constant and identical for all individuals.

There are 2 possible events :

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

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

depending on a parameter $\alpha$ and its initial size. Smaller individuals have a higher birth intensity. When a birth occurs, the new individual have the same size than his parent with a high 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: Due to competition between individuals, the death intensity of an individual depends on the size of other individuals in the population. Bigger individuals have a better chance of survival, and if a individual of size $x_0 +ga$ encounters an individual of size $x_0'+ ga'$, then it can die with the intensity

$$U (x_0 + g a, x_0' - g a'),$$

where the interaction function $U$ is defined by

$$U(x,y) = \beta \left( 1- \frac{1}{1+ c\exp(-4(x-y))}\right).$$

The death intensity of an individual of size $x_0 + ga$ at time $t$ is thus the result of interactions with all individuals in the population (including himself)

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

# Population creation

# 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. ) get_characteristics(pop_init) ## birth_size ## "double"

# Events and model creation

There are 2 possible events :

• Birth (with or without mutation)
• Death

Each event is characterized by its intensity and kernel code, described below.

## Birth event with individual intensity

An individual of size $x_0 \in [0,4]$ gives birth at the age independent rate given by $$b(x_0) = \alpha (4 - x_0)$$ Since the intensity only depends on the individual's characteristics, the event intensity is of type individual.

With probability $p = 0.03$ a mutation occurs, and with probability $1 - p$, the offspring inherits its parent’s trait, $x_0$. In the case of a mutation, the new trait is $x_0' = \min(\max(0, x_0 + G), 4)$, where $G$ is a Gaussian r.v. with expectation 0 and variance $\sigma^2=0.01$.

The birth event is then an individual event of type birth, created as follows:

### Parameters

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

### Event creation

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 with interaction

The death intensity of an individual with trait $x_0 \in [0, 4]$ and age $a \in [0, 2]$ is given by:

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

where

$$U(x,y) = \beta \left( 1- \frac{1}{1+c\exp(-4(x-y))}\right) \in \left[ 0, \beta\right]$$

This event intensity depends on the interaction kernel $U$, and is of type interaction.

### Parameters

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

### Event creation

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

Event bounds

Bounds for the birth intensity and the death interaction function $U$ have to be computed.

birth_intensity_max <- 4*params_birth$alpha interaction_fun_max <- params_death$beta
T = 500 # Multithreading is NOT possible due to interaction between individuals sim_out <- popsim(model = model, population = pop_init, events_bounds = c('birth'=birth_intensity_max, 'death'=interaction_fun_max), parameters = c(params_birth, params_death), age_max = 2, time = T) ## Simulation on [0, 500]
sim_out$logs["duration_main_algorithm"] ## duration_main_algorithm ## 0.131719 # Outputs sim_out$population is a data frame containing the date of birth, death, and characteristics of all individuals who lived in the population over the period [0,500].

str(sim_out$population) ## 'data.frame': 293742 obs. of 3 variables: ##$ birth : num 498 498 498 498 498 ... ## $death : num NA NA NA NA NA NA NA NA NA NA ... ##$ birth_size: num 2.95 2.47 2.76 2.43 2.47 ... pop_out <- sim_out$population Population size at$t=500$. pop_size <- nrow(population_alive(pop_out,t = 500)) pop_size ## [1] 338 Result from [@MR2562651] can be reproduce from the simulation. For each individual in the population, we draw below a line representing its birth size during its life time. ggplot(pop_out) + geom_segment( aes(x=birth, xend=death, y=birth_size, yend=birth_size), na.rm=TRUE, colour="blue", alpha=0.1) + xlab("Time") + ylab("Birth size") <img src="inter_output2-1.png")  # Simulation with different parameters The model can be simulated with different parameters without being recompiled. ## Impact of aging velocity The ageing velocity has an impact on the distribution of birth sizes of over time. r params_death$g <- 0.3


Events bounds are not modified since they do not depend on $g$.

sim_out <- popsim(model = model, population = pop_init, events_bounds = c('birth'=birth_intensity_max, 'death'=interaction_fun_max), parameters = c(params_birth, params_death), age_max = 2, time = T) ## Simulation on [0, 500]
pop_out <- sim_out$population ggplot(pop_out) + geom_segment(aes(x=birth, xend=death, y=birth_size, yend=birth_size), na.rm=TRUE, colour="blue", alpha=0.1) + xlab("Time") + ylab("Birth size") <img src="inter_output3-1.png")  **Evolution of age pyramid by birth size** ?age_pyramid returns the age pyramid data frame of the population, by birth size at a given time. r pyr <- age_pyramid(pop_out, ages = seq(0,2,by=0.2), time = 500) head(pyr) ## age birth_size value ## 1 0 - 0.2 1.735217 2 ## 2 0 - 0.2 1.756296 0 ## 3 0 - 0.2 1.844502 0 ## 4 0 - 0.2 1.857982 1 ## 5 0 - 0.2 1.858391 0 ## 6 0 - 0.2 1.859076 0  The age pyramid can be plotted, with a visualization of the individuals birth size, starting by defining discrete birth sizes subgroups, and by assigning a color to each subgroup. pyr$group_name <- as.character(cut(pyr$birth_size+1e-6, breaks = seq(0,4,by=0.25))) head(pyr) ## age birth_size value group_name ## 1 0 - 0.2 1.735217 2 (1.5,1.75] ## 2 0 - 0.2 1.756296 0 (1.75,2] ## 3 0 - 0.2 1.844502 0 (1.75,2] ## 4 0 - 0.2 1.857982 1 (1.75,2] ## 5 0 - 0.2 1.858391 0 (1.75,2] ## 6 0 - 0.2 1.859076 0 (1.75,2] library(colorspace) lbls <- sort(unique(pyr$group_name)) # Attribution of a color to each subgroup colors <- c(diverging_hcl(n=length(lbls), palette = "Red-Green")) names(colors) <- lbls

?plot_pyramid allows the user to plot the age pyramid at a given time of a population composed of several subgroups, given a population data frame with a column named group_name (only needed for displaying several subgroups).

plot_pyramid(pyr, group_colors = colors, group_legend = 'Birth size')

<img src="inter_pyramid-1.png")


Due to the interaction between individuals, only bigger individuals survive at higher ages.

Several age pyramids at different times can be computed similarly at different times by calling ?age_pyramids.

r
pyrs <- age_pyramids(pop_out, ages = seq(0,2,by=0.2), time = 50:500)
pyrs$group_name <- as.character(cut(pyrs$birth_size+1e-6, breaks = seq(0,4,by=0.25)))

lbls <- sort(unique(pyrs$group_name)) colors <- c(diverging_hcl(n=length(lbls), palette = "Red-Green")) names(colors) <- lbls # library(gganimate) # anim <- plot_pyramid(pyrs, group_colors = colors, group_legend = 'Birth size') + # transition_time(time) + # labs(title = "Time: {frame_time}") # animate(anim, nframes = 450, fps = 10) ## Increase in initial population size We can do the same simulation with a bigger initial population. In order for the population size to stay approximately constant, the birth (resp. death) intensity are increased (resp. decreased). 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 <- 4 params_birth$p <- 0.01 # Mutation probability params_death$beta <- 1/100 params_death$g <- 1 The birth intensity bound and interaction function bound must be updated before simulation. birth_intensity_max <- 4*params_birth$alpha interaction_fun_max <- params_death$beta sim_out <- popsim( model = model, population = pop_init_big, events_bounds = c('birth'=birth_intensity_max, 'death'=interaction_fun_max), parameters = c(params_birth, params_death), age_max = 2, time = T) ## Simulation on [0, 500] pop_size <- nrow(population_alive(sim_out$population, t = 500)) pop_size ## [1] 543
ggplot(sim_out$population) + geom_segment(aes(x=birth, xend=death, y=birth_size, yend=birth_size), na.rm=TRUE, colour="blue", alpha=0.1) + xlab("Time") + ylab("Birth size") <img src="inter_output-1.png")  # Model with "full" simulation algorithm In the presence of interactions, the randomized algorithm (activated by default in ?mk_event_interaction) is much faster than the standard algorithm (named full) which requires to iterate through the population vector at each candidate event time. r # Comparison full vs random death_event_full <- mk_event_interaction(type = "death", interaction_type= "full", 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))));" ) model_full <- mk_model(characteristics = get_characteristics(pop_init), events = list(birth_event, death_event_full), parameters = c(params_birth, params_death))  sim_out_full <- popsim(model = model_full, population = pop_init_big, events_bounds = c('birth' = birth_intensity_max, 'death' =interaction_fun_max), parameters = c(params_birth, params_death), age_max = 2, time = T) ## Simulation on [0, 500] sim_out_full$logs["duration_main_algorithm"]/sim_out\$logs["duration_main_algorithm"] ## duration_main_algorithm ## 39.58549