Code block 3.3.1.1

PSPMdimensions <- c(PopulationNr = 1, IStateDimension = 1, LifeHistoryStages = 2)

The software can simultaneously compute the population growth of more than a single population. The vector element PopulationNr of PSPMdimensions has to be defined equal to the number of structured populations accounted for in the model. For the Medfly example this is obviously equal to 1. The vector element IStateDimension of PSPMdimensions defines the dimension of the individual state. As only age characterizes the individuals in the Medfly model, this variable is defined equal to 1. Finally, the element LifeHistoryStages of PSPMdimensions has to be defined equal to the number of life stages that can be distinguished in the individual life history. While integrating the ODEs for the individual life history numerical problems may occur when the right hand side of the ODEs changes abruptly in value at a certain threshold value of the individual state, as a consequence of discontinuities in the development rate, the mortality rate or the fecundity. Each of such thresholds in the individual life history should be distinguished as a stage boundary. In the Medfly model the fecundity \(\beta(a)\) changes from 0 just before \(a=A_J\) to \(\beta_0\) at \(a=A_J\) and \(\beta_0\exp(-\beta_1(a-A_j))\) at larger ages. At \(a=A_j\) \(\beta(a)\) thus exhibits a discontinuity, which separates the juvenile and the adult stage from each other. The element LifeHistoryStages of PSPMdimensions is therefore set equal to 2.

Code block 3.3.1.2

NumericalOptions <- c(MIN_SURVIVAL  = 1.0E-9,   # Survival at which individual is considered dead
                      MAX_AGE       = 100000,   # Absolute maximum individual age
                      DYTOL         = 1.0E-7,   # Variable tolerance
                      RHSTOL        = 1.0E-8)   # Function tolerance

The specification of NumericalOptions is optional and can be left out, if default values are acceptable. A list of all possible vector elements that can be included in the NumericalOptions variable is provided in section 9.3.

The vector element MIN_SURVIVAL of NumericalOptions determines the threshold of the survival probability below which an individual is considered dead. The integration over the individual life history stops whenever the survival probability falls below this threshold value. In the code above the minimum survival is set to \(10^{-9}\), which is in fact the default value and is hence superfluous. Note that the value of MIN_SURVIVAL can not be set equal to 0. As an alternative to using 0 MIN_SURVIVAL can be set to a very small value like \(10^{-100}\).

The vector element MAX_AGE of NumericalOptions can be used as an alternative to determine the end of an individual life and to stop the integration over the individual life history. In the Medfly model there is no maximum individual age and hence the variable is set to a very high value (100000), which the individuals will never reach, because before that age their survival probability has already dropped below its threshold value (\(10^{-9}\)).

The remaining two vector elements DYTOL and RHSTOL of NumericalOptions determine whether a solution has been found. In general, both demographic analysis as well as equilibrium analysis of PSPMs boils down to solving a system of nonlinear equations that can be represented as \(G(y)=0\) for a set of unknowns \(y\) in an iterative manner. The subsequent estimates of the solution in the Newton iterations can be labeled as \(y_p\) and \(y_{p+1}\). A solution is now considered to be located if both of the following conditions hold:

\[\begin{align*} &\|{y_{p+1}-y_p}\| < \epsilon_y\\[2ex] &\|{G(y_{p+1})}\| < \epsilon_G \end{align*}\]

where \(\|.\|\) refers to the Euclidean norm. DYTOL and RHSTOL are the quantities \(\epsilon_y\) and \(\epsilon_G\), respectively. Increasing (decreasing) their value leads to easier (harder) acceptance of a set of unknowns as a solution to the system of equations \(G(y)=0\). The definition of these two accuracies in the code box above is in fact superfluous as they are defined equal to their default values (see section 9.3).

Code block 3.3.1.3

DefaultParameters <- c(Beta0 = 47.0, Beta1 = 0.04, Aj = 11.0, Mu0 = 0.00095, Mu1 = 0.0581)

The names of the vector elements (parameters) can be used in the programming of the life history functions of the model and are furthermore used to make the output files produced by the program more readable. These output files contain a small header text indicating among other details which parameter values were used for the computation of the results contained in the output file. In this report the parameter names are listed together with their value.

Code block 3.3.1.4

StateAtBirth <- function(E, pars)
{
  with(as.list(c(pars)),{
        # We model a single structured population with a single i-state variable (age)
        c(Age = 0.0)
      })
}

Notice that the arguments of the function StateAtBirth() contain a vector E in addition to the vector with parameter values pars. The vector E will contain the values of the environment variables during equilibrium computations of PSPMs (see sections 4.1 to 4.4). In demographic analysis of PSPMs this variable is non-functional and is best ignored, using it in a statement inside the routine may even cause the program to crash. The only reason for the presence of this variable among the function arguments is to keep the function declaration the same for both demographic and equilibrium analysis computations. In principle, the same model-specific file can hence be used for both types of analysis. The variable E will for the same reasons also be part of the headers of the next 2 routines.

Code block 3.3.1.5

LifeStageEndings <- function(lifestage, istate, birthstate, BirthStateNr, E, pars) {
  with(as.list(c(E, pars, istate)),{
        maturation  = switch(lifestage, Age - Aj, -1)
      })
}

The threshold value returned to the program in maturation will in general depend on the individual state variables, possibly on the individual’s state-at-birth and will be different for individuals in different life stages. For this reason, the function LifeStageEndings() has as arguments lifestage, specifying the life stage that the individual is currently in, istate, the individual state, and birthstate, the individual’s state-at-birth in addition to the arguments E and pars that have the same interpretation as discussed above for the function StateAtBirth().

This routine will be called as many times as there are possible states-at-birth, because the state-at-birth may influence the threshold between consecutive life stages. The same holds for the next 2 routines discussed, which define changes in the i-state variables, the fecundity and the mortality of individuals. In essence, individuals with different states-at-birth are treated as constituting subpopulations within the same structured population. Because of the possible dependence on the state-at-birth the variables birthstate and BirthStateNr are passed as arguments to the function LifeStageEndings() as well as to the 2 functions discussed below. These arguments contain the values of the i-state variables and the index in the set \(\left\{\boldsymbol{\phi}_1,\ldots,\boldsymbol{\phi}_m\right\}\), respectively, for which the routine is invoked and for which the threshold between consecutive life stages has to be evaluated.

Code block 3.3.1.6

LifeHistoryRates <- function(lifestage, istate, birthstate, BirthStateNr, E, pars) {
  with(as.list(c(pars, istate)),{
        list(
            # We model a single structured population (nrow=1) with a single i-state variable (age)
            development = 1.0,
            fecundity   = switch(lifestage, 0, Beta0*exp(-Beta1*(Age - Aj))),
            mortality   = Mu0*exp(Mu1*Age)
        )
      })
}

In the Medfly model the specification of the development rate is obviously trivial, as age is the only i-state variable. Furthermore, the code fragment above implements the function \(\beta(a)=\beta_0e^{-\beta_1(a-A_j)}\) for fecundity and the function \(\mu(a)=\mu_0e^{\mu_1a}\) for the mortality, which is not influenced by the life stage specifically and is only age-dependent.

Refer to the remarks in the discussion of the function LifeStageEndings() concerning the dependence on the individual’s state-at-birth.

Code block 3.3.1.7

DiscreteChanges <- function(lifestage, istate, birthstate, BirthStateNr, E, pars) {
  with(as.list(c(E, pars)),{
       # No discrete changes in this problem, function is commented out, which 
       # would be equivalent to returning a copy of the input argument 'istate'
       istate
     })
}

If defined, the function DiscreteChanges() is called whenever a transition between two consecutive life stages is reached during the integration over the individual life history. The function should return a vector of length equal to the number of i-state variables. Each element should specify the value of the particular i-state variable after the transition to the current state. In case the model accounts for multiple, structured populations this function should return a matrix with the number of rows equal to the number of structured populations in the problem and the number of columns equal to the number of i-state variables.

It should be noted that the value of the variable lifestage indicates the life stage that is entered, that is, following the current stage boundary. This routine will hence never be called with a value of one of the elements lifestage equal to 1. The discrete changes in the individual state variables have to be implemented by assigning new values to the variables istate. These assignments may as before depend on the life stage that is entered, as specified by the variable lifestage, the (old values) of the individual state variables, contained in the argument istate, and possibly on the individual’s state-at-birth, specified in the argument birthstate. If no assignment of a value to istate is implemented, that particular individual state variable will keep its current value.

Code block 3.3.2.1

// Dimension settings: Required
#define POPULATION_NR       1               // Structured consumer population
#define STAGES              2               // Juvenile & adult
#define I_STATE_DIM         1               // See below
#define PARAMETER_NR        5

// Numerical settings: Optional (default values adopted otherwise)
#define MIN_SURVIVAL        1.0E-9          // Survival at which individual is considered dead
#define MAX_AGE             100000          // Give some absolute maximum for individual age

#define DYTOL               1.0E-7          // Variable tolerance
#define RHSTOL              1.0E-8          // Function tolerance

The software can simultaneously compute the population growth of more than a single population. At the start of the problem file the variable POPULATION_NR has to be defined equal to the number of structured populations accounted for in the model. For the Medfly example this is obviously equal to 1 (line 2 in the code box above).

The variable STAGES has to be defined equal to the number of life stages that can be distinguished in the individual life history (line 3 in the code box above). While integrating the ODEs for the individual life history numerical problems may occur when the right hand side of the ODEs changes abruptly in value at a certain threshold value of the individual state, as a consequence of discontinuities in the development rate, the mortality rate or the fecundity. Each of such thresholds in the individual life history should be distinguished as a stage boundary. In the Medfly model the fecundity \(\beta(a)\) changes from 0 just before \(a=A_J\) to \(\beta_0\) at \(a=A_J\) and \(\beta_0\exp(-\beta_1(a-A_j))\) at larger ages. At \(a=A_j\) \(\beta(a)\) thus exhibits a discontinuity, which separates the juvenile and the adult stage from each other. The variable STAGES is therefore set equal to 2.

The variable I_STATE_DIM (line 4 in the code box above) defines the dimension of the individual state. As only age characterizes the individuals in the Medfly model, this variable is defined equal to 1.

The last required parameter that has to be specified is the number of parameters in the model (set in line 5 in the code box above). In the Medfly model this equals 5 (\(\beta_0\), \(\beta_1\), \(A_j\), \(\mu_0\) and \(\mu_1\)).

The remaining definitions in the code box are all optional and can be left away. A list of all possible variables that can be changed by a definition in this code section is provided in section 9.3. The variable MIN_SURVIVAL determines the threshold of the survival probability below which an individual is considered dead. The integration over the individual life history stops whenever the survival probability falls below this threshold value. In the code above (line 8) the minimum survival is set to \(10^{-9}\), which is in fact the default value and is hence superfluous. Note that the value of MIN_SURVIVAL can not be set equal to 0. As an alternative to using 0 MIN_SURVIVAL can be set to a very small value like \(10^{-100}\).

The variable MAX_AGE (line 9 in the code box above) can be used as an alternative to determine the end of an individual life and to stop the integration over the individual life history. In the Medfly model there is no maximum individual age and hence the variable is set to a very high value (100000), which the individuals will never reach, because before that age their survival probability has already dropped below its threshold value (\(10^{-9}\)).

The remaining two quantities DYTOL and RHSTOL determine whether a solution has been found. In general, both demographic analysis as well as equilibrium analysis of PSPMs boils down to solving a system of nonlinear equations that can be represented as \(G(y)=0\) for a set of unknowns \(y\) in an iterative manner. The subsequent estimates of the solution in the Newton iterations can be labeled as \(y_p\) and \(y_{p+1}\). A solution is now considered to be located if both of the following conditions hold:

\[\begin{align*} &\|{y_{p+1}-y_p}\| < \epsilon_y\\[2ex] &\|{G(y_{p+1})}\| < \epsilon_G \end{align*}\]

where \(\|.\|\) refers to the Euclidean norm. DYTOL and RHSTOL are the quantities \(\epsilon_y\) and \(\epsilon_G\), respectively. Increasing (decreasing) their value leads to easier (harder) acceptance of a set of unknowns as a solution to the system of equations \(G(y)=0\). The definition of these two accuracies in the code box above is in fact superfluous as they are defined equal to their default values (see section 9.3).

