Simulation Frameworks for Morphogenetic Problems

Abstract

Morphogenetic modelling and simulation help to understand the processes by which the form and shapes of organs (organogenesis) and organisms (embryogenesis) emerge. This requires two mutually coupled entities: the biomolecular signalling network and the tissue. Whereas the modelling of the signalling has been discussed and used in a multitude of works, the realistic modelling of the tissue has only started on a larger scale in the last decade. Here, common tissue modelling techniques are reviewed. Besides the continuum approach, the principles and main applications of the spheroid, vertex, Cellular Potts, Immersed Boundary and Subcellular Element models are discussed in detail. In recent years, many software frameworks, implementing the aforementioned methods, have been developed. The most widely used frameworks and modelling markup languages and standards are presented.

1. Introduction

Morphogenesis is the carefully orchestrated process by which the shapes and structures of organs and organisms form, such that they can reliably fulfill their function. Research in developmental biology cumulated in myriads of genetic, molecular biological and biophysical experiments in various model organisms over the last decades. However, many important core mechanisms, which may be based on self-organization, are still not understood.

Complementary to experimental approaches, mathematical modelling and computer simulations can be used to answer morphogenetic questions on a systems biology level. Simulations allow researchers to study a beforehand formulated hypothesis of a governing core mechanism in a very controlled way. Rather than just describing and characterizing the existing knowledge using phenomenological models, mechanistic modelling aims at creating explanatory models to gain confidence into the underlying mechanisms [1]. The emerging phenomena on the tissue or organ level, which are accessible to validation, act as strong indications for the validity or falsity of the hypothesis.

At the coarsest level, the tissue is described by a macroscopic continuous material. For instance, the tissue might be described by an incompressible viscous fluid [2]. However, many cellular processes such as cell migration and adhesion, cell polarity, directed division, monolayer structures and differentiation cannot be cast into a continuous formulation in a straightforward way. Although it has been shown for particular cases, e.g., anisotropic cell division [3], this exercise is not always possible without incisive simplifications and assumptions. For instance, Osborne et al. [4] compared two cell-based models (spheroid and vertex model) to a continuous model of healthy and cancerous crypt renewal. Although all models recover the fundamental behaviour, specific differences exist. The cell-based models require averaging over many simulation runs to account for stochasticity. On the other hand, stochasticity leads to phenomena (such as the monoclonality hypothesis in crypts) which cannot be predicted by deterministic continuous models. Cell-based models usually suffer from a larger number of parameters as compared to continuous models. Therefore, Murray et al. [5] derived a 1D continuous model from cell-based models, and found good agreement. Instead of either using cell-based or continuous approaches, they vote for the derivation of continuous formulations based on discrete models to understand a problem on different scales.

In this review, we discuss a selection of frequently and increasingly used tissue modelling approaches (illustrated in Figure 1), as well as software frameworks. In Figure 2, the discussed modelling techniques are related to the biological scales which they are able to represent, as well as to the computational complexity. At the highest level of detailism, molecules and their interactions can be simulated using molecular dynamics. The accessible length and time scales, however, are small due to the computational complexity. Continuous signalling models (such as deterministic reaction-diffusion equations and stochastic partial differential equations (SPDE’s)) allow to simulate large biochemical interaction networks. Both single cells and entire tissues can be modelled using continuous formulations. Resolving cells, and possibly subcellular structures explicitly, requires complex meshing and moving boundaries, thus leading to relatively high computational complexity. To simulate entire organs, however, the continuous approach has the lowest computational complexity. Currently, the vertex and spheroid models use the cell as the smallest resolved unit. In contrast to the vertex model, the spheroid model does not resolve the cell geometry explicitly and is therefore the most efficient cell-based approach among the presented models. The Cellular Potts, Subcellular Element and Immersed Boundary models are capable to explicitly represent subcellular structures (e.g., the nucleus), and they resolve arbitrarily deformable cell shapes explicitly. The high computational complexity limits these approaches to smaller tissues. These methods allow for coarse-graining, i.e., not every individual cells, but multiple cells are associated with a computational entity. In the case of the Immersed Boundary method, an entire tissue/organ can be represented by an Immersed Boundary “cell” [11], thus closing the loop to the continuous approach.

The organization is as follows. In the Section 2, we introduce the basic theory of the most widely used models for morphogenetic simulations, that is, the vertex model, the Subcellular Element model, the Cellular Potts model and the spheroid model, and highlight application examples. Existing major software frameworks are reviewed in Section 3. In Section 4, we conclude by comparing the modelling approaches and suggesting future directions.

2. Tissue Models

2.1. Continuum Tissue Models

At the highest level of abstraction, the tissue is modelled as a continuous material. The tissue can, amongst others, be described as a (visco-) elastic material, a viscous fluid, or reaction-diffusion processes. In the following, these approaches are discussed in the light of morphogenetictissue modelling.

A continuum theory for soft incompressible hyperelastic tissue, similar to the description of nonlinear elastic materials, has been developed [12], extended [13] and applied to various morphogenetic problems. To account for finite volumetric (anisotropic) tissue growth and remodelling, the total deformation gradient tensor F = F^★ · G is split into two tensors, namely the (passive) elastic deformation gradient tensor F^★, and the (active) growth tensor G. In order to generalize the theory of finite volumetric growth, it is necessary to find a growth tensor G to model the tissue properties. The dependency of the growth tensor G on the stress tensor σ is called “morphogenetic law” and is a constitutive law depending on the tissue of interest [14]. One possible hypothesis is Beloussov’s hyper-restoration hypothesis, which assumes that cells respond (and overshoot) to local stress changes [15]. To solve the equations, Finite Element methods are used.

