DEB-IBM for Abatus cordatus

1 collaborator

Anse du Halage version (with Delille1989 as resources)

ODD

1. Purpose

The model was developed from the individual mechanistic Dynamic Energy Budget (DEB) model to extrapolate the response of the populations of the sea urchin Abatus cordatus, endemic to the Kerguelen Plateau (sub-Antarctic region), to changes in environmental conditions (temperature and resources), through comparisons between sites and between predicted future scenarios.

2. Entities, State variables & Scales

DEB notation ..... equivalent in the code ..... (dimension) signification

U_H ...... UH ..................... (t·L²) scaled maturity ∂U_H ..... dUH ................... change of scaled maturity in time U_E ...... UE ..................... (t·L²) scaled reserves ∂U_E ..... dUE ................... change of scaled reserves in time e ........ escaled ................ (-) scaled reserve density per unit of structure U_R ...... UR ..................... (t·L²) scaled energy in reproduction buffer ∂U_R ..... dUR ................... change of energy in reproduction buffer L ......... L ......................... (L) structural length ∂L ....... dL ........................ change of structural length in time L_w ....... Lphy .................... (L), physical length (L/∂_M) L_b ....... Lb ...................... (L) structural length at birth

̈q ......... qacceleration ...... (t^-2) ageing acceleration ∂̈q ....... dqacceleration ..... change of ageing acceleration in time ḣ ......... hrate ................. (t^-1) specific death probability rate ∂ḣ ....... dhrate ................ change of hazard rate in time ḧ_a ....... h_a ...................... (t^-2) Weibull ageing acceleration s_G ....... sG ....................... (-) Gompertz stress coefficient

S_A ...... SA ...................... (L²) assimilation flux (scaled) S_C ...... SC ...................... (L²) mobilisation flux (scaled)

g ......... g ......................... (-) energy investment ratio ̇v ......... vrate .................. (L·t^-1) energy conductance κ ......... kap ...................... (-) fraction of mobilised reserve allocated to soma κ_R ....... kapR .................. (-) reproduction efficiency ̇k_M ....... kMrate .............. (t^-1) somatic maintenance rate coefficient ̇k_J ........ kJrate ............... (t^-1) maturity maintenance rate coefficient U_H^b ..... UH^b .................. (t·L²) scaled maturity at birth U_H^p ..... UH^p .................. (t·L²) scaled maturity at puberty

y_VE ..... yVE ..................... (mol/mol) yield of structure on reserve s_M ....... sM ..................... (-) acceleration factor L_m ....... zoom ................... (L) maximum volumetric length

Add-My-Pet parameters ('amp parameters' in interface) [̇p_M] .... pM ...................... (e·L^-3·t^-1) specific volume-linked somatic maintenance rate E_H^b .... EH^b ................... (e) maturity at birth
E_H^p .... EH^p ................... (e) maturity at puberty [E_G] .... EG ...................... (e·L^-3) volume-specific costs of structure

The model includes two types of entities: individuals and the environment.
Individuals are divided into 4 types of sub-agents, depending on life-stage and sex: embryos, juveniles, adult males and adult females. Life-stages are considered using DEB definitions. Individuals are characterized by four primary state variables, based on DEB theory (Kooijman 2010): scaled reserve (U_E), structural length (L), scaled maturity (U_H) and scaled reproduction buffer (U_R). Additional state variables for individuals are age (in months and in years), ageing acceleration (̈q) and death probability rate (ḣ). State variables in DEB are originally dependent on energy (unit of Joules), but in order to simplify the calculcations of differential equations so that they do not require measurements of energy, each state variable was discarded of their energy unit by dividing with the maximum surface-area specific assimilation rate {ṗ_Am} following Martin 2010 and Kooijman 2010. Females have supplementary attributes linked to reproduction processes for which the state variable is the gonado-somatic index (GSI). Finally, each individual has a variable called scatter-multiplier, used to implement a slight variation in three standard DEB parameters (U_H^b, U_H^p, g) and initial energy reserve at birth (UE0yo). This scatter-multiplier is the exponential of a number taken randomly on a normal distribution of mean 0 and standard deviation cv (set by the user in interface).

The environment in the model is characterized by two state variables, monthly average temperature (T, unit: °C) and monthly resources availability represented by the proportion of food an individual can intake on a scale of 0 to 1 (f, no dimension). Values for these variables are input into the model from external files, as time-series of monthly values. Temperature data was collected from existing thermo-recorders on the corresponding sites (implemented by the PROTEKER programme), while resources data comes from the publication by Delille & Bouvy 1989 for the site Anse du Halage. The variable f can be affected by the population density on site, due to intra-population competition for food.