Code block 3.3.2.2

// Descriptive names of parameters in parameter array (at least two parameters are required)
char  *parameternames[PARAMETER_NR] =
    { "Beta0", "Beta1", "AJ", "Mu0", "Mu1"};

// Default values of all parameters
double  parameter[PARAMETER_NR] =
    {47.0, 0.04, 11.0, 0.00095, 0.0581};

Model parameter values are stored by the program in the vector variable parameter. The lines 2-3 above assign each of the elements this vector a more meaningful, model-specific name. These name strings can not be used in the remaining parts of the model implementation, they only serve to make the output files produced by the program more readable. These output files contain a small header text indicating among other details which parameter values were used for the computation of the results contained in the output file. In this report the parameter names are listed together with their value. To adapt the above code to a different model, the code on line 2 of the code box above should remain the same, only change line 3 as needed (possibly extending it over multiple lines in case there are many parameters).

The default values to use for the model parameters are specified by the declaration of the vector parameter[PARAMETR_NR] on line 6-7 of the previous code box. The values should be specified as a comma-separated array of values within braces (don’t forget the closing semi-colon at the end of the statement!). To adapt the above code to a different model, the code on line 6 of the code box above should remain the same, only change line 7 as needed (possibly extending it over multiple lines in case there are many parameters).

Code block 3.3.2.3

// Aliases definitions for all istate variables
#define AGE                 istate[0][0]

// Aliases definitions for all parameters
#define BETA0               parameter[ 0]   // Default: 47.0
#define BETA1               parameter[ 1]   // Default: 0.04
#define AJ                  parameter[ 2]   // Default: 11.0
#define MU0                 parameter[ 3]   // Default: 0.00095
#define MU1                 parameter[ 4]   // Default: 0.0581

The developmental rates in individual state, fecundity and mortality in any model depend on the individual state itself, on the individual’s state at birth and on model parameters. The value of the individual state variables at a particular age are always referred to with the program variable istate[\(p\)][\(i\)], where the index \(p\) refers to the number of the population and the index \(i\) refers to the number of the individual state variables. Notice that in C array indices run from 0 (as opposed to 1 like in R)! Similarly, the value of the individual’s state variables at birth are always referred to with the program variable birthstate[\(p\)][\(i\)]. In case there are multiple populations and/or more than a single individual state variable, it is up to the user to keep track of which index pertains to which population or individual state variable. In the Medfly model there is only a single population and a single individual state variable, while the state at birth is rather irrelevant as it equals age 0. Therefore, istate[0][0] is the only program quantity to give a more meaningful name (line 2 in the code box above).

As discussed in the previous section all model parameters are contained in a vector named parameter in the code. Which element of this vector represents which model-specific parameter is up to the user. To prevent mixing up the interpretation of the different vector elements and hence to prevent mistakes, it is strongly advised to define meaningful, model-specific aliases for each of the elements of the vector parameter as is illustrated in lines 5-9 in the code box above. It is best to avoid completely the direct use of the program variable parameter in any part of the model specification and only use the models-specific aliases.

As can be seen in the code block above all aliases for program variables used in the C-implementation of the model are names in capitals. It is advisable to use only capitals when introducing these aliases (or global variables if they are needed) to avoid any conflict between these aliases and variables that defined elsewhere in any of the C files with numerical routines that are included in the package.

Code block 3.3.2.4

/*
 * Specify the number of states at birth for the individuals in all structured
 * populations in the problem in the vector BirthStates[].
 */

void SetBirthStates(int BirthStates[POPULATION_NR], double E[])
{
  BirthStates[0] = 1;

  return;
}

For each population with index \(p\) the variable BirthStates[\(p\)] has to be set to the number of possible states at birth. Because individual age is the only i-state variable the Medfly model, the state-at-birth is unique and hence BirthStates[\(0\)] is set to 1.

Note that different populations may have different numbers of states-at-birth. BirthStates[\(p\)] hence does not need to be the same for all \(p\).

Code block 3.3.2.5

/*
 * Specify all the possible states at birth for all individuals in all
 * structured populations in the problem. BirthStateNr represents the index of
 * the state of birth to be specified. Each state at birth should be a single,
 * constant value for each i-state variable.
 *
 * Notice that the first index of the variable 'istate[][]' refers to the
 * number of the structured population, the second index refers to the
 * number of the individual state variable. The interpretation of the latter
 * is up to the user.
 */

void StateAtBirth(double *istate[POPULATION_NR], int BirthStateNr, double E[])
{
  AGE    = 0.0;

  return;
}

For every population (\(p=0, 1,\ldots,\) POPULATION_NR-1) the value of each individual state variable istate[\(p\)][\(i\)] (\(i=0, 1,\ldots,\) I_STATE_DIM-1) has to be assigned a unique value, from which individual development will start at age 0. As shown in the example of the Medfly model, if the life history depends on the age of the individual, age should be explicitly included as one of the individual state variables. The program does not automatically include individual age in its characterization of the individual state, even though integration over the entire life history (as a function of age) is carried out. For the Medfly model age is the only individual state variables and set to 0 at birth.

The routine StateAtBirth() will be called as many times as there are possible states-at-birth. The variable BirthStateNr indicates the index \(j\) of the state-at-birth in the set \(\left\{\boldsymbol{\phi}_1,\ldots,\boldsymbol{\phi}_m\right\}\) for which the values have to be set in the current invocation of the routine. The routine will thus be called with BirthStateNr set equal to a value in \(0, 1, \dots, m-1\) (Remember the starting index 0 in C!). If there are multiple states-at-birth (BirthStates[\(p\)] $ > 1$) the definition of the values of the i-state variables has to depend explicitly on the index BirthStateNr to make the states-at-birth different from each other. Furthermore, if the problem involves multiple structured populations the number of possible states-at-birth can be different for each of them, which might lead to a situation that the routine above is called with a value of the index BirthStateNr that is larger than the maximum number of states-at-birth for a particular population (BirthStateNr \(\geq\) BirthStates[\(p\)]). The program safely ignores such inappropriate state-at-birth specifications.

Code block 3.3.2.6

/*
 * Specify the threshold determining the end point of each discrete life
 * stage in individual life history as function of the i-state variables and
 * the individual's state at birth for all populations in every life stage.
 *
 * Notice that the first index of the variable 'istate[][]' refers to the
 * number of the structured population, the second index refers to the
 * number of the individual state variable. The interpretation of the latter
 * is up to the user.
 */

void IntervalLimit(int lifestage[POPULATION_NR], double *istate[POPULATION_NR],
                   double *birthstate[POPULATION_NR], int BirthStateNr, double E[],
                   double limit[POPULATION_NR])
{
  if (lifestage[0] == 0)
    limit[0] = AGE - AJ;

  return;
}

In this routine the variable limit[\(p\)] has to be defined, which has as many elements as there are populations (\(p=0\ldots\) POPULATION_NR-1). The life stage that the individual is in at the moment this routine is called, is determined by the variable lifestage[\(p\)], which has a value of 0 if the individual is in the first life stage and a value of STAGES-1 if it is in the last life stage. The element limit[\(p\)] should now indicate when the current life stage as given in lifestage[\(p\)] ends. In particular, the program considers the current life stage to end when limit[\(p\)] turns from negative to positive. For the end of the last life stage the death of old age, either by reaching the maximum age MAX_AGE or by reaching the minimum survival threshold MIN_SURVIVAL, does not have to be specified separately, the program takes care of that automatically. In the Medfly model therefore only the end of the larval stage has to be specified, as expressed in lines 16-17 of the code box above.

The threshold value that has to be stored and returned to the program in limit[\(p\)] will depend on the individual state variables, possibly on the individual’s state-at-birth and will be different for individuals in different life stages. For this reason, the routine IntervalLimit() has as arguments lifestage[], specifying the life stage that the individual is currently in, istate[][], the individual state, and birthstate[][], the individual’s state-at-birth.

Like the previous routine, this routine will be called as many times as there are possible states-at-birth, because the state-at-birth may influence the threshold between consecutive life stages. The same holds for the routines discussed in sections 3.3.2.7-3.3.2.10 below, which define changes in the i-state variables, the fecundity and the mortality of individuals, respectively. In essence, individuals with different states-at-birth are treated as constituting subpopulations within the same structured population. Because of the possible dependence on the state-at-birth the variables birthstate[][] and BirthStateNr are passed as arguments to this routine and the once discussed in sections 3.3.2.7-3.3.2.10. These arguments contain the values of the i-state variables and the index in the set \(\left\{\boldsymbol{\phi}_1,\ldots,\boldsymbol{\phi}_m\right\}\), respectively, for which the routine is invoked and for which the threshold between consecutive life stages has to be evaluated.

If the problem involves multiple structured populations and the number of possible states-at-birth differs among them, the routine above may be called with a value of the index BirthStateNr that is larger than the maximum number of states-at-birth for a particular population (BirthStateNr \(\geq\) BirthStates[\(p\)]). Although this circumstance may seem confusing, the user does not have to worry about it, as the program is designed to safely ignore such assignments of thresholds between consecutive life stages, changes in the i-state variables, fecundity and mortality of individuals for states-at-birth with index BirthStateNr \(\geq\) BirthStates[\(p\)]that are inappropriate for the structured population with index \(p\).

Notice that the function header shown in the code box above also contains an array E[] as a variable. This array will contain the values of the environment variables during equilibrium computations of PSPMs (see sections 4.1 to 4.4). In demographic analysis of PSPMs this variable is non-functional and is best ignored, using it in a statement inside the routine may even cause the program to crash. The only reason for the presence of this variable in the function header is to keep the function declaration the same for both demographic and equilibrium analysis computations. In principle, the same model-specific file can hence be used for both types of analysis. The variable E[] will for the same reasons also be part of the headers of the next 4 routines.

Tip: The more advanced user who wants to perform both demographic and equilibrium analysis using the same model-specific file should notice that the array of environment variables E[] can in principle be used inside all the routines, if the dimension ENVIRON_DIM determining the number of environment variables has been set (see code box 4.3.2.1 for details). The appropriate value to use for the environment variables should be assigned to the elements E[\(e\)] in the routine StateAtBirth() (see code box 3.3.2.5), after which it will keep the same value throughout all the subsequent routines.

Code block 3.3.2.7

/*
 * Specify the development of individuals as a function of the i-state
 * variables and the individual's state at birth for all populations in every
 * life stage.
 *
 * Notice that the first index of the variables 'istate[][]' and 'development[][]'
 * refers to the number of the structured population, the second index refers
 * to the number of the individual state variable. The interpretation of the
 * latter is up to the user.
 */

void Development(int lifestage[POPULATION_NR], double *istate[POPULATION_NR],
                 double *birthstate[POPULATION_NR], int BirthStateNr, double E[],
                 double development[POPULATION_NR][I_STATE_DIM])
{
  development[0][0] = 1.0;

  return;
}

This routine specifies the right-hand side of the ODE: \[\dfrac{d\boldsymbol{\chi}}{da}\;=\;g(\boldsymbol{\chi}, \boldsymbol{\chi}_b)\] that determines the continuous development of the individual state variables during the life history. In the Medfly model the specification is obviously trivial. More generally, the value of development[\(p\)][\(i\)] determines for each structured population \(p\) the development in the individual state variable \(i\). Notice that these development rates may differ in different life stages, for example growth in body size may be different for juveniles and adults in case adults invest a lot of energy into reproduction. The development rates should then be specified dependent on the current life stage the individual is in. This current life stage at the moment the routine is evaluated is contained in the variable lifestage[\(p\)]. The development rate may furthermore depend on the individual state variables and possibly on the individual’s state-at-birth, which is the reason for istate[][], the individual state, and birthstate[][], the individual’s state-at-birth, as arguments to this routine.

Refer to the remarks in section 3.3.2.6 concerning the dependence on the individual’s state-at-birth.

Code block 3.3.2.8

/*
 * Specify the possible discrete changes (jumps) in the individual state
 * variables when ENTERING the stage specified by 'lifestage[]'.
 *
 * Notice that the first index of the variable 'istate[][]' refers to the
 * number of the structured population, the second index refers to the
 * number of the individual state variable. The interpretation of the latter
 * is up to the user.
 */

