Chapter 4 PSPMecodyn: ecological dynamics

The ecological dynamics of a PSPM can be simulated using the Escalator Boxcar Train (EBT in short). The EBT is a numerical method especially designed for integrating the partial differential equations that are the mathematical representations of PSPMs (De Roos 1988; De Roos, Diekmann, and Metz 1992). The EBT method differs from standard numerical integration methods for partial differential equations, as its design is inspired by the biological underpinning of PSPMs, rather than by their mathematical expressions in terms of partial differential equations. A recent comparison of 4 different methods for the numerical integration of size-structured population models therefore concluded that the EBT method performs overall best and is the only method that can be straightforwardly extended to PSPMs with more than a single individual state variable (Zhang, Dieckmann, and Brännström 2017).

The EBTtool is a software package that implements the Escalator Boxcar Train and provides a graphical user interface for its operation. This software package has been around already for quite some years. The EBTtool is itself implemented in C and also requires that a particular PSPM also is implemented in C. The current PSPManalysis package includes a function PSPMecodyn that is a trimmed-down version of the EBTtool. This trimmed-down version does not include the graphical user interface and is limited to the type of models that can be implemented using the model specification templates of the PSPManalysis package (see the discussion in sections 3.3.1 and 3.3.2). For example, whereas pulsed reproduction is easily implemented and dealt with in the EBTtool, the function PSPMecodyn can not handle it as reproduction is assumed to be continuous. Nevertheless, PSPMecodyn is useful to study the ecological dynamics of the models that can be analyzed with the PSPManalysis package and works with exactly the same model implementation file that is used for all other computations as well. Notice, however, that the PSPM should be non-linear and hence include environmental state variables as well as feedback of the population on these environmental state variables. The Medfly model that was discussed in chapter 2 in the context of demographic analysis of PSPMs can hence not be studied using PSPMecodyn.

Even though the EBTtool can not be used to numerically integrate a model that is implemented in the PSPManalysis template, it might nonetheless be a good idea to download and install this software as it provides a useful graphical user interface to explore the results of numerical integrations. In particular, it uses the same type of data files as produced by the PSPManalysis package (output files with extensions .out and .csb; see section 4.3 below) and is particularly useful to analyze the changes in the population distribution over time. The EBTtool moreover implements methods to output sequences of these population states to an animated GIF file that can be used for presentations.

The PSPManalysis package implements both the simplified EBT method proposed by Brännström, Carlsson and Simpson (2013) as well as the original EBT method (De Roos 1988; De Roos, Diekmann, and Metz 1992). The simplified EBT method is faster, but possibly less accurate than the original EBT method. Switching between the two methods can be done by defining the EBTMETHOD equal to \(0\) or 1 in the model implementation file. See also chapter 8.

4.1 Arguments and output of the PSPMecodyn function

The use of the PSPMecodyn function will be illustrated with the PNAS model as described in section 3.2 and analyzed in section 3.4. A demo script "deRoosPerssonDynamics" is included with the PSPManalysis package that carries out numerical integration of the PNAS model, starting from three different initial conditions.

The general call to the PSPMecodyn function is shown in the command box below.

> output <- PSPMecodyn(modelname = NULL, startstate = NULL, timepars = NULL, bifpars = NULL,
                       parameters = NULL, options = NULL, clean = FALSE, force = FALSE,
                       debug = FALSE, silent = FALSE)

