Morpheus implements explicit finite different solvers of numerical integration of ODEs:

• Euler's method: 1st order
• Heun's method (a.k.a. midpoint): 2nd order
• Runge-Kutta's method: 4th order

and stochastic differential equations (SDEs):

• Euler-Maruyama's method for 1st and 2nd order.

The latter is automatically used when a stochastic white noise term rand_norm(0,[amplitude]) is included as below, but is not available for Runge-Kutta.

 <DiffEqn solver="heun" time-step="0.01" symbol-ref="X">
   <Expression> rand_norm(0, noise) </Expression>
 </DiffEqn>
 <Constant symbol="noise" value="1e-3"/>

Solvers with adaptive time stepping, such as the Runge-Kutta-Fehlberg method, are not available. Currently, we use fixed time stepping methods, because this greatly simplifies temporal scheduling of updates for automated model integration.

Although the implemented explicit solvers are general and flexible, they are not applicable to stiff systems. Stiff systems are ODE systems that can change so rapidly that explicit solvers will result in numerical instability, unless the time step is extremely small. This situation may occur in systems where very large and small values are multiplied.

Morpheus is a lattice-based simulation platform. This means that the spatial models are discretized on a lattice.

The following lattices are available:

• linear 1D lattice
• square 2D (orthogonal) lattice
• hexagonal 2D lattice
• cubic 3D (orthogonal) lattice

For MembraneProperties, properties that are resolved on the surface of cells, Morpheus uses a special lattices with polar coordinate system

• circular 1D lattice (for 2D models)
• spherical 2D lattice (for 3D models)

See description in FAQ

Partial differential equations (PDEs) are numerically approximated by separating the reaction and diffusion steps using the method of lines. This method discretizes the PDE into a system of coupled ODEs. These ODEs are then solved using the numerical methods mentioned above.

For diffusion, Morpheus uses the simple and general forward Euler scheme. The time step $\delta t$ is automatically adjusted according to the CFL condition to guarantuee numerical stability. Alignment of multidimensional lattices to 1D memory (using valarrays) results in highly efficient implementation of diffusion.

The unconditionally stable method of alternate direction implicit (ADI) is also implemented (for constant (Dirichlet) boundary conditions), but is currently not used.

Cell shape and motility is implemented according to the cellular Potts model (CPM)(Graner and Glazier, 1992).

In this model formalism, biological cells are represented as a domains of lattices sites $\boldsymbol{x}$ with identical index $\sigma$. In its simplest form (ignoring differential adhesion), the changes in configuration of cells on the lattice are governed by the Hamiltonian:

$H=\sum_{\sigma > 0} \lambda_V (v_\sigma - V_t)^2 + \sum_{\sigma > 0} \lambda_P (p_\sigma - P_t)$

where $v_\sigma$ and $p_\sigma$ are the actual volume and perimeter of the cell with index $\sigma$ and $V_t$ and $P_t$ are the target volume and perimeter. Deviations from these target values increase the free energy $H$ according to the scalars $\lambda_V$ and $\lambda_P$.

The CPM is a Monte Carlo method in which the lattice is updated by randomly sampling lattice sites.

Within the CPM, a Monte Carlo step (MCS) is often taken as a discrete unit of time. A single Monte Carlo step is defined as the number of random sampled updates equal to the number of lattice sites, i.e. within one step, each lattice sites would have had chance to be updated.

In each update:

1. A lattice site $\boldsymbol{x}$ is chosen at random.

2. From the neighborhood $N$ of $\boldsymbol{x}$, a second lattice site $\boldsymbol{x'}$ is chosen.

3. Then, the change in free energy $\Delta H$ is calculated if the state $\sigma$ at $\boldsymbol{x'}$ ($\sigma_\boldsymbol{x'}$) would be copied to $\boldsymbol{x}$ according to the specific Hamiltonian $H$.

4. The proposed update is always accepted when $\Delta H<0$ and is accepted with a Boltzmann probability when $\Delta H>0$:

$P(\sigma_\boldsymbol{x'} \rightarrow \sigma_\boldsymbol{x})=\begin{cases} 1 &\mbox{if } \Delta H<0 \\ e^{{(\Delta H-Y)}/T} &\mbox{otherwise} \end{cases}$

where $T$ (for 'temperature') modulates the probability of unfavorable updates to be accepted and can be taken to represent local protrusions/retractions of the cell membrane. The parameter $Y$ (for 'yield') is sometimes used to avoid oscillations with energy-neutral updates. This can be said to represent cytoskeletal resistance to membrane fluctuations.

By default, Morpheus uses the Mersenne Twister 19937 pseudo-random number generator included in the CPP GNU compiler (TR1). Alternatively, the Boost RNG can be used (specified in CMake options).

Note: In multithreaded simulations, each thread gets its own RNG (the random seeds for which are based on the specified seed of the master thread RNG). Therefore, to reproduce simulation results, not only the random seed needs to be specified (Time/RandomSeed), but also the same number of threads (Settings → Local→ threads).

To be documented.

In order to reference variables in mathematical terms in a multi-model environment one requires to identify the respective sub-model a variable is defined in and to perform a mapping between different (spatial) contexts, i.e. determine the corresponding CPM cell from a spatial position to access a cell property.

This ensemble of a variable definition and contextual meta-data we call a symbol, which can be referenced all over the multi-model environment. (Also allows to iteration over the elements od a symbol context).

In case of spatial ambiguity Morpheus provides a set of predefined Reporters to take care for the spatial mapping. For example, a PDEReporter should be defined to compute the average concentration of a diffusible within the area occupied by a particular cell.

Example Autocrine Chemotaxis

<IMAGE>

extern> http://imc.zih.tu-dresden.de/morpheus/examples/Multiscale/AutocrineChemotaxis.xml

The time scheduler takes care (i) to evolve numerical schemes in time, (ii) to schedule and execute spatial data mappings, and (iii) to run event-based schemes, i.e. schemes which are run in certain time intervals. In order to combine different schemes in a single scheduler we exploit the method of fractional time steps. The basis is for our automatic scheduling is extensive knowledge of the data dependencies of the schemes and what kind of updated data they provide. This information is represented as a symbolic dependency graph based on the expressions, input and output symbols defined within the schemes.

We envisage the following criteria for the design of our scheduler implementation.

• Time stepping for all schemes shall be automatically adjusted, as far as possible, to the maximum permitted time step.
• In addition, time stepping can be overridden from in the MDL, allowing for user guided step sizes.
• The order of the specification of the schemes in the MDL shall not matter.
• The sequential order of updates within a temporal step shall be determined by the symbolic interdependency graph.
• Tightly coupled systems, often bearing circular dependencies, have to be defined in a System tag to be updated jointly.

A robust automatic scheduling methods necessarily needs reliable information about the internal limits of numerical schemes and its external dependencies. External dependencies are extracted as a dependency tree. Internal time step limits are reported by the schemes and can be adoptive in time. Solving Reaction-Diffusion Systems, for example, using a finite difference scheme and operator splitting allows us to schedule the diffusion solver depending on the CFL condition and running the reaction part independently.

1. Adjusting time steps:
   • Continuous Time Solvers: Ensure numerical stability and .
   • Mappers: Ensuring updated data is available as needed, but not computed more often than necessary. We let the time step size depend on the time steps of the downstream dependency tree (e.g a spatial Reporter shall run as often as it's output is needed by another scheme), and on stability criteria of the numerical schemes (e.g. CFL for Diffusion).
2. Sorting of schemes to be updated sequentially: The order is determined by the dependency tree. Circular dependencies cannot be sorted unambiguously and thus are rejected.
3. Initialisation of initial time state (e.g. mappings)

1. Phase I: Synchronous Time Evolution Schemes. The schemes are executed sequentially or in parallel, but the results are buffered as intermediates and only applied at the very end of that phase, such that any scheme exclusively depends on the previous time step solutions. Tightly coupled systems, defined in a System tag, are evolved by a system solver.
2. Time Step Progression
3. Phase II: Sequentially Updated Schemes (Mappers, Events), depending on data that has already been updated.
4. Phase III: Output generation in Analysis plug-ins. Can be performed in any order, just because there is nor interdependency.
5. Checkpointing of simulation state.

Morpheus provides a plugin interface to extend the feature set. Plugins are written in C++ and requires recompilation from source. Therefore, this type of extensibility is only available when building from source (currently limited to developers and collaborators).

Morpheus has interfaces for different types of plugins:

• Initialization: Plugins with methods to define an initial condition. Examples: TIFFReader and InitPDEExpression.
• Reporters: Plugins with customizable spatial mappings used in situations in which spatial mapping cannot be automated due to ambiguity. Examples are NeighborsReporter and PDEReporter.
• Motility: Plugins with methods for cell movement. Examples: Chemotaxis and Persistence.
• Shape: Plugins with methods to constrain cell shape. Examples: VolumeConstraint and SurfaceConstraint.
• Analysis: Plugins for data logging, plotting or processing. Examples: Gnuplotter and HistogramLogger.
• Miscellaneous: Set of plugins that do not fit in the categories above. Examples: Proliferation and InsertMedium.

Writing new plugins requires you to add three files:

• Header (MyPlugin.h): declares variables and functions
• Implementation (MyPlugin.cpp): implements the core functions
• XSD (MyPlugin.xsd): specifies the rules what constitutes valid XML input (used by GUI)

See an example for an Analysis plugin below.