simecol: An Object-Oriented Framework for Ecological Modeling in R

The simecol package provides an open structure to implement, simulate and share ecological models. A generalized object-oriented architecture improves readability and potential code re-use of models and makes simecol -models freely extendable and simple to use. The simecol package was implemented in the S4 class system of the programming language R . Reference applications, e.g. predator-prey models or grid models are provided which can be used as a starting point for own developments. Compact example applications and the complete code of an individual-based model of the water ﬂea Daphnia document the eﬃcient usage of simecol for various purposes in ecological modeling, e.g. scenario analysis, stochastic simulations and individual based population dynamics. Ecologists are encouraged to exploit the abilities of simecol to structure their work and to use R and object-oriented programming as a suitable medium for the distribution and share of ecological modeling code. Note


Introduction
The R system with the underlying S programming language is well suited for the development, implementation and analysis of dynamic models. It is, in addition to data analysis, increasingly used for model simulations in many disciplines like pharmacology (?), psychology (?), microbiology (?), epidemiology (?), ecology (???) or econometrics (??). Existing applications already cover a range from small conceptual process and teaching models up to coupled hydrodynamic-ecological models (?). Small models can be implemented easily in pure R (?) or by means of the XML-based Systems Biology Markup Language SBML and the corresponding Bioconductor package (?). For more complex or computation intensive simulations R is primarily used as an environment for data management, simulation control and data analysis, while the model core may be implemented in other languages like C/C++, FORTRAN or JAVA.
This works perfectly at the extremes, but problems appear with medium-sized models. While larger modeling projects usually start with an extensive planning and design phase carried out by experienced people, small models can be implemented ad-hoc without problems in R from scratch or by modification of online help examples. On the other hand, mediumsized applications often start by extension of small examples up to ever increasing size and complexity. An adequate design period is often skipped and at the end of a modeling project no time remains for re-design or appropriate documentation. The resulting programs are necessarily ill-structured in most cases or at least exhibit a very special, proprietary design.
The situation is even worse in ecological modeling, because this discipline is broad, modeling strategies vary substantially and ecological modelers are very creative. Different families of models (e.g. statistical, differential equations, discrete event, individual-based, spatially explicit) are applied alone or in mixtures to different ecological systems (terrestrial, limnetic, marine) and scales (individuals, laboratory systems, lakes/rivers/forests, oceans, biosphere). Not enough that there is a Babel of programming languages and simulation systems, there is also a Babel of approaches. There are often cases where it seems to be necessary to understand the whole source code when one tries to modify only one single parameter value or to introduce a new equation and it is not seldom easier to re-write code from scratch than to reuse an existing one. We aim to propose a possible way out of the dilemma, an open structure to implement and simulate ecological models. The R package simecol is provided to demonstrate the feasibility of this approach including a starter set of examples and utility functions to work with such models.
After giving a description of the design goals and the specification of the simObj class in Sections ?? and ?? we demonstrate basic use of the package in Section ??. The different mechanisms available to implement and simulate simObj models are explained in Section ??. A complete individual-based model is given in Section ?? to elucidate how to use and extend simecol in the modeling process. Finally, we discuss perspectives of R and simecol in ecological modeling as well as relations to other packages (Section ??).