The model as given here runs on a monthly temporal resolution over a temporal extent of 210 years which can be modified by the user in the interface. In other words, changes are applied to individuals and the environment on a monthly basis and thus each update corresponds to the state of the system at the end of the diplayed month: one timestep in the model represents one month and simulations are run for a total of 210 years. The first ten years are assumed to be the model initialisation phase and should be removed for the analysis of results. The model is non spatial, but can be adapted for spatialization. In this implementation, the model runs on one single patch of environment representing 1 square meter, and thus density of population is equal to number of individuals present in the model. Movements of individuals are not taken into account, and each individual born on the patch grows and dies on that same patch.

3. Process overview & Scheduling

At each timestep, the model runs the following commands in that order:

Reset the death counts
Update calendar

(For each patch:

    Update environmental variables
    If competition ON [
        Calculate competition 
    ]
    Calculate f 
)

(For each individual:

    Remove if marked as deceased 
    Convert relevant parameters with temperature correction factor
    Calculate change in reserve
    If not mature [
        Calculate change in maturity
    ]
    If mature [
        Calculate change in reproduction buffer 
    ]
    Calculate change in structural length
    If scaled reserve < scaled length [ 
        Starve
    ]
    Calculate ageing
)

Update individuals

(For each individual:

    Update reproduction timers
)

(For each females:
    Update birth timer
    If first month of reproduction period [
        If reproduction ON [    
            Mark GSI down 
            Prepare eggs
            ]
    ]
    Calculate GSI
    If GSI >= 0.07 [
        Turn ON reproduction
        If month before reproduction period [
            Mark U_R down 
            Launch reproduction (with birth_time)
            ]
    ]
    If GSI < 0.07 [ 
        Turn OFF reproduction
    ]
    If reproduction ON [
        If within reproduction period [ 
            reproduce
            ]
        If within birth-giving period [
            release offsprings
            ]
    ]
)

Background mortality
Check temperature
Monitoring of population
Update time

Individuals execute the same command one by one in a fixed order before going to the next command, but all have access to the same state of the environment since it updates once at the beginning of the timestep only.

In this model, time is represented continuously using ordinary differential equations (ODE) for the individual state variables, and all other variables are calculated in a discreet manner (every month, with one month rounded to 30.5 days).

Below is an illustration of the different calendars and timers used in the model, beginning at the point where the simulation starts (Gregorian calendar for reference):

repro_time:      7     8     9    10     11     0     1     2     3     4     5     6,     7     8 .......
GSI_time:        5     6     7     8     9     10    11    12     1     2     3     4,     5     6 .......
birth_time:      1     0     0     0     0      8     7     6     5     4     3     2,     1     0 .......
month_time:      10    11    12    01    02    03    04    05    06    07    08    09,    10    11 .......
Gregorian calendar: Oct   Nov   Dec   Jan   Feb   Mar   Apr   May   Jun   Jul   Aug   Sep,   Oct   Nov.......

'reprotime' follows the reproduction year, runs on a twelve steps loop, 'GSItime' follows the GSI cycle year, runs on a twelve steps loop, 'birthtime' is a countdown tracker for the time between the start of reproduction and the following release of offpsrings in a year. It is triggered by the launching of reproduction and counts backward, staying at 0 if not triggered. 'monthtime' follows the Gregorian calendar

4. Design concepts

Basic Principles

The IBM was built using the DEB-IBM model developed by B. Martin for Daphnia magna, along with its DEB-IBM user manual and model description (Martin et al. 2010). The underlying theory for the individual development in the model follows the Dynamic Energy Budget theory (Kooijman 2010). The population is studied following the IBM principles (Railsback & Grimm 2019), as a dynamic system composed of autonomous individuals affected by the environmental conditions throughout their life-cycle. Each individual undergoes a continuous development from birth till death, following the DEB principles with a slight variation between individuals at their initialisation, and represents a component of the IBM population, which is itself affected as a whole by variables such as population death rates and density-dependent processes. The emerging state of the population is then observed, and compared between different scenarios of environmental variations.

Emergence

The model illustrates the evolution of the population structure following the response of the individuals to the environmental conditions input. Metabolic responses, life-stages, ability to reproduce, starvation and ageing processes of the individuals, emerging from the mechanistic representation of their development, affect the population structure and average characteristics. A background mortality rate and a mortality caused by above normal temperatures are forced into the model, and the same reproductive output is imposed to all females that are able to reproduce.

Adaptation

Agents do not have an adaptive behavior. Individual traits vary among individuals in a population, but each individual carry the same traits along their entire lifespan and do not change nor learn from the events they experience or from each other. Consequently, the design concepts "objectives", "learning", "prediction", and "sensing" do not apply.

