Chapter 11 Reusable, Extensible Simulation Code
11.1 Tabular data
Biologists want intimate control of simulation details to incorporate realism. Typically this involves events, that is significant biological changes that can be marked and counted. The events themselves can often be laid out in tables, which we can then use as raw inputs to simulations. These tables not only describe the life events, but they also identify the organisms, units, and other features, allowing the simulation software itself to be free of specific reference to any one system. At present these files all reside in library(ewing) in the data folder.
11.1.1 Organism Features
The simulation code has no organism information per se. All organism information is designed by the scientist in tables as ordinary text files that can be incorporated into our internal structure, which is appropriately called Organism. This internal structure contains simulation features, schedules for future events, and process information for interactions among organisms. For our prototype simulation, the primary organisms are Red Scale, Aphytis, Encarsia, Comperiella and a greatly simplified orange tree as a substrate. Here, orange is actually a static substrate with certain properties that are important for the life histories of the other organisms. The “master” information on these organisms is in the organism.features.txt file:
Notice that the orange “organism” has missing values (NA) for the features since it is considered static. In addition, redscale (Californai Red Scale) has no attack feature, as it is only a host for aphytis (Aphytis sp.), Encarsia (Encarsia sp.), Comperiella and other parasitoids. Some features such as offspring may have either a number or the name another organism. The number signifies a mean value, while the name of an organism signifies that there is some interaction for this feature. For instance, Red Scale have on average 10 offspring, but the number of Aphytis offspring depends on the condition of its Red Scale host. Units of “biological time” may differ among organisms; Red Scale is temperature sensitive and grows by degree-days (DD), while Aphytis is diurnal, dependent on hours (hr). Other columns have to do with other features of the simulation:
Note that we may have to model species differently in different environments. This could be done in a variety of ways. For instance, we could include another column for environment in the above table, with corresponding changes in code.
11.1.2 Organism Future Events
An organism may have many potential future events, many of which correspond to stages of its life history. Other events, such as sex determination, feeding, starvation and death can be included as well. For instance, here are the future event schedules and other information by life stages for Red Scale, in object future.redscale.txt:
The times here are the mean number of degree-days (DD) to the event. The time to a scheduled future event, say from third.1 to third.2 would be a random draw T from the future event distribution, in this case F(t) = 1-exp(-M(t/55)). The mean value function M(t) is assumed to be the identity unless it is otherwise provided (see the next section). The future events for Aphytis, future.aphytis.txt, are similar with some notable exceptions:
Here the adult Aphytis may either feed or oviposit (ovip) on its host. Future events 5 and 6 are competing risks, which in this case are sorted out in the simulation code based on the size of the selected host. [This is one part that is still particular to this host-parasitoid system.] The event column identifies the type of event; either future event, birth or death (both immediate events); or attack (another immediate event that involves interaction with other individuals). The “fid” column simplifies some coding when relating current to future stages. The last column specifies the initial relative abundance of life stages: both species are weighted heavily toward young organisms.
The other columns for both Red Scale and Aphytis are used for plotting. The pch and color are used when plotting the spatial arrangement (temporarily disconnected), while the ageclass specifies how the life events are to be summarized. Recall that Red Scale “host” and Aphytis “adult” ageclasses are to be plotted by substrate according to the organism feature file.
11.1.3 Interactions among Organisms
The third type of spreadsheet concerns interactions among organisms. While this could be encoded directly into subroutines, such a table can provide more generality to the simulation modeling system. Here is redscale.aphytis.txt, including information about Aphytis can use various stages of Red Scale:
Most columns have a 0-5 scale which serve as individual host weights during sampling. That is, if there are 10 hosts in the area, Aphytis will chose one based on the weights associated with their life stage. Thus, Red Scale prefers to oviposit on second and third instars, and is more likely to produce male offspring for smaller Red Scale. It will feed if the Red Scale is very small. The mean number of Aphytis eggs that an emerging Aphytis would have ready to lay is dependent on Red Scale stage as well.
The following type of interaction, orange.aphytis.txt, indicates how a species interacts with its substrate. In this case, Aphytis can move around the orange. Our prototype orange consists of a two-sided leaf, a 20-sided orange (icosohedron), and a twig connecting the two. Again, a 0-5 scale is employed. This file really only has information about movement among fruit-twig-leaf, while a second file, not shown, has topological information about movement within each of these (e.g. movement on icosohedron is restricted to the three adjacent triangles).
Notice that Aphytis is more likely to find fruit than leaf top, and twig has the lowest probability. Also, Aphytis tends to move short distances if on a fruit, and long distances if on a twig. Perhaps we need only rely on probabilities to transfer between orange components. A similar file exists for orange.redscale.txt, with perhaps slightly different values.
We are not completely settled on this system, but it provides a convenient starting point for spatial simulation. Fine detail movement is based on triangular grid system and is introduced in Section 5.
11.2 Event subroutines
Once the various organism files discussed in the previous section have been constructed and placed in the appropriate location, the simulation is ready to be initialized. The simulation is initialized with the command:
In this example, this command reads in the appropriate information for both Red Scale and Aphytis, and any dependent organisms identified in the organism.features.txt file (in this case, orange), and it prompts the user for the number of organisms, and initializes the populations. It also initializes other simulation features, such as the temperature regime (see Section 5). Once the simulation is initialized, the simulation is executed by using the command
This command executes 1000 future events of reproduction, death, parasitization and aging through the various life stages. The summary plot changes with each 50 events. In this de novo use, the time to next event is drawn from an exponential with the mean time to event given by the “future.redscale.txt” file as explained in the previous section. During a simulation, as currently implemented, there are two graphical time summaries for each species. These show number of organisms by age class (left) and number of organisms by spatial location (right). For instance, Red Scale is counted by crawler, host, and gravid (host is all stages that can be parasitized by Aphytis), while Aphytis is counted by young and adult. The three spatial locations are leaf, twig and fruit.
<<Figure 1 here>>
Summary tables are provided at the end of the run indicating timing results and, if desired, age distributions before, during and after the simulation.
Changing the number of either species, or the relative timing of future events, can lead to an exponential growth or decay in the Red Scale population. The simulation displayed in Figure 1 is far from equilibrium. Because the system is inherently nonlinear, it is extremely hard to predict behavior. It is challenging to adjust the number of each species, time to events, movement and number of offspring to balance the population.
11.2.1 Competing Risk Structure
The simulation sorts all events such that the next future event has minimal time. This sort technique typically only involves rearranging one or two events, and is very fast. Essentially we are reordering a priority queue in such a way that the top of queue contains the next event to be processed (see Section 6.2 for details). That event is processed by setting the current event to the future event, processing any immediate events (e.g. scheduling birth of crawlers for gravid females), and then proceeding to the next event. Death is treated in this way as an immediate event. At the time of death, the individual is removed from the simulation. Competing risks arise when there are more than one possible future events following a given current event.
The competing risk structure has been described in a companion paper (Ewing et al. 2001) and is only briefly referred to here. The important steps in the current implementation are
efficient computation of M(?) and M-1(?) using forward and backward splines
static lists of individual characteristics that do not change over the span
interaction among species
In a future development of the simulation, we plan to incorporate more realistic settings and feedbacks. Currently the competing risk structure for Aphytis allows for a choice between feeding and ovipositing, as shown in the future event structure of the previous section. Aphytis proceeds from egg to larvae to prepupae to pupae to adult with mean future event times in DD as indicated earlier. However, in the adult stage, Aphytis may either feed or oviposit. After each feed or oviposit, the adult has a number of choices. Feeding sustains the insect, allowing it to potentially oviposit more offspring at a later time, while ovipositing depletes resources, but contributes to survival of the species. The choice Aphytis makes depends on a number of environmental factors including the number of eggs it has in reserve, the last time ovipositing took place, the availability of hosts, and the condition of the Aphytis in question. The decision to oviposit is also determined by the size of the scale and the number of Aphytis eggs previously oviposited. Many but not all of these features have been implemented.