Design goals
Our first goal is to provide a generalized architecture for the implementation of ecological models. Such a unified style, which can be considered as a template or prototype of model implementations, provides manifold advantages for a scientific community. The structured architecture will increase readability of the code by separating model equations from other code elements, e.g. for numerical techniques. This will enable ecological modellers to use R as a communication medium and allows to distribute model source code together with its documentation, e.g. as executable part of the "standard protocol for describing individualbased and agent-based models" suggested by ?. elsewhere e.g. at https://www.systemdynamics.org/ or Chapter 8 of ?. These simulation systems can be extremely effective for the specific class of applications they are intended for, but, they often lack the full power and flexibility of a programming language. In such cases, model frameworks or simulation libraries are commonly used to support one specific model family, e.g. PASCAL templates for ordinary differential and delay-differential equations (?), an object-oriented C++ framework like OSIRIS (?) for individual-based models or the Objective C framework SWARM 5 for agent-based simulations.
An alternative approach is the use of high-level programming environments and matrix oriented languages like MATLAB 6 or R (?). Such languages allow a more interactive development cycle, compared to compiled languages, and outweigh their performance handicap by efficient algorithms and compiled libraries for numerics and data management. Both openness and interactivity have made the R system a universal scripting interface for the free combination of a large diversity of applications in statistics and scientific computing.
The second design goal is to be as open as possible and to take advantage of the open philosophy of R. Users should be allowed to employ the full power of R's graphical, statistical and data management functions and to use arbitrary code written in R or compiled languages. The complete code of simecol-models should be published under a license that minimizes dependence from others and guarantees unrestricted use in science and education, including the possibility to be modified by others. Within this context, simecol is intended to provide a framework on the meta-level identifying structure-components of ecological simulation models.
Our third design goal is ease of use and simplicity. One of the main characteristics of programming languages like S and R is that users become programmers (?). Unfortunately, ecologists are commonly not well-trained in programming, which often hampers their application of models. Therefore, we aim to provide a software layer that bridges this gap and helps ecologists to work with models. In consequence, this means for simecol that simplicity of implementation is more important than efficiency. The system should support a broad level of user experience -in our case ecological models covering the whole range from teaching models to research applications.
From the perspective of a first time user it should be possible to run simulations without knowing too much about R and implementation details. A simulation of an ecological model should be as easy as fitting a linear model in R (see example in Section ??). A number of memorable "commands", i.e. a few essential but not overwhelmingly extensive generics for simulation, printing, plotting and data access, and utility functions accompany this package. Both the functions and also the simulation models should have meaningful defaults to enable new users to get almost immediate success and to enable experienced developers to structure their applications and to avoid unnecessary copy and paste (?).

Approach
The approach follows directly from the design goals to provide (i) a standardized structure, (ii) open and reusable code and (iii) ease of use of "the model". It is almost self-evident to apply an object-oriented design, consisting of: 1. A general and extensible class description suitable for ecological simulation models that allows sub-classes for different model families and multiple instances (objects of class simObj) which can be used simultaneously without interference, 2. Generic functions which work on objects of these classes and behave differently depending on the model family they work with.
All equations, constants and data needed for one particular simulation should be included in the model object, with the exception of general and widely needed functions, e.g. numerical algorithms. In the following sections we first analyse what is generally needed and then describe the particular approach.

State space approach
Most ecological models can be formulated by means of a state space representation, known from statistics and control theory ( Figure ??). This applies to dynamic (discrete resp. continuous) systems as well as to static, time independent systems when postulating that the latter case is a subset. A general description that is valid for both linear and nonlinear systems can be given as:ẋ where x is the state of the system andẋ its first derivative, t the time, u(t) is the input vector (or vector of boundary conditions), and y is the output vector. The functions f and g are the state transition function and the observation function, respectively, which rely on a vector p of constant parameters.
A simulation of a dynamic system is obtained by applying a suitable numerical algorithm to the function f . This algorithm can be a simple iteration or, when f is a system of ordinary differential equations, an appropriate ODE solver or a function giving an analytical solution.
Compared to the usual statistical models in R, ecological models are more diverse in their structure and exhibit tight relationships between procedural code (methods, equations) and data. Non-trivial ecological models are based on more or less modular building blocks (submodels), which are either base equations or complex models themselves.

The simObj specification
In essence, what do we need to implement a not too narrow class of ecological models? We need self-contained objects derived from classes with suitable properties and methods (data slots and function slots) resulting from the state space description: state variables, model equations and algorithms, model parameters (constants), input values, time steps, the name of an appropriate numerical algorithm (solver), and an optional set of possibly nested submodels (sub-equations). These parts are implemented as slots of the simObj class from which subclasses for different model families can inherit (Figure ??).
A small set of supporting functions is provided to work with these objects, namely:   ❼ Generic functions for simulation, printing, plotting, slot manipulation (accessor functions) and object creation (initialize functions), ❼ Utility functions, e.g. neighborhood relations for cellular automata.