The obligatory and optional arguments to the PSPMecodyn function are the following:

  1. The first, obligatory argument to the function PSPMecodyn is the name of the file specifying the PSPM, passed as a string argument. It is unnecessary to include the extension '.R' or '.h' as part of the file name, the PSPMecodyn function will automatically try to locate the appropriate file, checking first for a file implemented in C (with an extension '.h') and subsequently for a file implemented in R (with an extension '.R'). If both a file with an extension '.h' and a file with an extension '.R' are found, the program will use the first one. The program can be forced to use the file with an extension '.R' by including the extension explicitly as part of the filename. The R-command to simulate the dynamics of the PNAS model that will be used as illustrations below will therefore all take "PNAS2002" as their first argument. If the file specifying the PSPM can not be found in the current directory, the PSPMequi function will ask the user to search in the package directory for a model file with the specified name.

  2. The second, obligatory argument is the initial condition of the computation. This initial point should be a list that has the same structure as the list returned by the function csbread (see, for example, section 3.4.5). The list should at least include an element $Environment and elements $Pop00, $Pop01, $Pop02, etc., for as many structured populations are accounted for in the model. The element $Environment should be a vector with initial values of the environmental state variables for the simulation, while the elements $Pop00, $Pop01, $Pop02, etc., should be matrices with a number of columns that is exactly 1 larger than the number of individual state variables in the model. Each row of these matrices should specify the initial state of a cohort of individuals in the particular population. The first column of each of these matrices, i.e. $Pop00[i,1], should specify the number of individuals in cohort \(i\), while the average value of the individual state variable with index \(0\) and 1 (average age and average length in the PNAS model) have to be specified in the second and third column of each matrix, $Pop00[i,2] and $Pop00[i,3], respectively. Naturally, if individuals are characterized by more than two individual state variables, the values of these have to follow in additional columns. In general, there are 3 different methods to construct an initial state for the function PSPMecodyn: (1) the initial state can be a population state that results from an equilibrium computation with PSPMequi; (2) the initial state can be generated using the function PSPMind that computes the individual life history under a specific set of environmental conditions (see chapter 7); and, (3) the initial state can be constructed by the user with list-construction commands in R. These 3 different methods will be discussed in more detail in section 4.2.

    If the list in the second argument contains an element $Parameters and this vector is of the appropriate length, this vector will be used as values for the model parameter, for which to carry out the integration. This vector should have the same length as the number of parameters in the model (the length of the vector DefaultParameters in R, or the value of PARAMETER_NR in C). When of this length the values will replace the default values of the parameters that are listed in the model specification file. However, if a vector of parameter values is included as fifth argument in the call to the function PSPMecodyn, the latter supersede any values that may be specified in the initial state list.

    In principle, the list that is given as the second argument to the function PSPMecodyn should also include elements called $Pop00_BirthStates, $Pop01_BirthStates, $Pop02_BirthStates, etc., that specify the state at birth of the individuals in the different populations as the state at birth may influence the life history dynamics. These elements are automatically included when the initial state is produced by either a call to the function PSPMequi or to the function PSPMind. They consist of a vector (in case of a single state at birth) or a two-dimensional matrix of values (in case of multiple states at birth) with the values of the individual state variables with which individuals can be born, each row representing a state at birth and each column an individual state variable. If their are multiple states at birth the initial state list should contain elements $Pop00_Bstate00, $Pop00_Bstate01, $Pop00_Bstate02, etc., specifying the subpopulations originating from each of the states at birth (see command box 9.1.1.1.B in section 9.1). If these elements are missing, however, the function PSPMecodyn will automatically generate valid individual states at birth through calls to the functions SetBirthStates() and StateAtBirth() (see sections 3.3.2.4 and 3.3.2.5) or their appropriate counterparts when the model is implemented in R.

  3. The third, obligatory argument to the PSPMecodyn function is a row vector consisting of 4 elements:

    c(\(<\)cohort limit\(>\),\(<\)output interval\(>\),\(<\)state output interval\(>\),\(<\)maximum time\(>\))

    In this vector \(<\)cohort limit\(>\) refers to the time interval during which a new cohort of individuals is formed. The EBT method is based on the idea to collect all individuals that are born within a interval \(\Delta t\) into a single cohort of the population (De Roos 1988; De Roos, Diekmann, and Metz 1992). The element \(<\)cohort limit\(>\) sets the value of this time interval \(\Delta t\). Smaller values of \(\Delta t\) will mean that the computed model trajectories are more accurate, but will also slow down the computation as the populations will include more cohorts and hence the system of equations to integrate is larger. The vector element \(<\)output interval\(>\) defines the time interval between data output to the output file with extension .out (see section 4.3 below), while the vector element \(<\)state output interval\(>\) similarly defines the time interval between output of the entire population state to the output file with extension .csb (see section 4.3 below). Finally, the vector element \(<\)maximum time\(>\) sets the maximum time at which to halt the numerical integration.

  4. The fourth, optional argument to the PSPMecodyn function can be used to carry out bifurcation analysis of the model dynamics, in case the model predicts non-equilibrium dynamics. An easy way to create a bifurcation graph is to carry out a large number of numerical integrations of the model equations, each with a slightly different value for the bifurcation parameter. This approach, however, has its problems, as the computed outcome will depend on the choice of the initial values for the model variables. One way to circumvent the problems associated with the initial value of model variables is to carry out only a single numerical integration of the model equations but over a very long time period, in which the entire integration period is subdivided into intervals during which the value of the bifurcation parameter is constant, while from one interval to the next the value of the bifurcation parameter is increased or decreased by a small amount. In this way, the range of values of the bifurcation parameter is scanned either from low to high or vice versa. This stepwise increase or decrease of the bifurcation parameter implies that the final values of the model variables obtained for a particular value of the bifurcation parameter are used as initial values of the model variables for the subsequent parameter value.

    The advantage of this approach can best be explained in the context of stable model equilibria, but it works equally well in case of non-equilibrium dynamics. Consider that after a numerical integration over a sufficiently long time interval for a particular value of the bifurcation parameter all model variables have ended up close their equilibrium value. These final values of the model variables will also be close to their equilibrium value for a setting of the bifurcation parameter that is slightly larger or smaller, provided that the particular equilibrium still exists for this new parameter value. Hence, adopting these final values of the model variables as initial values for a subsequent integration with a slightly different bifurcation parameter value ensures that we continue to follow the curve of a particular model equilibrium as a function of the bifurcation parameter as long as the equilibrium exists and is stable. Only when the equilibrium becomes unstable or does not occur at all any more for the new value of the bifurcation parameter, will the model variables approach an entirely different equilibrium or a different type of dynamics, such as a limit cycle. By scanning a particular interval of the bifurcation parameter with increasing as well as decreasing parameter values in most cases also reveals the co-occurrence of alternative stable equilibria or alternative types of stable dynamics, such as the co-occurrence of different types of limit cycles, for the same value of the bifurcation parameter.

    If specified, this argument has to be a vector with 6 elements:

    c(\(<\)index\(>\),\(<\)start\(>\),\(<\)step\(>\),\(<\)stop\(>\),\(<\)output period\(>\),\(<\)state output period\(>\))

    In this vector \(<\)index\(>\) indicates the index of the bifurcation parameter to vary, \(<\)start\(>\) indicates the starting value of the bifurcation parameter, \(<\)step\(>\) the step size in the bifurcation parameter, \(<\)stop\(>\) the final value of the bifurcation parameter, \(<\)output period\(>\) the time interval at the end of each bifurcation interval during which data output is written to the file with .out (see section 4.3 below), and \(<\)state output period\(>\) determines the time interval at the end of each bifurcation interval during which output of the entire population state is written to the file with .csb (see section 4.3 below).

  5. The fifth, optional argument of the PSPMecodyn function is a (row) vector of model parameter values. When used, this array should have the same length as the number of parameters in the model (the length of the vector DefaultParameters in R, or the value of PARAMETER_NR in C). When of this length the values will replace the default values of the parameters that are listed in the model specification file. If the array used for this sixth argument is not of the correct length or when it is not specified at all, it will simply be ignored.

  6. The sixth, optional argument of the PSPMecodyn function is a (row) vector of string elements, containing possible options that modify the behavior of the computational module. Only two options are possible and both require a value and hence occur as a pair of option name and option value. The option argument is therefore either of the form:

    c("name 1", "value 1")

    or of the form:

    c("name 1", "value 1","name 2", "value 2")

    The options can be specified in any order, but the option value should always immediately follow after the option name. Possible options are:

    1. Option pair c("info", "i"): This option modifies the output of the DOPRI5 integration method to the output file with extension .err (see section 4.3 below). The value of i can be set equal to 1, 2, 3 or 4. The default behavior of the function PSPMecodyn is equivalent to i equal to \(0\), in which case no output at all is produced in the file with extension .err. With higher values of i more information is produced detailing the performance of the integration method (i.e. step sizes used, number of successful and unsuccessful integration steps, etc.). This output is only useful in case problems occur during the numerical integration of the model.

    2. Option pair c("report", "i"): This option, which should have a positive value i, determines the time interval between consecutive lines of data output to the console. The default behavior of the function PSPMecodyn is to write a line of output to the console whenever a new cohort of individuals is constructed. However, this may result in a lot of output. By choosing a higher value for this option the program can be forced to write to the console only every 10th time that a cohort is constructed. Notice that this option does not affect the frequency with which output is written to the output file with extension .out (see section 4.3 below).

