The first decision was that of programming language. We wanted a design driven by the immunological domain we are simulating, while maintaining the high performance necessary for such a complex simulation, which would have to deal with reaction-diffusion equations as well as the intracellular dynamics and interactions of thousands to millions of leukocytes and parenchymal cells. Our specific needs for a language included support for robust pseudorandom number libraries, multi-dimensional arrays and parallel computation. While languages such as Fortran and C fit these requirements well, they failed at our other requirement of a flexible high-level interface for description of the model domain. C++ includes C as a subset while providing high-level object-oriented facilities, and also has the advantage of being well-known and familiar to the scientific community. For better or worse, C++ has strong static typing, and this basic feature to a large extent dictates usage of the language. Strongly typed languages are comparatively inflexible and are poorly suited to rapid development, particularly in exploratory modes of programming. As a compromise, we opted to use a mixed language environment, where an interpreted dynamically typed language extends a statically typed compiled language. This hybrid system to some extent inherits the advantages of both languages.
For pragmatic reasons, we decided on a Python/C++ hybrid. Python is a popular object-oriented interpreted language in widespread use as a scripting language, but it also has good and ever-improving scientific support. It is also well-suited as a glue language, facilitating the connection of the base simulation code with databases and visualization systems. An important technical tool that made such a hybrid system feasible was the existence of the Boost Python library [6
], which offers a sophisticated C++ API to convert functions and data types automatically to and from the two languages. One approach to exploiting such a hybrid system is to code most of the simulation in Python, using C++ only for the time-intensive routines. We were constrained, however, by the need for high performance and the consequent need to develop for implementation on Beowulf cluster architectures; the Python MPI wrappers are, for this purpose, unacceptably slow. We chose instead to create two functionally identical systems, one in Python and one in C++. Synchrony between the two was enforced by writing a suite of Python unit tests, which depending on a switch, tests either the Python directly or C++ system via Boost Python wrappers. This seems to enhance productivity: it is simpler to code and test prototypes in Python and, at the appropriate point in their development, port them to C++ than to develop and test prototypes from scratch in C++. We also take advantage of flexible and efficient arrays in Fortran by calling individual Fortran routines from C++ for low level matrix and vector computations.
Just as fundamental were the choices we made for our software development environment. We chose to use the distributed version control system Mercurial [7
], because a distributed version control system offered superior support for multiple developers in terms of branching, merging and disconnected operations. For software builds, we chose Scons [8
]. Among other benefits, Scons automatically manages file dependencies, and allows us to customize our build scripts using the full power of Python. We elected to use Trac for bug and issue management, as it allows us to manage version control, bug tracking and documentation within the same web-based package. Not coincidentally, all the tools we chose are Python based, which makes it possible for us to extend these tools if necessary.
One of our design constraints was that the primary domain for our simulation was molecular and cellular immunology, not mathematics, physics or computer science. At the least, this meant that we had to be able to specify and describe a simulation run using the vocabulary of immunology. Another constraint was that we were constantly researching new prototypes, adding more functionality to the system and investigating more efficient data structures and algorithms for performance bottlenecks, so the system had to be extensible; ideally, we would be able to swap in new functionality for old with minimal additional effort. To meet these constraints, we decided on a hierarchical modular system. The implementation of such a system was naturally object-oriented, a paradigm both C++ and Python support. At the top level, we split the project into Catalog, Core and Graphics/Analysis modules. Within each of these were separate sub-modules and classes that actually implemented the functionality of the system. This organization allowed us to develop individual modules independently, so long as there was a standard API that allowed communication between modules.
The Catalog module is a repository of biological and physical information. Its entries include, for example, parameters governing the locomotion and behavior of specific leukocyte types and the diffusion coefficients and reaction rates of cytokines. The information in the Catalog is stored in a PostgreSQL relational database, and simulation classes are mapped to database tables using the strategy of mapping each inheritance tree to a table. To facilitate data entry into the Catalog, the Python library
SQLAlchemy citesqlalchemy was used to provide an object-relational mapping, making it possible for developers to create new class specifications without knowing SQL. Our plans call for the development of a web interface to the Catalog, so that experimentalists and other non-developers can contribute items to the database. After describing the Experiment module, we will describe how the Catalog is used to construct class entities in a simulation at run-time.
The Core module is where the code for running a simulation actually resides. To allow the description of a computational experiment using immunological terminology, there are only a few base classes in the public interface that a typical simulation is required to specify, namely Simulation, Environment, Cell, Vessel and Soluble_factor (Figure ). We have designed our classes to have a relatively flat inheritance structure and to use abstract base classes where possible. We find that this design promotes loose coupling and facilitates maintenance.
Figure 1 Class composites and aggregates in the Core module. Top level classes (Environment, Vasculature, Cell, Solfac) are shown as diamonds and component classes that provide specific functionality as ovals. Instances of component classes are passed in as parameters (more ...)
The Simulation class manages administrative tasks, with the ability to run or step through a simulation and log simulation details. It will also provide hooks to interface with run-time logging, computational steering and visualization systems when these are developed.
The Environment class represents a physical spatial volume, which may be simple, such as when modeling an in vitro
environment, or complex, such as when modeling physiological tissue. Since our task allocation is done in a parallel environment using a spatial decomposition, each processor in a parallel run will contain its own Environment instance. Most of the work done by the Environment class is delegated to three private classes which handle the reaction-diffusion equations (Diffusion), collision detection (Occupation) and cellular immigration and emigration (Trafic). Since these classes are an implementation of the well-known Strategy design pattern, it is simple to upgrade the functionality provided by simply pointing to a new class. For example, we recently replaced the original forward Euler scheme with a Multigrid reaction-diffusion algorithm [9
], resulting in an order of magnitude speedup. This was accomplished by simply replacing the Diffusion class in Environment. The Environment also serves as a container for various cell types (Cell), blood and lymphatic vessels (Vessel), and various cytokines and chemokines (Soluble_factor), and manages the interactions among these components.
Similar to Environment, the Cell base class is also a composite class, with distinct behaviors delegated to their own classes, and hence also extensible. A cell's functional state is given by, among other things, the local concentrations of various proteins in each of their post-translationally modified states (e.g., phosphorylated on specific residues). These characteristics of a protein are themselves dynamically regulated, and typically modeled using nonlinear ordinary differential equations. Importantly, these proteins determine the functional behavior of the cell, and it is sufficient to couple cell behavior modules to these protein concentrations to model cell stimulus-response behavior. The protein readouts are encapsulated in a State class, which tracks and exposes the fluctuating protein levels. In our current implementation, there is a mechanistic model that governs the protein levels, in which the rates of protein synthesis and/or degradation are governed by an explicit cascade involving model classes for Receptors, Inducers, Promoters, Genes and regulatory Proteins. It would be straightforward, however, to replace this mechanistic model with a probabilistic model, or even a much more detailed mechanistic model depending on the goals of the simulation. By having a flexible granularity, we can study an immune response at different levels of spatiotemporal resolution. When studying the behavior of individual cells, for example, we can increase the resolving power, while for studying large population behaviors, a coarser approximation of internal cell dynamics might be sufficient. Other behaviors like Motility and Secretion can then base the quality/quantity of their response at any time step on the internal environment represented by the State class and/or the external environment (e.g. concentration or gradient of various cytokines or cell-cell contacts). This design allows an extremely flexible stimulus-response coupling between the internal and external cellular environments and cellular behavior. Because of the loose coupling and the ability to change cellular behavior by "plug-in", an inheritance hierarchy to represent the different cell types is not necessary. Instead, we just "plug-in" different behavior modules for the various cell types to implement their observed behaviors. In practice, this means that specific behaviors for a particular cell type are implemented as instances of their respective classes that are given to the cell class constructor as parameters.
The Vessel class is used for the representation of the blood vessels that permeate host tissue. Leukocytes typically enter tissue by migrating across the endothelial barrier at capillaries, and egress at post-capillary venules to enter the lymphatic system. The rate of entry and egress of various cell types is regulated by pro-inflammatory cytokines, which regulate the adhesiveness and permeability of the blood vessels, as well as by chemokines, which influence leukocyte migratory properties. Our implementation of blood vessels is simple; the entry and egress vessels are stored in matrices representing positions in space, and the probability of cell trafficking at a point depends on the permeability and adhesiveness of the corresponding matrix entry, which in turn is sensitive to the local concentration of the relevant soluble factors. We currently assume that the blood vessels are homogeneously distributed; it is trivial to relax this assumption to reflect tissue histology more accurately by masking certain matrix entries.
The Soluble_factor class is essentially a container for a array of concentrations, with a few accessor and mutator methods for convenience and some basic attributes including its diffusion coefficient. An instantiation of the Soluble_factor class is made for every cytokine, chemokine and other soluble factor of concern in the simulation. The array is updated by the Diffusion class which takes into account secretion from cellular sources, adsorption by cells and extracellular matrix, binding reactions between soluble factors, and diffusion. There are several classes derived from Diffusion that represent different methods for solving the diffusion problem, including a forward Euler method and a backward Euler method which incorporates multigrid.
While the bundling of distinct behaviors into separate classes certainly simplifies the dynamical configuration of classes, their instantiation can be clumsy, accounting for a significant proportion of the code required to run a simulation. Biological cells, however once adequately characterized, can be reused in many different simulations, just as an experimentalist would use the same cultured cell lines (possibly ordered from a catalog!) for many different experiments. This was the idea driving the creation of the Catalog module, and here we will describe how it is used for instantiation of Experiment classes.
We started by writing generic code that could construct a class instance without advance specification of the instance type or even its constructor parameters. This design provided the flexibility to easily incorporate new classes in the database, or new subclasses of currently existing types. To do so, we adopted the object factories scheme, which allows flexible creation of object instances using C++ templates [10
We also extended the factory to take arbitrary constructor parameters with partial template specialization, and wrote a Python script to generate templates for up to 20 arbitrary constructor arguments. Our specific implementation of the factory pattern follows [11
]. To create a factory for a particular class, we instantiate the factory and register it with the appropriate constructor arguments passed in as a string. While this process is a somewhat involved, it has to be done only once for each class; the typical user does not have to be exposed to any of it, since all that is required once the infrastructure is in place is to call a function such as create_cell("macrophage"), which will call the appropriate constructor for a macrophage. As described earlier, a cell such as a macrophage is a hierarchically organized composite class, and each subclass in the hierarchy will require the appropriate parameters to be passed in to its factory. We store the parameters for instances of each class in the open source relational database PostgreSQL, using foreign key constraints to relate the correct subclass instances in any given object hierarchy (Figure ). Foreign key constraints allow one field in a database table to refer to fields in another table, allowing information retrieval from tables with a matching entry. The component parameters can then be retrieved via such foreign key references when a specific cell type like macrophage is constructed, and this ensures that the appropriate parameters are used to instantiate all levels of the composite class. In practice, we connect to a local or remote database containing the data using the libpqxx C++ API for PostgreSQL at the beginning of the simulation, and also use the libpqxx utility functions to execute the necessary SQL statements to find the correct parameters when the constructor is called.
Database schema for the Cell class in the Catalog.
For visualization and analysis, we export simulation variables using the standard HDF5 format [12
], which can then be visualized using our own wxPython and OpenGL based visualization application or imported into a statistical package such as R or Matlab for further analysis.
The standard equations that describe the dynamics of the soluble factors as well as the cells are as follows. Suppose ci
) is the local concentration of soluble factor i
at the position x
and time t
. Let xμ
be the position of the center of cell μ
in the bounded domain Ω
. Then the reaction-diffusion equations that describe the behavior of the soluble factors from [13
is a smoothly cut-off Gaussian with support over the volume of the cell. Here Di
denotes the diffusion coefficient for soluble factor i
is the rate at which soluble factor i
is removed by interaction with soluble factor j
, and λi
is the rate of removal of soluble factor i
by other processes. These positive constants are stored in the Soluble_factor class. The secretion of soluble factor i
through the surface of cell μ
is approximated by a source term centered at the cell position xμ
with secretion rate Jiμ
). Although the model in [13
] contains point sources, the existence of a solution to (1) requires a smooth source term as discussed in [9
]. In the simulation, the rate Jiμ
) is determined by the Secretion member of each cell. The source term for each soluble factor is a sum over the M
individual cells indexed by μ
The system (1) is coupled with M
systems of stochastic differential equations that describe the motions of each cell individually. For each cell, indexed by μ
, denote by xt
the position and by
the velocity. The cell motion is modeled by the Langevin process
is a standard Wiener process in
. The cell velocity is stochastic but biased toward the direction of the gradients of the relevant soluble factors by the relation
Here χi is the chemotactic constant that controls how much the drift is influenced by the gradient of soluble factor i and h0 is the length of the gradient of i at which h attains half its maximum value. The magnitude of the Wiener process depends on the soluble factor concentration by
The positive constants γ, c0, h0, σ0, q and χi are stored in the Motility class, which is a member of Cell.
The reaction-diffusion-stochastic system (1), (2) contains a set of coupled partial differential equations that interacts with a set of stochastic differential equations. Instead of solving this complex system, we consider the following three problems for which there are known numerical methods:
1. The diffusion of the soluble factors,
2. The reactions of the soluble factors,
3. The motion of the cells (2).
Given initial concentrations c0
), initial positions x0
and initial velocities v0
, we can approximate ci
) and v
) by the following splitting scheme [9
1. Solve the diffusion, (5), for Δt using the initial data ui(x, t), and x0 and v0.
2. Solve the reaction, (6), for Δt using the solution from the previous step as initial data.
3. Move the cells according to (2) for one time step Δt using the soluble factors from the previous step and the initial cell positions and velocities, x0 and v0, as initial data.
The error due to this splitting is
). This result is well-known for ordinary and partial differential equations; we have extended it to the reaction-diffusion-stochastic system [9
]. Furthermore, if the numerical schemes for (5), (6) and (2) are at least
), then the resulting error is
). Here we see that operator splitting allows us to solve three tractable problems rather than one complex system. Splitting the reaction and the diffusion reduces a system of nonlinear partial differential equations to a system of linear partial differential equations for the diffusion, and a system of ordinary differential equations for the reaction. The partial differential equations can be solved using a backward Euler scheme with multigrid for the linear solve as described in [15
] and [16
]. The multigrid scheme is implemented as a member function of the Multigrid_diffusion class which is derived from Diffusion. This class contains a necessary set of parameters, as well as functions to compute the relevant sinks and sources from the list of cells. The reaction equations are solved numerically using a first order semi-implicit Euler scheme. A separate scheme for the reaction is also a member function of the Diffusion class; it updates the entire list of soluble factors for one time step. Finally, the Langevin process for the cell motion can be simulated exactly as described in [17
]. The Langevin motility class stores a set of parameters that are relevant to the equations and numerical scheme. It also contains member functions that use the concentrations of the soluble factors to update the velocity and return the displacement of the cell for one time step. Although the soluble factors and the cells depend on each other, operator splitting allows us to update them sequentially instead of simultaneously and further modularizes the simulation.
We have also implemented multiple schemes for diffusion, reaction and cell motility, which are derived classes from Diffusion and Motility (Figure ). The derived classes may represent a change in the dynamical equations themselves or in the numerical methods used to solve them. For each simulation, the user can choose a set of equations and schemes that best models the desired behavior, approximates the mathematical equations most accurately or produces results most efficiently. In particular, different cell types can display different motility behaviors within the same simulation. As researchers develop more realistic models and more accurate and efficient numerical methods, the system facilitates simple swapping of these new classes for old ones and simplifies comparison and testing.
Motility and Diffusion inheritance trees.
We run unit tests using the Python unittest module, with the only novelty being that since all our C++ classes are compiled to a shared library readable in Python, we can drive both Python and C++ unit tests with the same test code. At a slightly higher level, the user may compare simulated behavior with analytical predictions where closed form solutions (either deterministic or statistical) are available. Examples include testing cell motility modeled as Langevin processes or simple diffusion from a single point-source secreting at a constant rate. There are also informal eyeball tests, in which simulations with simple configurations are run and visually checked to ensure that the resulting behavior is qualitatively consistent with experimental results. We have not incorporated code coverage analysis [18
] or continuous build systems [19
] yet, but these would certainly be sensible additions.
Due to the highly modular nature of our code, it is generally simple to replace an inefficient module with a more efficient one. Identification of bottlenecks is done using the callgrind tool widely available on Unix platforms, and the kcachegrind GUI to visualize the resulting profile. Examples of such high level optimizations include the replacement of a forward Euler diffusion scheme with a Multigrid scheme, and replacing collision detection by looping over all cells with a grid-based method implemented as a hash table. At a lower level, we have removed unnecessary copies by using references where possible for large data structures, and rearranged code to minimize unnecessary looping. So far we have avoided doing micro-optimizations (e.g. loop unrolling) believing that such tasks are best left to the compiler. One simple optimization we did incorporate was to specify the appropriate compiler flags for each machine architecture, allowing the compiler to make optimal use of pipelining and vectorization.