Morpheus

This page is read only. You can view the source, but not change it. Ask your administrator if you think this is wrong.
====== Technical specification ======

===== Solvers =====

==== Ordinary differential equations (ODE) ====

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

  * [[http://en.wikipedia.org/wiki/Euler_method|Euler's method]]: 1st order
  * [[http://en.wikipedia.org/wiki/Heun's_method|Heun's method]] (a.k.a. midpoint): 2nd order
  * [[http://en.wikipedia.org/wiki/Runge%E2%80%93Kutta_method|Runge-Kutta's method]]: 4th order

and stochastic differential equations (SDEs):

  * [[http://en.wikipedia.org/wiki/Euler%E2%80%93Maruyama_method|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.

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

=== Fixed time stepping ===

Solvers with adaptive time stepping, such as the [[http://en.wikipedia.org/wiki/Runge%E2%80%93Kutta%E2%80%93Fehlberg_method|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.

=== Stiff systems ===

Although the implemented explicit solvers are general and flexible, they are not applicable to stiff systems. [[http://en.wikipedia.org/wiki/Stiff_equation|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.
===== Lattice =====

==== Lattice structure ====
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

==== Boundary conditions ====

See description in [[documentation:faq#Modeling|FAQ]]

==== Membrane Properties ====
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)


===== Diffusion =====

=== Method of lines ===

Partial differential equations (PDEs) are numerically approximated by separating the reaction and diffusion steps using the [[http://www.scholarpedia.org/article/Method_of_lines|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.

=== Diffusion equation ===

For diffusion, Morpheus uses the simple and general forward Euler scheme. The time step $\delta t$ is automatically adjusted according to the [[http://en.wikipedia.org/wiki/Courant%E2%80%93Friedrichs%E2%80%93Lewy_condition|CFL condition]] to guarantuee numerical stability. Alignment of multidimensional lattices to 1D memory (using [[http://www.cplusplus.com/reference/valarray/|valarrays]]) results in highly efficient implementation of diffusion.  

The unconditionally stable method of [[http://en.wikipedia.org/wiki/Alternating_direction_implicit_method|alternate direction implicit (ADI)]] is also implemented (for constant (Dirichlet) boundary conditions), but is currently not used.

===== Motility =====

=== Cellular Potts model ===

Cell shape and motility is implemented according to the [[http://en.wikipedia.org/wiki/Cellular_Potts_model|cellular Potts model]] (CPM)[[http://dx.doi.org/10.1103%2FPhysRevLett.69.2013|(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$.
 
=== Monte Carlo ===

The CPM is a [[http://en.wikipedia.org/wiki/Monte_Carlo_method|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.

=== Metropolis kinetics ===

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.

===== Random numbers =====

By default, Morpheus uses the [[http://en.wikipedia.org/wiki/Mersenne_twister|Mersenne Twister 19937]] pseudo-random number generator included in the Cpp GNU compiler (TR1). Alternatively, the [[http://www.boost.org/doc/libs/1_53_0/doc/html/boost_random.html|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'').




===== Plugins =====

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.
=== Plugin interfaces === 
[{{:documentation:tech_specs:plugin1.png?direct&400|Plugin interface provides interfaces for various types of plugins.}}]


=== Header file (MyPlugin.h) === 
[{{:documentation:tech_specs:plugin2.png?direct&400|example header file}}]

=== Implementation (MyPlugin.cpp) === 
[{{:documentation:tech_specs:plugin3.png?direct&400|example implementation file}}]

=== XML Schema (MyPlugin.xsd) === 
[{{:documentation:tech_specs:plugin4.png?direct&400|example XML schema file}}]