Four other optional arguments can be passed to the PSPMecodyn function: clean, force, debug and silent. These are all boolean arguments that hence have to be passed to the PSPMecodyn function as \(<\)option name\(>\)=TRUE or \(<\)option name\(>\)=FALSE, the latter being the default value of all options (Specifying these options as argument is hence only useful when setting them equal to TRUE). Unlike the previous arguments, which all modify the computations to be performed, these options modify the behavior of the PSPMecodyn function itself, in particular the compilation of the model specific file into a dynamic library module that can be executed from R. Also unlike all the previous arguments that can be passed, these arguments can be passed in any order and at any position, the PSPMecodyn function will filter these 3 optional arguments from the argument list before passing the filtered argument list to the computational routine.

  • Option clean: When clean=TRUE is passed as argument, this argument instructs the PSPMecodyn function to delete all result files that have been generated during previous calculations with the model. These result files have names of the form <Modelname>-<Type>-<NNNN>.err, <Modelname>-<Type>-<NNNN>.csb and <Modelname>-<Type>-<NNNN>.out, in which <Modelname> refers to the name of the model, <Type> refers to the type of computation that has been performed, which in the case of PSPMecodyn equals ECODYN, and <NNNN> is a unique number that distinguishes consecutive computations of the same type of curve with the same model. Deleting all the output files from previous computations and/or the compiled program executables that the package has generated can also be done separately. The package implements a function PSPMclean() to delete all .bif, .err, .csb and .out files and/or all executable files that are present in the current working directory.

  • Option force: When force=TRUE is passed as argument, it instructs the PSPMecodyn function to force re-compilation of the model specific file into a dynamic library module that can be executed by R. This option will usually not be needed by normal users, as the PSPMecodyn function automatically recompiles the computational module when the model specific file with an '.h' or '.h' extension is more recently changed than the compiled dynamic library file. However, if for some unclear reason this automatic recompilation fails, the force option can be used to initiate re-compilation.

  • Option debug: When debug=TRUE is passed as argument, it instructs the PSPMecodyn function to turn on debugging flags while compiling the model specific file into a dynamic library module. This option can be useful to detect programming mistakes in the model-specific file that are otherwise hard to track down. The downside is that depending on the version of R that is used, turning on debugging flags during compilation may generate a lot of output, including warnings about standard files of the operating system that are perfectly correct. It is hence not so easy to spot among all these messages the warnings that relate to the model-specific code that has been implemented.

  • Option silent: When silent=TRUE is passed as argument, it instructs the PSPMecodyn function to suppress all messages from the compilation of the model specific file into a dynamic library module. This option is useful to prevent cluttering the console with superfluous messages once a model specific file has been tested sufficiently and functions without problems.