Generic functions
In the S4 class model of the S language methods are based on generic functions. They specify the behavior of a particular function, depending on the class of the function arguments (?). All generic functions in simecol are defined as default methods for the class simObj and specific methods exist if necessary for subclasses. If new subclasses are defined for additional model families by the user it may be necessary to create new methods that work with these user-defined data types and provide the required functionality.

Simulation
The core function to work with simObjects is the generic function sim(simObj, ...), which, for dynamic systems, simulates an initial value problem using the initial state, boundary conditions, equations and parameters stored in one particular simObj instance by calling the numerical algorithm referred by its name in the solver slot of simObj. Common for all versions of sim is the pass-back modification behavior, i.e. a modified version of the original simObj is returned with a newly added or updated slot out holding the simulation results: The functionality of sim can vary for different subclasses of simObj e.g. odeModel, gridModel, rwalkModel. This behavior results mainly from a different data structure of the state variables and the set of numerical algorithms that are adequate for a given family of ecological models. Whereas ODE models have a vector for state and a data frame for outputs, grid models may have a grid matrix for state and a list of grids (one grid per time step) as output, and finally, random walk models may have a list for the initial state of the particles and a list of lists for the output.
The returned simObj can be printed and plotted directly with appropriate functions, the simulation results can be extracted with out or the resulting simObj can be used in another simulation with modified data or functionality.

Accessor functions
Similar to the out function other accessor functions are available for all slots with (in opposite to out) not only read but also write access. These functions are used similarly like the base function names and work with the appropriate data structures, see help files for details. The functions allow to change either the whole content of the respective slot or to change single elements, e.g. parameter values, time steps or equations. For example, the following will change only the parameter value of k1: An entirely new parameter is added to the parameter vector via: Elements can be deleted when a modified version of the parameter vector is assigned: The behavior is analogous for all other slots with the exception of out, given that the correct data type for the respective slot (vector, list or matrix) is used.
In addition to the command line accessor functions, graphical Tcl/Tk versions exist (editParms, editTimes, editInit) 7 , however, more complex data types cannot be handled yet by these functions.

Numerical algorithms
In order to simulate ecological models of various types, the appropriate numerical algorithms can be plugged into the sim function either by using an existing function, e.g. from this package, by imported solvers of package deSolve or by user-defined algorithms.
The algorithm used for one particular simObj is stored as character string in the solver slot of the object. User-defined algorithms have to provide interfaces (parameter line, output structure) and functionality (see below) that fit into the respective object specification and are compatible to the data structures of the particular class.

Utility functions
A few utility functions are provided for overcoming frequently occuring problems. However, it is not planned to overload simecol with numerous utilities as most of them are applicationspecific. Additional supporting functions should be written in the user workspace when they are needed or may be included in optional packages.

Interpolation
Dynamic systems often require interpolation of input data. This is particularly important for ODE solvers with automatic step size adjustment and there are cases where excessive interpolation outweighs the advantages of automatic step size determination.
The performance of linear approximation is crucial and we found that the performance of the respective functions from the base package can be increased if approxfun is used instead of approx, if matrices are used instead of data frames and if the number of data (nodes) in the inputs is limited to the essential minimum. In addition to this, two special versions approxTime and approxTime1 provided by simecol may be useful, see the help file for details.

Neighborhood functions
The computation of neighborhood is time critical for cellular automata. Two C++ functions, eightneighbors (direct neighbors) and neighbors (generalized neighborhood with arbitrary weight matrices) are provided for rectangular grids. The implementation of these functions is straightforward and may serve as a starting point for even more efficient solutions or other grid types, for example hexagonal or 3D grids.
Neighborhood functions can also be used for spatially explicit models. Models of this family commonly include both, an explicit spatial representation of organisms (in most cases with real-valued locations) and a grid-based representation of environmental factors (??).