The finite volumetric growth theory has been used to study many morphogenetic processes (reviewed in [16]), including stress-dependent bone growth [12], the gastrulation process in sea urchin embryos [14], drosophila ventral furrow formation [17,18], avian head fold formation [19], brain morphology development [20], initialization of budding processes in branching morphogenesis [21], and bending morphogenesis of the primitive heart tube [22].

To study mechanical cell-substrate interactions, the Oster-Murray-Harris model assumes that cells, which are described by a scalar cell concentration, exert active forces generated by the actomyosin and move within an extracellular matrix (ECM). The latter is modelled as a linear viscoelastic medium. By modelling the motion of the cells being biased into the direction of the highest strain, the cells’ forces deform the substrate, and cells many cell diameters away are influenced. This mechanism has the capability to form patterns in a self-organizing fashion, thus being an interesting biomechanical mechanisms for morphogenesis, e.g., mesenchymal cells embedded in an ECM [23].

Oftentimes, embryonic tissue can be modelled as a viscous fluid when being interested in the long timescale behaviour [24,25]. Furthermore, it has been shown that high proliferative activity leads to a fluidization of the tissue [26]. Dillon et al. [2] used the Navier-Stokes equation to describe the embryonic limb as fluidic tissue:

$$\rho ({\partial }_t\mathit{u}+(\nabla \cdot \mathit{u})\mathit{u})=-{\nabla }_p+\mu \left(\mathrm{\Delta }\mathit{u}+\frac{1}{3}\nabla (\nabla \cdot \mathit{u})\right)+\mathit{f}$$

(1a)

$$\rho \nabla \cdot \mathit{u}=\mathcal{S}$$

(1b)

where u denotes the velocity field, p the pressure field, and f an external force field. The local mass source $\mathcal{S}$ can be used to model proliferation, which may be controlled by biomolecular signalling. Usually, the Reynolds number $\mathit{Re}=\frac{UL}{\nu }$, with U being a characteristic velocity, L a characteristic length scale, and ν the kinematic viscosity, is very small. Assuming L = 10⁻³ [m], U = 10⁻⁸ [m/s] and ν = 10¹ … 10² [m²/s] [25], then Re = 10⁻¹³ … 10⁻¹² can be estimated. By non-dimensionalizing Equation (1), and assuming Re → 0, the Stokes equations can be written as:

$$0=-\nabla p+\mathrm{\Delta }\mathit{u}+\frac{1}{3}\nabla (\nabla \cdot \mathit{u})$$

(2a)

$$\nabla \cdot \mathit{u}=\mathcal{S}$$

(2b)

where all variables (and differential operators) now denote dimensionless quantities. Instead of computing the velocity field based on a tissue model, it can be acquired from experimental data. To this end, tissue shapes are taken at regular developmental stages using microscopy and subsequent image processing. Morphing and interpolation are needed for the continuity of the velocity fields [27,28].

The movement of chemotactic cells of the slime mold (Acrasidae, e.g., Dictyostelium) in a chemical gradient has been described by Keller and Segel [29]. In their original work, they derived a system of reaction-diffusion partial differential equations (PDE’s) for two chemical species, their complex, and a mass density of the amoebae. The interplay between diffusion, complex formation and chemotaxis leads to complex pattern formation, as observed in nature.

If cell mechanics is ignored, tissue growth can be achieved by formulating laws how to move the tissue boundaries as a function of the biomolecular signalling. This approach has been used to study the ability of a Turing-type signalling mechanism to control the branching of embryonic lungs [6]; an example is given in Figure 1A.

2.2. Spheroid Models

A wide range of models have been developed which model cells as particle-like (colloidal) objects (cf. Figure 1B). The cells are assumed to have spherical, isotropic shape, which is represented by soft sphere interaction potentials. It has been shown that the Johnson-Kendall-Roberts (JKR) theory for adhesive spheres can be applied to cells [30]. The JKR model features a hysteresis effect, i.e., the distance at which cells lose contact when being pulled apart is larger than the distance at which they formed contact when having approached each other earlier [31]. Alternatively, Hertz contact models [32], harmonic interaction potentials [33] and dashpot-spring elements [34] have been applied. Because the used interaction potentials are very diverse and purely phenomenological, the interested reader may find details elsewhere [32–35].

To move the cells, the over-damped (i.e., acceleration and thus inertial effects are ignored) equations of motion can be solved deterministically [36]. Drasdo et al. [32,33] used a Metropolis algorithm to solve the friction dominated stochastic (Langevin) equations since cells perform random walks when being isolated. The Langevin equation for a cell i with velocity u_i at position x_i is given as [32]:

$$\left(\mathrm{\gamma }+{\underset{¯}{\underset{¯}{\mathrm{\Gamma }}}}_{is}^f\right){\mathit{u}}^i+{\displaystyle \sum _{j\ne i}{\underset{¯}{\underset{¯}{\mathrm{\Gamma }}}}_{ij}^f({\mathit{u}}^i-{\mathit{u}}^j)={\displaystyle \sum _{j=0}^N{\mathit{F}}^{ij}+{\mathit{\eta }}^i(t)}}$$

(3)

where γ is an effective friction coefficient. The tensor ${\underset{¯}{\underset{¯}{\mathrm{\Gamma }}}}_{iy}^f={\mathrm{\gamma }}_{||}^{(iy)}{\mathit{n}}_{iy}\otimes {\mathit{n}}_{iy}+{\mathrm{\gamma }}_{\perp }^{(iy)}(\underset{¯}{\underset{¯}{E}}-{\mathit{n}}_{iy}\otimes {\mathit{n}}_{iy})$ (with y = s for substrate- and y = j for cell-cell interactions) is composed of a parallel and perpendicular contribution with friction coefficients ${\mathrm{\gamma }}_{||}^{(ij)}$ and ${\mathrm{\gamma }}_{\perp }^{(ij)}$, respectively. ${\mathit{n}}_{iy}=\frac{{\mathit{x}}_y-{\mathit{x}}_i}{|{\mathit{x}}_y-{\mathit{x}}_i|}$ is the unit vector from cell i to entity y. The repulsive and attractive forces F^ij are derived from the interaction potentials, and ηⁱ (t) is a noise term.