The computational module generates on execution a single list object as output with 2 member elements (see the help page on PSPMecodyn using ?PSPMecodyn). The first element of the output list, output$curvepoints contains the numerical information of the points along the computed curve. This variable output$curvepoints is a matrix, in which each row represents one solution point along the curve. The columns contain the time value, the value of all environment variables, the value for the birth rate of all structured populations in the problem, followed by the value of all interaction variables defined in the routine Impact().

The second member element of the output list, output$curvedesc, contains the description of the executed calculation, which includes the command-line that is used for the invocation of the computational routine, the values of all parameters used for the current computation and a header line indicating the meaning of all the output variables produced by the computational module. This textual information is also printed to the R console at the end of calculations. In fact, the PSPMecodyn function prints its report on the calculations by execution of the statement cat(output$curvedesc, sep='\n').

4.2 An example session using the PSPMecodyn function

The demo script "deRoosPerssonDynamics" illustrates the use of the PSPMecodyn function by simulating the ecological dynamics of the PNAS model as presented in section 3.2, starting from 3 different initial conditions. The first numerical integration starts from an equilibrium state computed with the function PSPMequi. The particular calls to the PSPMequi and PSPMecodyn functions are shown in the following R command box:

Command box 4.2.A
> output <- PSPMequi("PNAS2002", "EQ", c(3.0E-04, 1.561E-04, 1.270E-04, 4.008E-06, 2.761E-04), 0.1, 
                     c(1, 0, 1), options = c("single"), clean = TRUE, force = TRUE)

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

  3.00000000E-04, 1.56127570E-04, 1.27032702E-04, 4.00801603E-06, 2.76129173E-04

> initstate <- csbread("PNAS2002-EQ-0000.csb", 1)
> output1 <- PSPMecodyn("PNAS2002", initstate, c(1, 1, 10, 1000), options = c("report", "50"),
                        clean = TRUE, force = TRUE)

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

        0.00, 1.56127570E-04, 1.27032702E-04, 3.26244320E-06
       50.00, 1.57798334E-04, 1.26716011E-04, 3.96498271E-06
      100.00, 1.59593562E-04, 1.26024893E-04, 3.92329007E-06
<...output lines suppressed in this box...>
      900.00, 1.56843463E-04, 1.27713985E-04, 3.95446742E-06
      950.00, 1.58701063E-04, 1.26750658E-04, 3.93952452E-06
     1000.00, 1.59249607E-04, 1.25798353E-04, 3.95558086E-06


RUN PNAS2002-ECODYN-0000 COMPLETED at T = 1000.00:
** Program terminated. Normal closure of output files succeeded.          **

> cat(output1$curvedesc)
 #
 # Executing : PSPMecodyn("PNAS2002", <Initial state>, c(1, 1, 10, 1000), NULL, NULL, c("report", "50"))
 #
 # Parameter values  : 
 #
 #  Rho       :  0.1            Rmax      :  0.0003         Lb        :  7            
 #  Lv        :  27             Lj        :  110            Lm        :  300          
 #  Omega     :  9E-06          Imax      :  0.0001         Rh        :  1.5E-05      
 #  Gamma     :  0.006          Rm        :  0.003          Mub       :  0.01         
 #  A         :  5000           Th        :  0.1            Epsilon   :  0.5          
 #  Delta     :  0.01         
 #
 #  Cohort cycle time interval                                         : 1.0
 #  Output time interval                                               : 1.0
 #  Complete state output interval                                     : 10.0
 #  Maximum integration time                                           : 1000.0
 #
 #
 #   1:Time            2:E[0]          3:E[1]          4:E[2]          5:b[0]       6:I[0][0] ....

> output1$curvepoints
        Time         E[0]         E[1]         E[2]         b[0] ..      I[0][1]      I[0][2]      I[0][3]
   [1,]    0 0.0001561276 0.0001270327 3.262443e-06 0.0000000000 .. 3.262443e-06 1.414324e-05 0.0001710999
   [2,]    1 0.0001562952 0.0001268490 3.571646e-06 0.0002817937 .. 3.571646e-06 9.659792e-06 0.0001755738
   [3,]    2 0.0001563517 0.0001267466 3.775143e-06 0.0002807738 .. 3.775143e-06 1.000305e-05 0.0001752195
<...output lines suppressed in this box...>
 [999,]  998 0.0001592588 0.0001258301 3.954380e-06 0.0002686477 .. 3.954380e-06 1.157543e-05 0.0001699083
[1000,]  999 0.0001592545 0.0001258141 3.954975e-06 0.0002686223 .. 3.954975e-06 1.158832e-05 0.0001698792
[1001,] 1000 0.0001592496 0.0001257984 3.955581e-06 0.0002685985 .. 3.955581e-06 1.160123e-05 0.0001698510