void DiscreteChanges(int lifestage[POPULATION_NR], double *istate[POPULATION_NR],
                     double *birthstate[POPULATION_NR], int BirthStateNr, double E[])
{
  return;
}

This routine is not relevant in case of the Medfly model and hence its contents are empty (apart for the necessary return; statement).

This routine is called whenever a transition between two consecutive life stages is reached during the integration over the individual life history. It should be noted that the value of the variable lifestage[\(p\)] indicates the life stage that is entered, that is, following the current stage boundary. This routine will hence never be called with a value of one of the elements lifestage[\(p\)] equal to 0. The discrete changes in the individual state variables have to be implemented by assigning new values to the variables istate[\(p\)][\(i\)]. These assignments may as before depend on the life stage that is entered, as specified by the variable lifestage[], the (old values) of the individual state variables, contained in the argument istate[][], and possibly on the individual’s state-at-birth, specified in the argument birthstate[][]. If no assignment of a value to istate[\(p\)][\(i\)] is implemented, that particular individual state variable will keep its current value.

Refer also to the remarks in section 3.3.2.6 concerning the dependence on the individual’s state-at-birth.

Code block 3.3.2.9

/*
 * Specify the fecundity of individuals as a function of the i-state
 * variables and the individual's state at birth for all populations in every
 * life stage.
 *
 * The number of offspring produced has to be specified for every possible
 * state at birth in the variable 'fecundity[][]'. The first index of this
 * variable refers to the number of the structured population, the second
 * index refers to the number of the birth state.
 *
 * Notice that the first index of the variable 'istate[][]' refers to the
 * number of the structured population, the second index refers to the
 * number of the individual state variable. The interpretation of the latter
 * is up to the user.
 */

void Fecundity(int lifestage[POPULATION_NR], double *istate[POPULATION_NR],
               double *birthstate[POPULATION_NR], int BirthStateNr, double E[],
               double *fecundity[POPULATION_NR])
{
 if (lifestage[0] == 1)             // Only for adults
    {
     fecundity[0][0]  = BETA0*exp(-BETA1*(AGE - AJ));
    }
  else
    fecundity[0][0] = 0;

  return;
}

In this routine not only the fecundity (i.e. the number of offspring produced per unit time) has to be specified, but also the state-at-birth of the produced offspring. Therefore, this routine has to assign values to the matrix fecundity[\(p\)][\(j\)], which determines for the population with index \(p\) the number of offspring produced per unit time with state-at-birth with index \(j\) in the set \(\left\{\boldsymbol{\phi}_1,\ldots,\boldsymbol{\phi}_m\right\}\). This fecundity will certainly depend on the life stage that the individual is in (only adults reproduce), which is contained in the argument lifestage[], and on the individual state variables, i.e. the values of the argument istate[][]), but possibly also on the individual’s state-at-birth, the values and index of which are specified by the arguments birthstate[][] and BirthStateNr, respectively.

In the most common case of a unique state-at-birth and a single structured population, like in the Medfly model, the only valid indices are \(p=0\) and \(j=0\) and hence only the variable fecundity[0][0] has to be assigned.

For more detailed remarks about models with multiple states-at-birth consult section 3.3.2.6.

Code block 3.3.2.10

/*
 * Specify the mortality of individuals as a function of the i-state
 * variables and the individual's state at birth for all populations in every
 * life stage.
 *
 * Notice that the first index of the variable 'istate[][]' refers to the
 * number of the structured population, the second index refers to the
 * number of the individual state variable. The interpretation of the latter
 * is up to the user.
 */

void Mortality(int lifestage[POPULATION_NR], double *istate[POPULATION_NR],
               double *birthstate[POPULATION_NR], int BirthStateNr, double E[],
               double mortality[POPULATION_NR])
{
  mortality[0] = MU0*exp(MU1*AGE);

  return;
}

For each population the corresponding element of the array mortality[\(p\)] should be assigned the mortality rate, possibly dependent on the life stage the individual is in at the moment this routine is called (given in argument lifestage[]), the current i-state of the individual (given in argument istate[][]) and the individual’s state-at-birth (current values and index in the set \(\left\{\boldsymbol{\phi}_1,\ldots,\boldsymbol{\phi}_m\right\}\) given by birthstate[][] and BirthStateNr, respectively).

In the Medfly model the mortality is not influenced by the life stage specifically and is only age-dependent. The line 16 in the code box above implements the function \(\mu(a)=\mu_0e^{\mu_1a}\).

Refer to the remarks in section 3.3.2.6 concerning the dependence on the individual’s state-at-birth.

Code block 4.3.1.1

PSPMdimensions <- c(PopulationNr = 1, IStateDimension = 2, LifeHistoryStages = 3, ImpactDimension = 4)

The software allows for the analysis of models with multiple structured populations, each of which consists of individuals that are characterized by a finite number of individual state variables. The number of state variables characterizing an individual should, however, be the same for each of the structured populations in the model. The vector element PopulationNr of PSPMdimensions has to be defined equal to the number of structured populations accounted for in the model. For the PNAS2002 example this is obviously equal to 1. The vector element IStateDimension of PSPMdimensions defines the dimension of the individual state. In the PNAS2002 model this variable is defined equal to 2 as the individual age is included among the individual state variable in addition to the individual length.

The element LifeHistoryStages of PSPMdimensions has to be defined equal to the number of life stages that can be distinguished in the individual life history. While integrating the ODEs for the individual life history numerical problems may occur when the right hand side of the ODEs changes abruptly in value at a certain threshold value of the individual state, as a consequence of discontinuities in the development rate, the mortality rate or the fecundity. Each of such thresholds in the life history should be distinguished as a stage boundary. In the PNAS model the mortality changes discontinuously at \(\ell=\ell_v\), while the fecundity changes discontinuously at \(\ell=\ell_j\). Three life stages can hence be distinguished: vulnerable juveniles, invulnerable juveniles and adults, and the changes in the life history rates are indeed abrupt at the transition boundaries between these stages. The element LifeHistoryStages of PSPMdimensions is therefore set equal to 3.

Finally, the ImpactDimension of PSPMdimensions has to be defined equal to the number of functions that represent the impact of an individual on its environment. In the PNAS model this feedback of an individual consumer on its environment consists of its grazing rate \(I(\ell,R)\) and the biomass-length relation \(\omega\ell^3\) determining the biomass of vulnerable consumers, as it represents food for predators. Therefore, the value of ImpactDimension should be at least set equal to 2. However, all interaction functions are also saved to the output file generated during a computation. The interaction functions can hence be conveniently used to produce arbitrary output quantities of the form

\[\int_\Omega h(\boldsymbol{\chi}, \boldsymbol{\chi}_b, E)\,\tilde{n}(\boldsymbol{\chi})\,d\boldsymbol{\chi}\]

where \(h(\boldsymbol{\chi}, \boldsymbol{\chi}_b, E)\) is an interaction (weighing) function that can depend on the values of the individual state \(\boldsymbol{\chi}\), the state-at-birth of individuals \(\boldsymbol{\chi}_b\) and the values of the environment variables \(E\), and \(\tilde{n}(\boldsymbol{\chi})\) is the stable population distribution in equilibrium. Such quantities can therefore represent the total population density in equilibrium (when \(h(\boldsymbol{\chi}, \boldsymbol{\chi}_b, E)=1\)), the total population biomass (when \(h(\boldsymbol{\chi}, \boldsymbol{\chi}_b, E)\) equals the biomass of an individual with individual state \(\boldsymbol{\chi}\) and state-at-birth \(\boldsymbol{\chi}_b\)) or the total population birth rate in equilibrium (when \(h(\boldsymbol{\chi}, \boldsymbol{\chi}_b, E)\) is the fecundity of an individual with individual state \(\boldsymbol{\chi}\) and state-at-birth \(\boldsymbol{\chi}_b\)). In the PNAS model I want in addition to the biomass of vulnerable consumers, also the biomass of non-vulnerable juvenile consumers and the biomass of adult consumers as output of the model and hence have set the value of the element ImpactDimension in PSPMdimensions equal to 4.

Code block 4.3.1.2

NumericalOptions <- c(MIN_SURVIVAL  = 1.0E-9,   # Survival at which individual is considered dead
                      MAX_AGE       = 100000,   # Give some absolute maximum for individual age
                      DYTOL         = 1.0E-7,   # Variable tolerance
                      RHSTOL        = 1.0E-6,   # Function tolerance
                      ALLOWNEGATIVE = 0,        # Negative solution values allowed?
                      COHORT_NR     = 100)      # Number of cohorts in population state output

The specification of NumericalOptions is optional and can be left out, if default values are acceptable. A list of all possible vector elements that can be included in the NumericalOptions variable is provided in section 9.3.

The vector element MIN_SURVIVAL of NumericalOptions determines the threshold of the survival probability below which an individual is considered dead. The integration over the individual life history stops whenever the survival probability falls below this threshold value. In the code above the minimum survival is set to \(10^{-9}\), which is in fact the default value and is hence superfluous. Note that the value of MIN_SURVIVAL can not be set equal to 0. As an alternative to using 0 MIN_SURVIVAL can be set to a very small value like \(10^{-100}\).

The vector element MAX_AGE of NumericalOptions can be used as an alternative to determine the end of an individual life and to stop the integration over the individual life history. In the PNAS2002 model there is no maximum individual age and hence the variable is set to a very high value (100000), which the individuals will never reach, because before that age their survival probability has already dropped below its threshold value (\(10^{-9}\)).

The two vector elements DYTOL and RHSTOL of NumericalOptions determine whether a solution has been found. In general, both demographic analysis as well as equilibrium analysis of PSPMs boils down to solving a system of nonlinear equations that can be represented as \(G(y)=0\) for a set of unknowns \(y\) in iterative manner. The subsequent estimates of the solution in the Newton iterations can be labeled as \(y_p\) and \(y_{p+1}\). A solution is now considered to be located if both of the following conditions hold:

\[\begin{align*} &\|{y_{p+1}-y_p}\| < \epsilon_y\\[2ex] &\|{G(y_{p+1})}\| < \epsilon_G \end{align*}\]

where \(\|.\|\) refers to the Euclidean norm. DYTOL and RHSTOL are the quantities \(\epsilon_y\) and \(\epsilon_G\), respectively. Increasing (decreasing) their value leads to easier (harder) acceptance of a set of unknowns as a solution to the system of equations \(G(y)=0\). The definition of these two accuracies in the code box is in fact superfluous as they are defined equal to their default values (see section 9.3).

The vector element ALLOWNEGATIVE of NumericalOptions is a flag that can only have a value of 0 or 1 and determines whether or not computations should stop when one of the variables to solve for reaches a negative value. In most population models negative solution values are biologically not relevant and ALLOWNEGATIVE is hence set to 0 by default. Line 16 in the code box above is only included to illustrate the use of ALLOWNEGATIVE and does not change the value of this variable from its default value. Most likely, setting ALLOWNEGATIVE equal to 1 as opposed to 0 will only be useful in specific cases.

The last vector element COHORT_NR of NumericalOptions defines the number of cohorts making up the equilibrium population output. During computations of the equilibrium a number of output files will be generated (see section 4.4.5), one of which is a file containing the population equilibrium state for each parameter that the equilibrium values are computed for. The vector element COHORT_NR specifies how many cohorts should be used to represent these equilibrium population states. Larger values will generate more detailed representations of the equilibrium population state at the expense of larger file sizes.

Code block 4.3.1.3

EnvironmentState <- c(R = "GENERALODE", P = "PERCAPITARATE", Bv = "POPULATIONINTEGRAL")

The first environment variable represents the resource density, which follows semi-chemostat growth in the absence of consumers. As a consequence, the resource \(R\) does NOT have an equilibrium value \(\tilde{R}=0\). The name of the first environment variable is therefore defined as R and its value is set to "GENERALODE". The second environment variable represents the predator density, the dynamics of which is described by the ODE

\[\dfrac{dP}{dt}\;=\;\left(\epsilon\dfrac{aB}{1+T_hB}\:-\:\delta\right)\,P\]

Because \(P=0\) is a regular fixed point of this ODE, the type of the second environment variable is set to "PERCAPITARATE" in the code box above. This environmental variable is given the name ‘P’. Finally, the third environment variable in the PNAS model is the total biomass of small juvenile consumers \(B\), which exerts a direct density-dependent effect on the mortality rate of small juvenile consumer individuals themselves, as it influences the predator functional response. The type of the third environment variable is therefore set to "POPULATIONINTEGRAL" and its name is defined as Bv in the code box above.