Various implementations for cell cycle, cell growth, division, migration and signalling have been developed, e.g., [32,37–39]. In order to introduce geometrical anisotropy and cell polarity, the cells can be represented by ellipsoids (or ellipses in 2D). Dallon, Othmer and Palsson [36,40] applied an ellipsoid cell model to study the chemotactical behaviour of Dictyostelium discoideum (slime mold) cells.

The spheroid model origins from a spheroid model variant where each cell is represented by a computational particle, but explicit cell boundaries are represented by Voronoi tesselation (also known as Dirichlet domains) [41]. It has been used to study pattern formation of the epidermis in butterfly [42], and intestinal crypt dynamics [43]. Gevertz et al. [44] used a spheroid-Voronoi model in combination with a reaction-diffusion based signalling model to study the interplay between vascular endothelial growth factor (VEGF) signalling, blood vessel formation and tumorogenesis. The epithelial renewal mechanisms in intestinal crypts have been studied with a spheroid model, where the cells are interconnected by springs, but the cell shape represented by Voronoi tesselation [45]. Finally, hybrid (multi-scale) models, which combine the cell representation of the spheroid model and the efficiency of continuous tissue description have been proposed [46,47]. Panthamanathan et al. [48] show how the mechanical properties relate to the model parameters both in the spheroid and vertex model.

Due to the simple geometric representation, spheroid models are predestined for problems with large numbers of cells where cellular and subcellular details (such as arbitrarily deformable shapes) can be omitted or modelled into the interaction potentials. More detailed tissue modelling approaches are discussed in the following.

2.3. Vertex Models

In two dimensions, the vertex model approximates each cell by a polygon. The interface between each pair of neighbouring cells is represented by a common edge, and the points of intersection by vertices. Since two neighbouring cells share a common edge which means that gaps (interstitial space) cannot be natively represented. On the other hand, undesired overlaps are avoided. Hence the vertex model is most appropriate to represent densely packed epithelial tissue. The vertex model originates from the work introduced by Honda and co-workers [49]. They use Voronoi tesselation to construct polyhedral geometries around cell centers. Extensive reviews can be found in [50,51]. Different approaches have been developed to move the vertices; however, they have in common that the forces solely act on the vertices itself. The force F_i acting on a vertex i can be formulated explicitly (e.g., [52,53]), or by using an energy potential E [8,54]:

$${\mathit{F}}_i=\frac{\partial E}{\partial {\mathit{R}}_i}$$

(4)

with R_i being a junctional direction of vertex i. A typical energy function includes area elasticity, line tension and contractility of the cell perimeter [8]:

$$E({\mathit{R}}_i)={\displaystyle \sum _{\alpha }\frac{K_{\alpha }}{2}{(A_{\alpha }-A_{\alpha }^0)}^2+}{\displaystyle \sum _{<i,j>}{\Lambda }_{ij}{l}_{ij}}+{\displaystyle \sum _{\alpha }\frac{{\mathrm{\Gamma }}_{\alpha }}{2}}{L}_{\alpha }^2$$

(5)

The first term on the right hand side describes the area elasticity, with K_α being the area elasticity coefficient, A_α the current area and $A_{\alpha }^0$ the resting area. The second term denotes line tension, which, for vertex i, considers all neighbouring vertices j. Λ_ij is the line tension coefficient, and l_ij the edge length, and <i, j> denotes the summation over all edges. The last term represents the perimeter contractility, with Γ_α being the cell perimeter contractility coefficient and L_α the cell’s perimeter.

Once the forces are defined, different approaches have been used to move the vertices. Either, at each time step, the entire morphology is relaxed to its local steady state, e.g., using a conjugate gradient method [8]. In this case, the vertex model describes stable and stationary topologies, because force balance is assumed, and the local minimum of the energy function is computed. Or, in the simplest deterministic way, the cells are moved according to the overdamped Newtonian mechanics, where inertia can be neglected [54]:

$$\eta \frac{d{\mathit{x}}_i}{dt}={\mathit{F}}_i$$

(6)

where x_i denotes the position of vertex i, η the mobility coefficient, and F_i the force acting on vertex i. As a third approach, a Monte Carlo algorithm to shorten the cell boundaries and thus relaxing the tissue (termed Boundary Shortening procedure in [41]) can be applied.

More related models can be found for plant biology. Brodland et al. [55,56] use a Finite Element formulation. Alternatively, an explicit spring network can be relaxed to equilibrium at each time step [57]. A similar formulation is used in the open source platform VirtualLeaf [58], but there a Monte Carlo energy minimization approach, similar to the Potts models [59], is used for relaxation.

Different biologically motivated processes can alter the connectivity of the topology. An number of basic mesh restructuring processes have been identified (nicely summarized in [50]): cells can change neighbours (T1 transition), delaminate or die (element removal, T2 transition), or clot together (element intersection, T3 transition). Furthermore, node switching has to be applied to prevent self-intersection. Cell division is modelled by splitting the mother cell and introducing a new edge.

Several vertex models have been combined with biomolecular signalling. Schilling et al. [60] used the cell polygons to solve the partial differential equations, representing reaction-diffusion systems, in a Finite Volume framework. However, the chemical species can only diffuse freely across the entire tissue, or are bound to cells without diffusion, thus being described by ordinary differential equations (ODE’s). Advection has been ignored. A similar approach has been taken by Smith et al. [61], where the polygonal elements serve as a mesh for a finite element method.