(In the output shown in the box above a large number of output lines have been suppressed and some intermediate columns have been deleted, because the page width does not allow them to be shown completely.)

The first R command in the box above calls the function PSPMequi with the option argument c("single"). This instructs the PSPMequi function to compute the equilibrium state for only a single parameter value (in this case the value \(3.0\cdot10^{-4}\) for the parameter with index 1, which is the parameter \(R_{max}\) in the PNAS model; see code block 3.3.2.2). This call to the function PSPMequi produces a result file called PNAS2002-EQ-0000.csb containing a single set of data for the environmental state variables and the population, which is imported into the R workspace using the function csbread. This initial state, contained in the variable initstate, is subsequently used as initial condition for the numerical integration of the dynamics using the function PSPMecodyn.

This starting point of the numerical integration is passed to the function PSPMecodyn as its second argument, as shown above, whereas the first argument defines the basename of the file with the model implementation "PNAS2002". The third argument to the function PSPMecodyn sets the interval during which a new cohort is formed equal to 1.0, the time interval between output written to the file with extension .out also equal to 1.0 and the time interval between complete state output written to the file with extension .csb equal to 10.0. As last element of the third argument the maximum integration is set equal to 1000.0. The final argument to the function PSPMecodyn, the option argument c("report", "50"), forces the function PSPMecodyn to only write output to the console every 50 time units. Every 50 time units the function PSPMecodyn hence reports the current values of the time and the environmental state variables to the console.

In the R command box above the values of the computed points along the trajectory of the model dynamics are saved in the output list element output1$curvepoints, the contents of which are inspected after the PSPMecodyn function finishes and it has printed out the textual information about the computation. The output list element output1$curvepoints contains in consecutive columns the time value, the values of all environmental state variables, the current population birth rate of all structured populations and the values of all interaction variables defined in the function Impact() (see section 3.3.2.11). The demo script "deRoosPerssonDynamics" uses the data contained in output1$curvepoints to plot the juvenile biomass (the sum of column 7 and 8 of output1$curvepoints) as well as the adult biomass (column 9 of output1$curvepoints) as a function of time.

In the following R command box, another trajectory of the dynamics is computed, but now starting from an initial state generated by the function PSPMind. This function computes the individual life history at a specific set of values for the environmental state variables. The details about the function PSPMind are discussed in chapter 7 and hence will not be covered here. The only thing to point out here is that the second argument in the call to the function PSPMind is a vector c(1.561276e-04, 1.270327e-04, 4.008016e-06, 0.01) with as its first 3 elements the values of the environmental state variables at which to compute the individual life history (see section 3.3.2.1 for the interpretation of these environmental state variables). The last element of this vector (0.01) defines the population birth rate, which value will be used to scale the number of individuals in all cohorts of the population. Larger values of this birth rate will hence imply that the size-structured consumer population is initially larger. The output of the function PSPMind is a list with the same structure as produced by the function PSPMequi in its state output file with extension .csb. This state produced as output by the function PSPMind is assigned to the variable initstate, which is subsequently used as starting point for the numerical integration. The subsequent call to the function PSPMecodyn is identical to the one shown in R command box 4.2.A, while also its output is similar to the output shown in command box 4.2.A and discussed above. These will hence not be further discussed here. The demo script "deRoosPerssonDynamics" uses also the data contained in output2$curvepoints to plot the juvenile biomass (the sum of column 7 and 8 of output2$curvepoints) as well as the adult biomass (column 9 of output2$curvepoints) of this trajectory as a function of time.

Command box 4.2.B
> initstate <- PSPMind("PNAS2002", c(1.561276e-04, 1.270327e-04, 4.008016e-06, 0.01),
                       options = c("isort", "1"))

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

                                Istate[ 0]  Istate[ 1]    Survival        R0    Impact[ 0] ....
Pop. # 0 - Bstate  0 - (Final):    1248.79    273.555       1E-09          1     0.0521033 ....


> output2 <- PSPMecodyn("PNAS2002", initstate, c(1, 1, 10, 1000), options = c("report", "50"))
Dynamic library file PNAS2002ecodyn.so is up-to-date

        0.00, 1.56127600E-04, 1.27032700E-04, 1.00801215E-04
       50.00, 1.32044964E-05, 1.68821129E-04, 2.82764765E-06
      100.00, 2.10361075E-05, 1.52692631E-04, 3.64796500E-06
<...output lines suppressed in this box...>
      900.00, 1.90701042E-04, 1.25787901E-04, 3.13606041E-06
      950.00, 1.96875492E-04, 1.13526555E-04, 3.29816659E-06
     1000.00, 1.81028462E-04, 1.07370790E-04, 3.90578812E-06


RUN PNAS2002-ECODYN-0001 COMPLETED at T = 1000.00:
** Program terminated. Normal closure of output files succeeded.          **