Code block 4.3.1.4

DefaultParameters <- c(Rho = 0.1, Rmax = 3.0E-4, Lb = 7.0, Lv = 27.0, Lj = 110.0, Lm = 300.0, Omega = 9.0E-6, Imax = 1.0E-4, Rh = 1.5E-5, Gamma = 0.006, Rm = 0.003, Mub = 0.01, A = 5000.0, Th = 0.1, Epsilon = 0.5, Delta = 0.01)

The names of the vector elements (parameters) can be used in the programming of the life history functions of the model and are furthermore used to make the output files produced by the program more readable. These output files contain a small header text indicating among other details which parameter values were used for the computation of the results contained in the output file. In this report the parameter names are listed together with their value.

Code block 4.3.1.5

StateAtBirth <- function(E, pars)
{
  with(as.list(c(E, pars)),{
        # We model a single structured population with two i-state variables:
        # 1: age (initial value 0); 2: length (initial value equal to parameter Lb)
        c(Age = 0.0, Length = Lb)
      })
}

Code block 4.3.1.6

LifeStageEndings <- function(lifestage, istate, birthstate, BirthStateNr, E, pars) {
  with(as.list(c(E, pars, istate)),{
        maturation  = switch(lifestage, Length - Lv, Length - Lj, -1)
      })
}

In the model there is a discontinuous change at the length threshold \(\ell=\ell_v\) when individuals turn from vulnerable to completely invulnerable to predation. The value that indicates the end of the first life stage (when lifestage equals 1) is hence set to \(\ell-\ell_v\). Furthermore, individuals mature at \(\ell=\ell_j\), which changes their fecundity discontinuously from a 0 value just before maturation to a positive value just after maturation. The value that indicates the end of the second (juvenile) stage (when lifestage equals 2) is hence set to \(\ell-\ell_j\).

The threshold value returned to the program in maturation will in general depend on the individual state variables, possibly on the individual’s state-at-birth and will be different for individuals in different life stages. For this reason, the function LifeStageEndings() has as arguments lifestage, specifying the life stage that the individual is currently in, istate, the individual state, and birthstate, the individual’s state-at-birth. In addition, the threshold value marking the end of a particular stage may depend on the value of the environment variables contained in E and the parameter values contained in pars, which are hence also passed arguments to the function LifeStageEndings().

This routine will be called as many times as there are possible states-at-birth, because the state-at-birth may influence the threshold between consecutive life stages. The same holds for the next 2 routines discussed, which define changes in the i-state variables, the fecundity, mortality and the impact of individuals. In essence, individuals with different states-at-birth are treated as constituting subpopulations within the same structured population. Because of the possible dependence on the state-at-birth the variables birthstate and BirthStateNr are passed as arguments to the function LifeStageEndings() as well as to the 2 functions discussed below. These arguments contain the values of the i-state variables and the index in the set \(\left\{\boldsymbol{\phi}_1,\ldots,\boldsymbol{\phi}_m\right\}\), respectively, for which the routine is invoked and for which the threshold between consecutive life stages has to be evaluated.

Code block 4.3.1.7

LifeHistoryRates <- function(lifestage, istate, birthstate, BirthStateNr, E, pars) {
  with(as.list(c(E, pars, istate)),{
    list(
      # We model a single structured population (nrow=1) with two i-state variables:
      # 1: age (developmental rate 1.0); 2: Length (vonBertalanffy growth rate)
      development = c(1.0, Gamma*(Lm*R/(R + Rh) - Length)),
      fecundity   = switch(lifestage, 0, 0, Rm*R/(R + Rh)*Length^2),
      mortality   = switch(lifestage, Mub + A*P/(1+A*Th*Bv), Mub, Mub),
      impact      = switch(lifestage, c(Imax*R/(R + Rh)*Length^2, Omega*Length^3, 0, 0),
                                      c(Imax*R/(R + Rh)*Length^2, 0, Omega*Length^3, 0),
                                      c(Imax*R/(R + Rh)*Length^2, 0, 0, Omega*Length^3))
    )
  })
}

The first individual state variable in the model corresponds to the individual age, which obviously has a rate of development equal to 1. The rate of development in individual length, the second individual state variable, follows the vonBertalanffy growth function, as specified by the function \(g(\ell, R)\) (refer to the model formulation in section 4.2).

The code fragment above implements a non-zero fecundity for individuals in the third life stage (when lifestage equals 3), which corresponds to the adult individuals with \(\ell>\ell_j\). The implemented expression corresponds to the function \(\beta(\ell, R)=r_mR/(R_h+R)\ell^2\) as assumed in the PNAS model (see section 4.2).

In the PNAS model all individuals experience a background mortality rate \(\mu_b\), while small juvenile individuals, which are in the first distinguished life stage (when lifestage equals 1), experience on top of the background mortality a predation mortality equal to \(aP/(1+T_hB)\) as expressed by the function \(\mu(\ell,P)\) in section 4.2.

In the PNAS model the impact of an individual consumer on its environment consists of its grazing rate \(I(\ell,R)\) and the biomass-length relation \(\omega\ell^3\) determining the biomass of vulnerable consumers, as it represents food for predators. The grazing rate of an individual consumer in the PNAS model is independent of the life stage it is in. As is shown in the code box above, this impact is assigned to the first interaction variable. The code box above furthermore shows that biomass of juvenile consumers that are vulnerable to predation is assigned to the second interaction variable, whereas the biomass of the invulnerable juvenile and adult consumers is assigned to the third and fourth interaction variable, respectively. These last two interaction variables are obviously not needed for the specification of the equilibrium, but are only included as additional output.

Refer to the remarks in the discussion of the function LifeStageEndings() concerning the dependence on the individual’s state-at-birth.

Code block 4.3.1.8

DiscreteChanges <- function(lifestage, istate, birthstate, BirthStateNr, E, pars) {
  with(as.list(c(E, pars)),{
       # No discrete changes in this problem, function is commented out, which 
       # would be equivalent to returning a copy of the input argument 'istate'
       istate
     })
}

If defined, the function DiscreteChanges() is called whenever a transition between two consecutive life stages is reached during the integration over the individual life history. The function should return a vector of length equal to the number of i-state variables. Each element should specify the value of the particular i-state variable after the transition to the current state. In case the model accounts for multiple, structured populations this function should return a matrix with the number of rows equal to the number of structured populations in the problem and the number of columns equal to the number of i-state variables.

It should be noted that the value of the variable lifestage indicates the life stage that is entered, that is, following the current stage boundary. This routine will hence never be called with a value of one of the elements lifestage equal to 1. The discrete changes in the individual state variables have to be implemented by assigning new values to the variables istate. These assignments may as before depend on the life stage that is entered, as specified by the variable lifestage, the (old values) of the individual state variables, contained in the argument istate, and possibly on the individual’s state-at-birth, specified in the argument birthstate. If no assignment of a value to istate is implemented, that particular individual state variable will keep its current value.

Code block 4.3.1.9

EnvEqui <- function(I, E, pars) {
  with(as.list(c(E, pars)),{
        c(Rho*(Rmax - R) - I[1], Epsilon*A*I[2]/(1+A*Th*I[2]) - Delta, I[2])
      })
}

The first environment variable in the model represents the resource density, which follows semi-chemostat growth in the absence of consumers. As a consequence, the resource \(R\) does NOT have an equilibrium value \(\tilde{R}=0\). Its equilibrium condition is therefore specified in the code box above as Rho*(Rmax - R) - I[1]. The first term in this expression implements the semi-chemostat dynamics in the absence of consumers and the quantity I[1] represents the grazing rate by the total population of consumers. These two quantities should cancel each other for the resource density to be in equilibrium.

The second environment variable represents the predator density, the dynamics of which is described by the ODE \[\dfrac{dP}{dt}\;=\;\left(\epsilon\dfrac{aB}{1+T_hB}\:-\:\delta\right)\,P\] Because \(P=0\) is a regular fixed point of this ODE, its per capita growth rate \(\epsilon aB/(1+T_hB)-\delta\) is used to define its equilibrium condition. Notice in this respect that the population feedback quantity I[2] represents \(B\), the total biomass of small juvenile consumers that are vulnerable to predation, which is calculated from the individual-level impact quantity impact[2] defined in code box 4.3.1.7.

Finally, the third environment variable in the PNAS model is the total biomass of small juvenile consumers \(B\), which exerts a direct density-dependent effect on the mortality rate of small juvenile consumer individuals themselves, as it influences the predator functional response. The equilibrium condition of this third environment variable is specified by identifying it with the population feedback quantity I[2], which the program computes on the basis of the individual-level impact quantity impact[2] that represents the biomass of an individual consumer that is in the first life stage, where it is vulnerable to predation (see the definition of impact in code box 4.3.1.7).

In case of multiple states-at-birth, individuals with different states-at-birth may have different impacts on their environment. The expected impact of an individual during its entire life on its environment is then given by an integral of the form: \[\Gamma_j\;=\;\int_0^\infty\gamma(\boldsymbol{\chi}(a, \boldsymbol{\phi}_j), \boldsymbol{\phi}_j, E)\,\mathcal{F}_j(a)\,da\] in which \(\gamma(\boldsymbol{\chi}(a, \boldsymbol{\phi}_j), \boldsymbol{\phi}_j, E)\) quantifies the impact of an individual that is born with state \(\boldsymbol{\phi}_j\) state dependent on its individual state at age \(a\), \(\boldsymbol{\chi}(a, \boldsymbol{\phi}_j)\), its state-at-birth \(\boldsymbol{\phi}_j\) and the environment variables. \(\mathcal{F}_j(a)\) now represents the probability that an individual born with state-at-birth \(\boldsymbol{\phi}_j\) survives until age \(a\), which is defined as: \[\mathcal{F}_j(a)\;=\;\exp\left(-\int_0^a\mu(\boldsymbol{\chi}(a, \boldsymbol{\phi}_j), \boldsymbol{\phi}_j, E)\,da\right)\] with \(\mu(\boldsymbol{\chi}(a, \boldsymbol{\phi}_j), \boldsymbol{\phi}_j, E)\) the individual’s instantaneous mortality rate.

If the possible states-at-birth are given by the set \(\left\{\boldsymbol{\phi}_1,\ldots,\boldsymbol{\phi}_m\right\}\), the individual-level impact is an \(m\)-dimensional vector \(\Gamma=\left(\Gamma_1,\ldots,\Gamma_m\right)\). The population-level impact on the environment in this case equals the dot product of this vector \(\Gamma\) with the \(m\)-dimensional vector \(\tilde{b}\), representing the equilibrium distribution of produced offspring over the possible states-at-birth \(\left\{\boldsymbol{\phi}_1,\ldots,\boldsymbol{\phi}_m\right\}\) (refer to Diekmann, Gyllenberg, and Metz (2003) for details).

The program automatically computes the equilibrium distribution of produced offspring over the possible states-at-birth and uses it to compute the population-level impacts contained in I from the individual-level impacts that have been specified in impact.

Code block 4.3.2.1

// Dimension settings: Required
#define POPULATION_NR       1
#define STAGES              3
#define I_STATE_DIM         2
#define ENVIRON_DIM         3
#define INTERACT_DIM        4
#define PARAMETER_NR        16

// Numerical settings: Optional (default values adopted otherwise)
#define MIN_SURVIVAL        1.0E-9          // Survival at which individual is considered dead
#define MAX_AGE             100000          // Give some absolute maximum for individual age

#define DYTOL               1.0E-7          // Variable tolerance
#define RHSTOL              1.0E-6          // Function tolerance

#define ALLOWNEGATIVE       0               // Negative solution values allowed?
#define COHORT_NR           100             // Number of cohorts in state output

The variable STAGES has to be defined equal to the number of life stages that can be distinguished in the individual life history (line 3 in the code box above). While integrating the ODEs for the individual life history numerical problems may occur when the right hand side of the ODEs changes abruptly in value at a certain threshold value of the individual state, as a consequence of discontinuities in the development rate, the mortality rate or the fecundity. Each of such thresholds in the life history should be distinguished as a stage boundary. In the PNAS model the mortality changes discontinuously at \(\ell=\ell_v\), while the fecundity changes discontinuously at \(\ell=\ell_j\). Three life stages can hence be distinguished: vulnerable juveniles, invulnerable juveniles and adults, and the changes in the life history rates are indeed abrupt at the transition boundaries between these stages. The variable STAGES is therefore set equal to 3.