Example models
A set of small ecological models is supplied with the package. These models are intended as a starting point for testing the package and for own developments. The models are provided in two versions, as binary objects in the data directory and in full source code in the directory "examples". The number of example models is intentionally limited and will grow only moderately in the future. In addition to this, ecological models which follow the simObj specification are well suited to be published and shared between scientists either as single code objects or in domain specific packages. At the first level of experience, users can simply explore example models supplied with the package or provided by other users without carrying too much on implementation details. They can be loaded with source from harddisk or the Internet, for example the stochastic cellular automaton shown in Figure ??: Note, that the sim function uses pass-back modification, i.e. the result is the complete simObj with the model outputs inserted. The advantage is that the resulting simObj is consistent, i.e. the model output corresponds to the equations, parameters and other settings of the simObj. Now, the settings may be inspected and changed, e.g. the number of time steps:

Predator-prey model
A second built-in demonstration example of simecol, is the elementary Lotka-Volterra predatorprey model, which can be given by two ordinary differential equations: In order to reproduce a schoolbook example two scenarios may be created by modifying two copies (clones) of lv: We now inspect default settings of initial values and parameters, modify them as required for lv2 and simulate both scenarios: The outputs of lv1 and lv2 can be compared visually using the plotting method of the odeModel class (plot(lv1)) or with regular plotting functions after extracting the outputs (Figure ??). It is quite obvious that scenario 1 produces stable cycles and that scenario 2 is at equilibrium for the given initial values and parametrization, because of: It is a particular advantage of R, that the complete set of statistical functions is immediately available, e.g. to inspect summary statistics like the range: The identity of the lower and upper limits for scenario 2 confirm the equilibrium state. Moreover, the period length of the cycles of scenario 1 can be analysed by means of spectral analysis: R> tlv <-times(lv1) R> ots <-ts(o1[c("predator", "prey") (dx1, dx2)) }, parms = c(k1=0.2, k2=0.2, k3=0.2), times = c(from=0, to=100, by=0.5), init = c(prey=0.5, predator=1), solver = "rk4" ) 5. Implementation of simecol models

Lotka-Volterra model
The implementation of the Lotka-Volterra equations is straightforward and results in a compact S4 object (Table ??). The two equations (Eq. ??, ??) can easily be put into the main function and there is no need for sub-equations. The code can be made even simpler without the two assignments at the beginning of main, but with respect to more structured models we found it generally advantageous to keep the default values of the names in the parameter line and on the other hand to use common symbols in the equations.

Models with nested subequations
For large models with numerous equations or for models with alternative (i.e. exchangeable) submodels it may be preferable to use a separate structure. Although simecol principally allows implementing subroutines as local functions of the main slot or even directly in the user workspace such a strategy would not be in line with our design goals. Instead, the equation-slot of the simObj class definition provides the structure where relevant submodels and model equations are stored. Consequent usage of the equation slot helps to increase the readability of the main function, leads to more structurized code and complies with the objectoriented paradigm. Moreover, the equation slot can be used to store alternative submodels, see Table ?? for a small example.
In this example, two versions of the functional response can be enabled alternatively by assigning one of f1 or f2 to f via equations (last line of the Table ??) and with the same mechanism it is possible to introduce further functional response curves.

Input data
The simulation models presented so far are autonomous, i.e. they have no external forcing data (matrix u in Figure ??). Such time dependent data, e.g. food availability or meteorological conditions, which are required in many practical cases can be provided in the inputs slot. In order to give a minimal example we may create a new odeModel by modifying a clone lv_ef of the elementary predator-prey model. To enable external forcing a modified version of the main slot is introduced, that simulates substrate (S) dependent growth of the prey population: Note, that inputs are converted into a matrix for performance reasons because otherwise repeated conversions were performed by approxTime1, or similarly by approx, which would be time consuming, especially for larger input data sets.
The resulting model can then be easily simulated and plotted and results in Figure ??:

Initializing
Sometimes, it may be required to perform computations while initializing a simObj. This may be either required to ensure consistency between different slots (e.g. parameters, inputs and initial values) to perform error checking or to create non-deterministic variants. Initializing methods, which exist in R as class methods of the generic initialize, are called either explicitly or implicitly during object creation. The syntax allows, in principal, two different modes of use. One can either provide all slots of the object in creation as named arguments to new or one can provide an existing simObj as the first un-named argument to initialize in order to get a re-initialized clone. 9 In the case of simObj this mechanism is extended by an optionally existing function slot initfunc, which is executed during the object creation process. Object creation is then as follows: in the first step an incomplete object is created internally via new according to the slots given and in the second step this object in creation is passed to the obj argument of initfunc which performs the final steps and returns the complete object.  Table 3: Predator-prey simulation with stochastic input variables. The example is derived from the externally forced object lv_ef. An initialisation function initfunc is provided which is called by initialize and returns a re-initialized obj with a new random sample of input values. The utility function fromtoby is used to expand the time vector from its compact form c(from=, to=, by=) into a sequence.  Table ??.
In the example shown in Table ?? new instances with different stochastic realizations of the input variables are created and simulated (see Figure ??). Note that initfunc is called automatically every time, when new instances are created via new or initialize.
6. Creating own models 6.1. The modeling cycle The modeling process is an iterative cycle of tasks (see ?). It begins with the formulation of (i) questions and (ii) hypotheses, (iii) the translation of these questions into a specific model structure, (iv) the implementation of the model in form of computer software, (v) the analysis, test and (in most cases) revision of the model, and (vi) communication of the model and its results to the scientific community. Another view is given by ?, who with respect to software modeling suggested to distiguish three different perspectives: 1. Conceptional perspective, 2. Specification perspective, 3. Implementation perspective.
These perspectives are complementary to the tasks defined by ? when tasks (i)-(ii) are regarded as conceptional and task (iii) as specification. In the following we concentrate on task (iv) to explain by means of a real, but still simple example, how a specified model can be implemented using R and the simecol software.

Conceptional perspective: an individual-based model of Daphnia
The scientific purpose of the Daphnia model given here was the analysis of demographic effects of Daphnia (water flea) populations. Two main hypotheses should be tested: ❼ Size-selective predation leads to an increased population mortality rate, compared to non-selective predation by fish (?).
❼ In comparison to predictions from the conventional Lotka-Volterra approach the inclusion of demographic effects results in a delayed but then inexpectedly rapid decline of abundance during periods of food limitation due to ageing effects (?).
Due to a multiple number of important features, the genus Daphnia (water flea) is an outstanding model organism in limnology, toxicology and population ecology, so results derived on this example may be of general interest to other areas as well.
The Daphnia model consists of three general parts: 1. A semi-empirical model of temperature and food dependent somatic growth and egg production derived from laboratory experiments (TeFI = temperature-food-interaction model) according to ?, 2. An empirical function of egg development time after ?,

A non-spatial individual-based simulation of population dynamics.
Individual-based models (IMBs) are a popular technique in ecological modeling (???). It is our aim to demonstrate how such models can be implemented with simecol.

Model specification
The state of the system is defined as a sample of individuals, each having four states: age, size, number of eggs and age of eggs. Population development is dependend on two environmental variables, food (phytoplankton, given in mg L −1 carbon) and temperature (in ❽). The model is simulated in fixed time steps (usually 0.1 day) over a period of several days up to a few months. The time scales are selected with respect to the egg development time, which is about 4.4 days at 15 ❽(?).
The life cycle of Daphnia starts with the release of neonate individuals of a given size (L 0 ) from the brood chamber of the mother into the water. Somatic growth follows the von Bertalanffy growth equation (?), depending on several empirical parameters. As soon as an individual reaches a fixed size at maturity (SAM) a clutch of eggs is produced (spawned), whereby the clutch size (number of eggs) is controlled by food availablility. After a temperature dependend egg development time the individuals from this clutch are released (hatched) and the cycle is started again (parthenogenetic, i.e. asexual reproduction).
Mortality can be modelled with arbitrary deterministic or stochastic mortality functions, e.g. size dependent mortality due to fish predation, but for the first simulation a deterministic fixed maximum age is used. All equations and parameters are given in detail in ? and although more elaborate bioenergetic Daphnia models are available in the meantime (??), the relatively simple model given here should be sufficient for the intended purpose.

Definition of a user-defined subclass
A subclass for non-spatial individual-based models was not available in past versions of simecol, but could be easily derived from simObj as class with appropriate data types, in particular with a data.frame for the table of the individuals stored in init: R> setClass("indbasedModel", + representation( + parms = "list", + init = "data.frame" + ), + contains = "simObj" + ) Since simecol version 0.8-4 class indbasedModel is a built-in class. The code snippet above is left in this document as an example how to derive user-defined subclasses.