> cat(output2$curvedesc)
 #
 # Executing : PSPMecodyn("PNAS2002", <Initial state>, c(1, 1, 10, 1000), NULL, NULL, c("report", "50"))
 #
 # Parameter values  : 
 #
 #  Rho       :  0.1            Rmax      :  0.0003         Lb        :  7            
 #  Lv        :  27             Lj        :  110            Lm        :  300          
 #  Omega     :  9E-06          Imax      :  0.0001         Rh        :  1.5E-05      
 #  Gamma     :  0.006          Rm        :  0.003          Mub       :  0.01         
 #  A         :  5000           Th        :  0.1            Epsilon   :  0.5          
 #  Delta     :  0.01         
 #
 #  Cohort cycle time interval                                         : 1.0
 #  Output time interval                                               : 1.0
 #  Complete state output interval                                     : 10.0
 #  Maximum integration time                                           : 1000.0
 #
 #
 #   1:Time            2:E[0]          3:E[1]          4:E[2]          5:b[0]       6:I[0][0] ....
 
> output2$curvepoints
        Time         E[0]         E[1]         E[2]         b[0] ..      I[0][1]      I[0][2]      I[0][3]
   [1,]    0 1.561276e-04 0.0001270327 1.008012e-04 0.0000000000 .. 1.008012e-04 1.769969e-04 0.0010387677
   [2,]    1 2.133163e-05 0.0001569627 8.017346e-05 0.0011930961 .. 8.017346e-05 1.817787e-04 0.0010364129
   [3,]    2 4.103512e-06 0.0001795904 4.068710e-05 0.0004298024 .. 4.068710e-05 1.806636e-04 0.0010172865
<...output lines suppressed in this box...>
 [999,]  998 0.0001821318 0.0001074357 3.873354e-06 0.0001903570 .. 3.873354e-06 1.988401e-05 0.0001212593
[1000,]  999 0.0001815848 0.0001074010 3.889513e-06 0.0001911060 .. 3.889513e-06 2.009461e-05 0.0001214775
[1001,] 1000 0.0001810285 0.0001073708 3.905788e-06 0.0001918771 .. 3.905788e-06 2.030538e-05 0.0001217083

(In the output shown in the box above a large number of output lines have been suppressed and some intermediate columns have been deleted, because the page width does not allow them to be shown completely.)

In the following R command box, another trajectory of the dynamics is computed, but now starting from an initial state that is constructed using list construction commands from the R command line. The first command in the box below constructs a list with elements Environment and Pop00. The first element Environment contains the initial values for the 3 environmental state variables in the PNAS model (see section 3.3.2.1 for the interpretation of these environmental state variables). The last element Pop00 is a matrix with cohort data for the initial population to start the numerical integration with. This matrix should have 3 columns, specifying the number of individuals in the cohort, their age and their body length, given that in the PNAS model age and length are the two individual state variables. If more columns are specified in the element Pop00, they will be ignored. Each row of the matrix contained in Pop00 specifies one cohort. The first command below shows that the initial population for the following integration consists of a cohort of newborn individuals with age \(0\) and length \(\ell=\ell_b=7.0\). The density of these newborn individuals is \(0.001\). In addition, the initial population includes a cohort of adult individuals with age 300 and length \(\ell=111\) with a cohort density equal to \(1.0\cdot10^{-5}\). The list initstate that is thus produced in the first command in the box below is subsequently used as starting point for the numerical integration. The subsequent call to the function PSPMecodyn is identical to the one shown in R command box 4.2.A, while also its output is similar to the output shown in command box 4.2.A and discussed above. These will hence not be further discussed here. The demo script "deRoosPerssonDynamics" uses also the data contained in output3$curvepoints to plot the juvenile biomass (the sum of column 7 and 8 of output3$curvepoints) as well as the adult biomass (column 9 of output3$curvepoints) of this trajectory as a function of time.

Command box 4.2.C
> initstate <- list(Environment = c(1.561276e-04, 1.270327e-04, 4.008016e-06),
                    Pop00 = matrix(c(0.001, 0, 7.0, 1.0E-5, 300, 111), ncol = 3, byrow = TRUE))
> output3 <- PSPMecodyn("PNAS2002", initstate, c(1, 1, 10, 1000), options = c("report", "50"))
Dynamic library file PNAS2002ecodyn.so is up-to-date

        0.00, 1.56127600E-04, 1.27032700E-04, 3.08700000E-06
       50.00, 1.27379538E-04, 1.39392354E-04, 4.48274388E-06
      100.00, 1.34147338E-04, 1.42360458E-04, 3.99924658E-06
<...output lines suppressed in this box...>
      900.00, 1.29856560E-04, 1.27516287E-04, 4.56714543E-06
      950.00, 1.27055262E-04, 1.35483257E-04, 4.37381850E-06
     1000.00, 1.36470065E-04, 1.39125331E-04, 4.05246328E-06


RUN PNAS2002-ECODYN-0002 COMPLETED at T = 1000.00:
** Program terminated. Normal closure of output files succeeded.          **