The variable I_STATE_DIM (line 4 in the code box above) defines the dimension of the individual state. For the PNAS model this is defined equal to 2 to account for both individual age and individual length.

The variable ENVIRON_DIM is required in nonlinear PSPM, whereas it is optional for demographic analysis of linear PSPMs. It represents the number of environment variables that determine the life history of an individual. In the PNAS model the growth and fecundity of individual consumers are functions of the resource density \(R\), whereas the mortality is a function of the predator density \(P\) and the biomass of vulnerable consumers \(B\). The latter only influences the mortality of the vulnerable consumers, because it determines the value of the predator functional response. ENVIRON_DIM hence equals 3 in the PNAS model.

The variable INTERACT_DIM defines the number of functions that represent the impact of an individual on its environment. In the PNAS model this feedback of an individual consumer on its environment consists of its grazing rate \(I(\ell,R)\) and the biomass-length relation \(\omega\ell^3\) determining the biomass of vulnerable consumers, as it represents food for predators. Therefore, the variable INTERACT_DIM should be at least set equal to 2. However, all interaction functions are also saved to the output file generated during a computation. The interaction functions can hence be conveniently used to produce arbitrary output quantities of the form

\[\int_\Omega h(\boldsymbol{\chi}, \boldsymbol{\chi}_b, E)\,\tilde{n}(\boldsymbol{\chi})\,d\boldsymbol{\chi}\]

where \(h(\boldsymbol{\chi}, \boldsymbol{\chi}_b, E)\) is an interaction (weighing) function that can depend on the values of the individual state \(\boldsymbol{\chi}\), the state-at-birth of individuals \(\boldsymbol{\chi}_b\) and the values of the environment variables \(E\), and \(\tilde{n}(\boldsymbol{\chi})\) is the stable population distribution in equilibrium. Such quantities can therefore represent the total population density in equilibrium (when \(h(\boldsymbol{\chi}, \boldsymbol{\chi}_b, E)=1\)), the total population biomass (when \(h(\boldsymbol{\chi}, \boldsymbol{\chi}_b, E)\) equals the biomass of an individual with individual state \(\boldsymbol{\chi}\) and state-at-birth \(\boldsymbol{\chi}_b\)) or the total population birth rate in equilibrium (when \(h(\boldsymbol{\chi}, \boldsymbol{\chi}_b, E)\) is the fecundity of an individual with individual state \(\boldsymbol{\chi}\) and state-at-birth \(\boldsymbol{\chi}_b\)). In the PNAS model I want in addition to the biomass of vulnerable consumers, also the biomass of non-vulnerable juvenile consumers and the biomass of adult consumers as output of the model and hence have set the variable INTERACT_DIM equal to 4.

The last required parameter that has to be specified is the number of parameters in the model (set in line 7 in the code box above). In the PNAS model this equals 16 (\(\rho\), \(R_{max}\), \(\ell_b\), \(\ell_v\), \(\ell_j\), \(\ell_m\), \(\omega\), \(I_{max}\), \(R_h\), \(\gamma\), \(r_m\), \(\mu_b\), \(a\), \(T_h\), \(\epsilon\) and \(\delta\)).

The remaining definitions in the code box are all optional and can be left away. A list of all possible variables that can be changed by a definition in this code section is provided in section 9.3. The variable MIN_SURVIVAL determines the threshold of the survival probability below which an individual is considered dead. The integration over the individual life history stops whenever the survival probability falls below this threshold value. In the code box above (line 10) the minimum survival is set to \(10^{-9}\), which is in fact the default value and is hence superfluous. Note that the value of MIN_SURVIVAL can not be set equal to 0. As an alternative to using 0 MIN_SURVIVAL can be set to a very small value like \(10^{-100}\).

The variable MAX_AGE (line 11 in code box above) can be used as an alternative to determine the end of an individual life and to stop the integration over the individual life history. In the PNAS model there is no maximum individual age and hence the variable is set to a very high value (100000), which the individuals will never reach, because before that age their survival probability has already dropped below its threshold value (\(10^{-9}\)).

The two quantities DYTOL and RHSTOL determine whether a solution has been found. In general, both demographic analysis as well as equilibrium analysis of PSPMs boils down to solving a system of nonlinear equations that can be represented as \(G(y)=0\) for a set of unknowns \(y\) in iterative manner. The subsequent estimates of the solution in the Newton iterations can be labeled as \(y_p\) and \(y_{p+1}\). A solution is now considered to be located if both of the following conditions hold:

\[\begin{align*} &\|{y_{p+1}-y_p}\| < \epsilon_y\\[2ex] &\|{G(y_{p+1})}\| < \epsilon_G \end{align*}\]

where \(\|.\|\) refers to the Euclidean norm. DYTOL and RHSTOL are the quantities \(\epsilon_y\) and \(\epsilon_G\), respectively. Increasing (decreasing) their value leads to easier (harder) acceptance of a set of unknowns as a solution to the system of equations \(G(y)=0\). The definition of these two accuracies in the code box is in fact superfluous as they are defined equal to their default values (see section 9.3).

The quantity ALLOWNEGATIVE is a flag that can only have a value of 0 or 1 and determines whether or not computations should stop when one of the variables to solve for reaches a negative value. In most population models negative solution values are biologically not relevant and ALLOWNEGATIVE is hence set to 0 by default. Line 16 in the code box above is only included to illustrate the use of ALLOWNEGATIVE and does not change the value of this variable from its default value. Most likely, setting ALLOWNEGATIVE equal to 1 as opposed to 0 will only be useful in specific cases.

The last quantity COHORT_NR defines the number of cohorts making up the equilibrium population output. During computations of the equilibrium a number of output files will be generated (see section 4.4.5), one of which is a file containing the population equilibrium state for each parameter that the equilibrium values are computed for. COHORT_NR specifies how many cohorts should be used to represent these equilibrium population states. Larger values will generate more detailed representations of the equilibrium population state at the expense of larger file sizes.

Code block 4.3.2.2

// Descriptive names of parameters in parameter array (at least two parameters are required)
char  *parameternames[PARAMETER_NR] =
    { "Rho", "Rmax", "Lb", "Lv", "Lj", "Lm", "Omega", "Imax", "Rh", "Gamma", "Rm", "Mub",
      "A", "Th", "Epsilon", "Delta"};

// Default values of all parameters
double  parameter[PARAMETER_NR] =
    { 0.1, 3.0E-4, 7.0, 27.0, 110.0, 300.0, 9.0E-6, 1.0E-4, 1.5E-5, 0.006, 0.003, 0.01,
      5000.0, 0.1, 0.5, 0.01};

Model parameter values are stored by the program in the vector variable parameter. The lines 2-4 above assign each of the elements this vector a more meaningful, model-specific name. These name strings can not be used in the remaining parts of the model implementation, they only serve to make the output files produced by the program more readable. These output files contain a small header text indicating among other details which parameter values were used for the computation of the results contained in the output file (see section 4.4.5). In the output file parameters are referred to with their names as defined in the array of strings *parameternames[\(k\)]. To adapt the above code to a different model, the code on line 2 of the code box above should remain the same, only change lines 3-4 as needed (possibly extending it over more lines in case there are many parameters).

The default values to use for the model parameters are specified by the declaration of the vector parameter[PARAMETR_NR] on line 7-9 of the previous code box. The values should be specified as a comma-separated array of values within braces (don’t forget the closing semi-colon at the end of the statement that is required in the C-language!). To adapt the above code to a different model, the code on line 7 of the code box above should remain the same, only change lines 8-9 as needed (possibly extending it over more lines in case there are many parameters).

Code block 4.3.2.3

// Aliases definitions for all istate variables
#define AGE                 istate[0][0]
#define LENGTH              istate[0][1]

// Aliases definitions for all environment variables
#define R                   E[0]
#define P                   E[1]
#define B                   E[2]

// Aliases definitions for all parameters
#define RHO                 parameter[ 0]   // Default: 0.1
#define RMAX                parameter[ 1]   // Default: 3.0E-4

#define LB                  parameter[ 2]   // Default: 7
#define LV                  parameter[ 3]   // Default: 27
#define LJ                  parameter[ 4]   // Default: 110
#define LM                  parameter[ 5]   // Default: 300

#define OMEGA               parameter[ 6]   // Default: 9.0E-6

#define IMAX                parameter[ 7]   // Default: 1.0E-4
#define RH                  parameter[ 8]   // Default: 1.5E-5

#define GAMMA               parameter[ 9]   // Default: 0.006
#define RM                  parameter[10]   // Default: 0.003

#define MUB                 parameter[11]   // Default: 0.01

#define A                   parameter[12]   // Default: 5000.0
#define TH                  parameter[13]   // Default: 0.1
#define EPSILON             parameter[14]   // Default: 0.5
#define DELTA               parameter[15]   // Default: 0.01

The value of the environment variables are contained in an array E[\(e\)] with \(e\) an index in the range \(0\ldots\) ENVIRON_DIM-1. Again, it is up to the user to keep track of which index pertains to which environmental state variable. The use of aliases is really beneficial for this purpose. As defined in code box 4.3.2.1 three environment variables are identified in the PNAS model: the resource density, the density of predators and the biomass of juvenile consumers that are vulnerable to predation. Lines 6-8 in the code box above identifies these with the first, second and third element of the array E[\(e\)], respectively, and introduces the aliases R, P and B for these quantities. All code shown below will make use of these aliases as opposed to their real names in the program (E[0], E[1] and E[2]).

Similar arguments as given above for the individual state variables contained in the array istate[\(p\)][\(i\)] and the environment variables contained in the array E[\(e\)] hold for the model parameters. All model parameters are contained in a vector named parameter[\(k\)] in the code (see the previous section) with \(k\) an index in the range \(0\ldots\) PARAMETER_NR-1. Which element of this vector represents which model-specific parameter is up to the user. To prevent mixing up the interpretation of the different vector elements and hence to prevent mistakes, it is strongly advised to define meaningful, model-specific aliases for each of the elements of the vector parameter[\(k\)] as is illustrated in lines 11-32 in the code box above. It is best to avoid completely the direct use of the program variable parameter[\(k\)] in any part of the model specification and only use the models-specific aliases.

As can be seen in the code block above all aliases for program variables used in the C-implementation of the model are names in capitals. It is advisable to use only capitals when introducing these aliases (or global variables if they are needed) to avoid any conflict between these aliases and variables that defined elsewhere in any of the C files with numerical routines that are included in the package.

Code block 4.3.2.4

/*
 * Specify the number of states at birth for the individuals in all structured
 * populations in the problem in the vector BirthStates[].
 */

void SetBirthStates(int BirthStates[POPULATION_NR], double E[])
{
  BirthStates[0] = 1;

  return;
}

For each population with index \(p\) the variable BirthStates[\(p\)] has to be set to the number of possible states at birth. Because in the PNAS model all individuals are born with age 0 and length \(\ell=\ell_b\), all individuals have the same, unique state at birth and hence BirthStates[\(0\)] is set to 1.

Note that different populations may have different numbers of states-at-birth. BirthStates[\(p\)] hence does not need to be the same for all \(p\).

Code block 4.3.2.5

/*
 * Specify all the possible states at birth for all individuals in all
 * structured populations in the problem. BirthStateNr represents the index of
 * the state of birth to be specified. Each state at birth should be a single,
 * constant value for each i-state variable.
 *
 * Notice that the first index of the variable 'istate[][]' refers to the
 * number of the structured population, the second index refers to the
 * number of the individual state variable. The interpretation of the latter
 * is up to the user.
 */

void StateAtBirth(double *istate[POPULATION_NR], int BirthStateNr, double E[])
{
  AGE = 0.0;
  LENGTH = LB;

  return;
}

For every population (\(p=0, 1,\ldots,\) POPULATION_NR-1) the value of each individual state variable with index \(i\), istate[\(p\)][\(i\)] (\(i=0, 1,\ldots,\) I_STATE_DIM-1), has to be assigned a unique value, from which individual development will start at age 0. Notice that the program does not automatically include individual age in its characterization of the individual state, even though integration over the entire life history (as a function of age) is carried out. For the PNAS model age at birth is (obviously) set to 0, while the length at birth is given by the parameter \(\ell_b\) (line 15 and 16, respectively in the code box below).