Objectives

Learning

Prediction

Sensing

Interaction

Individuals do not have any direct interaction. They only affect each other indirectly, as the size of the population influences the resources availability and thus the capacity of each individual to ingest a maximum proportion of food.

Stochasticity

In the model, stochasticity is used in the ageing submodel: there is a 50% chance that the individual will look into their probability of dying. This stochastic element can be modified in the code by changing the numbers x and y in the 'update individual' procedure. Stochasticity is also implemented in four of the initial variables for each individual (scaled maturity at birth, scaled maturity at puberty, energy investment ratio, energy reserve at birth), using the scatter-multiplier, a number randomly chosen on a log-normal distribution (taken from Martin 2010, Kooijman 2010).

Collectives

The individuals are grouped under a particular type of entity depending on their life-stage and sex, and will move between groups in one direction only. Thus, individuals are born into the "juveniles" group and move onto either the "males" or "females" group after reaching puberty. The age at which a juvenile reaches puberty is an emergent property of its development, but the sex is arbitrarily and randomly imposed on the individual that becomes an adult so that the sex-ratio of the population is around 0.99. Depending on which group they belong to, some variables are different: juveniles do not modify their reproduction buffers, while males and female do not modify their maturity compartment, and females possess some proper variables such as GSI (Gonado-somatic index), eggs (number of eggs produced) and Ri (reproductive output). These collectives do not emerge from individual behaviour, but instead are implemented by the modeller in order to distinguish the life stages and sex of the individuals.

Observations

The main output of the model are plots of the population structure with densities of population in the different life-stages, of the cumulative counts of individuals dead of different causes, as well as of mean values in state variables U_R, U_E and L and change in those state variables (∂U_R, ∂U_E and ∂L) for the different type of agents. These plots allow to observe the response of the population to the different environmental conditions and the metabolic responses of the individuals in the population in relation to those environmental conditions. Additionally, plots of the mean age at death of individuals dying due to the ageing submodel were used to calibrate the submodel itself.

NB: The plots are built on the lists compiled during the calc-pop-varb procedure (called by the update-pop procedure at the end of each month), which are built so that in case none of the individuals possessing the variable is present, an outlier mock value is applied instead. This is used to facilitate processing of the outputs afterwards.

Details

5. Initialisation

For the standard model, the following elements must be selected in the interface:

Sites: 'Anse du Halage'
projection: 'present'
future: 'mixed temp & food' (not used under 'present' projection)
sensitivity: 'resistant'
competition: 'ON'
run_time: '210'
cv: '0.1'

The initial DEB parameters can be calculated by the model if the 'add-my-pet?' switch is set to ON in the interface and the basic DEB parameters [̇p_M], E_H^b, E_H^p, [E_G] and L_m for the species as taken from the Add-my-Pet database are input into the relevant boxes (respectively: pM, EH^b, EH^p, EG, zoom).

The standard model is run for 210 years in total for the site Anse du Halage under present-day conditions, with a population sensitivity set to 'resistant' and competition affecting resources availability. At setup, lists are made out of the input files of temperature and f data for the site. The model is initialized for the month of October (monthtime 10) and the environment is setup with the corresponding environmental variables: if the model is used for a simulation under present-day conditions, the values are taken as is. If the model is set for future projections, the values are modified according to the chosen scenario (i.e. either one of RCP 2.6 and RCP 8.5 with food only, temperature only or food and temperature combined). Values for calculations of temperature correction factor are given as Tref = 293.15 and T_A = 9000. The carrying capacity is set at 200 ind./m² and the proportion of females at 0.5.