Warning messages:
1: In PSPMecodyn("PNAS2002", initstate, c(1, 1, 10, 1000), options = c("report",  :
  
Initial state list does not contain an element "Parameters" with parameter values
Parameter values taken from model file

2: In PSPMecodyn("PNAS2002", initstate, c(1, 1, 10, 1000), options = c("report",  :
  
Birth states of population Pop00 not specified, will be obtained through calls to SetBirthStates() and StateAtBirth()

> cat(output3$curvedesc)
 #
 # Executing : PSPMecodyn("PNAS2002", <Initial state>, c(1, 1, 10, 1000), NULL, NULL, c("report", "50"))
 #
 # Parameter values  : 
 #
 #  Rho       :  0.1            Rmax      :  0.0003         Lb        :  7            
 #  Lv        :  27             Lj        :  110            Lm        :  300          
 #  Omega     :  9E-06          Imax      :  0.0001         Rh        :  1.5E-05      
 #  Gamma     :  0.006          Rm        :  0.003          Mub       :  0.01         
 #  A         :  5000           Th        :  0.1            Epsilon   :  0.5          
 #  Delta     :  0.01         
 #
 #  Cohort cycle time interval                                         : 1.0
 #  Output time interval                                               : 1.0
 #  Complete state output interval                                     : 10.0
 #  Maximum integration time                                           : 1000.0
 #
 #
 #   1:Time            2:E[0]          3:E[1]          4:E[2]          5:b[0]       6:I[0][0] ....
> output3$curvepoints
        Time         E[0]         E[1]         E[2]         b[0] ..      I[0][1]      I[0][2]      I[0][3]
   [1,]    0 0.0001561276 0.0001270327 3.087000e-06 0.0000000000 .. 3.087000e-06            0 0.0001230868
   [2,]    1 0.0001545920 0.0001269028 4.032747e-06 0.0003394561 .. 4.032747e-06            0 0.0001250936
   [3,]    2 0.0001529027 0.0001270109 4.607412e-06 0.0003415640 .. 4.607412e-06            0 0.0001270795
<...output lines suppressed in this box...>
 [999,]  998 0.0001359284 0.0001390871 4.065588e-06 0.0003399111 .. 4.065588e-06 7.546175e-06 0.0002128082
[1000,]  999 0.0001361979 0.0001391073 4.059020e-06 0.0003393768 .. 4.059020e-06 7.490502e-06 0.0002126396
[1001,] 1000 0.0001364701 0.0001391253 4.052463e-06 0.0003388324 .. 4.052463e-06 7.436078e-06 0.0002124645

(In the output shown in the box above a large number of output lines have been suppressed and some intermediate columns have been deleted, because the page width does not allow them to be shown completely.)

4.3 Output files generated by the PSPMecodyn function

The computational module that is produced by the PSPMecodyn function generates 3 output files. The name of these files is always of the form <Modelname>-ECODYN-<NNNN>.<ext>, in which <Modelname> is the same as the name of the file specifying the model excluding its '.R' or '.h' extension, <NNNN> is a 4-digit number that is unique for the current computation and .<ext> is the extension, which can be either .csb, .err or .out. The unique number distinguishes the same types of curve computations for the same model from each other. The number is obtained by considering increasing values of <NNNN> (i.e., 0000, 0001, 0002 and so forth) and testing whether result files with the particular index are already present. The program uses the first value of <NNNN> that is not in use.

The file called <Modelname>-ECODYN-<NNNN>.err that is generated during the computations contains information about the numerical progress of the computations. It reports details on the steps taken during the numerical integration, such as step sizes used, number of successful and failed integration steps and information about the detection of stage transitions. The amount of detail is dependent on the value of the option "report". The default behavior is to produce no output at all (c("report", "0"); see section 4.1). This file can be informative in case the computation of a particular curve stops for unknown reasons, but is otherwise of little use.

The output file called <Modelname>-ECODYN-<NNNN>.out holds the same information as is contained in the two elements of the output list returned by the PSPMecodyn function, output$curvepoints and output$curvedesc (see the help page on PSPMecodyn using ?PSPMecodyn). The first lines of this file all start with a '#' sign and contain the information about the run performed, which is also contained in output$curvedesc and can be listed by the statement cat(output$curvedesc, sep='\n'). Following this descriptive header the file contains columns with computational results that are also contained in the variable output$curvepoints (see, for example, R command box 4.2.A). In fact, the two elements of the output list, output$curvepoints and output$curvedesc, are generated by reading the contents of the file <Modelname>-ECODYN-<NNNN>.out from disk after the computations have ended, storing all lines that start with a '#' sign into a single string variable output$curvedesc, while storing the information on all other lines into the data matrix output$curvepoints.

The file called <Modelname>-ECODYN-<NNNN>.csb is an output file containing the values of all environment variables and the distribution of all structured populations in the model. This is a binary file, the content of which can be accessed from R using the function csbread. Output is written to this output file at regular time intervals, where the interval between consecutive output times is specified by the third element of the obligatory argument timepars to the function PSPMecodyn (see section 4.1). For example, the file PNAS2002-ECODYN-0002.csb is generated by the invocation of the PSPMecodyn function in R command box 4.2.C. Its contents can be listed by:

Command box 4.3.A
> csbread("PNAS2002-ECODYN-0002.csb")

States in file PNAS2002-ECODYN-0002.csb:

    1: State-0.000000E+00
    2: State-1.000000E+01
    3: State-2.000000E+01
<...output lines suppressed in this box...>
   99: State-9.800000E+02
  100: State-9.900000E+02
  101: State-1.000000E+03

The structure called State-2.000000E+01 contains the population state at time point \(t=20.0\) during the simulation of the ecological dynamics, as its name suggests. Its contents can be read into the workspace by issuing the command csbread("PNAS2002-ECODYN-0002.csb",3) or csbread("PNAS2002-ECODYN-0002.csb","State-2.000000E+01").

Loading this state into the R workspace reveals it to be a list containing various arrays of numbers, as shown in the following box:

Command box 4.3.B
> popstate <- csbread("PNAS2002-ECODYN-0002.csb","State-2.000000E+01")
> popstate
$Time
[1] 20

$Parameters
 [1] 1.0e-01 3.0e-04 7.0e+00 2.7e+01 1.1e+02 3.0e+02 9.0e-06 1.0e-04 1.5e-05 6.0e-03 3.0e-03 1.0e-02 5.0e+03 1.0e-01 5.0e-01 1.0e-02

$Environment
[1] 1.364491e-04 1.324733e-04 4.008016e-06

$Pop00
           Density    Istate00   Istate01     Istate02 Istate03 Istate04 Istate05 Istate06 Istate07
 [1,] 2.687298e-04   0.4443476   7.700723 0.0002687298        0        7        0        0       19
 [2,] 1.371993e-04   1.4444489   9.272342 0.0002681358        0        7        0        0       18
 [3,] 7.013400e-05   2.4445505  10.835272 0.0002674988        0        7        0        0       17
<...output lines suppressed in this box...>
[20,] 6.246577e-08  19.4462036  36.144616 0.0002496189        0        7        1        0        0
[21,] 1.892823e-07  20.0000000  36.932027 0.0010000000        0        7        1        0        0
[22,] 8.187308e-06 320.0000000 129.171752 0.0010000000        0        7        2        0        0

The first element $Time in the list popstate contains the value of the integration time \(t\) at which the population state is stored. The second element of the list, a vector called $Parameters, contains the values of all the model parameters used in the numerical integration. The third element of the list is a vector called $Environment containing the values of the environmental state variables at time \(t\). The last element in the list called $Pop00 is a two-dimensional array characterizing the population distribution at time \(t\) with the first column $Pop00[,1] representing the number of individuals in a particular cohort in the population and the subsequent columns $Pop00[,2] and $Pop00[,3] representing the average values of their individual state variables with index \(0\) and 1 (corresponding to individual age and body size in the PNAS model, see section 3.2), as shown in the R command box above. If individuals would be characterized by more than two individual state variables, the values of these would follow in additional columns of the two-dimensional array $Pop00. The additional columns with labels Istate02 to Istate07 are used by the PSPMecodyn function internally for bookkeeping purposes. The first column labeled Istate02 represents the initial number of individuals in the cohort at the moment it was formed. The following two columns, labeled Istate03 and Istate04, specify the state at birth of the individuals in the cohort, which in case of the PNAS model equals age \(0\) and length \(\ell=\ell_b=7.0\), respectively. In case the number of individual state variables is larger, the number of these columns representing the state at birth will increase accordingly. The next column, labeled Istate05, indicates the life stage that the individuals in the cohort are currently in. In the PNAS model life stage \(0\), 1 and 2 refers to the juveniles that are vulnerable to predation, juveniles that are invulnerable to predation and adult individuals, respectively. The next column, labeled Istate06, contains the index of the state at birth, with which individuals are born. Finally, the last column, labeled Istate07, contains the time value at which the cohort was formed.

The R command box above also illustrates that the dimension of the array $Pop00 indicates that the population at time \(t\) consists of 22 cohorts. However, this number of cohorts will vary over time and is here still rather low because the initial population consisted of only two cohorts of individuals (see command box 4.2.C).

4.4 Additional remarks

The PSPMecodyn function allows for studying the dynamics of models can that can be analyzed with the other functions of the PSPManalysis package as well. In these models reproduction should be continuous in time. The package does allow to carry out demographic analysis for models in which reproduction is pulsed in time (see section 9.1.1.2), but pulsed reproduction in combination with feedback of the environment on individual life history can not be dealt with by PSPManalysis. Moreover, the computations of the equilibrium state of a model, as carried out by the PSPMequi function, rely heavily on the assumption that the environmental variables are constant, which would exclude that the environmental variables exhibit seasonality or other time-dependent dynamics. However, the PSPMecodyn function can be used to study time-dependent dynamics. To allow for this a global variable called Time is defined that holds the current time value. This global variable can for example be used in the function specifying the dynamics of the environmental variables (EnvEqui(); see section 3.3.2.12). Notice that this global variable Time is always equal to \(0\) when using all other functions in the PSPManalysis package.