This routine will be called as many times as there are possible states-at-birth. The variable BirthStateNr indicates the index \(j\) of the state-at-birth in the set \(\left\{\boldsymbol{\phi}_1,\ldots,\boldsymbol{\phi}_m\right\}\) for which the values have to be set in the current invocation of the routine. The routine will thus be called with BirthStateNr set equal to a value in \(0, 1, \dots, m-1\) (Remember the starting index 0 in C!). If there are multiple states-at-birth (BirthStates[\(p\)] > 1) the definition of the values of the i-state variables has to depend explicitly on the index BirthStateNr to make the states-at-birth different from each other. Furthermore, if the problem involves multiple structured populations the number of possible states-at-birth can be different for each of them, which might lead to a situation that the routine above is called with a value of the index BirthStateNr that is larger than the maximum number of states-at-birth for a particular population (BirthStateNr \(\geq\) BirthStates[\(p\)]). The program safely ignores such inappropriate state-at-birth specifications.

Code block 4.3.2.6

/*
 * Specify the threshold determining the end point of each discrete life
 * stage in individual life history as function of the i-state variables and
 * the individual's state at birth for all populations in every life stage.
 *
 * Notice that the first index of the variable 'istate[][]' refers to the
 * number of the structured population, the second index refers to the
 * number of the individual state variable. The interpretation of the latter
 * is up to the user.
 */

void IntervalLimit(int lifestage[POPULATION_NR], double *istate[POPULATION_NR],
                   double *birthstate[POPULATION_NR], int BirthStateNr, double E[],
                   double limit[POPULATION_NR])
{
  switch (lifestage[0])
  {
    case 0:
      limit[0] = LENGTH - LV;
      break;
    case 1:
      limit[0] = LENGTH - LJ;
      break;
  }

  return;
}

In this routine the variable limit[\(p\)] has to be defined, which has as many elements as there are populations (\(p=0\ldots\) POPULATION_NR-1). The life stage that the individual is in at the moment this routine is called, is determined by the variable lifestage[\(p\)], which has a value of 0 if the individual is in the first life stage and a value of STAGES-1 if it is in the last life stage. The element limit[\(p\)] should now indicate when the current life stage as given in lifestage[\(p\)] ends. In particular, the program considers the current life stage to end when limit[\(p\)] turns from negative to positive. For the end of the last life stage the death of old age, either by reaching the maximum age MAX_AGE or by reaching the minimum survival threshold MIN_SURVIVAL, does not have to be specified separately, the program takes care of that automatically.

The threshold value that has to be stored and returned to the program in limit[p] will depend on the individual state variables, possibly on the individual’s state-at-birth and will be different for individuals in different life stages. For this reason, the routine IntervalLimit() has as arguments lifestage[], specifying the life stage that the individual is currently in, istate[][], the individual state, birthstate[][] and BirthStateNr, the value and index of the individual’s state-at-birth, respectively. The threshold value marking the end of a particular stage may, however, in addition depend on the value of the environment variables (E[]).

In the PNAS model there is a discontinuous change at the length threshold \(\ell=\ell_v\) when individuals turn from vulnerable to completely invulnerable to predation. The value that indicates the end of the first life stage (when lifestage[\(p\)] \(=0\)) is hence set to \(\ell-\ell_v\) (line 18-20 in the code box above). Furthermore, individuals mature at \(\ell=\ell_j\), which changes their fecundity discontinuously from a 0 value just before maturation to a positive value just after maturation. The value that indicates the end of the second (juvenile) stage (when lifestage[\(p\)] \(=1\)) is hence set to \(\ell-\ell_j\) (line 21-23 in the code box above).

Like the previous routine, this routine will be called as many times as there are possible states-at-birth, because the state-at-birth may influence the threshold between consecutive life stages. The same holds for the routines discussed in sections 4.3.2.7-4.3.2.11 below, which define changes in the i-state variables, the fecundity and the mortality of individuals and their impact on the environment, respectively. In essence, individuals with different states-at-birth are treated as subpopulations within the same structured population. Because of the possible dependence on the state-at-birth the variables birthstate[][] and BirthStateNr are passed as arguments to this routine and the once discussed in sections 4.3.2.7-4.3.2.11. These arguments contain the values of the i-state variables and the index in the set \(\left\{\boldsymbol{\phi}_1,\ldots,\boldsymbol{\phi}_m\right\}\), respectively, for which the routine is invoked and for which the threshold between consecutive life stages has to be evaluated.

If the problem involves multiple structured populations and the number of possible states-at-birth differs among them, the routine above may be called with a value of the index BirthStateNr that is larger than the maximum number of states-at-birth for a particular population (BirthStateNr \(\geq\) BirthStates[\(p\)]). Although this circumstance may seem confusing, the user does not have to worry about it, as the program is designed to safely ignore such assignments of thresholds between consecutive life stages, changes in the i-state variables, fecundity, mortality and impact on the environment of individuals for states-at-birth with index BirthStateNr \(\geq\) BirthStates[\(p\)]that are inappropriate for the structured population with index \(p\).

Code block 4.3.2.7

/*
 * Specify the development of individuals as a function of the i-state
 * variables and the individual's state at birth for all populations in every
 * life stage.
 *
 * Notice that the first index of the variables 'istate[][]' and 'development[][]'
 * refers to the number of the structured population, the second index refers
 * to the number of the individual state variable. The interpretation of the
 * latter is up to the user.
 */

void Development(int lifestage[POPULATION_NR], double *istate[POPULATION_NR],
                 double *birthstate[POPULATION_NR], int BirthStateNr, double E[],
                 double development[POPULATION_NR][I_STATE_DIM])
{
  development[0][0] = 1.0;
  development[0][1] = GAMMA*(LM*R/(R + RH) - LENGTH);

  return;
}

For each individual state variable \(i\) of every structured population \(p\) that is part of the individual state istate[\(p\)][\(i\)], its rate of development during the life history has to be specified in development[\(p\)][\(i\)]. Notice that these development rates may differ in different life stages, for example growth in body size may be different for juveniles and adults in case adults invest a lot of energy into reproduction. The development rates should then be specified dependent on the current life stage the individual is in. This current life stage at the moment the routine is evaluated is contained in the variable lifestage[\(p\)]. The development rate may furthermore depend on the individual state variables, on the individual’s state-at-birth and on the value of the environment variables, which is the reason for istate[][], the individual state, birthstate[][] and BirthStateNr, the value and index of the individual’s state-at-birth, respectively, and E[], the environment variables, as arguments to this routine.

In the PNAS model the first individual state variable corresponds to the individual age, which obviously has a rate of development equal to 1. The rate of development in individual length, the second individual state variable, follows the vonBertalanffy growth function, as specified by the function \(g(\ell, R)\) (refer to the model formulation at the start of this chapter).

Refer to the remarks in section 4.3.2.6 concerning the dependence on the individual’s state-at-birth.

Code block 4.3.2.8

/*
 * Specify the possible discrete changes (jumps) in the individual state
 * variables when ENTERING the stage specified by 'lifestage[]'.
 *
 * Notice that the first index of the variable 'istate[][]' refers to the
 * number of the structured population, the second index refers to the
 * number of the individual state variable. The interpretation of the latter
 * is up to the user.
 */

void DiscreteChanges(int lifestage[POPULATION_NR], double *istate[POPULATION_NR],
                     double *birthstate[POPULATION_NR], int BirthStateNr, double E[])
{
  return;
}

This routine is not relevant in case of the PNAS model and hence its contents are empty (apart for the necessary return; statement).

This routine is called whenever a transition between two consecutive life stages is reached during the integration over the individual life history. It should be noted that the value of the variable lifestage[\(p\)] indicates the life stage that is entered, that is, following the current stage boundary. This routine will hence never be called with a value of one of the elements lifestage[\(p\)] equal to 0. The discrete changes in the individual state variables have to be implemented by assigning new values to the variables istate[\(p\)][\(i\)]. These assignments may as before depend on the life stage that is entered, as specified by the variable lifestage[], the (old values) of the individual state variables, contained in the argument istate[][], the individual’s state-at-birth, determined by the arguments birthstate[][] and BirthStateNr, and on the environment variables E[]. If no assignment of a value to istate[\(p\)][\(i\)] is implemented, that particular individual state variable will keep its current value.

Refer to the remarks in section 4.3.2.6 concerning the dependence on the individual’s state-at-birth.

Code block 4.3.2.9

/*
 * Specify the fecundity of individuals as a function of the i-state
 * variables and the individual's state at birth for all populations in every
 * life stage.
 *
 * The number of offspring produced has to be specified for every possible
 * state at birth in the variable 'fecundity[][]'. The first index of this
 * variable refers to the number of the structured population, the second
 * index refers to the number of the birth state.
 *
 * Notice that the first index of the variable 'istate[][]' refers to the
 * number of the structured population, the second index refers to the
 * number of the individual state variable. The interpretation of the latter
 * is up to the user.
 */

void Fecundity(int lifestage[POPULATION_NR], double *istate[POPULATION_NR],
               double *birthstate[POPULATION_NR], int BirthStateNr, double E[],
               double *fecundity[POPULATION_NR])
{
  fecundity[0][0] = 0.0;
  if (lifestage[0] == 2)
    fecundity[0][0] = RM*R/(R + RH)*LENGTH*LENGTH;

  return;
}

In this routine not only the fecundity (i.e. the number of offspring produced per unit time) has to be specified, but also the state-at-birth of the produced offspring. Therefore, this routine has to assign values to the matrix fecundity[\(p\)][\(j\)], which determines for the population with index \(p\) the number of offspring produced per unit time with state-at-birth with index \(j\) in the set \(\left\{\boldsymbol{\phi}_1,\ldots,\boldsymbol{\phi}_m\right\}\). This fecundity will certainly depend on the life stage that the individual is in (only adults reproduce), which is contained in the argument lifestage[], and on the individual state variables, i.e. the values of the argument istate[][]), but possibly also on the individual’s state-at-birth, the values and index of which are specified by the arguments birthstate[][] and BirthStateNr, respectively, and on the value of environment variables, provided by the argument E[], at the moment this routine is called.

In the most common case of a unique state-at-birth and a single structured population, like in the PNAS model, the only valid indices are \(p=0\) and \(j=0\) and hence only the variable fecundity[0][0] has to be assigned. The code fragment above implements a non-zero fecundity for individuals in the third life stage (lifestage[0] == 2), which corresponds to the adult individuals with \(\ell>\ell_j\). The implemented expression corresponds to the function \(\beta(\ell, R)=r_mR/(R_h+R)\ell^2\) as assumed in the PNAS model (see section 4.2). The code provides a good example of how to assign a different value for a particular life history rate dependent on the life stage that an individual is in. The same approach can also be used in the other routines specifying the life history rates of individuals.

For more detailed remarks about models with multiple states-at-birth consult section 4.3.2.6.

Code block 4.3.2.10

/*
 * Specify the mortality of individuals as a function of the i-state
 * variables and the individual's state at birth for all populations in every
 * life stage.
 *
 * Notice that the first index of the variable 'istate[][]' refers to the
 * number of the structured population, the second index refers to the
 * number of the individual state variable. The interpretation of the latter
 * is up to the user.
 */

void Mortality(int lifestage[POPULATION_NR], double *istate[POPULATION_NR],
               double *birthstate[POPULATION_NR], int BirthStateNr, double E[],
               double mortality[POPULATION_NR])
{
  if (lifestage[0] == 0)
    mortality[0] = MUB + A*P/(1+A*TH*B);
  else
    mortality[0] = MUB;

  return;
}

For each population the corresponding element of the array mortality[\(p\)] should be assigned the mortality rate, possibly dependent on the life stage the individual is in at the moment this routine is called (given in argument lifestage[]), the current i-state of the individual (given in argument istate[][]), the individual’s state-at-birth (values and index in the set \(\left\{\boldsymbol{\phi}_1,\ldots,\boldsymbol{\phi}_m\right\}\) given by birthstate[][] and BirthStateNr, respectively) and the value of environment variables (E[]).