Initial parameters are based on DEB parameters, taken from the individual mechanistic DEB model for A. cordatus developped by C. Guillaumot in 2019 and published on the add-my-pet database. Two simulations procedures are run for the initialisation of individuals: a simulation of embryonic development to determine the initial reserve at birth, and a simulation of the development of one individual from 0 till 5 years old at constant f and temperature. The first simulation uses a bisection method to go through possible embryonic developments depending on the initial energy invested into the egg at conception. The simulation runs with the aim to obtain an embryo with a scaled reserve e close (within ±5%) to the reference scaled reserve when it reaches birth. The initial cost is taken as the mean of an upper and lower bound, which are arbitrarily set in the first loop and then determined by the result of the simulated development at the end of the loop: if the scaled reserve e is within the range of the reference e before the embryo reaches birth, the value that was taken for the initial cost is used as the lower bound in the next loop; if the scaled reserve e reaches the needed value after the embryo has reached birth, the value that was taken for initial cost is used as the upper bound in the next loop. When the embryo reaches birth with a reserve density within the aimed range, the simulation stops and the reserve density is set aside. The second simulation starts off where the first one finished, using the resulting reserve density at birth. It runs a loop for the development of the individual from birth till five years old, with standard parameters and a constant functional response f = 1. The simulation keeps track of the age of the individual, and for each year the values for the state variables are set aside and the simulation continues until the next year. Once these two simulations are run, the model creates the initial population of 120 individuals and first gives all of them initial parameters for individuals at birth (0 years old) with initial values from the DEB model slightly altered by a scatter-multiplier. Then the total number of individuals is divided in equal proportions into six classes corresponding to six age classes from 0 to 5 years old. Initial parameters are modified for individuals from classes 1 to 5 years old, by putting the initial parameters corresponding to their age that were set aside from the second simulation. All individuals are given an age in months and years, and a breed corresponding to their life-stage or sex (juveniles from 0 to 2 years old, female or male for 3 to 6 years old). Each individual sets their calendars with GSItime at 5 and reprotime at 7, and females are given initial values of 0.03 for their GSI and set their birth_time timer at 0.

Initialisation differs from run to run in one aspect only, the individuals' own parameter values: each individual created at initialisation gets its parameters slightly altered by a stochastic element, the scatter-multiplier. This scatter-multiplier is the exponential of a number taken randomly on a normal distribution of mean 0 and standard deviation cv (set by the user in the interface of the model, at 0.1 for the standard model).

Temperature data for two other sites is available: Port Couvreux and Ile Haute. In order to use the model with those sites, they must be selected in 'Sites' on the interface. NB: Keep in mind that even if these sites are selected, the food resources data will be that of Anse du Halage, as detailed data on food availability for other sites are not available. This is merely to test the model with a different kind of temperature data.

6. Input data

The model reads environmental variables from input files containing time-series data on monthly resources availability (from Delille & Bouvy 1989) and monthly average temperatures (PROTEKER programme IPEV n°1044). The two files are text files containing an ordered list of values (see below for Anse du Halage data). The temperature files contains 72 values of monthly average temperatures corresponding to temperature records from october 2012 to september 2018. The resources file contains 12 values, taken from the measures of organic carbon content in sediments published in Delille & Bouvy 1989. From each files, two ordered lists are compiled by the model during the setup procedure: one list with values from 0 to the total number of values, used as a timestamp, and one list with the corresponding value. Then the lists are used in parallel to take out the value found at the same level as the timestamp corresponding to the model's timing.