Implementation of the model equations
The implementation which is provided in Table ?? and ?? starts with the selection of an appropriate data structure for the state variables. A table of individuals with four columns: age, size, number of eggs and age of eggs (eggage) is realized as data frame with one row for each individual. The data frame is initialized with an arbitrary number of start individuals (e.g. one single juvenile animal in the init slot).
The main function simulates the life cycle of Daphnia and calls the sub-equations live, survive and hatch which implement the following processes: live: The age of all individuals and the egg-age for individuals with eggs is incremented by the actual time step DELTAT. Then, the empirical function tefi is called to estimate length and potential number of eggs as a function of age, food and temperature. The data frame of individuals is then updated and for all adult individuals (size > size at maturity, SAM) which actually have no eggs the appropriate number of eggs is initialized.
survive: The survival function returns the subset of surviving individuals. Note that it is particularly easy in R to implement survival with the subset function by simply applying a logical rule to the table of individuals and returning only those rows which match the condition.
hatch: In a first step the the actual egg age is compared with the egg development time.
Then the total number of mature eggs is counted and a data frame (newinds) with an appropriate number of individuals is created (function newdaphnia. Population growth occurs by appending the data frame with newborns (newinds) to the data frame of the surviving (inds).
All functions of the life cycle receive the actual table of individuals (init) as their first argument and return an updated table inds which is then passed back to init.
The model is simulated by iteration over the time vector. Note that in contrast to ODE models the main function explicitly returns the new state and not derivatives. To account for this, the iteration algorithm is to be used here and not one of the ODE solvers like euler, rk4 or lsoda.
A number of constant parameters is needed by the empirical model (see ?, for details), which are represented as list within the parms slot.

Model simulation
With the ibm_Daphnia object derived from the indbasedModel class and given as complete source code in Tables ?? and ?? it is now possible to perform a first simulation:

R> solver(ibm_daphnia) <-"iteration" R> ibm_daphnia <-sim(ibm_daphnia)
This already works with the iteration method provided by the package, but the default behavior may not be optimal for the concrete subclass.
One disadvantage here is the fact that the default iteration algorithm stores the complete state data structure (i.e. the complete data frame) for each time step as list in the out slot. This behavior is rather memory consuming for individual-based simulations with several hundred or thousand individuals. Moreover, no adequate plotting method is currently available for such a list of data frames and therefore the default plot(simObj) method simply returns a warning message.

Class-specific functions and methods
Depending on the complexity of the model it may be necessary to supply either an own solver function or a complete sim method. The difference is that only one sim method is available for one particular class, but several solver functions may be provided as alternatives for one class, either as S4 methods with different names or as ordinary (non generic) functions.
In most cases it should be sufficient to write class specific solvers, but for complicated data structures or hybrid model architectures it may be necessary to provide class specific sim methods. In case of the individual-based Daphnia model, the solver should be of type "iterator" but with additional functionality to reduce the amount of data stored in out. To do this, mysolver given in Table ?? has a local function observer, which for each time step returns one line of summary statistics for all individuals. Additionally, it would be also possible to write data to logfiles, to print selected data to screen or to display animated graphics during simulation.
The argument naming of the solver functions is compatible with the ODE solvers of the deSolve-package with respect to the first four arguments. Moreover, some essential functionality must be provided by all solvers: 1. Extraction of slots of the simObj (argument y) to local variables, expansion of y@times via fromtoby, attach and detach for the list of equations as given in the example, 2. Iteration loop or any other appropriate numeric algorithm, 3. Assignment of the special parameter DELTAT (optional, if needed by the model), 4. Accumulation of essential simulation results to the outputs (out-slot) and assignment of explanatory variable names, in this case done by observer.

R> solver(ibm_daphnia) <-"myiteration" R> ibm_daphnia <-sim(ibm_daphnia) R> plot(ibm_daphnia)
It is beyond the scope of this paper to provide an overview over simulation techniques or to answer domain specific questions about Daphnia population dynamics; however, the following example is intended to give an impression how simecol models can be used in practice. The example deals with the effect of size-selective predation, similar to the more extensive analysis of ?. Four scenarios will be compared: At the first step we create one clone of the daphnia_ibm-object, assign settings common to all scenarios and an initial sample population:

))
Then we replace the default survive-function with a more general one which depends on a user-specified mortality function fmort: Copies of object Sc0 are created and modified according to the scenario specification. In the example below we have two functions with constant mortality and two other functions where per capita mortality is higher for the larger or smaller individuals, respectively:  The result shows very clearly the influence of demography on population growth (Figure ??). Given that the population growth rate r without any mortality (i.e. equal to the birth rate b) is approximately 0.11d −1 in Sc0 and the mortality rate d is set to 0.1d −1 , it is plausible that the population growth rate in Sc1 is: In case of size-selective predation, demography has to be taken into account in order to get realistic estimations of r. The simulation shows an increased population loss in case of fish predation (Sc3, r = −0.05) and a lower effect in case of Chaoborus (Sc2, r = 0.07). Please see ?? for details and how the results may depend on fecundity of the prey, the shape of the selection function and the dynamics of predator and prey.