Although the vast majority of models have been formulated in two spatial dimensions, 3D versions have been proposed. In the simplest case, the vertices are allowed to move in all spatial dimensions. This approach has been used to study buckling behaviour of an infinitely thin single layer tissue [62]. Honda and co-workers [63,64] used Voronoi tesselation to represent polyhedra, the three-dimensional counterpart of the vertex model’s polygons. In recent work, viscosity has been added to the model to study epithelial vesicle growth in 3D [65].

The vertex model has been used to study numerous morphogenetic problems, such as to explore the hypothesis of a mechano-sensing mechanism on Drosophila wing disc growth control [66]. Brodland et al. [67] developed a model variant and achieved realistic cell rearrangement and sorting behaviour as observed in embryonic epithelia. However, they find that cells are artifactually stiffening when cell edges shorten towards zero length. Also in the wing disc, Schilling et al. [60] studied the effect of signalling on the mechanical tension between compartments of different cells, namely cell sorting. Farhadifar et al. [8] analyzed the topological transition of the irregular polygons of the Drosophila wing disc into regular, largely hexagonal patterns (cf. Figure 1C). They compare simulated network morphologies to experimental data, and they found parametric spaces which lead to good agreement. Furthermore, they estimated the parameter values from laser ablation experiments.

In general, for many vertex model variants, it is unclear how they can be extended to include the mechanical effects of the cytoskeleton (e.g., the viscous behaviour). The Finite Element formulation of Brodland et al. [67] is an interesting approach to combine cell-based vertex models with a continuous description of the viscoelastic material.

2.4. Cellular Potts Model

The Cellular Potts model (CPM), first introduced by Graner and Glazier [59,68], originates from the Ising model [69] and on its later generalization, the Potts model [70], which show phase transition phenomena in ferromagnetic materials. The CPM is also known—according to its founders and main contributors—as the Glazier-Graner-Hogeweg (GGH) model [71]. The key element of the CPM is the lattice. Each lattice site x carries a spin σ (x) ∈ ℤ ^+,0, where each spin value represents a cell identity (with the exception of σ = 0, which usually denotes extracellular space). The central element of the CPM, similar to the Ising model, is the Metropolis algorithm [72]. The update algorithm seeks to minimize a Hamiltonian energy function H that is defined over the entire domain. In the original form [59,68], the Hamiltonian H is composed of two contributions, a volume constriction term H_v and a cell-cell adhesion term H_a:

$$H=H_v+H_a={\displaystyle \sum _{\mathrm{\sigma }}{\mathrm{\lambda }}_v}{(V_{\mathrm{\sigma }}-V_{\mathrm{\sigma }}^T)}^2+{\displaystyle \sum _{(\mathit{x},{\mathit{x}}^{\prime })}J(\tau (\mathrm{\sigma }(\mathit{x})),\tau (\mathrm{\sigma }({\mathit{x}}^{\prime })))}\cdot (1-\delta (\mathrm{\sigma }(\mathit{x}),\mathrm{\sigma }({\mathit{x}}^{\prime })))$$

(7)

where V_σ and $V_{\mathrm{\sigma }}^T$ are the actual and target volume of cell σ, respectively. The coefficient λ_v controls the energy penalization. The term J (τ, τ′) is the surface energy between two cell types τ and τ′. Here, τ (σ (x)) denotes the cell type of the cell σ at position x. The inverse Kronecker delta 1 − δ (·, ·) equals one if σ ≠ σ′, which is true if the two neighbouring lattice sites are of different cell type. The summation ${\displaystyle {\sum }_{(\mathit{x},{\mathit{x}}^{\prime })}}$ runs over all lattice sites x and its neighbouring lattice sites x′. To summarize, the cell-cell adhesion contribution is only non-zero across cell-cell interfaces, and only if the two cells are of different cell types. The modified Metropolis algorithm is implemented as follows: a lattice site x and a neighbouring lattice site x′ are chosen randomly. The probability of converting the spin σ (x′) of the lattice site x′ into the spin value σ (x) is set to be:

$$p(\mathrm{\sigma }(\mathit{x})\to {\mathit{x}}^{\prime })=\{\begin{array}{ll}1& \mathrm{if}\mathrm{\Delta }H(\mathrm{\sigma }(\mathit{x})\to {\mathit{x}}^{\prime })<0\\ \exp (-\mathrm{\Delta }H(\mathrm{\sigma }(\mathit{x})\to {\mathit{x}}^{\prime })/T)& \text{otherwise}\end{array}$$

(8)

Thus a lattice site conversion is accepted if it decreases the global Hamiltonian. An increased Hamiltonian is accepted with a Boltzmann-distributed probability, which exponentially decreases with an increasing energy jump. If the lattice consists of N lattice sites, then a series of N spin-copy attempts is called a Monte Carlo step (MCS). The factor T corresponds to the “temperature” in the antecedent Ising model. Similar to the Ising model, the CPM may show phase transition behaviour as a result of changing temperature. Thus, restrictions on the temperature have been proposed [71] to guarantee domain connectivity. In general, it is difficult to interpret the biophysical interpretation of the “temperature” [71,73].

Since the entire computational domain has to be finely resolved by the lattice (and not only the cells of interest), the memory requirements and computational costs may be high. However, as for many lattice-based algorithms, high-performance-computing and parallelization techniques can be applied efficiently [74]. As it is the case for every inherently stochastic approach, multiple repetitions are required for statistical reasons, thus increasing the computing time.

The CPM has been coupled to reaction-diffusion equations by applying Finite Difference discretization and explicit time integration [75–77], and has been applied to study blood vessel formation, where it was found that cell elongation is required to explain the vascular network formation [78], and where a vessel branching mechanism was shown to reproduce experimental observations [9] (cf. Figure 1D). Szabo et al. [79] studied the tip growth of blood vessel sprouts, and Bauer et al. [80] the inducing effect of tumors on angiogenesis. Furthermore, the CPM was used to simulate cell sorting [59,68], kidney branching [81], somitogenesis [82] and chicken limb morphogenesis [83].

2.5. Immersed Boundary Cell Model

In the 2D IBCell (Immersed Boundary Cell) model, the shapes (i.e., the membranes) of the cells are approximated by finely resolved polygons [84–86]. However, in contrast to the vertex model, cells do not share edges, but each cell rather owns its own edge in the IBCell model. An illustration is given in Figure 1E. These polygons are immersed in a fluid, and, in a first step, the vertices are passively advected by the fluid velocity field u. Hence, the membrane is impermeable to the fluid. In other words, the curvilinear coordinates of the boundary points (q, r, s) are attached to a material fluid particle, and the Eulerian coordinates of that particle at time t are denoted as X (q, r, s, t). The material derivative $\frac{D\mathit{u}}{Dt}(\mathit{x},t)=\frac{{\partial }^2\mathit{X}}{\partial t^2}(q,r,s,t)$ is the acceleration of whatever material point is at position x at time t. On the cell boundary, various force generating processes F can be modeled, including membrane tension, cell-cell and cell-substrate junctions, and forces exerted by the cytoskeleton. For example, membrane tension can be modelled as a set of Hookean springs connecting each pair of vertices of the polygon. The cellular geometries thus represent the elastic material properties of the tissue. The fluid, in which the polygons describing the cell shapes are immersed in, is divided into compartments: the fluid inside the cells represents the viscous properties of the cytoplasm, whereas the fluid in the inter-cellular space represents the viscous properties of the extracellular matrix and interstitial fluid. The fluid is, as already given in Equation (1), described by the Navier-Stokes equation. The local mass source $\mathcal{S}$ may be used to model cell growth. In contrast to the Lagrangian polygons, the fluid equations are solved on an Eulerian mesh.

To solve the mechanical interaction between the elastic structures and the fluid, the Immersed Boundary method is used [87]. The conversion from Lagrangian variables (e.g., force density F (q, r, s, t)) to their Eulerian counterparts is executed by integrating over the material coordinates:

$$\mathit{f}(\mathit{x},t)=\int \mathit{F}(q,r,s,t)\mathrm{\delta }(\mathit{x}-\mathit{X}(q,r,s,t))dqdrds$$

(9)

where δ (·) denotes the delta Dirac function. For the numerical implementation, the delta Dirac function is approximated by a kernel function in such a way that it covers multiple grid points on which the conversions in Equations (9) and (10) are evaluated. Examples of widely used kernel functions can be found in [88]. The equation of motion of the Lagrangian particles reads:

$$\frac{\partial \mathit{X}}{\partial t}(q,r,s,t)=\mathit{u}(\mathit{X}(q,r,s,t),t)=\int \mathit{u}(\mathit{x},t)\mathrm{\delta }(\mathit{x}-\mathit{X}(q,r,s,t))d\mathit{x}$$

(10)

where the same kernel function is used for the interpolation of the velocity field to the Lagrangian vertices of the polygon.

By iteratively solving the fluid equations, interpolating the velocity field to the cell geometries, moving the cells, recomputing the forces acting on the cell geometries, distributing the forces to the local fluid neighbourhood and restarting the process, fully coupled fluid-structure interaction is achieved and fluid-structure interaction, an inherently difficult problem, can be handled.

The IBCell model has been independently developed by Rejniak [85] and Dillon et al. [86]. While the model of Rejniak uses a single membrane with a single mass source in the middle of the cell and compensating mass sinks distributed in the interstitial neighbourhood of the membrane, the model of Dillon et al., has two interconnected closed polygons to represent membrane of finite thickness and thus non-zero bending stiffness, and cell growth is mimicked by water channels along the membrane. The different versions of the model have been applied to ductal tumorogenesis, and no fundamental difference of the results was found [89].

As compared to the Subcellular Element model and the vertex model, the mechanical processes (elasticity and viscosity) are explicitly represented and can be directly associated with experimental measurements. Although the membrane and cell-cell junction forces are phenomenological (like the potentials in the spheroid and Subcellular Element model), the IBCell model’s explicit representation of the cell membranes, the interstitial space and the cell-cell junctions allows to study cell intercalations and topological changes in a much more intuitive way than in the vertex model, and cytoskeletal viscosity is natively included. Additionally, the polygonal shapes of cells in a densely packed tissue is an emergent property rather than a model assumption (vertex model).

A drawback of the IBCell model is, however, its inherent computational complexity. Even for 2D models, the need for Navier-Stokes solvers, the finely discretized cell shape and the iterative Immersed Boundary method currently limit the number of cells to a few thousand. In contrast to the Subcellular Element model and the the vertex model, the generalization of the IBCell model into 3D is an expensive exercise, requiring complex triangulation of the cell surface. 3D simulations similar to the IBCell model to study the membrane of red blood cells in blood flow are limited to a few cells [90].

On the other hand, the IBCell model could be used to represent clusters of cells or even entire tissues instead of individual cells, thus closing the gap to the continuous approach to model tissue as fluid (cf. Section 2.1). Dillon et al. [11] modeled the limb bud tissue as a viscous fluid, and used the Immersed Boundary method to model the tissue boundary.

2.6. Subcellular Element Models

The Subcellular Element model (SEM), first presented by Newman [91], assumes that the cytoskeleton of a cell is divided into a number of subcellular elements. Each element is mathematically and computationally represented by a point particle, and biophysically, they interact via forces derived from interaction potentials. In the simplest version, two potentials have to be specified: an intracellular interaction potential V_intra for the interaction of elements within an individual cell, and an intercellular interaction potential V_extra for the interaction of elements between different cells. Besides representing the cytoskeleton, dedicated elements can represent the extracellular matrix, the cell membrane, the nucleus, or other organelles or cellular structures. To that end, the modeller has to model the interaction potentials between the element types accordingly. Since the Lagrangian elements move freely in space, there is no need for an underlying grid. Neglecting inertia effects, the equation of motion for the position ${\mathit{y}}_{{\alpha }_i}$ of a subcellular element α_i of cell i reads [91]:

$$\eta \frac{\partial {\mathit{y}}_{{\alpha }_i}}{\partial t}={\zeta }_{{\alpha }_i}-{\mathrm{\Delta }}_{{\alpha }_i}{\displaystyle \sum _{{\beta }_i\ne {\alpha }_i}{V}_{\mathrm{int}ra}}(|{\mathit{y}}_{{\alpha }_i}-{\mathit{y}}_{{\mathrm{\beta }}_i}|)-{\nabla }_{{\alpha }_i}{\displaystyle \sum _{j\ne i}{\displaystyle \sum _{{\mathrm{\beta }}_j}{V}_{\mathrm{int}er}(|{\mathit{y}}_{{\alpha }_i}-{\mathit{y}}_{{\mathrm{\beta }}_i}|)}}$$

(11)

where ${\zeta }_{{\alpha }_i}$ is Gaussian noise, and η the viscous damping coefficient. The first term represents intra-cellular interactions, and the summation runs over all remaining elements β_i of cell i. The inter-cellular interaction described by the second term takes all pair-interactions between the elements β_j of other cells j and the elements α_i of cell i into account. Numerically, the forward Euler and a two-step Runge-Kutta scheme have been used to integrate Equation (11) [92].

Similar to the inter-atomic interaction potentials in molecular dynamics simulations, the phenomenological potentials V_intra, V_inter are chosen such that the subcellular elements prefer to be in an equilibrium distance to their neighbouring elements, i.e., that the subcellular elements repel each other when the distance between two cells is low, and attract each other when the distance is large. However, the attractive force rapidly decreases to zero at large distances. The authors of [91] used a Morse potential for both the inter- and intracellular potentials V_intra, V_inter:

$$V(r)=U_0\exp (-r/{\xi }_1)-V_0\exp (-r/{\xi }_2)$$

(12)

where r denotes the distance between two subcellular elements. The energy scale parameters U₀, V₀ and the length scale parameters ξ₁, ξ₂ control the magnitude and shape of the potential, as well as the location of the equilibrium distance.

It has been shown that the simplest implementations of the SEM successfully recover biophysical cell behaviour on intermediate time scales and for low strains. An illustration of a stretching simulation is given in Figure 1F. The high-frequency response of the cytoskeleton, as well as the long time and high strain behaviour could not be recovered. An extension has been proposed to capture active cytoskeleton restructuring (active cytoskeleton polymerization and depolimerization of actin filaments and microtubules) [93]. This is realized by adding and removing subcellular elements in a smooth way in order to polarize the cell. In this way, cell polarization and migration could be realized.

Using mean field theory, it has been shown that the coarse-grained mechanical behaviour can be written as an analytically derived partial differential equation [91], and that the macroscopic properties do not depend on the number of discretizing elements [94]. Indeed, it has been found that macroscopic tissue viscosity is an emerging property of cellular processes [10].

Due to the explicit representation of arbitrary cell shapes, the SEM requires a high number of computational particles to represent a cell, making it a computationally expensive method as compared to the 3D vertex model or the 3D spheroid model. However, due to its algorithmic simplicity borrowed from molecular dynamics, efficient library implementations can be used [92,93]. It has been reported that the SEM is suitable to study problems with up to 10,000 cells [91].

To date, the SEM has not been widely used to study morpghogenetic or related problems. To study the dynamics of layers of epidermal cells, the SEM has been efficiently implemented for graphics processing units (GPU’s) and coupled to an ODE-based model of Notch-signalling [92]. The authors of [93] improved the cytoskeleton remodelling algorithms, and studied cell migration, growth and proliferation, and coupled the model to simplified juxtacrine signalling models. Cell motility and invasion of tissue, which are for instance important processes in primitive streak formation, have been demonstrated by Sandersius et al. [10].

3. Software Frameworks and Standards

3.1. Chaste

The original aim of Chaste ( http://www.cs.ox.ac.uk/chaste/) [95,96] was to create a comprehensive simulation environment for physiological simulations (cardio electrophysiology), and later for a wider range of computational biology. Chaste, which is developed by a large team from the University of Oxford (Oxford, UK), makes extensive use of modern software engineering methods such as agile and test-driven software development. The framework offers core functionality for meshing, solving ODE’s, PDE’s and I/O. Furthermore, it makes use of highly parallelized libraries such as PETSc (Portable, Extensible Toolkit for Scientific Computation, http://www.mcs.anl.gov/petsc/) to provide parallelized data structures and (FEM) solvers.

Cell-based Chaste denotes a collection of cell-based tissue model implementations. This includes Cellular Automata models (as reviewed in [97]), Cellular Potts models (see Section 2.4), spheroid models (cf. Section 2.2) and vertex models (cf. Section 2.3, [50]). Currently, cell-based Chaste does not offer the possibility to use the markup languages CellML (Cell Markup Language) [98] or SBML (Systems Biology Markup Language, cf. Section 3.6).

In recent years, cell-based Chaste has been used to study several biological problems, most notably the renewal dynamics of intestinal crypts in the gut. Van Leeuwen et al. [45] developed a spheroid-Voronoi model [43] (cf. Section 2.2) in combination with a signalling and cell cycle model. They perform virtual experiments to predict the mechanisms and conditions of clonal expansion and niche succession. Furthermore, the vertex model has been used to simulate healthy and mutageneous crypt dynamics. Fletcher et al. [99] use a vertex model to study different hypotheses on monoclonal conversion, and find that the stem cell niche hypothesis is more favourable than the immortal stem cell hypothesis. To take mechanical properties, such as cell adhesion, into account to simulate the behaviour of mutageneous cells in the crypt, particle-Voronoi models [100] and spheroid models [101] have been used. Dunn et al. [102,103] study the mechanical interplay between the basement membrane and an epithelial monolayer using a particle-Voronoi cell model. In particular, they study the homeostasis between cell migration and death. Figueredo et al. [104] developed an on-lattice agent-based model which incorporates a multitude of cell behaviour such as proliferation, cell death and motility. Furthermore, they coupled it to cell-cycle and reaction-diffusion-signalling models.

3.2. CompuCell3D

The open-source modelling framework CompuCell3D ( http://www.compucell3d.org/), whose development is mainly driven by researchers from the Indiana University and led by James Glazier, implements the Cellular Potts (Glazier-Graner-Hogeweg) model [74]. The framework is highly modular, and customized modules—written either in C++ or Python—can be developed and added as plugins. This way, the user can add, for instance, own Hamiltonian energy functions as given in Equation (7). The core is parallelized using the shared-memory paradigm, thus limiting the software to desktop computers and servers. CompuCell3D uses its own XML-based format to store and load models. To solve PDE-based signalling models, a built-in facility is offered. Most noteworthy, CompuCell3D offers a mature user interfacing infrastructure: a model editor, a graphical tool to draw and configure cells, import and segment experimental microscopy images, and for visualizing and analysing the simulation results. A thorough collection of manuals and tutorials are supplementing the framework.

CompuCell3D has been deployed to study a broad range of developmental and biological problems, and a few of them shall be highlighted here. To explain the regular pattern formation in somitogenesis, Glazier et al. [105] studied a clock based mechanism. A clock- and wavefront-free mechanisms for somitogenesis has been proposed [106]. They propose a self-organizing patterning mechanism based on cell-cell interaction and polarization. Computer simulations using CompuCell3D were indeed successful in achieving self-organizing patterns. Many authors studied angiogenesis and vasculogenesis [107,108], lumen formation during blood vessel formation [109], and the interplay between angiogenesis and tumorogenesis [110]. CompuCell3D has also been used to study tumor and cancer growth [110–112]. In addition, classical morphogenetic problems have been explored, including primitive streak formation in chick embryos [113] and salivary gland branching morphogenesis [114]. Popawski et al. [83] studied the influence of the AER on limb bud shape formation.

3.3. CellSys

The freely available, but closed-source 3D framework CellSys ( http://msysbio.com/) [115], mainly developed by Stefan Hoehme and Dirk Drasdo and co-workers, is built around the spheroid model (as described in Section 2.2). The software, written in C++, offers graphical user interfaces and real-time rendering of the results. CellSys offers already parametrized cell models, cell migration, cell-cycle and signalling models. The authors of CellSys developed the software in the context of liver regeneration. They studied the influence of parameters and mechanisms on liver regeneration when hepatocytes were damaged by drugs [116], and found that, besides hepatocyte proliferation, cell polarization plays a crucial role. The model has been compared with experimental data [7] to obtain a so-called virtual liver that may help to understand the mechanisms of liver regeneration and diseases [117,118].

3.4. EPISIM

The modelling platform EPISIM, which was initially applied to epidermal tissue homoeostasis and which is available as free closed-source software [119,120], is developed by researchers from the Heidelberg University and is an attempt to standardize models in systems biology, in particular multi-scale agent-based tissue models, cell behaviour and signalling models. To this end, it is inherently using the Systems Biology Markup Language SBML (cf. Section 3.6). The platform, offering graphical user interfaces (SBML model importer, signalling and cell behaviour model editor), aims to be easy and intuitive to operate, such that it may be used by biologists without computer and computational science background. The compilation into executable code, and the multi-scale time step algorithms are hidden from the user. However, the simple and entirely graphical usability comes at the price of less flexibility and extensibility. Currently, an ellipsoid (agent-based) and a lattice-based tissue model are implemented. Although EPISIM exploits multicore processors by using multithreading, it is not distributed-memory parallelized. However, the possibility to run parameter screens by launching independent instances on cluster computers is offered [119].

EPISIM has been used by the authors to model epidermal homeostasis in 2D and 3D [121] and wound healing [122]. A range of functionalities such as cell behavioural, cell-cycle, differentiation, signalling and mitosis have been implemented.

3.5. Other Frameworks

The open-source framework VirtualLeaf ( https://code.google.com/p/virtualleaf/) [58] implements a vertex model [54], which is appropriate for plant tissue modelling since it restricts relative cell movements. The Hamiltonian energy function is minimized by a Monte Carlo algorithm. The framework allows for solving ODE signalling models for each cell, as well as PDE based reaction-diffusion signalling models.

The aim of the freely available, but closed-source software Biocellion ( http://biocellion.com) [123] is to allow for high-performance computing on (massively parallelized super-) computers with millions to billions of cells. To that end, a parallelized data structures, domain decomposition, load balancing and PDE solvers are implemented, but hidden from the user. Biocellion offers an agent-based (spheroid) cell model. The authors used the software to reproduce cell sorting and microbial pattern formation.

The free 2D/3D software Morpheus ( http://imc.zih.tu-dresden.de/wiki/morpheus) [124] offers a Cellular Potts tissue model, coupled ODE lattices (for intracellular processes) and deterministic and stochastic reaction-diffusion finite difference solvers. Signalling models can be loaded in SBML format. Morpheus has been applied to pancreatic [125,126] and vascular pattern formation [127,128].

The open-source software OpenCMISS ( http://physiomeproject.org/software/opencmiss) [129] focuses on organ-level modelling by means of continuum tissue models.

Recently, the software framework LBIBCell (Lattice Boltzmann Immersed Boundary Cell), which couples the IBCell model with signalling solvers, has been developed and used to study the behaviour of Turing-type pattern mechanisms on growing cellular tissues [130].

3.6. Mark-Up Languages, Standards and Guidelines

In order to efficiently and unambiguously describe and store systems biology models (as for example in the database BioModels [131]), multiple standards have been developed recently. MIRIAM (Minimum Information Requested in the Annotation of Biochemical Models) [132] and MIASE (Minimum Information about a Simulation Experiment) [133] are guidelines to describe biological models and their simulation execution. The exchange format SED-ML (Simulation Experiment Description Markup Language) [134] uses MIASE to define simulation experiments. The Systems Biology Markup Language SBML ( http://www.sbml.org) [135] is a machine-readable standard to describe biological reactions. The standard supports the description of, amongst others, mathematical functions, species, kinetic parameters, reactions and compartments. To date, more than 260 software packages, including CompuCell3D, make use of SBML. Similarly to SBML, CellML [98] is an XML-based exchange format, but aims at the general description of the structure and mathematics of cells. Apart from compartment (lumped-parameter) models, CellML does currently not support spatial processes. The Cell Behaviour Ontology ( http://cbo.biocomplexity.indiana.edu/) [136] describes (spatial) physical and biological processes of cells and tissues. It is, however, not yet realised in SBML. The standard language (ontology) BioPAX (Biological Pathway Exchange) [137] focuses more on reaction network visualization and qualitative analysis. Conversion tools, such as from BioPAX to SBML, are available [138]. FieldML (Field Markup Language) [139] is a general format for spatial models including partial differential equations.

4. Discussion and Future Directions

The choice of an appropriate tissue model depends on many factors: the cellular processes which can be modelled, computational cost and 3D capability, and many more. In the following, we discuss these aspects for the presented tissue models.

The coarsest of the presented cell-based models, the spheroid model, is a good choice when the specific model requires discrete cells, but sophisticated biomechanics can be ignored. This may be the case for bulk tissue, 2D cell cultures on substrates, or epithelia being mechanically connected to a basal membrane. For instance, cell growth and division in random direction can be well reproduced, but anisotropic cell division cannot because the shapes of the individual cells impact the behaviour. As a generalization of the spheroid model, the Subcellular Element model allows for arbitrarily deformable shapes. It suffers, however from the lack of hydrodynamic interactions of the elements [94], which may be improved in the future—it will be interesting to see more biological applications of this model. The vertex model is particularly suitable to model densely packed 2D epithelia due to the polygonal structure of the apical surface of the cells. Since the forces are usually computed from energy potentials, which are subject to modelling, the vertex model is ideal to study mechanical processes where the membrane is involved, such as interface tension and cell sorting. In most of the formulations, however, the mechanical influence of the cytoskeleton is missing. The Immersed Boundary Cell model represents the viscous behaviour of the cytoskeleton and the interstitial fluid and extracellular matrix. Furthermore, cell-cell junctions are explicitly resolved and their dynamics can be studied in great detail, which renders it an excellent approach to study, e.g., directed cell division, cell intercalation and migration. As the vertex model, the Cellular Potts model is therefore an appropriate framework to study mechanisms and hypotheses which can be formulated on the basis of energy potentials. If the modellers succeed in deriving continuous equations of the cellular processes, then the problem of interest is ideally analysed both on a cellular or subcellular scale and on the macroscopic scale. Unfortunately, this exercise may be difficult, and it has to be decided in each case whether the cellular mechanisms do manifest on the macroscopic scale, or if they can be neglected. Studies, where it has been shown how certain cell-based models can be described by means of macroscopic continuous description, may be helpful for modellers to choose the appropriate cell model [91].

Unsurprisingly, the computational cost of the different tissue models depends on its resolution and on its mechanical complexity. The spheroid model, followed by the vertex model, has, in its basic form, a comparably low mathematical and algorithmic complexity, and is therefore capable of simulating larger tissues. Although the Subcellular Element model can make use of very efficient techniques borrowed from molecular dynamics, the large number of elements to represent a single cell limits this approach to smaller tissues. The Immersed Boundary cell model, finally, buys its very high mechanical accuracy and structural complexity with very high computational costs, and is currently only able to handle small tissues.

Whenever the modeller aims at explaining phenomena which rely intrinsically on three-dimensional effects, a 3D cell model might be required. Oftentimes, if morphogenetic problems are intrinsically three-dimensional in nature, then 2D models require incisive simplifications or compensations to account for the missing third dimension. For example, the mechanics of branching morphogenesis cannot be studied in 2D cross-sectional models without compensation of the circumferential forces which are pivotal for budding and bud lengthening processes. For the same reason, Dillon et al. [2] introduced forces from the anterior to the posterior boundary in their 2D limb development model. Although many 3D models exist (namely the spheroid model, the Subcellular Element model and the Cellular Potts model), 3D versions of the vertex model and the Immersed Boundary model [140] are not yet implemented in software frameworks. The implementation of highly detailed 3D cell models (i.e., with arbitrarily deformable cell geometries) requires high-performance computing (HPC) considerations in order to be able to study large numbers of cells. To our knowledge, currently no computational framework for arbitrary cell shapes is capable to exploit large supercomputers with hundreds to thousands of cores. In order to integrate and standardize cell, signalling, and behavioural model descriptions, promising steps have been taken to develop representation formats, namely SBML [135]. The extension to the description of cell-based models and cellular processes, such as cellular force generation and cell division, is highly desirable in the future.