Data for resources and temperatures at Anse du Halage (to use, copy onto two separate txt files with the same name as in the code, only copy the values (no header, no text, make two columns):

   Resources:
month  f value
01  0.748
02  0.775
03  0.756
04  0.712
05  0.559
06  0.477
07  0.432
08  0.432
09  0.477
10  0.648
11  0.909
12  1.000

   Temperatures:
time   mean temperature
0   3.576178523
1   4.272769444
2   5.693903226
3   6.680989247
4   7.540075893
5   7.95569852
6   7.359872222
7   5.870103495
8   4.220568056
9   3.152399194
10  2.534522849
11  2.528454167
12  2.880237584
13  3.636294643
14  5.96272449
15  6.222271505
16  7.38149256
17  7.409368775
18  6.966427778
19  5.492860215
20  4.282763889
21  3.241544355
22  2.01715457
23  1.521733333
24  2.291408054
25  3.276548936
26  5.257983871
27  6.709844086
28  7.32331994
29  7.503302826
30  6.7152625
31  5.828805108
32  4.465727778
33  3.162452957
34  2.993697581
35  2.861191667
36  2.898080537
37  3.910002778
38  5.598921622
39  6.695869792
40  7.683181034
41  7.759943472
42  7.293447222
43  6.056247312
44  4.463266667
45  3.685244624
46  2.871793011
47  2.45
48  3.922534228
49  4.934752793
50  6.197569892
51  7.121584677
52  8.077349702
53  8.012516824
54  7.298430556
55  5.771642473
56  4.241619444
57  3.590569892
58  3.162024194
59  3.287647222
60  4.061998658
61  4.965858136
62  6.445442204
63  7.738768817
64  8.714534226
65  8.409689098
66  7.613561111
67  5.821275538
68  4.157559722
69  3.148745968
70  2.482346774
71  2.836618056

7. Submodels

Update calendar The model advances by one month on the gregorian calendar.

Update environmental variables The model takes the relevant temperature and f values for the current month in the input list.

Competition and f The model calculates the current population density and quantifies the competition effect on the food availability depending on how far from the carrying capacity the population density is, and updates the f value in accordance.

If popdensity < 1.9 * carcap, then foodcompet = ((1 - f-scaled) * (1 - (popdensity / (2 * carcap - popdensity)))) If popdensity ≥ 1.9 * carcap, then foodcompet = ((1 - f-scaled) * (1 - (popdensity / (car_cap / 10))))

The use of two types of formula is because if popdensity = 2 * carcap, the first formulat gives an error due to a division by 0, and if popdensity > 2 * carcap, then the formula gives the untrue result of less competition with a bigger population.

Explanation of the procedure: Competition is only effective if the food availability is less than the maximum (hence the use of '1 - f-scaled'). In which case, proportionnally to how much the f is lessened compared to maximum, the size of the current population has an influence on how important the competition is: If the population is below the carrying capacity (carcap), then food is more available for the individuals present, but if the population is over the carrying capacity, the availability of food is lessened. Therefore, if pop > carcap, foodcompet < 0 <=> f decreases if pop = carcap, foodcompet = 0 <=> f constant if pop < carcap, food_compet > 0 <=> f increases

Meaning, the competition is actually calculated depending on how far from the carrying capacity is the population density, and how far from maximum f is the food availability. The less food there is and the bigger the population, the higher the competition.

The regulatory effect of this competition lies in the starvation of individuals at lower food availability, which leads to a reduction of the population size with higher competition (combination of low food availability and big population size).

f is always contained between 0 and 1. If the competition is turned off, f is the direct value taken from the input list.

Convert parameters with TC A temperature correction factor (TC) is calculated using the Arrhenius temperature (TA) and applied to conductance ̇v, somatic maintenance rate ̇k_M, maturity maintenance rate ̇k_J and Weibull ageing acceleration ḧa, which are all affected in the same way by the correction factor. TC = exp((TA/Tref - TA/T) A given metabolic rate X at temperature T is thus modified with: X(T) = X * exp(TA/Tref - TA/T) with Tref being the reference temperature (20°C) and T the actual temperature (in Kelvin) of the organism’s life environment.

Change in reserve The reserve is supplied from ingested food, represented in the model by the functional response f (from 0 to 1) proportional to food availability. Scaled assimilation rate S_A is found with ̇p_A /{ṗ_Am}, where ̇p_A is the flux of assimilation (in energy per time) and ṗ_Am the maximal assimilation flux per surface area of structure (with the same dimension). Since ̇p_A = {ṗ_Am} * f * L², with L the structural length, then S_A = f * L². A flux of mobilized energy goes outside of the reserve compartment: the scaled mobilisation flux S_C is the scaled equivalent of ̇p_C / {ṗ_Am} therefore equal to: L² * (ge/(g+e)) * (1+L * ̇k_M/̇v) = S_C, where g is the energy investment ratio (the cost of the added volume for this timestep relative to the maximum potentially available energy for growth and maintenance), e is the scaled reserve density (reserve density relative to maximum reserve density) and L the structural length, and with ̇k_M the rate of mobilisation of the κ fraction of S_C for somatic maintenance, proportional to structural length, and ̇v the energy conductance. The reserve dynamics calculated at each time step correspond to ∂U_E = S_A - S_C.

Change in maturity or reproduction buffer Before puberty (ability to reproduce), changes in maturity level are calculated as the flux of energy going into the maturity compartment, that is the fraction 1-κ of the mobilisation flux after paying for the maintenance costs of the maturity compartment U_H: ṗ_R = (1-κ) * ṗ_C - ṗ_J.

The maturity level of the compartment U_H changes each month through the scaled formula for ∂U_H: S_R = (1-κ) * S_C - ̇k_J * U_H = ∂U_H,

while the reproduction buffer U_R receives no energy, thus ∂U_R is set to 0.

Juveniles keep growing until they reach puberty, when the maturity level U_H is equivalent to U_H^p. At this point, they are able to reproduce, thus the energy flux S_R is redirected entirely to the reproduction buffer U_R and the maturity compartment does not increase anymore: U_H stays set as U_H^p. Therefore, after puberty, and except for females undergoing reproduction: S_R = (1-kap) * S_C - ̇k_J * U_H^p = ∂U_R .

Change in structural length The structural length L changes with the flow of energy that remains from the fraction κ of mobilisation flux S_C after somatic maintenance has been paid for: the structural length at the end of the month is different from the month before by ∂L, calculated with the scaled formula ∂L = (1/3) * ((̇v/gL²) * S_C - ̇k_M * L).

If scaled reserve falls below the scaled structural length, starvation is launched.

Starvation When scaled reserve falls below the scaled structural length l (length relative to maximum length), that is when e < l, it is assumed that the individual is confronted to starving conditions: the kappa rule is then altered as energy is entirely redirected to somatic maintenance and all other fluxes (growth, reproduction or maturation) are set to 0, the following calculations are made:

Mobilisation flux S_C = ([ṗ_M] / L³)/({ṗ_Am} ) . Since [ṗ_M] = [E_G] * ̇k_M and [E_G] = g * kappa{ṗ_Am} / ̇v ,

we can rewrite S_C = (L³ * ̇k_M * g * kappa{ṗ_Am} / ̇v ) / ({ṗ_Am}) = (L³ * ̇k_M * g * kappa)/ ̇v

then recalculate ∂U_E = S_A – S_C .

When e < 0, the organism doesn’t have enough energy to pay somatic maintenance and dies.

The starvation strategy used in the population model was chosen among the ones presented in Kooijman 2010 based on Magniez 1983 research on A. cordatus reproduction and development.

Ageing Two ordinary differential equations are calculated: changes in ageing acceleration ̈q (also called the scaled density of damage inducing compounds) and changes in hazard mortality rate ḣ. ∂̈q = (̈q (L/L_m)³ * s_G + ḧ_a) * e * (̇v/L- ̇r)- ̇r * ̈q

∂ḣ = ̈q - ̇r * ḣ , with ̇r = (3/L)*∂L

These calculations are used for the simulation of the accumulation of damage inducing compounds and their effect, following the DEB theory for ageing (Kooijman 2010). Damage inducing compounds density is proportional to reserve mobilisation S_C and influences the hazard mortality rate ḣ which is a function of the damage accumulated in the body. Damage inducing compounds are diluted via growth ̇r, and additionally ageing is calculated with two other parameters, the Weibull ageing acceleration ḧ_a and the Gompertz stress coefficient s_G . In other words, the hazard mortality rate is the simulation of the vulnerability of the individual towards damage, such as the risk of dying from an illness increasing as the individual ages. Additionally, in our model, the ageing submodel relies on a stochastic element, where the individual has a 50% chance of looking into its death probability rate ḣ.

Update individuals The calculated changes are applied to each state variables of the individual: For a state variable X, X = X + ∂X * 30.5 . The temporal resolution is a monthly interval: each ODE is calculated then the resulting ∂ is applied *30.5. If the individual has reached a maturity level corresponding to a threshold, it updates its life stage (i.e. its breed) accordingly. The individual also updates its age.

Update reproduction and birth timers The reproduction calendar (reprotime) and the GSI calendar (GSItime) advance by one month each timestep, and fall back to the start in a twelve months cycle. The starting date of the two calendars is not the same (March for reprotime and June for GSItime). The birth timer (birth_time) is updated if it has been set off, by counting down (e.g. if it was at 7, the timer will be set to 6). If the timer has not been set off (i.e. is still set to 0), it will stay at 0. This timer is borne by females only and will set off if the female has launched reproduction.

Reproduction Only females are considered in the reproduction processes. The value of the Gonado-somatic index (GSI = 100 * ((ash-free gonads dry weight) / (ash-free body dry weight)) when the reproduction period comes is what allows a female to reproduce. Whenever the level of GSI is high enough (at least 0.07%), the female signals being ready for reproduction. If the GSI falls back below the threshold, the signal is turned off.

Following is the explanation of the process in the order of execution (see pseudo-code): If the signal is on when the reproduction period comes (i.e. it has been turned on the month before the reproduction period), females note down their value of GSI and prepare their eggs. Females calculate their GSI every month (see submodel 'Calculate GSI'). If it is the month before the reproduction period and the GSI is above the threshold, females note down the value of energy density in their reproduction buffer and in addition to turning their reproduction signal on, they start their birth timer (by setting the 'birthtime' countdown to 9) which serves as a tracker in the reproduction process from conception of offsprings till their birth. Once the reproduction period starts, the birth timer is what allows to verify if the individual is undergoing reproduction and to adapt its state variables accordingly: at birthtime 8, 7 and 6, females are reproducing (i.e. conceiving offsprings by decreasing the energy in their reproduction buffer, see below); at birth_time 3, 2 and 1, females release offsprings (i.e. a number of new juveniles proportional to the number of females having reproduced is initiated into the model, see below).

When females are reproducing, conception of offsprings causes a decrease in energy in their reproduction buffer: their usual ∂U_R (see submodel 'Change in maturity or reproduction buffer') is set to 0 for the three months, while U_R is forced to decrease: for each month of the reproduction period the female decreases its buffer by a third of 52% of the energy stored: ∂₂U_R = (U_Rstart - 0.52 * U_Rstart) / 3, with U_{R_start} the reproduction buffer at the start of the period. The GSI follows a similar pattern (see submodel 'Calculate GSI').

When females release offsprings (i.e. give birth), five months after conception, 65% of the eggs are assumed to have survived until birth. The reproductive output (Ri) is: Ri = 0.65 * eggs For each of the three months of offspring release, the number of juveniles (Ri / 3) are intiated into the model with parameters from the calculations of the initial simulation.

Calculate GSI The GSI is estimated for each month according to the time of accumulation of energy into the reproduction buffer from the end of the reproduction, with this equation: GSI = ((timeofaccumulation * ̇k_M * g) / (f³ * (f + κ * g * y_VE ))) * ((1-κ) * f³ - ((̇k_J * U_H^p) / (L_m² * s_M³)) ) , where the time of accumulation is the number of days since the end of the reproduction period, ̇k_M the somatic maintenance rate coefficient and ̇k_J the maturity maintenance rate coefficient, g the energy investment ratio, f the scaled functional response, y_VE is the parameter for the yield of structure on reserve, that is the number of moles of structure that can be produced with one mole of reserve, s_M the acceleration factor, U_H^p the scaled energy in the complexity compartment at puberty, κ the fraction of energy directed towards structure and 1-κ the fraction of energy directed towards complexity and L_m the maximum volumetric length (see "Entities, State variables and scales" for the dimensions and in-code notations).

The GSI of a reproducing female will decrease by 52% of its initial value over the 3 months period of reproduction (i.e. a decrease of one third of 52% per month): ∂GSI = (GSI_start - 0.52 * GSI_start) / 3, where GSI_start is the level of gonadal index at the onset of reproduction.

Background mortality The barckground mortality rate is applied to the overall population: 3.42% of juveniles and 2% of adults (males and females) die each month. Depending on the cause of death, the individuals set on a certain flag (deceasedbg or deceasedold) and a 'deceased' flag. They are then removed only at the start of the next tick (for monitoring purpose).

Check temperature Depending on the temperature for the current and prior month and on the type of sensitivity to temperatures chosen for the model, a mortality rate is applied to the population for temperatures from 8 to 12°C. For a "vulnerable" setting, temperatures exceeding thresholds of 8, 9.5, 11 and 12°C for two consecutive months cause a mortality rate of 25%, 35%, 45% and 100% respectively. For an "intermediate" setting, temperatures exceeding thresholds of 8, 9.5, 11 and 12°C for only one month cause a mortality rate of 10%, 20%, 30% and 100% respectively. For a "resistant" setting, temperatures exceeding thresholds of 8, 9.5, 11 and 12°C for two consecutive months cause a mortality rate of 10%, 20%, 30% and 100% respectively. When the cause of death is temperature, individuals die on the spot (contrary to background mortality).

Monitoring of population At the end of each timestep, the population density is calculated and the data collected for monitoring and plotting mean values of state variables. Plots are built on the lists compiled out of all individuals state variables values. To facilitate processing of the outputs afterwards, mock values are implemented in case none of the individuals possessing the variable is present at the timestep. These mock values are values purposely exceedingly higher than the normal range for the variable, so that they can then be easily spotted and taken out during processing of results.

Update time The counters for the number of years and months since the beginning of the simulation advance by respectively one twelfth and one, each timestep.

References (short list)

Delille, D., Bouvy, M., 1989. Bacterial responses to natural organic inputs in a marine subantarctic area. Hydrobiologia, 182(3), 225-238. doi:10.1007/BF00007517. Guillaumot, C., 2019. AmP Abatus cordatus, version 2019/01/17. https://www.bio.vu.nl/thb/deb/deblab/add_my_pet/entries_web/Abatus_cordatus/Abatus_cordatus_res.html Kooijman, S.A.L.M., 2010. Dynamic energy budget theory for metabolic organisation, third ed. Cambridge university press. Martin, B.T., Zimmer, E.I., Grimm, V., Jager, T., 2010. DEB-IBM User Manual: Dynamic Energy Budget theory meets individual‐based modelling: a generic and accessible implementation. Available at https://www.bio.vu.nl/thb/deb/deblab/debibm/DEB_IBM_manual.pdf.
Railsback, S. F., Grimm, V., 2019. Agent-based and Individual-based Modeling: A Practical Introduction. Princeton university press.

Supplementary information:

For further details and information such as sensitivity analysis, detailed explanations of processes in submodels, the ecophysiological data used and full references, see our research article on the construction of this model: Arnould-Pétré Margot, Guillaumot Charlène, Danis Bruno, Féral Jean-Pierre and Saucède Thomas "Modelling individual-based population dynamics for an endemic species of the Southern Ocean, the sea urchin Abatus cordatus, under contrasting environmental conditions”, submitted to the Ecological Modelling Special Issue of the 2019 Dynamic Energy Budget Symposium.

Main things to keep in mind in this implementation:

For the model to work, the files for the environmental data (temperatures: "temptimemonthavgHalage.txt"_ and resources: "inputRScesMfDelille.txt"_) need to exist in the same folder as where the model is stored.

For now, the model only fully works if the site Anse du Halage is selected. Sites Port Couvreux and Ile Haute can also be selected, however the site-specific data is only available for temperature and not for resources. Thus, if one of those two sites is selected, the file for the resources data will be that of Anse du Halage. This is merely to test the model for different temperature data using real time-series records rather than future projections. NB: Here again make sure that the relevant temperature files ("temptimemonthavgCouvreux.txt"_ or "temptimemonthavgHaute.txt"_) are present in the same folder as the model.
The temporal resolution is a monthly interval : each ODE is calculated then the resulting ∂ is applied *30.5
Temperature and functional response f are in the form of time-series. (To test with constant and DEB standard values, uncomment the corresponding lines at the beginning of the go procedure in the code.)
This species lives in water at ~5°C average, however we are missing the upper limit of the arrhenius relationship. Because of this, the mortality due to temperature is mainly randomly applied among individuals in the population without distinction.
Depending on the cause of death, the individuals set on a certain flag (deceasedbg or deceasedold) + a 'deceased' flag and are removed only at the start of the next tick (for monitoring purpose). If the cause is temperature, they die on the spot.
State variables plots are built on the lists compiled during the calc-pop-varb procedure (called by the update-pop procedure at the end of each run/month), which are built so that in case none of the individuals possessing the variable is present, an outlier mock value is applied instead. This is used to facilitate processing of the outputs afterwards. Death plot is non cumulative: death counts reset at the beginning of each tick.

Comments and Questions

Margot Minju Arnould-Pétré

How to run the model:

A few simple steps are necessary to run this model under its basic implementation (you can find the same information in the file ODD.pdf - "How to run the model + Model description"): 1/ The files for the environmental data need to exist in the same folder as where the model is stored. To do this, after downloading the model, download the data files from the ‘’Files’’ tab in the NetLogo modeling commons (http://modelingcommons.org/browse/one_model/6201#model_tabs_browse_files). For the basic implementation, the two files needed are “temp_time_monthavg_Halage.txt” for temperatures and “inputRSces_M_f_Delille.txt” for resources. Make sure the model (.nlogo) and the data (.txt) files are stored in the same folder in the computer. 2/ Once you open the model (.nlogo), the interface is generally the first thing that is visible. Navigation between the interface, the information page (this page), and the code of the model is done through the three tabs at the top of the software (‘Interface’, ‘Info’, ‘Code’). Make sure that the following elements are selected in the interface: • Sites: ‘Anse du Halage’ • projection: ‘present’ • future: ‘mixed temp & food’ • sensitivity: ‘resistant’ • competition: ’On’ • run_time: ‘210’ • cv: ‘0.1’ • add-my-pet?: ‘Off’ Also ensure that none of the green boxes (‘input paramaters’) is empty. If any of them is empty, switch the add-my-pet? button ‘On’ and fill the relevant boxes with the basic DEB parameters for A. cordatus as taken from the Add-my-Pet database (https://www.bio.vu.nl/thb/deb/deblab/add_my_pet/entries_web/Abatus_cordatus/Abatus_cordatus_res.html): [ṗM], EHb, EHp, [EG] and Lm, which correspond to these boxes respectively: p_M, E_H^b, E_H^p, E_G, zoom. Except in the aforementioned case, do not modify any of the parameters in the green boxes placed under the line « Input parameters » on the interface. 3/ Click on the purple setup button. This initializes the model, and should barely take a second on an average computer. A sure way of knowing the model has finished setup, is color shapes appearing in the small black square that is on the bottom-right of the purple buttons. 4/ Once setup is finished, click on the purple go button. This will run the model for the simulated duration input in the run_time box (number of years). Clicking on the go button again before the end of the simulation will pause the model, clicking on it after the end of the simulation will continue the simulation without a temporal limit. The go once button will only run the model for a single loop, that is a simulation of one month. 5/ The model can be run for future projections with different combinations of food and/or temperature scenarios. For this, select in the interface the desired RCP scenario under projection and the wanted combination under future. Three types of sensitivity to high temperature are also available under sensitivity (see « Check temperature » submodel in the ODD below for a short explanation of the difference between the three types).

Posted over 5 years ago