Discussion
The main contribution of the simecol package is the proposal of a generalized, declarative structure to describe ecological models. This structure was inspired by the state space representation used in control theory and is intended as a pragmatic solution to unify the upcoming diversity of R implementations of ecological models. The object-oriented simObj structure may be useful also in other areas and for other models like continuous-time Markov processes and stochastic differential equations.
With the set of examples presented and some additional models developed in our workgroup, the matrix-oriented R language was found to be well suited for model development (rapid prototyping) and model evaluation. According to our experience, a structured OOP style is more efficient compared to a purely functional style, or even worse, ad-hoc programming. The functional OOP system of R is different from languages like JAVA or C++ and the approach of generic functions for common tasks seems to be more appropriate for statistical data analysis than for ecological simulation models which have not only variable data but also variable code. Moreover, the lack of references and the invisibility of member variables in slot functions of the same object was seemingly inconvenient and needed re-thinking. However, the R language with its Scheme heritage (?) is a "programmable programming language". Lexical scoping and local environments (?) allow to change its default behavior if needed. There were temptations to apply an alternative OOP paradigm that allows for references e.g. R.oo (?) or proto (?), but it was decided to stay with the default behavior as much as possible. Similarly, we used only flat object hierarchies and abandoned delegation-based approaches and instead suggest cloning (creation time sharing) as a standard technique to create derived objects.
At a first look R seems to be less suited for large applications, e.g. turbulence models, where C and FORTRAN are standard or for complex individual-based simulations with large numbers of interacting individuals, where class-based OOP in the flavor of C++ or JAVA is regarded as more natural (?). However, even such applications can take advantage of simecol, either because of vectorization in R (subset is in fact highly efficient) or due to the possibility to embed compiled code as shared library. For large applications or external simulation programs, simecol objects can be constructed as an interface provided that the external program is open enough to be linked or at least is callable in batch mode.
The package is designed to be open for local extensions and further evolution of the package itself. A limited number of classes will follow, e.g. for individual-based models similar to the Daphnia example or for purely statistical models like neural networks. An integrated parameter estimation functionality may follow as well as an interface to quantitative and qualitative model evaluation criteria (?). Moreover, interfaces to other promising approaches to solve simulation models in R may be worth to be established, e.g. to the XML based description language of the bioconductor package SBMLR (?), or to the nonlinear mixed effects modelling package nlmeODE of ? who independently developed a similar list-based object structure for a class of ordinary differential equation models.
Another appealing approach is the stoichiometry-matrix based approach for ODE models of aquatic systems (wastewater treatment, biofilm, rivers and lakes, https://www.eawag.ch/ organisation/abteilungen/siam/lehre/Modaqecosys/). These R scripts, developed by a prominent group in water modeling (e.g. ??), are currently used for teaching of aquatic modeling together with model assessment like sensitivity and uncertainty analysis, optimization, and frequentist or Bayesian model tests.