In the PNAS model all individuals experience a background mortality rate \(\mu_b\), while small juvenile individuals, which are in the first distinguished life stage (lifestage[0] == 0), experience on top of the background mortality a predation mortality equal to \(aP/(1+T_hB)\) as expressed by the function \(\mu(\ell,P)\) in section 4.2.

Refer to the remarks in section 4.3.2.6 concerning the dependence on the individual’s state-at-birth.

Code block 4.3.2.11

/*
 * For all the integrals (measures) that occur in interactions of the
 * structured populations with their environments and for all the integrals
 * that should be computed for output purposes (e.g. total juvenile or adult
 * biomass), specify appropriate weighing function dependent on the i-state
 * variables, the individual's state at birth, the environment variables and
 * the current life stage of the individuals. These weighing functions should
 * be specified for all structured populations in the problem. The number of
 * weighing functions is the same for all of them.
 *
 * Notice that the first index of the variables 'istate[][]' and 'impact[][]'
 * refers to the number of the structured population, the second index of the
 * variable 'istate[][]' refers to the number of the individual state variable,
 * while the second index of the variable 'impact[][]' refers to the number of
 * the interaction variable. The interpretation of these second indices is up
 * to the user.
 */

void Impact(int lifestage[POPULATION_NR], double *istate[POPULATION_NR],
            double *birthstate[POPULATION_NR], int BirthStateNr, double E[],
            double impact[POPULATION_NR][INTERACT_DIM])
{
  impact[0][0] = IMAX*R/(R + RH)*LENGTH*LENGTH;

  switch (lifestage[0])
    {
    case 0:
      impact[0][1] = OMEGA*LENGTH*LENGTH*LENGTH;
      impact[0][2] = 0;
      impact[0][3] = 0;
      break;
    case 1:
      impact[0][1] = 0;
      impact[0][2] = OMEGA*LENGTH*LENGTH*LENGTH;
      impact[0][3] = 0;
      break;
    case 2:
      impact[0][1] = 0;
      impact[0][2] = 0;
      impact[0][3] = OMEGA*LENGTH*LENGTH*LENGTH;
      break;
    }

  return;
}

As explained in section 4.3.2.1 interaction functions are of the form

\[\int_\Omega h(\boldsymbol{\chi}, \boldsymbol{\chi}_b, E)\,\tilde{n}(\boldsymbol{\chi})\,d\boldsymbol{\chi}\]

where \(h(\boldsymbol{\chi}, \boldsymbol{\chi}_b, E)\) is an interaction (weighing) function that can depend on the values of the individual state \(\boldsymbol{\chi}\), the state-at-birth of individuals \(\boldsymbol{\chi}_b\) and the values of the environment variables \(E\), and \(\tilde{n}(\boldsymbol{\chi})\) is the stable population distribution in equilibrium. Such quantities can therefore represent the total population density in equilibrium (when \(h(\boldsymbol{\chi}, \boldsymbol{\chi}_b, E)=1\)), the total population biomass (when \(h(\boldsymbol{\chi}, \boldsymbol{\chi}_b, E)\) equals the biomass of an individual with individual state \(\boldsymbol{\chi}\) and state-at-birth \(\boldsymbol{\chi}_b\)) or the total population birth rate in equilibrium (when \(h(\boldsymbol{\chi}, \boldsymbol{\chi}_b, E)\) is the fecundity of an individual with individual state \(\boldsymbol{\chi}\) and state-at-birth \(\boldsymbol{\chi}_b\)). These interaction variables, in fact, determine the equilibrium of the model, because the nonlinearities that make an equilibrium possible arise through the impact of an individual on its environment.

As also explained in section 4.3.2.1 all interaction functions are saved to the output file when an equilibrium has been computed. Interaction functions are hence not only used to compute the density-dependent feedback in the model, but also to obtain model output quantities of the form shown above.

In the PNAS model this feedback of an individual consumer on its environment consists of its grazing rate \(I(\ell,R)\) and the biomass-length relation \(\omega\ell^3\) determining the biomass of vulnerable consumers, as it represents food for predators. The grazing rate of an individual consumer in the PNAS model is independent of the life stage it is in. As is shown in line 23 of the code box above, this impact is assigned to the first interaction variable (the one with index 0). Line 25-42 of the code box above show that biomass of juvenile consumers that are vulnerable to predation is assigned to the second interaction variable (lines 27-31), whereas the biomass of the invulnerable juvenile and adult consumers is assigned to the third (lines 32-36) and fourth (lines 37-41) interaction variable. These last two interaction variables are obviously not needed for the specification of the equilibrium, but are only included as additional output.

It should be pointed out that the routine Impact() should only specify the impact of an individual on its environment, given its current life stage (function argument lifestage[]), its individual state (argument istate[][]), its state-at-birth (values and index in the set \(\left\{\boldsymbol{\phi}_1,\ldots,\boldsymbol{\phi}_m\right\}\) given by birthstate[][] and BirthStateNr, respectively) and the value of environment variables (E[]). In other words, the routine should only specify the weighing function \(h(\boldsymbol{\chi}, \boldsymbol{\chi}_b, E)\). The program automatically translates this individual-level impact function to the feedback of the total population on its environment, as explained in the next section.

Refer to the remarks in section 4.3.2.6 concerning the dependence on the individual’s state-at-birth.

Code block 4.3.2.12

/*
 * Specify the type of each of the environment variables by setting
 * the entries in EnvironmentType[ENVIRON_DIM] to PERCAPITARATE, GENERALODE
 * or POPULATIONINTEGRAL based on the classification below:
 *
 * Set an entry to PERCAPITARATE if the dynamics of E[j] follow an ODE and 0
 * is a possible equilibrium state of E[j]. The ODE is then of the form
 * dE[j]/dt = P(E,I)*E[j], with P(E,I) the per capita growth rate of E[j].
 * Specify the equilibrium condition as condition[j] = P(E,I), do not include
 * the multiplication with E[j] to allow for detecting and continuing the
 * transcritical bifurcation between the trivial and non-trivial equilibrium.
 *
 * Set an entry to GENERALODE if the dynamics of E[j] follow an ODE and 0 is
 * NOT an equilibrium state of E. The ODE then has a form dE[j]/dt = G(E,I).
 * Specify the equilibrium condition as condition[j] = G(E,I).
 *
 * Set an entry to POPULATIONINTEGRAL if E[j] is a (weighted) integral of the
 * population distribution, representing for example the total population
 * biomass. E[j] then can be expressed as E[j] = I[p][i]. Specify the
 * equilibrium condition in this case as condition[j] = I[p][i].
 *
 * Notice that the first index of the variable 'I[][]' refers to the
 * number of the structured population, the second index refers to the
 * number of the interaction variable. The interpretation of the latter
 * is up to the user. Also notice that the variable 'condition[j]' should
 * specify the equilibrium condition of environment variable 'E[j]'.
 */

const int EnvironmentType[ENVIRON_DIM] = {GENERALODE, PERCAPITARATE, POPULATIONINTEGRAL};

void EnvEqui(double E[], double I[POPULATION_NR][INTERACT_DIM],
             double condition[ENVIRON_DIM])
{
  condition[0] = RHO*(RMAX - R) - I[0][0];
  condition[1] = EPSILON*A*I[0][1]/(1+A*TH*I[0][1]) - DELTA;
  condition[2] = I[0][1];

  return;
}

In the code box above, the array EnvironmentType[ENVIRON_DIM] has to define for each environment variable separately its type (PERCAPITARATE, GENERALODE or POPULATIONINTEGRAL). This array hence has as many elements as there are environment variables. Secondly, in the routine that follows the specification of EnvironmentType[ENVIRON_DIM] the equilibrium conditions have to be implemented as a function of the values of the environment variables itself E[] and the value of the population feedback functions I[][]. These latter two arrays are passed as arguments to the routine EnvEqui, while the equilibrium conditions have to be specified in the array condition[]. Notice that the ordering of the elements in the arrays condition[] and E[] are the same, meaning for example that the equilibrium condition for the second environment variable E[1] has to be returned in condition[1].

Notice that the feedback functions I[][] are the population-level representations of the individual-level impact functions impact[][] as they are defined in section 4.3.2.11. In case of a unique state-at-birth the expected life history is the same for all individuals in the population and hence not dependent on a state-at-birth \(\boldsymbol{\chi}_b\). In this case, such an expected impact of an individual during its entire life on its environment is given by an integral of the form:

\[\Gamma\;=\;\int_0^\infty\gamma(\boldsymbol{\chi}(a), E)\,\mathcal{F}(a)\,da\]

in which \(\gamma(\boldsymbol{\chi}(a), E)\) quantifies the impact dependent on the individual state and the environment variables, \(\boldsymbol{\chi}(a)\) is the individual state at age \(a\) and \(\mathcal{F}(a)\) represents the probability that the individual survives until age \(a\), which is defined as:

\[\mathcal{F}(a)\;=\;\exp\left(-\int_0^a\mu(\boldsymbol{\chi}(a), E)\,da\right)\]

with \(\mu(\boldsymbol{\chi}(a), E)\) the individual’s instantaneous mortality rate. The population-level impact on the environment now simply equals the product of the total population birth rate in equilibrium \(\tilde{b}\) and the individual-level impact:

\[\tilde{I}=\tilde{b}\,\Gamma\]

The elements of I[][] and impact[][] therefore correspond one-to-one, such that for example in the PNAS model, in which impact[0][0] is defined as the individual feeding rate on the resource, the quantity I[0][0] equals the grazing rate of the entire population on the resource.

In the PNAS model, the first environment variable represents the resource density, which follows semi-chemostat growth in the absence of consumers. As a consequence, the resource \(R\) does NOT have an equilibrium value \(\tilde{R}=0\). The type of the first environment variable is therefore set to GENERALODE and its equilibrium condition is specified in line 34 of the code box above as RHO*(RMAX - R) - I[0][0]. The first term in this line of C-code implements the semi-chemostat dynamics in the absence of consumers and the quantity I[0][0] represents the grazing rate by the total population of consumers as mentioned above.

The second environment variable represents the predator density, the dynamics of which is described by the ODE

\[\dfrac{dP}{dt}\;=\;\left(\epsilon\dfrac{aB}{1+T_hB}\:-\:\delta\right)\,P\]

Because \(P=0\) is a regular fixed point of this ODE, the type of the second environment variable is set to PERCAPITARATE on line 29 of the code box above, while its per capita growth rate \(\epsilon aB/(1+T_hB)-\delta\) is used to define its equilibrium condition on line 35. Notice in this respect that the population feedback quantity I[0][1] represents \(B\), the total biomass of small juvenile consumers that are vulnerable to predation, which is calculated from the individual-level impact quantity impact[0][1] defined on lines 27-31 in code box 4.3.2.11.

Finally, the third environment variable in the PNAS model is the total biomass of small juvenile consumers \(B\), which exerts a direct density-dependent effect on the mortality rate of small juvenile consumer individuals themselves, as it influences the predator functional response. The type of the third environment variable is therefore set to POPULATIONINTEGRAL on line 29 of the code box above. On line 36 of the code box above the equilibrium condition of this third environment variable is specified by identifying it with the population feedback quantity I[0][1], which the program computes on the basis of the individual-level impact quantity impact[0][1] that represents the biomass of an individual consumer that is in the first life stage, where it is vulnerable to predation (lines 27-31 in the code box above).

In case of multiple states-at-birth, individuals with different states-at-birth may have different impacts on their environment. The expected impact of an individual during its entire life on its environment is then given by an integral of the form:

\[\Gamma_j\;=\;\int_0^\infty\gamma(\boldsymbol{\chi}(a, \boldsymbol{\phi}_j), \boldsymbol{\phi}_j, E)\,\mathcal{F}_j(a)\,da\]

in which \(\gamma(\boldsymbol{\chi}(a, \boldsymbol{\phi}_j), \boldsymbol{\phi}_j, E)\) quantifies the impact of an individual that is born with state \(\boldsymbol{\phi}_j\) state dependent on its individual state at age \(a\), \(\boldsymbol{\chi}(a, \boldsymbol{\phi}_j)\), its state-at-birth \(\boldsymbol{\phi}_j\) and the environment variables. \(\mathcal{F}_j(a)\) now represents the probability that an individual born with state-at-birth \(\boldsymbol{\phi}_j\) survives until age \(a\), which is defined as:

\[\mathcal{F}_j(a)\;=\;\exp\left(-\int_0^a\mu(\boldsymbol{\chi}(a, \boldsymbol{\phi}_j), \boldsymbol{\phi}_j, E)\,da\right)\]

with \(\mu(\boldsymbol{\chi}(a, \boldsymbol{\phi}_j), \boldsymbol{\phi}_j, E)\) the individual’s instantaneous mortality rate. > > If the possible states-at-birth are given by the set \(\left\{\boldsymbol{\phi}_1,\ldots,\boldsymbol{\phi}_m\right\}\), the individual-level impact is an \(m\)-dimensional vector \(\Gamma=\left(\Gamma_1,\ldots,\Gamma_m\right)\). The population-level impact on the environment in this case equals the dot product of this vector \(\Gamma\) with the \(m\)-dimensional vector \(\tilde{b}\), representing the equilibrium distribution of produced offspring over the possible states-at-birth \(\left\{\boldsymbol{\phi}_1,\ldots,\boldsymbol{\phi}_m\right\}\) (refer to Diekmann, Gyllenberg, and Metz (2003) for details). > > The program automatically computes the equilibrium distribution of produced offspring over the possible states-at-birth and uses it to compute the population-level impacts contained in I[][] from the individual-level impacts that have been specified in impact[][].

Code block 9.1.1.1.A

/*
 *===========================================================================
 *      DEFINITION OF THE LIFE HISTORY MODELS FOLLOWS BELOW
 *===========================================================================
 * Specify the number of states at birth for the individuals in all structured
 * populations in the problem in the vector BirthStates[].
 *===========================================================================
 */

void SetBirthStates(int BirthStates[POPULATION_NR], double E[])
{
  BirthStates[0] = 2;

  return;
}

The values of the i-state variables in the different states-at-birth have to be defined separately. This means that in the routine StateAtBirth() the assignment of the values to the variables istate[][] has to be made conditional on the value of the index of the state-at-birth BirthStateNr. The model implemented in the file KlanjscekDEB2.h has 2 states-at-birth \(\left\{\boldsymbol{\phi}_1,\boldsymbol{\phi}_2\right\}\). In the code box below it is shown that the state-at-birth \(\boldsymbol{\phi}_1\) (with index BirthStateNr \(= 0\)) corresponds to the small-sized offspring with body size \(0.7\cdot V_b\), while the state-at-birth \(\boldsymbol{\phi}_2\) (BirthStateNr \(= 1\)) corresponds to the large-sized offspring with body size \(1.3\cdot V_b\).

Code block 9.1.1.1.B

/*
 *===========================================================================
 * Specify all the possible states at birth for all individuals in all
 * structured populations in the problem. BirthStateNr represents the index of
 * the state of birth to be specified. Each state at birth should be a single,
 * constant value for each i-state variable.
 *
 * Notice that the first index of the variable 'istate[][]' refers to the
 * number of the structured population, the second index refers to the
 * number of the individual state variable. The interpretation of the latter
 * is up to the user.
 *===========================================================================
 */

void StateAtBirth(double *istate[POPULATION_NR], int BirthStateNr, double E[])
{
  if (BirthStateNr == 0)
    {
      AGE    = 0.0;
      VOLUME = 0.7*VB;
      Q      = 0.0;
      H      = 0.0;
    }
  else
    {
      AGE    = 0.0;
      VOLUME = 1.3*VB;
      Q      = 0.0;
      H      = 0.0;
    }

  return;
}

Finally, the last part of the code that has to be changed in case of multiple states-at-birth is the assignment of fecundity to different states-at-birth. The model implemented in the file KlanjscekDEB2.h assumes that individuals with a different state-at-birth differ in their offspring production. More specifically, individuals that are born with a small size (\(V=0.7\cdot V_b\)) are assumed to invest 2/3 of the energy that they have available for reproduction on producing small offspring (i.e. with \(V=0.7\cdot V_b\)) and 1/3 of the reproductive energy producing large-sized offspring with \(V=1.3\cdot V_b\). Vice versa, individuals that are born with a large body size (\(V=1.3\cdot V_b\)) are assumed to invest 2/3 of the energy that they have available for reproduction on producing large offspring (i.e. with \(V=1.3\cdot V_b\)) and 1/3 of the reproductive energy producing small-sized offspring with \(V=0.7\cdot V_b\). Parents bias their energetic investment into reproduction therefore toward producing offspring with the same size at birth as they were born with themselves. These different energetic investments into the two types of offspring are subsequently converted into a number of offspring by dividing them by the energy costs to produce a single offspring. For small- and large-sized offspring these costs are proportional to \(0.7\cdot V_b\) and \(1.3\cdot V_b\), respectively. For mothers born with a small size-at-birth this implies that the biased energetic investment in producing offspring with small sizes-at-birth is even more pronounced when considered in terms of number of offspring produced, whereas for mothers born with a large size-at-birth the bias is dampened by the conversion to number of offspring produced.

The size-at-birth of the offspring produced is therefore on average smaller in the model implemented in KlanjscekDEB2.h, compared to the model with a single state-at-birth, which is implemented in KlanjscekDEB.h. As a consequence, the number of offspring produced is larger, which is the most likely reason for the finding that the population growth rate of the model with 2 states-at-birth is consistently larger than in case of a single state-at-birth (see the graphical output of the demo script KlanjscekDEB).

Code block 9.1.1.1.C

/*
 *===========================================================================
 * Specify the fecundity of individuals as a function of the i-state
 * variables and the individual's state at birth for all populations in every
 * life stage.
 *
 * The number of offspring produced has to be specified for every possible
 * state at birth in the variable 'fecundity[][]'. The first index of this
 * variable refers to the number of the structured population, the second
 * index refers to the number of the birth state.
 *
 * Notice that the first index of the variable 'istate[][]' refers to the
 * number of the structured population, the second index refers to the
 * number of the individual state variable. The interpretation of the latter
 * is up to the user.
 *===========================================================================
 */

void Fecundity(int lifestage[POPULATION_NR], double *istate[POPULATION_NR],
               double *birthstate[POPULATION_NR], int BirthStateNr, double E[],
               double *fecundity[POPULATION_NR])
{
  double        Er;

  if (lifestage[0] == 1)                // Only for adults
    {
      Er = (1-KAPPA)*EM*FOOD*G*(NU*pow(VOLUME, 2.0/3.0) + M*VOLUME)/(FOOD + G);
      fecundity[0][0]  = fecundity[0][1]  = max(Er - (1-KAPPA)*EM*M*G*VP,0);
      if (BirthStateNr == 0)
        {
          fecundity[0][0] *= (2.0/3.0)*KAPPA_R/(EM*(KAPPA*G + FOOD)*0.7*VB);
          fecundity[0][1] *= (1.0/3.0)*KAPPA_R/(EM*(KAPPA*G + FOOD)*1.3*VB);
        }
      else
        {
          fecundity[0][0] *= (1.0/3.0)*KAPPA_R/(EM*(KAPPA*G + FOOD)*0.7*VB);
          fecundity[0][1] *= (2.0/3.0)*KAPPA_R/(EM*(KAPPA*G + FOOD)*1.3*VB);
        }
    }
  else
    {
      fecundity[0][0] = 0;
      fecundity[0][1] = 0;
    }

  return;
}

In case of multiple states at birth the structures in the output file containing the stable population states is more complex. Consider for example the computation with the file KlanjscekDEB2.h that is executed when running the demo script KlanjscekDEB:

Command box 9.1.1.1.A

> output2 <- PSPMdemo("KlanjscekDEB2", c(0, 1.0, -0.02, 0.4, 1.0), clean=TRUE, force=TRUE)

Building executable KlanjscekDEB2demo.so ...

<...compilation output lines suppressed in this box...>

  1.00000000E+00, 6.95508116E-01
  9.80000000E-01, 6.82539595E-01
  9.60000000E-01, 6.69183490E-01
<...output lines suppressed in this box...>
  4.40000000E-01, 9.73367353E-02
  4.20000000E-01, 5.99646757E-02
  4.00000000E-01, 2.03431517E-02

> cat(output2$curvedesc)
 #
 # Executing : PSPMdemo("KlanjscekDEB2", c(0, 1, -0.02, 0.4, 1), NULL, NULL)
 #
 # Parameter values  : 
 #
 #  Food      :  1              Kappa     :  0.8            Kappa_R   :  0.001        
 #  Nu        :  0.075          m         :  0.583          g         :  1.286        
 #  Vb        :  1E-09          Vp        :  1.73E-06       [Em]      :  0.7          
 #  ha        :  0.15         
 #
 # Index and name of bifurcation parameter #1                   :  0 (Food)
 #
 #         1:Food        2:PGR[0]         3:Tc[0]       4:S[0][0]       5:S[0][1]       6:S[0][2] ....

Obviously, this model contains more parameters and hence there are many more columns in the output representing sensitivities of the population growth rate with respect to model parameters. Loading the first structure from the output file KlanjscekDEB2-PGR-0000.csb containing the stable population states and displaying its contents reveals the additional elements due to the multiple states at births:

Command box 9.1.1.1.B

> csbread("KlanjscekDEB2-PGR-0000.csb", 1)
$BifPars
[1] 1

$Parameters
 [1] 1.000e+00 8.000e-01 1.000e-03 7.500e-02 5.830e-01 1.286e+00 1.000e-09 1.730e-06 7.000e-01 1.500e-01

$PGR
[1] 0.6955081

$Pop00_StableBirthDist
      Bstate00  Bstate01
[1,] 0.7097291 0.2902709

$Pop00_BirthStates
     Istate00 Istate01 Istate02 Istate03
[1,]        0  7.0e-10        0        0
[2,]        0  1.3e-09        0        0

$Pop00_Bstate00
         StableDist   Istate00     Istate01     Istate02   Istate03  ReproVal
  [1,] 7.097291e-01  0.0000000 7.000000e-10 0.000000e+00 0.00000000  1.069154
  [2,] 6.570771e-01  0.1101390 8.928356e-09 7.625625e-09 0.01058176  1.154825
  [3,] 6.074248e-01  0.2202781 3.423292e-08 3.154797e-08 0.02537620  1.249224
<...output lines suppressed in this box...>
 [98,] 1.551992e-09 10.6834855 3.313724e-04 9.494863e-04 2.86850013  5.715531
 [99,] 1.045388e-09 10.7936245 3.367025e-04 9.735928e-04 2.91605413  3.380026
[100,] 7.097291e-10 10.9001686 3.418510e-04 9.971985e-04 2.96246947  0.000000

$Pop00_Bstate01
         StableDist   Istate00     Istate01     Istate02    Istate03  ReproVal
  [1,] 2.902709e-01  0.0000000 1.300000e-09 0.000000e+00 0.000000000 0.8309155
  [2,] 2.687551e-01  0.1101019 1.177775e-08 9.742355e-09 0.009863215 0.8974364
  [3,] 2.484770e-01  0.2202038 4.088023e-08 3.735042e-08 0.024398778 0.9706759
<...output lines suppressed in this box...>
 [98,] 6.349473e-10 10.6798867 3.321082e-04 9.527934e-04 2.871120854 4.3466391
 [99,] 4.276192e-10 10.7899886 3.374354e-04 9.769332e-04 2.918718535 2.5706185
[100,] 2.902709e-10 10.8964976 3.425811e-04 1.000571e-03 2.965176457 0.0000000

The first additional element of the list representing the equilibrium population state is Pop00_StableBirthDist, which specifies the stable distribution of offspring produced with the 2 possible states at birth that are defined in the model. Each of the rows of the element Pop00_BirthStates of the structure specifies a different state at birth with its columns specifying the value of the 4 individual state variables in that particular state. For each state at birth, a stable population distribution for individuals born in that particular state is stored in two-dimensional arrays, called Pop00_Bstate00 and Pop00_Bstate01 respectively. As before, these two-dimensional arrays contain as the first and last column the stable population density and the reproductive value, respectively, while the intervening columns contain the values of the individual state variables. Notice, however, that the first column containing the stable population density does not start at an initial value of 1.0, as for each of the possible states at birth the stable density is multiplied by the corresponding element of the stable population distribution at birth, contained in the two-dimensional array Pop00_Bstate00.