Simulating tissue mechanics with agent-based models: concepts, perspectives and some novel results

P. Van Liedekerke¹,
M. M. Palm¹,
N. Jagiella¹ &
…
D. Drasdo¹

15k Accesses
166 Citations
9 Altmetric
Explore all metrics

Abstract

In this paper we present an overview of agent-based models that are used to simulate mechanical and physiological phenomena in cells and tissues, and we discuss underlying concepts, limitations, and future perspectives of these models. As the interest in cell and tissue mechanics increase, agent-based models are becoming more common the modeling community. We overview the physical aspects, complexity, shortcomings, and capabilities of the major agent-based model categories: lattice-based models (cellular automata, lattice gas cellular automata, cellular Potts models), off-lattice models (center-based models, deformable cell models, vertex models), and hybrid discrete-continuum models. In this way, we hope to assist future researchers in choosing a model for the phenomenon they want to model and understand. The article also contains some novel results.

Free and Interfacial Boundaries in Individual-Based Models of Multicellular Biological systems

Article Open access 08 October 2023

Introduction to Computational Modeling of Multicellular Tissues

Dynamical Rheological Properties of In-Silico Epithelial Tissue by Vertex Models

Use our pre-submission checklist

Avoid common mistakes on your manuscript.

1 Introduction

Biological cells are constantly exposed to mechanical stress. In response to this stress they may change their shape, migrate, or change their state in another way. In the past decades, research majority focused on the molecular processes and their impact on biological organization processes. Molecular biology was led by the expectation that unveiling the genome, identifying genes and their function, could be sufficient to understand functioning and failure of living organisms. Despite impressive insights triggered by molecular biology, these expectations were not fulfilled. Since then, stepwise higher levels of organization have been linked to the lower levels, arriving at gene expression, gene regulation, gene translation, post translational protein modifications, and signal transduction and metabolism. It is becoming increasingly obvious that molecular events can impact cell behavior and subsequently modify multi-cellular organization, which in turn can feed back to the molecular control inside the cell [1]. Beside feedback by signaling molecules outside of the cell, cells can sense mechanical stress by mechanotransduction, which is the mechanism by which cells transduce an external mechanical stimulus into an internal molecular signal. The dynamics of the cytoskeleton (CSK) and focal adhesion complexes plays a major role herein [2, 3]. Mechanotransduction can materialize itself in different types of cell responses: (i) primary mechanoreceptors like ion channels and integrins. Integrins are transmembrane adhesion proteins that instantaneously mediate coupling between the plasma membrane or the extracellular matrix (ECM) and the cystoskeleton; (ii) mechanosignaling complexes, like caveolae [4] or focal adhesions [5], which respond rapidly (seconds to minutes range) to forces. The protein talin seems to play a crucial role herein since it undergoes force-induced conformational changes and reveals cryptic sites to actin bundles (via vinculin) binding [265]; (iii) signal integrators like the actin CSK [6] which is essential to propagate extracellular forces to the nucleus. The nucleus has also been proposed to act as a mechanosensor [7], where integration of the changes in nuclear shape may cause conformational changes in chromatin organization; (iv) nuclear shuttling proteins [8] (like zyxin, YAP, Yes-associated protein, MRTFA, myocardin-related transcription factor), which transit from the cytoplasm to the nucleus where they affect gene expression. In turn, these signals lead to cell organization and function adjustments. In particular, they alter the mechanical properties and response of the cell. Mechanotransduction is therefore intrinsically a feedback mechanism. For these reasons the understanding of biological organization requires to ultimately include the contribution of mechanics, as well as the level of cells and tissues.

Mechanical stimuli play a prominent role in the development and constitution of an organism. Already in the early stages of development mechanosensing plays a role in stem cell differentiation [9]. This role is maintained in later stages of differentiation, e.g., in chondrocytes in bone remodeling. The tissue engineering community now acknowledges that physical influences, such as ECM geometry, ECM elasticity, or mechanical signals transmitted from the ECM to the cells (i.e., mechanotransduction) are of great importance to explain biological observations [10]. In order to get more insight and quantify the relations between mechanical stress and biological processes, people have started to study the mechanics of tissues such as bone, heart, and the artery system using mathematical models. Yet, experimental studies continue to raise important questions and reveal new observations.

In cancer, cells experience augmented stresses as the tumor develops. Helmlinger et al. [11] set up an experiment where cells were grown in an elastic gel and they found that the growth of the tumor was gradually stalled (reaching compressive stresses of several kPa) although no indication of apoptosis was reported. In similar experiments, however, Cheng et al. [12] concluded induced apoptosis in regions of high compressive stress and allowing proliferation in regions of low stress. Basan et al. [13] studied the influence of mechanical stress by compressing tissue in a cylindrical chamber. This led to concept of homeostatic pressure, i.e., the pressure at which the net result of growth and apoptosis is zero. The same group also studied spheroids in an osmotic solution to mimic external pressure. They observed diminished proliferation rates in the spheroid except at the periphery. Similarly, Alessandri et al. [14] showed that spheroids in a confined environment exhibited a declined growth rate but an increased migrative activity at the periphery of the tumor. The origin of this “skin” effect is not yet clarified. There is still discussion whether increased stress triggers migrative activity in the tumor [14–17].

In wound healing, cells that attempt to repair the tissue are steered by both passive and active (motility) forces. Traction force microscopy studies using a variety of in vitro monolayers have revealed that active motility forces are not only generated by the boundary cells but also by cells in the middle of the epithelia [18]. The leading edge of the epithelial tissue often does not move uniformly but forms finger-like protrusions while cells in the bulk can exhibit spontaneous swirl-like flow patterns, which depend on the cell density and viscous forces between cells and substrate.

Various experiments illustrate the importance of cell–matrix mechanics. Matrix mechanics can foster tissue formation by correlating the relative motions of cells, promoting the formation of cell–cell contacts. For instance, it has been reported that cardiac cells can synchronize their beating through substrate deformations [19]. In blood vessel formation and growth, diverse mechanical forces originating from key cell processes, such as lamellipodium and filopodium formation, regulate the direction and rate of formation of the vessels [20]. Also in this case there is a prominent role for cell–matrix mechanical interactions, which has been recognized in several modeling studies using a continuum approach [21–24] and was recently included in and agent-based angiogenesis model [25]. During metastasis in cancer, cells detach from the tumor and migrate through the stroma. Studies have revealed that cells remodel the matrix around them [26, 27], thereby possibly creating trails for other cells [28, 29].

Modeling cells and tissue mechanics is an emerging field in biomedical sciences. The abovementioned phenomena have created plenty of challenges for engineers, mathematicians, physicists, and biologists, and endowed a whole new subject in the biomechanics and mechanobiology modeling community. Clearly, modeling these phenomena is not a trivial task due to the biological complexity and various scales involved. Agent-based models were developed to understand tissue dynamics as a result of interplay of its individual cells. In these models, cells are by definition regarded as separate units. This approach is contrary to continuum methods, where the individual character of the cells is discarded and tissue dynamics is derived from mesoscopic or macroscopic conservation and constitutive laws. Agent-based systems provide an ideal framework for the integration of intracellular molecular processes as their directly represent the cell itself as individual which captures natural spatial inhomogeneities and variability between cells.

Inhomogeneities on the level of individual cells can be fundamental to understand the organization on the tissue level. For example, it has been realized that understanding of tumor progression and resistance to treatment requires to take into account the cell-to-cell variability [30, 31]. These inhomogeneities can readily be represented within an agent-based approach (e.g., [32]). Tissue architecture in vivo can closely be related to the function of the organ in which case it needs to be represented explicitly. For example hepatocytes, the parenchyma of liver, that detoxify blood from food toxins and other toxins, are arranged within a complex architecture around a network of micro-vessels in such a way that the exchange area between hepatocytes and blood is very large [33]. Another example are intestinal crypts forming one-cell-thick pockets in the intestinal wall hence optimizing resorption from the intestine (e.g., [34–37]). Such inhomogeneities and small scale effects favor agent-based approaches as they are often difficult to be represented in continuum models which usually deal with quantities locally averaged over a group of many cells. Agent-based models can also probe properties of in vitro multi-cellular systems such as monolayer cultures or multicellular spheroids as these systems can, due to their moderate cell population size, be modeled in a 1:1 manner, yielding an abstraction of the “wet” experiment on the computer.

Recent developments of experimental imaging techniques facilitate validation of simulation models at histological scales. Confocal laser scanning micrographs can be used to reconstruct the three-dimensional tissue architecture at histological scales [33]. For liver tissue, standard imaging and image processing protocols have been developed indicating that in the future more and more automated or half automated well established procedures will be established to speed up analysis of tissue organization processes at histological scales [38]. Staining for tissue components such as cell types or cell properties can be used to collect information about the multicellular state. The information can directly be fed into a pipeline of image analysis and modeling, whereby the modeling can be regarded as in silico experiment [39]. Modeling can be either started directly out of reconstructed 3D volume data sets, or be performed in virtual tissue architectures that represent the image information in a statistical sense. Live imaging permits direct observation of cellular arrangement processes in real time over a period of several hours. Hence, combining different image modalities operating on the histological scale can provide complementary information on different time and length scales that can be used to construct and validate agent-based tissue models.

Over the last decades a number of different agent-based approaches have been developed to mimic multicellular organization. Our focus here is on agent-based models in space. “Space-free” agent-based models have been considered for example in hematology and are reported elsewhere [40]. Spatial agent-based models can roughly been distinguished by those operating on space fixed lattices [41–60] and those operating without a lattice (“lattice-free”, “off-lattice”; [52, 61–66]). Among lattice models, either (i) a lattice site may be occupied by many cells (e.g., [67]) permitting modeling of large cell population sizes, (ii) cells may occupy at most a single lattice site (e.g., [68]), or (iii) many lattice sites may be occupied by one biological cell (e.g., [53]). In each of the aforementioned lattice approaches, division, death, and migration are modeled as stochastic processes. The degree to which biomechanical aspects can be captured depends much on which of the three types (i)–(iii) is chosen. Model types (i) and (ii) can represent excluded volume effects, (iii) can capture qualitatively cell deformation and compression, while each of the three approaches to some extend can describe the effect of mechanical forces of one cell on its neighbor, or on a group of neighbor cells. A last type of agent-based model defined on a lattice is the lattice gas cellular automaton (LGCA) [69] in which each lattice site also contains velocity channels. Among lattice-free models one can roughly distinguish between center-based models and deformable cell models. The first class models the interaction of cells based on forces or energies between cell centers and does not resolve cell shape precisely. The second class resolves cell shape in detail and permits for representation of complex cell shapes.

As agent-based models are used by an expanding modeling community, our goal is to provide an overview of agent-based modeling approaches on the modeling aspects with emphasis on the cell mechanics that they can capture, rather than to give an overview of mathematical models developed for the modeling of well-defined phenomena such as tissue growth [70] or tackling a specific clinical problem [71]. For each of the approaches we make the trial to list their specific capabilities, advantages and drawbacks (knowing this is partially worthy of discussion and therefore a bit a matter of taste as so far no absolute quantitative measures for it exist), and provide as such an short guide to people entering or already established community of biomechanical modeling. In addition, we present some complementary so far unpublished (novel) results.

The paper is organized as follows: in Sect. 2 we introduce the lattice models addressing different spatial resolutions. This is followed by off-lattice models, encompassing center-based models (CBM) and deformable cell models (Sect. 3), considering the cell as a object being able to move continuously in space. Finally, we give a short overview on how discrete models can be coupled to continuum models (Sect. 4). An overview of existing software tools is given in Sect. 5.

Where respective simulation results are available we discuss the different model approaches referring to monolayer and multicellular spheroid growth as they may represent a good reference (“traveling salesman”) example for in how far the models are able to capture tissue mechanics of growth processes. By monolayer we mean cell populations growing in vitro on a flat substrate without shortage of nutrients, growth factors etc.. Multicellular spheroids are their three-dimensional counterpart, cell populations growing in vitro either in liquid suspension or collagen-like gel.

2 Lattice models

Among lattice models, either only positions of cells are considered, or, with the class of lattice gas cellular automata (LGCA), in addition the velocity of the cells is denoted (e.g., [68]. We consider first those where the velocity is no explicit variable. In this case either (type A) a lattice site may be occupied by many cells [67] which permits simulation of large cell population sizes, (type B) cells may occupy at most a single lattice site [72], or (type C) many lattice sites may be occupied by one biological cell [53]. Model type C is derived from the Potts model and therefore usually called “Cellular Potts model” [73]. Our choice of the enumeration by letters A, B, C is led by the relative size of a cell compared to the size of the lattice spacing (constant): type A: cell size is much smaller than the lattice spacing, type B: cell size is the same or about the same as the lattice spacing, type C: cell size is much bigger than the lattice spacing. LGCA (denoted as model type D) usually also have several cells on the same lattice site, but in addition carry a velocity variable.

For didactic reasons we first discuss model type B, then model types A, C and D in subsequent subsections.

2.1 Cellular automaton with one cell per lattice site (type B)

We start with type B as reference case because type A can be described as a generalization of type B.

Cellular automaton models of type B are used quite a lot in cancer modeling to study how intrinsic cell mechanisms contribute to tumor growth, death, and morphology (e.g., [74–77]), and have been intensively used to model monolayers and multicellular spheroids (e.g., [67] and refs. therein).

The growth, death and migration dynamics of cells on a lattice is usually modeled by a stochastic process. In type B automata, the volume of a cell is usually identified with the volume associated with that lattice site. On a cubic lattice with lattice spacing a this is $a^d$, with d being the spatial dimension. A commonly used alternative to the structures lattice is the Delaunay lattice [72, 78], in which the number of cell neighbors corresponds to the value found in epithelial tissues [79, 80]. The Delaunay lattice is generated by seeding a set of constructions points and then generating the corresponding Voronoi mesh and Delaunay lattice (Fig. 1).

Block et al. [72] proposed generating a 2D Delaunay lattice by placing the construction points into each square of a square lattice to ensure that cell sizes only slightly vary around $a^2$. This way of constructing the cellular automaton lattice was later extended to three dimension by Radszuweit et al. [67]. For illustration, Fig. 2 represent an recent example for a simulation of population growth in a CA model on an unstructured lattice and within a graphic scheme illustrates the key processes that can occur. There are growth, migration, division, death, or pushing [72, 81]. The processes are explained below.

The cell cycle time distribution found experimentally is peaked and does not correspond to the assumption that division is a Poissonian process, as such a process generates a cell cycle time distribution decreasing exponentially with time. By introducing intermediate states with Poissonian transitions, the experimentally observed cell cycle distribution can be reproduce and the emerging cell cycle time distribution is an Erlang distribution:

$$\begin{aligned} f(\tau ) = \frac{\lambda (\lambda _m\tau )^{m-1}}{(m-1)!}e^{-\lambda _m\tau }, \end{aligned}$$

(1)

where $\lambda _m = m/\tau _\mathrm{e}$ is the transition rate between subsequent states, $\tau $ the cell cycle time, and m the total number of states. This yields $\langle \tau \rangle = \tau _\mathrm{e}$.

The magnitude of the expansion speed of dense one-cell-thick monolayers is incompatible with the assumption that only those cells at the monolayer border grow and divide [82]. Furthermore, in monolayers both glucose and oxygen are abundant, suggesting that the expansion speed of the monolayer may be controlled by physical forces [72, 82, 83]. Different from a monolayer in a multicellular spheroid, the three-dimensional counterpart of a monolayer, nutrients become limiting in the interior of the multicellular spheroid if the size of the multicellular spheroid grows beyond a certain size. A careful analysis of EMT6/Ro growing multicellular spheroids shows significant central necrosis above a certain size which is associated to the lack of glucose and oxygen (see [82] and refs. therein). The central necrotic core and the diameter of these multicellular spheroids grow linear in time [84] suggesting the existence of a proliferating rim of constant size. Moreover, change of glucose does not change the growth speed, while it does change the size of the necrotic core.

All these observations suggest, that biomechanical forces, presumably by the same mechanisms, play an important role in the control of the growth speed in both monolayer and multicellular spheroids, and that the proliferating rim is constant in time.

A growing and dividing cell can exert forces on neighboring cells pushing them aside to generate free space for its division. In type B cellular automaton this was taken into account by defining a certain length $\Delta L$ over which a cell is able to push neighbor cells away [72, 78, 83] (Fig. 2a). This would generate a constant proliferating rim and comply with the observations in growing monolayers and multicellular spheroids. A cell can only divide if at least one lattice site within radius $\Delta L$ is free. To make space for the daughter cells, a line is drawn linking the cell’s site to the closest free lattice site and all intersecting neighbors are shifted one position towards the free lattice site. Next, one daughter cell is placed at the site of the mother cell and the other one is placed at the lattice site that was freed by the shifting procedure (Fig. 2). This procedure ensures that the daughter cells of the same mother cell are neighbors directly after division. Instead of using for $\Delta L$ a fixed length, one may also use an upper bound for the number of neighbor cells a dividing cell can push aside. The algorithm mimics that a dividing cell exerts forces on its environment, leading to a rearrangement of cells in its neighborhood such that the total energy needed to push the neighbor cells away is minimal. Most models used $\Delta L$ as a threshold parameter with cells having the same chance of dividing if their distance to the next free lattice site does not exceed $\Delta L$ [72, 78, 83], but direct comparison with experimental data on multicellular spheroids suggest that the probability of division depends on the distance of the cell to the spheroid border leading to a slightly different dynamics [81].

Using the shifting procedure in a type B cellular automaton, the grow kinetics of one-cell-thick monolayers can be reproduced (Fig. 3).

The processes migration, growth, birth and death are usually modeled by a master equation providing an equation for the time development of the multivariate probability to find the entire system in a particular state. The master equation is a balance equation and permits to calculate transitions into accessible neighboring states. The transitions in a type B cellular automaton are described as:

$$\begin{aligned} \frac{\partial p(Z,t)}{\partial t} \!=\! \sum _{Z'} W(Z'\rightarrow Z) p(Z',t) - W(Z\rightarrow Z') p(Z,t), \end{aligned}$$

(2)

where the first term denotes all transition from neighbor states $Z'$ into the “current” state Z and the second term denotes all transitions from the current state Z into any accessible neighbor state $Z'$. This master equation may either be solved with a fixed time-step algorithm choosing the time step $\Delta t$ so small that only a single event is likely to occur within $\left[ t \ldots t+\Delta t\right] $, or by the Gillespie algorithm (often also called Kinetic Monte-Carlo) [85, 86]. In the Gillespie algorithm a variable time step is calculated by $\Delta t=-\left( 1/\sum _{Z'}W(Z\rightarrow Z')\right) ln(1-\xi )$, where $\xi \in [0,1)$ is a uniformly distributed random variable. $W(Z\rightarrow Z')$ denotes the rate of a transition from Z to $Z'$. Different from chemical particles, hopping of cells should not depend on the number of free neighbor sites as long as at least one free neighbor site exists. Accordingly, one may only sum over $\tilde{Z}(Z')$ states assuming that the rate at which a cell changes its state by a hop, a progress in the cell cycle, a division, or death process is independent of the number of accessible states as long as at least one state is accessible, that is, one free adjacent lattice site in case of a hop and one free site within a circle of radius $\Delta L$ in case of a division. For example, a cell in $d=3$ having three free neighbor sites would then be chosen for a move with rate $(\lambda _1+\lambda _2+\lambda _3)/3$ (instead of with $\lambda _1 + \lambda _2 + \lambda _3$) if $\lambda _i$ with $i=1, 2, 3$ denoting the individual hoping rates. Once it has been selected for a move, the move to lattice site j is performed with probability $\lambda _j/(\sum _{j=1}^3\lambda _j)$. An equivalent line of argument holds for division. Growth can be modeled by allowing a cell to occupy two adjacent lattice sites simultaneously without considering the cell as having divided [81]. This means, that in this case the same cell would occupy more than one lattice site. This principle can be extended to occupation of many lattice sites of a cell then changing into an occupation of at most double of many lattice sites during growth until division takes place. A migrating cells may also be able to exert forces to a neighbor cell which can be implemented in the same way as division [72]. Cell–cell or cell–extra-cellular matrix interactions can be included in the transition rate for each move, which then becomes $\propto e^{-\Delta E/F_\mathrm{T}}$ where $\Delta E$ denotes the total energy of change for a particular cell move, $F_\mathrm{T}$ a reference energy that was proposed to represent the membrane fluctuation energy as a formal equivalent to $k_\mathrm{B}T$ in fluids ($k_\mathrm{B}$: Boltzmann constant, T: temperature) [87].

By systematic simulations it was possible to infer the growth speed of a growing population with hopping rate $\phi $, shift radius $\Delta L$, cell cycle time $\tau $, and Erlang-number m in the linear growth regime for monolayers (Fig. 3a) to

$$\begin{aligned} v^2 \approx B^2\left( \underbrace{[\Delta L'(\Delta L)]^2/\tau _\mathrm{eff}^2}_{(1)}+\underbrace{\phi /\tau _\mathrm{eff}}_{(2)}\right) , \end{aligned}$$

(3)

with $\tau _\mathrm{eff}^{-1}=m(2^{1/m}-1)\tau ^{-1}$, $\tau = \tau _\mathrm{e}$ (with $\lambda _m=m/\tau _\mathrm{e}$ being the rate with which the internal state of the cell is increased by $+1$), $\Delta L'(\Delta L)=A(\Delta L-1)+1$, with $A\approx 0.68$, $B\approx 1.4$ for the 2D case, and $A\approx 0.708$, $B\approx 1.236$ for the 3D case. $\phi $ is here expressed in units of $a^2/\tau _\mathrm{e}$, a is the lattice spacing and $\tau _\mathrm{e}$ the expected cycle time. The term (1) expresses the front movement due to proliferation, (2) the front movement due to random migration. In case the proliferation length is much larger than the diffusion length, $v\propto \Delta L/\tau _\mathrm{eff}$ i.e., the mechanical pushing movement controls the expansion of the moving population surface. In the opposite case, where the diffusion length is much larger than the proliferation length, the limiting continuum equation of the Type-B cellular automaton can be shown by a statistical field theoretical approach to be a stochastic Fisher–KPP equation with and intrinsic multiplicative noise term [83]. The velocity change with increasing proliferation length on large scale can be included in the Fisher–KPP equation by some heuristic redefinition of the diffusion constant but the short range behavior of the cellular automaton type B model cannot correctly be captured by this approach. The reason is, that in the derivation of the Fisher–KPP equation, local homogeneity on small scales had to be assumed which is only valid if the diffusion length overcomes the proliferation length [83].

Recent direct comparison of model simulations with a type B cellular automaton model with KI 67 staining in multicellular spheroids of non-small cell lung cancer (NSCLC) cells shows that the assumption that a cells divide with equal probability as long as the dividing cell has to push at most $\Delta L/a$ (a: lattice spacing) cells away in order to divide does not suffice for the quantitative reproduction of the experimentally observed spatial proliferation pattern [81]. The experimentally found proliferation profiles can be correctly reproduced if the transition rate for division is assumed to be $\propto e^{-d_{i}/\Delta L}$. Here $d_{ij}$ is the distance between the dividing cell and the lattice site that will be occupied as a consequence of the division and the shift of neighbor cells coming with it. The distance is proportional to the number of neighbors that has to be shifted, and might be proportional to the energy to shift the neighbors. Also other model variants have been studied addressing cells pushing during migration, cell–cell adhesion etc. [72].

In the CA type B approach, cells can neither be compressed nor deformed. Only rigid cell body movements are possible. On a regular lattice, this corresponds to rigid body movements into discrete directions. However, in case of noise reduction achieved by choosing many intermediate steps between subsequent cell divisions and suppressing migration, the emerging growth pattern reflect the lattice symmetry [72, 83] (see Fig. 4). On an irregular, unstructured lattice, cell size of a moving cell is only conserved on the average. The advantage of unstructured lattice is that the emerging growth figures do not reflect the lattice symmetry even in the presence of noise reduction. This has great advantages if quantities sensitive to lattice artifacts should be studied by modeling, such as the surface scaling properties of growing monolayers that by using a type B cellular automaton model on an unstructured lattice could be shown by Block et al. [72] to clearly reveal the scaling properties of Kardar–Parisi–Zhang (KPZ) universality class [88, 89]. This model prediction contradicted claims emerging from interpretations of experimental findings at that time [90], saying that the surface scaling of growing monolayers and multicellular spheroids (generally: solid tumors) should be in Molecular-Beam Epitaxy (MBE) universality class. The model prediction was subsequently confirmed [91, 92] indicating that despite their simplicity the CA models can provide valid predictions.

Interestingly, introducing $\Delta L$ on a regular lattice reduces lattice effects [83], and may upon a slight modification of the above algorithm eliminate them [93].

Below we briefly summarize advantages and disadvantages of type B cellular automaton models from our perspective. Evaluations on the computation times are done with regard to a standard PC.

2.2 Cellular automaton with many cells per lattice site (type A)

In type A automata multiple cells $(N_{\max }>1)$ can occupy a lattice site. Consequently, each lattice site represents a spatial compartment much larger than the size of the individual cell. Birth, death and migration processes can occur within each compartment and between neighboring compartments. As with the type B cellular automata, the changes of the number of cells in a certain state can be formulated as chemical reactions and turned into a master equation for the multivariate probability of finding a certain multicellular configuration. However, the position of a cell cannot be resolved in this approach, only its compartment is known. Simulations with type A cellular automata can be speed up by tracking the number of cells in a certain state instead of tracking each cell individually. An excluded volume assumption is implemented by only permitting cells to divide within a compartment or into a neighbor compartment, or to hop into a neighbor compartment if this compartment is not already full. Similar to the type B cellular automata, mechanical pushing is included by defining a maximum length $\Delta L > b$ over which cells can be pushed, with $b > a \approx V^{1/3}$ denoting the lattice spacing, V the cell volume [67]. When a cell divides in a compartment that is not completely filled (i.e., $N < N_{\max }$), the number of cells in that compartment increases by 1. When the compartment is full, the number of cells in the closest non-filled compartment, within radius $\Delta L$, increases by one.

Type A and type B cellular automata can be matched very well when the compartment size is set to $b=\Delta L$ (Fig. 5). Also migration rate $\phi $ and internal state m as a representation of the cell cycle can be matched to produce the same outcome with type A and type B cellular automata models (Fig. 5). Therefore, a type A cellular automaton can be used as a coarse grained model of a type B cellular automaton model given the lattice spacing of the type B cellular automaton is properly chosen. This is advantageous if model parameters are calibrated with in vitro experiments (monolayer or multicellular spheroids) in which the cell population sizes do usually not exceed 300, 000–400, 000 cells, and this calibrated model should be used to predict possible growth scenarios of in vivo tumors with $10^9{-}10^{10}$ cells. As an example, Fig. 3 shows such a simulation for a Xenograft [67]. However, growth of tumors in vivo are largely controlled by the availability of oxygen, glucose and growth factors. Hence the spatial profiles of these molecules have to be equally coarse grained. Preliminary attempts at developing numerical schemes for this indicated that for high concentrations of molecules outside the tumor, coarse graining schemes yielding to the same growth and death of tumor cell populations can be found while for small concentrations of molecules growth and death of the tumor cells deviated strongly. As a result, the growth dynamics in the coarse grain type A CA can deviate significantly from the growth dynamics of the type B CA where each single cell is represented. These deviations may occur because at high consumption rates, caused by a large number of cells in a compartment, the approximation schemes fail [94].

A way to deal with this problem is to construct a hybrid lattice which represents the cell population at high resolution with a type B cellular automaton where strong gradients occur, and otherwise consider the coarse grained automaton type A. Such a hybrid model is presented in Fig. 6. In order to maximize the resolution where necessary and coarse grain where no accuracy is needed, ideally the hybrid model would perform automated switches of the resolution depending on some criterion. For the example shown in Fig. 6 the small compartments have the size of one cell while the large compartments contain many cells. During the simulation all compartments which are occupied by an heterogeneous mix of cells are refined to the single-cell-resolution. Whereas, when all lattice sites associated to a large compartment are occupied exclusively by cells of the same type, a group of lattice sites are represented by a single lattice site on the coarser scale. This way only those compartments with a heterogeneous mix of cells or not completely filled with cells are resolved at single-cell scale. This algorithm especially resolves cells individually in border regions as at the tumor-medium interface.

Below we briefly summarize advantages and disadvantages of type A cellular automaton models from our perspective. Evaluations on the computation times are done with regard to a standard PC.

2.3 Lattice gas cellular automata (type D)

LGCA are a special kind of cellular automata. In addition to the lattice model types discussed above, LGCA models include velocity channels i.e., particles are characterized by their position and velocity, which are both are discrete. Frisch, Hasslacher and Pomeau demonstrated in 1986 that by defining a set of rules for the interaction (in their case collision) and propagation of particles on a regular lattice with discrete time steps and a small set of velocities, one can generate dynamics that on large scales reproduce the behavior of the incompressible Navier Stokes equation [95–97]. A simple intuitive derivation can be found in Ref. [98]. In the simulations of fluids, the velocity channels represented particle collisions. Deutsch and co-workers extended the framework to model growth and migration of cells in multi-cellular environments [46, 69, 99, 100].

In our description we very closely follow the line of argument presented by Hatzikirou and Deutsch in [101], which we consider as an excellent presentation addressing the underlying concepts of LGCA for modeling of tissue organization and growth.

LGCAs are defined on a regular lattice, which, in two dimensions, are usually either hexagonal or square lattices. Each node possesses b velocity channels, where b is the number of neighbors of a lattice site, and $m\in {0, 1, 2, \ldots }$ rest channels (Fig. 7). Assuming the coordinate system is fixed at the center of the node, the rest channel(s) is (are) placed at (0, 0), and the four velocity channels of a rectangular lattice site are (with $v_k>0$) at $(1,0), (0,1), (-1,0), (0,-1)$. A node’s state is given by $\eta (\underline{r})=(\eta _1(\underline{r}), \eta _2(\underline{r}), \ldots , \eta _{b+m}(\underline{r}))$. The number of particles on a site represent a type of microscopic density (called node density) and can be calculated by summing up all occupation numbers $\eta _i(\underline{r}) \in {0, 1}$ at a node $\underline{r}$:

$$\begin{aligned} n(\underline{r},t)=\sum _i^{b+m}\eta _i(\underline{r},t). \end{aligned}$$

(4)

The dynamics of an LGCA emerges from applying superpositions of local probabilistic interactions and deterministic transport steps. The definition of these steps have to satisfy the exclusion principle, i.e., each channel can at most be occupied by one particle (here: cell). Cells can move to a neighboring lattice site, and divide or die at discrete time steps with an simultaneous update at all nodes. The temporal evolution of a state $\eta (\underline{r},t)\in \{0,1\}^{b+m}$ is determined by the temporal evolution of the occupation numbers $\eta _i(\underline{r}, t)$ for each $i\in \{1, \ldots , m+b\}$. With probability $P(\eta \longrightarrow \eta ^\mathrm{C})$, the pre-interaction state $\eta _i(\underline{r},t)$ is replaced by the post-interaction state $\eta _i^\mathrm{C}(\underline{r},t)\in (0,1)$ according to

$$\begin{aligned} \eta _i^\mathrm{C}(\underline{r},t)= & {} R_i^\mathrm{C}\left( \{\eta (\underline{r}, t)|\underline{r}\in \mathcal {N}_\mathrm{b}(\underline{r})\}\right) , \end{aligned}$$

(5)

$$\begin{aligned} \eta ^\mathrm{C}(\underline{r},t)= & {} \left( R_i^\mathrm{C}\left( \{\eta (\underline{r}, t)|\underline{r}\in \mathcal {N}_\mathrm{b}(\underline{r})\}\right) \right) ^{b+m}, \end{aligned}$$

(6)

$\mathcal {N}_\mathrm{b}(\underline{r}) := \{\underline{r}+\underline{v}_i : \underline{v}_i\in \mathcal {N}_\mathrm{b}), i=1,2, \ldots , b\}$ denotes a finite list of the node localized at position $\underline{r}$ on the lattice, $R_i^\mathrm{C}$ denotes the interaction rule. This is the time-independent probability $P(\eta \rightarrow \eta ^\mathrm{C})$ for transition from the pre-interaction to the post-interaction node state. In the deterministic propagation step, all particles are moved simultaneously to nodes in the direction of their velocity i.e., a particle residing in channel $(\underline{r}, \underline{v}_i)$ at time t moves to another channel $(\underline{r}+\tilde{m}\underline{v}_i, \underline{v}_i)$ ($\tilde{m}\in \mathbb {N}_0$) during one time step. This is represented by

$$\begin{aligned} \eta _i(\underline{r} +\tilde{m}\underline{v}_i, t+\tau ) = \eta _i^{\tilde{P}}(\underline{r},t), \end{aligned}$$

(7)

where, $\tilde{m}\in \mathbb {N}_0$ denotes the particle speed, $\tilde{m}\underline{v}_i$ the translocation of the particle, $\tau $ denotes the time step. Index $\tilde{P}$ indicates that it is a streaming step.

Propagation according to the last equation respects mass and momentum conservation. As a consequence all particles (cells) in the same velocity channel move the same number of $\tilde{m}$ lattice units. Thus, in contrast with the CA models of type A and B, cells in LGCAs can move in one time step over many lattice sites. Combining interaction (C) and propagation $(\tilde{P})$ leads to

$$\begin{aligned} \eta _i(\underline{r}+\tilde{m}\underline{v}_i, t+\tau ) = \eta _i^{C\tilde{P}}(\underline{r},t). \end{aligned}$$

(8)

The change of the occupation number due to the combined interaction and propagation is

$$\begin{aligned}&\eta _i(\underline{r}+\tilde{m}\underline{v}_i, t+\tau ) - \eta _i(\underline{r},t) = \eta _i^{CP}(\underline{r},t) - \eta _i(\underline{r},t) \nonumber \\&\quad = C_i\left( \eta _{\mathcal {N}_\mathrm{b}(\underline{r})}(t)\right) ,\quad i=1, \ldots , b+m, \end{aligned}$$

(9)

with

$$\begin{aligned}&C_i\left( \eta _{\mathcal {N}(\underline{r})}(t))\right) \nonumber \\&\quad = \left\{ \begin{array}{cl} 1 &{} \text {creation of a particle in channel}\, (\underline{r}, \underline{v}_i)\\ 0 &{} \text {no change in channel}\, (\underline{r}, \underline{v}_i) \\ -1 &{} \text {annihilation of a particle in channel}\, (\underline{r}, \underline{v}_i). \end{array} \right. ,\nonumber \\ \end{aligned}$$

(10)

where $C_i$ denotes the change of the occupation number.

So far, the upper outline summarizes local rules applied on each site of the lattice. However, the dynamics can also be formalized on scales larger than individual site of the lattice following the following conceptual ideas. The overall dynamics at this stage can be described by microdynamical difference equations that include the probabilistic interaction as well as the deterministic propagation. Ensemble averaging leads to an equation for the single-particle distribution functions, $f_{i}(\underline{r}, t) = \langle \eta _i(\underline{r},t)\rangle $ which typically contains a collision term. In a continuum expansion, the difference equation for the single-particle distribution function can be turned into a differential equation. The collision term generally involves higher-order particle distribution functions hence requiring to solve a hierarchy of coupled equations for the one, two, three, ...particle distribution functions. This is the BBGKY (Bogolyubov, Born, Green, Kirkwood, Yvon) hierarchy. Breaking up after the first order, and approximating the two particle distribution function by the one particle distribution function leads to the Boltzmann equation from which by integration over the velocity the density, momentum density and energy density can be calculated. As these operations can be equally performed for the LGCA, it permits to define a pressure tensor. The kinetic pressure in the LGCA is defined as $p_K=(1/V_0)\sum _i\underline{v}_{i}^2f_i$, with $V_0$ being the elementary volume associated to a node, $f_i$ the single-particle distribution function, and $\underline{v}_i$ with $i_1, \ldots , b$ denotes the velocity channels [97]. The thermodynamic pressure is defined by $p=(1/(V_0\beta )) \left( \sum _i f_i +(1/2)\sum _if_i^2 + \cdots \right) $ [97] with $\beta = 1/(k_\mathrm{B}T)$, T being temperature, $k_\mathrm{B}$ the Boltzmann constant. Hence relations for the kinetic and thermodynamic pressure can be derived for the LGCA. The local density is defined as $\varrho (\underline{r}, t)=\sum _i f_{i}(\underline{r}, t)$. Equivalently, a momentum density can be defined by summing $\varrho \underline{v} = \sum _i \tilde{m}\underline{v}_if_i(\underline{r},t)$. While the rules of classical LGCA are chosen in order to reproduce classical partial-differential equations describing fluid dynamics such as the mass conservation and the Navier–Stokes equation for an incompressible fluid as the momentum conservation, for cell systems momentum conservation is not assumed.

Deutsch and co-workers established extended the framework to include birth and death. In a most basic model they permit cells to move according to a propagation operator, reorient, and grow and divide. In the simplest case the rules are chosen such that isolated cells perform a random walk with transition probability

$$\begin{aligned} {\mathcal P(\eta \rightarrow \eta ^0)(\underline{r}, t)} = \frac{1}{Z}\delta (n(\underline{r},t), n^0(\underline{r},t)) \end{aligned}$$

(11)

with a normalization factor $Z\!=\!\sum _{\eta ^0(\underline{r},t))}\delta (n(\underline{r},t), n^0(\underline{r},t))$. The index O here denotes re-orientation.

The rule for birth is:

$$\begin{aligned} \mathcal {R}_i(\underline{r},t) = \xi _i(\underline{r},t)(1-\eta _i(\underline{r},t)), \end{aligned}$$

(12)

where $\xi _i$ are random Boolean variables with $\sum _{i=1}^{b+m} \xi _i(\underline{r},t)=1$. This rule takes into account that occupied velocity channels cannot be occupied by a second cell hence mimics contact-inhibited growth. The corresponding probabilities are:

$$\begin{aligned} \mathcal {P}(\xi _i(\underline{r},t)=1)=r_M\frac{\sum _{i=1}^{b+m} \eta _i(\underline{r},t)}{b+m}, \end{aligned}$$

(13)

with $r_M$ being a proportionality constant. Notice that if no cell is on the node, division cannot occur, and the probability is proportional to the number of occupied channels.

Death occurs with rate $r_d=(b+m-C)/(b+m)r_M$, with $C\le b+m$ being the maximum node occupancy. Results are shown in Fig. 8. As in the CA-types A and B models above without a velocity channel, growth is linear.

The above dynamics is fully specified by the following microdynamical equations, that form a useful starting point to derive the large-scale behavior (see below):

$$\begin{aligned} \eta _i^\mathrm{R}(\underline{r}, t)= & {} \eta _i(\underline{r}, t)+\mathcal {R}_i(\underline{r}, t) \end{aligned}$$

(14)

$$\begin{aligned} \eta _i^\mathrm{R}(\underline{r}+m\underline{v}_i, t +\tau )= & {} \sum _{j=1}^{b+m}\mu _j(\underline{r},t)\eta _j^\mathrm{R}(\underline{r},t). \end{aligned}$$

(15)

The first of the two equations refer to the growth operator $\mathcal {R}$, the second to the redistribution of cells on the velocity channels, and the propagation to the neighboring nodes. The growth operator associates a new occupation number for a given channels by a stochastic growth process. $\mu _j(\underline{r}, t) \in \{0, 1\}$ are Boolean random variables selecting only one of the $b+m$-terms on the rhs of the equation 15. The random walk is implemented as a simple reshuffling of cells within the node channels leading to the probability of choosing a channel: $\langle \mu _j\rangle = 1/(b+m)$ for $j=1,\ldots , b+m$. The terms $\mathcal {R}_i(\underline{r}, t)\in \{0,1\}$ for $i=0,1,\ldots ,b+m$ represent birth/death processes of cells in channel i and are applied to each channel independently.

By local ensemble averaging, the lattice Boltzmann equation can be obtained from Eqs. (14), (15):

$$\begin{aligned} f_i(\underline{r}+m\underline{v}_i, t+\tau ) - f_i(\underline{r},t)= & {} \sum _{j=1}^{b+m}\varOmega _{ij}f_j(\underline{r}, t) \nonumber \\&+ \sum _{j=1}^{b+m}(\delta _{ij} + \varOmega _{ij}) \tilde{\mathcal {R}}_j(\underline{r},t).\nonumber \\ \end{aligned}$$

(16)

Here $\varOmega _{ij} = 1/(b+m) - \delta _{ij}$ is the transition matrix of the shuffling process, $\tilde{\mathcal {R}}_i=F(\varrho )/(b+m)$ the averaged reaction term, that is assumed to be independent of the particle direction. $F(\varrho )$ is the cell reaction term for a single node. Inserting the upper equations, one obtains for the reaction term

$$\begin{aligned} \tilde{\mathcal {R}}_j(\underline{r},t) = r_Mf_i(\underline{r}, t)\left( 1 - \frac{r_D}{r_M}-f_i(\underline{r},t)\right) . \end{aligned}$$

(17)

Using a Chapman–Enskog expansion, the macroscopic dynamics can be obtained from the upper equation (14, 15) by using diffusive scaling according to $\underline{x}=\epsilon \underline{r}$, $\tilde{t}=\epsilon ^2t$. The emerging equation for slow growth rates is of the Fisher–KPP type:

$$\begin{aligned} \frac{\partial \varrho }{\partial \tilde{t}} = \frac{\tilde{m}^2}{(b+m)\tau }\nabla ^2\varrho + \frac{1}{\tau }F(\varrho ), \end{aligned}$$

(18)

with $F(\varrho )=r_M\varrho (C-\varrho (\underline{x}, \tilde{t}))$. The mean field equation shows qualitatively the same behavior as a simulation with the LGCA. The minimum growth speed for this equation, well-known and Fisher–KPP equation, is $v_{\min }\ge 2\frac{\tilde{m}}{\tau }\sqrt{\frac{r_M}{(b+m)}}$. This qualitatively agrees with the estimate for the type B automaton made in Ref. [83], for which in case where the proliferating rim is much larger than the diffusion length, $v\propto \Delta L/\tau _\mathrm{eff}$ has been found (see above), if $\Delta L\propto \tilde{m}$ is chosen. I.e., the introduction of the velocity channels indicate the same function as the introduction of the distance over which dividing cells are able to shift their neighbor cells.

However, to reach quantitative agreement between the above Fisher–KPP equation and the LGCA, the long tail of the traveling front needs to be cut, which may occur by introducing a cut-off. In this case the behavior is as those of the cellular automaton without velocity channel, for which on long time scales a stochastic Fisher–KPP equation could be obtained [83], where an intrinsic multiplicative noise term cuts off the long small density tail of the traveling front.

Below we briefly summarize advantages and disadvantages of type D cellular automaton models from our perspective. Evaluations on the computation times are done with regard to a standard PC.

2.4 Cellular automaton with many lattice sites per cell (type C; cellular Potts models)

The cellular Potts model (CPM) is an energy-based, lattice-based modeling method for cell-based modeling. It uses an energy functional generalized from the Potts model to evaluate a multi-cellular state, and it is therefore named Cellular Potts model. The Potts model is a generalization of the celebrated and extensively studied Ising model, both are used to describe phenomena in solid state physics as e.g., ferromagnets.

In CPM, one cell occupies many lattice sites. Different from the CA models discussed before, the CPM explicitly represents the cell shape which has made the CPM a popular tool to model morphogenic processes such as cell sorting [53, 103], cancer and tumor growth [104–110], and angiogenesis [25, 111–117].

In the CPM cells are positioned on a 2D or 3D lattice. Each site $\underline{x}$ on the lattice is assigned an identifier $\sigma \in \mathbb {N}^0$. All sites with the same, positive identity $\sigma $ are part of the same cell, and all remaining lattice sites, with $\sigma =0$, belong to the medium. Cell migration, cell growth and cell shape changes are modeled using a Markov chain Monte Carlo method. The CPM iteratively attempts to change the state of a randomly chosen lattice site $\underline{x}$ to that of a randomly chosen neighboring lattice site $\underline{x}^\prime $. To determine whether such a copy attempt is successful the change in effective energy of the configuration ($\Delta E$) is evaluated using the Metropolis algorithm: $p(\sigma (\underline{x}) \rightarrow \sigma (\underline{x}^\prime )) = \min (1,e^\frac{-\Delta E}{F_\mathrm{T}})$, where $F_\mathrm{T}$ denotes an effective fluctuation energy that controls the cell motility. Each time step, this sampling process is repeated until the number of copy attempts is equal to the number of lattice sites. Afterwards, the model parameters may be updated and non-energy driven process such as cell division may occur.

The effective energy E summarizes the effects of the modeled cell behavior. The standard effective energy of the CPM consists of two terms:

$$\begin{aligned} E= & {} \sum _{(\underline{x},\underline{x}^\prime )} J(\sigma (\underline{x}),\sigma (\underline{x}^\prime ))(1-\delta (\sigma (\underline{x}),\sigma (\underline{x}^\prime ))\nonumber \\&+\, \lambda _\mathrm{V}\sum _\sigma (V_{0}-V(\sigma ))^2 \end{aligned}$$

(19)

The first term represents both cell–cell adhesion and cell–ECM adhesion that is computed for all sets of neighboring lattice sites $(\underline{x},\underline{x}^\prime )$ where $J(\sigma (\underline{x}),\sigma (\underline{x}^\prime ))$ represents the energy cost of an interface between two cells. Because only interfaces at the cell membrane are associated with energy costs, the Kronecker delta $(\delta (x,y) = \{1,x=y;0,x\not =y\})$ is used such that the contact energy becomes zero for $\sigma (\underline{x}) = \sigma (\underline{x}^\prime )$. The second term is a volume constraint, with cell volume $V(\sigma )$ and target volume $V_{0}$, and the penalty parameter $\lambda $. For 2D simulations volume should be replaced with area.

In case the cell is assumed to be homogeneous isotropic elastic, one would expect $\lambda _\mathrm{V} = K_\mathrm{V} = E/(3(1-2\nu ))$ with E being the Young modulus of the cell, $\nu $ being the Poisson ratio, and $K_\mathrm{V}$ the bulk modulus. To model growth, the target volume V may be updated at each simulation step. Once the target volume has adopted twice its value after mitosis, and the ratio of actual volume and target volume $\varTheta = V_{0}/V$ overcomes a critical threshold, a cell splits into two daughter cells. A scenario is shown in Fig. 9 [118]. It can be shown that the dependency of the growth kinetics on the model parameters is the same for the CPM as for the lattice type A and B models. A larger Erlang parameter m leads to slower growth and, as in the center-based off-lattice models discussed below, growth is faster if $\varTheta $ is decreased (Fig. 10). Non-shear contributions to the elastic energy are controlled by the bulk modulus. The change of hydrostatic pressure dp and the volume change dV are related by

$$\begin{aligned} \frac{dp}{dV}=-K_\mathrm{V}\frac{1}{V}, \end{aligned}$$

(20)

leading to $p=p_{0}-K_\mathrm{V} \log (V/V_{0})$, with $p_{0}$ being the hydrostatic pressure of the cell at volume $V_{0}$. For small volume changes, Taylor expansion for small $(V-V_{0})/V_{0}$ yields the approximate pressure

$$\begin{aligned} p=p_{0} - K_\mathrm{V}\frac{V-V_{0}}{V_{0}}, \end{aligned}$$

(21)

with $K_\mathrm{V}=E/(3(1-2\nu ))$ as defined above, hence the volume threshold can be related to a pressure. With increasing E the pressure increment $p-p_{0}$ necessary to generate the same volume deviation increases or, vice versa, in order to generate the same pressure increment for a large Young modulus E than for a small Young modulus E, the volume difference $(V-V_{0})$ must be smaller. With an increasing Young modulus and a fixed $\varTheta $ the growth speed increases, as also found for the center-based model [82].

The volume constraint also demonstrates a shortcoming of the CPM. For incompressible cells, volume changes during migration are associated with infinite energy increase. However, a cell in the CPM moves by flipping lattice sites, which, according to the Metropolis scheme, has probability zero in case a flip is associated with an infinite energy increase. This means that an incompressible cell cannot move in the standard CPM described by Eqs. (19) and the choice of $\lambda _\mathrm{V}$ with $\nu = 0.5$. The reason is that cell movement in the standard CPM is inherently linked to volume changes because the CPM does not distinguish between the main cell body and filopodia. The relation to energy and temperature has to be properly chosen to permit reasonable movement.

Another commonly used term in the effective energy is the surface area conservation: $\lambda _S \sum _\sigma (S_{0}-S(\sigma ))^2$, with surface area $S_{0}$ and target surface area S [119]. In $d=2$ dimensions it becomes a line energy. This term is either added as an extra term, or replaces the volume conservation term. When the surface area or the perimeter constraint energy contribution are combined with contact energies (J) that are permitted be negative, then the roughness of the cell border can be controlled and may be linked to experimental observations [120]. Note that the magnitude of both the surface area conservation term and the volume conservation term depend on the fineness of the lattice. For example, when the number of lattice sites on a 2D lattice doubles, a deviation of 1 lattice site yields the same energy while the actual deviation is twice as small. To fix this issue, both the deviation of the volume or surface area may be divided by the target volume or target surface area [25, 121].

The effective energy only includes static cell properties such as size and neighbors. To include dynamic cell behavior, such as chemotaxis, extra terms may be added to the energy change $\Delta E$. The most common mechanism that is included in this way is chemotaxis [55]:

$$\begin{aligned} \Delta E_\text {chemotaxis} = -\chi \left( c(\underline{x}^\prime )-c(\underline{x}) \right) , \end{aligned}$$

(22)

where $\chi $ denotes the chemotactic sensitivity and $c(\underline{x})$ is the concentration at position $\underline{x}$. In Ref. [122], external forces have been defined for any external potential. To illustrate modeling external forces in the CPM we set up a model with a row of three cells and define a potential for right to left on the rightmost cell. We expect that the cell on the right and the cell in the middle will elongate because of this force. Figure 11 shows the evolution of the cells and the cell length during the simulation. As expected, the pushing force towards the left causes the cells to elongate. After $\sim $600 steps the cell the right starts crawling over the middle cell and migrates further to the left. As a result, the rightmost cell no longer pushes on the other two cells and they relax. Note that the exact moment of at which the rightmost cell starts crawling over its neighbor differs per simulation, therefore we stopped measuring the cell length after 600 steps.

Recently, several efforts have been made to extend the CPM with a heterogeneous ECM that interacts with cells and that is remodeled by the cells. Overall, three different approaches can be recognized. In the first approach the ECM is modeled as a density field that affects the behavior or the cells and is remodeled by the cells. Daub and coworkers [117] extended the CPM with both haptotaxis, migration towards higher ECM concentrations, and haptokinesis, preferential migration on intermediate levels of ECM, using an ECM density field. The cells affected the field by depositing matrix metalloproteinases (MMPs) that degrade the ECM: $\frac{\partial c(\underline{x},t)}{\partial t} = - \delta (\sigma (\underline{x}),0)\varepsilon c_\text {MMP}(\underline{x},t)c(\underline{x},t)$, with $c_\text {MMP}$ the MMP concentration and $\varepsilon $ the ECM degradation rate. Both haptotaxis and haptokinesis affect the dynamics of the cells, therefore they are added as changes in effective energy. The energy change due to haptotaxis is described as:

$$\begin{aligned} \Delta E_\text {haptotaxis} = \chi \left( \frac{c(\underline{x}^\prime )}{1+s\;c(\underline{x}^\prime )} -\frac{c(\underline{x})}{1+s\;c(\underline{x})}\right) , \end{aligned}$$

(23)

where $\chi $ denotes the sensitivity, $c(\underline{x})$ the ECM concentration ($c \in [0,1]$), and s the saturation coefficient. For haptokinesis a reverse Gaussian function is used to model the response of a cell to the local ECM concentration:

$$\begin{aligned} \Delta E_\text {haptokinesis} = \eta \delta (\sigma (\underline{x}),0)\left( 1-\frac{1}{\rho \sqrt{2\pi }} e^{- \frac{(\mu (\underline{x}^\prime )- \bar{\mu })^2}{2\rho ^2}} \right) , \end{aligned}$$

(24)

with haptokinesis strength $\eta $, and $\mu $ and $\rho $ the mean an standard deviation of the Gaussian function. The second approach for modeling cell–ECM interactions in the CPM is by including immobile objects that represent fibers [105, 115, 116, 124–126]. The cells adhere to the fibers and thereby the fibers facilitate migration. However, because the fibers are modeled as immobile objects, they may also block migration. In some of the models, cells are also able to degrade the fibers by overtaking their lattice sites [126]. The third method of modeling cell–ECM interactions is by mechanically linking the ECM to the cells. For this the ECM is modeled as a 2D compliant substrate using the finite element method (FEM) [25] and each CPM step is followed by a FEM step that resolves the ECM deformation. The CPM lattice is placed on top of the FEM mesh such that one CPM lattice site corresponds with one FEM node. Each FEM node that is covered by a cell pulls on all other nodes that are part of the same cell [25, 127], causing a resultant force $\underline{F}_i$ on each node i:

$$\begin{aligned} \underline{F}_{i} = \mu \sum ^{j \ne i}_{\forall j : \sigma (j) = \sigma (i)} \underline{d}_{ij}, \end{aligned}$$

(25)

with $\mu $ the tension per unit length, and $\underline{d}_{ij}$ the Euclidean distance between node i and j. Then, during the FEM step, the substrate deformation is computed according to: $K\underline{u} = \underline{F}$ with stiffness matrix K, displacement $\underline{u}$ and forces $\underline{F}$. The deformation causes strain stiffening, which is expressed in the Young’s modulus $E(\epsilon )$. Strain stiffening is fed back to the CPM via durotaxis: the preferential migration of cells towards areas of higher stiffness:

$$\begin{aligned} \Delta E_\text {durotaxis}= & {} -g(\underline{x},\underline{x}^\prime )\lambda _\text {durotaxis}\left( h(E(\epsilon _1))(\underline{v}_1\cdot \underline{v}_m)\right. \nonumber \\&\left. +\,h(E(\epsilon _2))(\underline{v}_2\cdot \underline{v}_m) \right) \end{aligned}$$

(26)

with $g(\underline{x},\underline{x}^\prime ) = 1$ for extensions and $g(\underline{x},\underline{x}^\prime ) = -1$ for retractions, $\lambda _\text {durotaxis}$ the magnitude of the response to stiffness, h(E) a sigmoid representing the shape of response of the cells to the stiffness, $\epsilon _1$ and $\epsilon _2$, and $\underline{v}_1$ and $\underline{v}_2$ the principal strains and strain orientation, and $\underline{v}_m = \widehat{\underline{x}-\underline{x}^\prime }$ is the unit vector in direction of $\underline{x}-\underline{x}^\prime $.

As described, cell-based models have been linked to the ECM in several ways using the CPM. The main reasons that this has happened with the CPM is that the method is computationally inexpensive in 2D and quasi-3D. Therefore, it is possible to perform large-scale parameter sweeps and sensitivity analyses with CPM-based models [128, 129]. Furthermore, open source software CPM modeling is readily available [130–132] which makes it easy to implement and simulate models using the CPM. However, there are some major disadvantages to the CPM. First, whereas the CPM performs well in 2D and quasi-3D, real 3D simulations take long and are therefore uncommon. Second, because the CPM minimizes the effective energy, all cell behavior is highly linked to cell motility. It is, for example, not possible to simulate a immobile, flexible cell using the CPM (see also the discussion above). Finally, some parameters of the CPM cannot be linked directly to physical parameters. Therefore, physical parameters measured in experiments cannot always be used to set model parameters in the CPM. One could derive the Poisson ratio and Young’s modulus from single-cell simulations in which the cell is compressed or stretched [133]. By calibrating the model such that the Poisson ratio and Young’s modulus measured in simulations correspond with those measured in experiments, the single-cell parameters can be estimated. This leaves only the contact energies between cells as an unknown parameter, which can only be estimated by calibrating the CPM simulations against experimental observations. Altogether, whereas the CPM offers a flexible framework for cell-based simulations that include interactions between cells and the ECM, the method is limited to simulating motile cells and is mainly applicable for qualitative modeling. Below we briefly summarize advantages and disadvantages of type C cellular automaton models (Cellular Potts models) from our perspective. Evaluations on the computation times are done with regard to a standard PC.

3 Off-lattice models

In off-lattice models, cells are represented by a single or a clusters of particles and interactions between them can then be described by forces or potentials. As in the lattice models, cells have the ability to grow, migrate, divide, and die. Position changes can be obtained by solving an equation of motion for each cell. Alternatively, the dynamics of a system of cells can be mimicked using energy-based methods using numerical procedures such as Monte Carlo sampling and the Metropolis algorithm [61, 82, 134], which is used also for the CPM (see previous section). The advantages of force-based models are a well-defined time scale, and a more intuitive way of taking into account complex interactions of cells with other cells or their environment which is why they became the standard approach. For this reason, we mainly limit ourselves to force-based models. Energy-based modeling with CBM are discussed in detail in Ref. [73].

3.1 Center-based models

3.1.1 Basic concepts

In CBM cells are represented by simple geometrical objects that can be described by one or a small number of centers. The basic assumption is that each trajectory of a cell in space can be described by an equation of motion in formal analogy to physical particles. Usually inertia effects are neglected as the Reynolds numbers are very small meaning friction dominates the system [135]. The basic equation considering forces between cells and between cells and substrate then reads for cell i

$$\begin{aligned} \varGamma ^\mathrm{cs}_{i} \underline{v_i} \!+\! \sum _{j}\varGamma ^\mathrm{cc}_{ij} (\underline{v_i}-\underline{v_j}) \!=\! \sum _{j}\underline{F}^\mathrm{int}_{ij} + \underline{F}^\mathrm{sub}_{i} + \underline{F}^\mathrm{mig}_{i}. \end{aligned}$$

(27)

On the left side, we have the friction forces of cells with substrate (first term) and among cells (second term), on the right side of this equation, we find the forces on the cells emerging from repulsion/adhesion with other cells j, $\underline{F}^\mathrm{int}_{ij}$, as well as with the substrate $\underline{F}^\mathrm{sub}_{i}$, as well as the cell migration force $\underline{F}^\mathrm{mig}_i$. The friction forces involve tensors for the cell–cell friction ($\varGamma ^\mathrm{cc}_{ij}$) and cell–substrate friction ($\varGamma ^\mathrm{cs}_i$). The latter can also capture capillaries or membranes [33]. If friction coefficients parallel and normal to the movement direction are different, the friction tensor becomes

$$\begin{aligned} \varGamma _{ij} = \gamma _{\perp } (\underline{u}_{ij}\otimes \underline{u}_{ij}) + \gamma _{||} (I - \underline{u}_{ij}\otimes \underline{u}_{ij}), \end{aligned}$$

(28)

with $\underline{u}_{ij}=(\underline{r}_j-\underline{r}_i)/|\underline{r}_j-\underline{r}_i|$, where $r_i$, $r_j$ denote the position of the centers (see Fig. 12) of cell i and object j. ($\otimes $ denotes the dyadic product. If object j denotes another cell, this would correspond to the superscript “cc”, otherwise to the superscript ”cs”. For example if cells move on a flat substrate (s), the substrate can be viewed as a sphere with infinite radius.) Here $\gamma _{\perp }$ denotes perpendicular, $\gamma _{||}$ parallel friction. This can be seen after multiplication of the friction tensor with the cell velocity:

$$\begin{aligned} \varGamma _{ij}\underline{v}_i= & {} \left( \gamma _{\perp } (\underline{u}_{ij}\otimes \underline{u}_{ij}) + \gamma _{||} (I - \underline{u}_{ij}\otimes \underline{u}_{ij})\right) \underline{v}_i\\= & {} \gamma _{\perp } \underline{u}_{ij}(\underline{u}_{ij}\underline{v}_i) + \gamma _{||} (I\underline{v}_i - \underline{u}_{ij}(\underline{u}_{ij}\underline{v}_i)). \end{aligned}$$

Friction occurring in the perpendicular direction of the movement may be associated with internal friction in particular if two cells are pushed against each other which leads to deformation and reorganization of the CSK, while friction in the parallel direction is a description of the dissipative forces when sliding the cell membranes along each other.

The migration force $\underline{F}^\mathrm{mig}_i$ is usually an active force. In absence of influences that impose a direction it is common to assume that the migration force is stochastic, with zero mean and uncorrelated at different points in time, $\langle \underline{F}^\mathrm{mig}_i\rangle = 0$, $\langle \underline{F}^\mathrm{mig}_i(t)\otimes \underline{F}^\mathrm{mig}_i(t')\rangle = A_i \delta (t-t')$. Here $A_i$ is a $d \times d$ matrix that may be specific for each individual cell i, d being the spatial dimension. Using the formal analogy to colloidal particles in suspension, one obtains for the fluctuation strength

$$\begin{aligned} \langle \underline{F}^\mathrm{mig}_{||;i}(t)\underline{F}^\mathrm{mig}_{||;i}(t')\rangle= & {} G_{||;i} \delta (t-t'), \end{aligned}$$

(29)

$$\begin{aligned} \langle \underline{F}^\mathrm{mig}_{\perp ;i}(t)\underline{F}^\mathrm{mig}_{\perp ;i}(t')\rangle= & {} G_{\perp ;i} \delta (t-t'), \end{aligned}$$

(30)

with $\underline{F}^\mathrm{mig}_{||;i}(t)=(I - \underline{u}_i\otimes \underline{u}_i)\underline{F}^\mathrm{mig}_i$, $\underline{F}^\mathrm{mig}_{\perp ;i}(t)=(\underline{u}_i\otimes \underline{u}_i)\underline{F}^\mathrm{mig}_i$. $G_{||;i}$ and $G_{\perp ;i}$ are scalars. $\underline{u}_i$ is the unit vector pointing into the direction of the cell movement. In colloidal particle physics, the autocorrelation amplitude of the noise is controlled by the friction coefficient and the thermal energy $k_\mathrm{B}T$, $k_\mathrm{B}$ being the Boltzmann constant, T temperature, according to $G_{||;i}=4\gamma _{||}k_\mathrm{B}T$, $G_{\perp ;i}=2\gamma _{\perp }k_\mathrm{B}T$. A formal analogy would lead to replacing $k_\mathrm{B}T$ by the membrane fluctuation strength $F_\mathrm{T}$ as proposed by Beysens et al. [87]. However, the mechanisms of cell movement cannot be expected to follow the equipartition theorem of statistical mechanics. Cell migration is active, depending on the local matrix density and orientation hence we cannot claim the autocorrelation amplitude matrix A is properly described by the aforementioned analogy to colloidal particle movement. Instead, by measuring $\langle ((\underline{r}_i(t+\tau )-\underline{r_i}(t))\otimes (\underline{r}_i(t+\tau )-\underline{r}_i(t))\rangle = 2D\tau $, one might determine the diffusion tensor D of the cell where again diffusion is not a consequence of collisions with smaller particles, but a parameter phenomenologically describing the cells’ spread. For a simple algorithm to simulate a force-based Brownian motion of a cell in a isotropic medium, see Appendix 2.

Some authors use the inertia term $m_i\text {d}\underline{v}_i/\text {d}t$ to mimic the effect of persistence [28] while in principle it can be directly simulated using a memory term. In this case, $\varGamma _{ij}^{cx} \underline{v}_k$ $(x\in \{c, s\}$, $k \in \{i, j\})$ has to be replaced by $\int _{-\infty }^{\infty } \varGamma _{ij}^{cx}(t-t') \underline{v}_k(t') \text {d}t'$, so that if $\varGamma _{ij}^{cx}(t-t') = \varGamma _{ij}^{cx}\delta (t-t')$, the basic equation of motion is recovered. Another way is to assume that the active movement component of a cell has a memory i.e., a cell has the tendency towards keeping its (active) migration direction. This can be modeled by replacing the uncorrelated noise by noise for which the auto-correlation amplitude decreases over a time period over which the cell has migrated of a distance of the order of its diameter or even larger. Sepulveda et al. [136] modeled this by assuming the active random migration force $\underline{F}^\mathrm{mig}_i$ follows an Ornstein-Uhlenbeck process modeled by the solution of the equation $(1/\beta )\text {d}\underline{F}^\mathrm{mig}_i/\text {d}t = -\underline{F}^\mathrm{mig}_i + \underline{f}^\mathrm{uc}_i$ were $\underline{f}^\mathrm{uc}_i$ again has zero mean and is uncorrelated, $\beta $ is the inverse of the autocorrelation time.

Directed cell migration was successfully modeled in Hoehme et al. [33] considering the closure of a necrotic lesion after drug induced damage in liver. Firstly, by a chemotaxis force $\underline{F}_i^\mathrm{chem} = \chi \nabla c$ was assumed with $\chi $ being the chemotaxis coefficient, and $c(\underline{r},t)$ the morphogen concentration assuming a morphogen being secreted by the necrotic cells. A second assumption was that cells attempt to move away from cells of the same type, modeled by the term $\underline{F}^\mathrm{mig}_i = (1-\theta (\nabla p_i\underline{f}^\mathrm{uc}_i))\underline{f}^\mathrm{uc}_i$. Here $\theta (x)=1$ if $x\ge 0$, else $\theta (x)=0$, is the Heaviside step function, $\underline{f}^\mathrm{uc}_i$ again uncorrelated white noise. $p_i$ is a pressure-like measure defined below (Eq. 36). The effect of this formula is, that it inhibits all active moves that would increase cell density as for those the Theta-function becomes one so the migration force becomes zero. Only active moves leading to a decrease of density can occur. Hence the expression reinforces cells to move towards empty spaces.

The cells interact by pairwise potentials having a repulsive and adhesive part, and are described by a function of the geometrical overlap $\delta _{ij}$. A number of different approaches have been used for the interaction force, including linear springs [34, 134, 137], a force derived from Lennard–Jones like potentials [138], and a Hertz force approximating cells by isotropic homogeneous elastic bodies that are moderately deformed if pressed against each other [63, 82]. We here consider the latter approach and extensions based upon it. Consider two spherical cells i and j which are in each others neighborhood (see Fig. 12). If the distance $d_{ij}$ between the cells becomes smaller than the sum of their radii, we assume a Hertzian contact force develops calculated by:

$$\begin{aligned} F^\mathrm{rep}_{ij}= {4/3 \hat{E}}\sqrt{\hat{R}}\delta ^{3/2}. \end{aligned}$$

(31)

in which $\delta $ is the overlap between the cells and $\hat{E}$ and $\hat{R}$ are defined as

$$\begin{aligned} \hat{E} = {\left( \frac{1-\nu _i^2}{E_i} + \frac{1-\nu _j^2}{E_j} \right) }^{-1} \; \text {and}\; \hat{R} = {\left( \frac{1}{R_i} + \frac{1}{R_j} \right) }^{-1} \end{aligned}$$

with $E_i$ and $E_j$ being the Young’s moduli, $\nu _i$ and $\nu _j$ the Poisson numbers and $R_i$ and $R_j$ the radii of the cells i and j, respectively. Galle et al. [63] and Ramis-Conde et al. [65] considered adhesion of cells to other cells or an underlying substrate by an additive term in the interaction energy which is proportional to the Hertzian contact area of the cells without adhesion and the strength of the adhesion formed by bonds $\sigma $:

$$\begin{aligned} F^\mathrm{adh}_{ij}= - \pi \sigma \hat{R}, \end{aligned}$$

(32)

and named the force “extended Hertz model”. Note that the total force $F^\mathrm{int}=F^\mathrm{rep}+F^\mathrm{adh}$ acts in the direction of the vector connection the centers of the spheres. This simplified approach is convenient as it permits to give an explicit analytical expression for the force. However, it neglects that adhesion is modifying the contact area, and disregards that the force distribution within the contact area are is inhomogeneous. A more realistic adhesive model taking this into account results from the Johnson–Kendal–Roberts (JKR) theory, which also takes into account a hysteresis effect if the cells are separated from each other. In this case the force is computed by

$$\begin{aligned} F^\mathrm{adh}_{ij}= \frac{4 \hat{E}}{3 \hat{R}} \left[ a(\delta _{ij})\right] ^3 - \sqrt{ 8 \pi \sigma \hat{E} \left[ a(\delta _{ij})\right] ^3}. \end{aligned}$$

(33)

The contact radius a in Eq. 33 must be obtained (from [139]):

$$\begin{aligned} \delta _{ij} = \frac{a^2}{\hat{R}} - \sqrt{\frac{2 \pi \sigma }{\hat{E}} a}. \end{aligned}$$

(34)

Compared to the JKR model, the modified Hertz model requires less computational effort, because Eq. 34 needs to be solved iteratively. On the other hand, the JKR model is a more realistic adhesion model for solid adhesive soft spheres for the reasons explained above. A study of multi-cellular dynamics in monolayers reported an effect of the adhesion model choice [62].

Once all forces are determined in Eq. 27, we arrive at a linear problem described by a sparse symmetric matrix, which can be solved efficiently by a Conjugate Gradient method [140, 141]. Note that the above system of equations of motion (Eq. 27) does not conserve total momentum. It assumes a substrate on which momentum can be transfered on is a petri dish in case of growing monolayers, or a stiff but not too dense ECM. If cells migrate over each other or if the ECM is not rigid enough, this assumption might be problematic. In the latter case of a soft ECM long-range communication between the cells mediated by the ECM could become important and should be taken into account explicitly. Models taking into account ECM are briefly discussed below.

In the cycle a cell increases its volume, which can be described by

$$\begin{aligned} \frac{dV_{i}}{\text {d}t} = \int ^{T_{i}}_{0}{\alpha _{i}(V_i, \ldots ) \text {d}t}. \end{aligned}$$

(35)

with a volume $V_{0;i}$ right after mitosis has occured, $T_{i}$ the cell cycle time, and $\alpha _{i}(V_i, \ldots )$ the volume growth rate. The latter one cannot be assumed to be constant over the cell cycle and maybe altered under the conditions of inhibited growth, and for example depend on the actual volume $V_i$ and many other parameters. If $\alpha _i$ is a constant, cell volume increase is linear. If the cell passed a critical volume, critical time point, or, in case a cell is modeled by a dumb-bells, a critical dumb-bell axis length, the cell undergoes mitosis, the process during which a cell splits into two daughter cell bodies, and two new cells with equal volume are created. During the mitosis process several internal forces develop in the mother cell, involving an actin contractile ring at the cell equator [142]. Different models of cell division have been studied [62]. The interphase adding up the $G_1$, S and $G_2$ phase and mitosis phase may be described by an increase the radius of a spherical cell, eventually transforming it into a dumbbell during mitosis. As dumbbells lack spherical symmetry an additional equation of motion for rotation is required. The equation for the angular momentum is complex [83], and therefore rotations are often simulated by a Monte-Carlo dynamics using the Metropolis algorithm [33]). A simplification of the division algorithm consists in skipping the transformation phase and placing two smaller daughter cells in the space originally filled by the mother cell at the end of the interphase (e.g., [63, 143]). With this algorithm no dumbbells are generated, hence rotations of to dumb-bell shapes does not occur. When the two daughter cells are created, the two daughter cells will be pushed from each other away until mechanically in equilibrium. The direction in which the cell divides may be chosen randomly, according to biological principles or the cells’ environment, or according to stress principles (see further). However, the simplified algorithm has two major shortcomings. Firstly, if the space filled by the mother cell is small already, which is often the case for cells in the interior of a cell population, the local forces occurring after replacing the mother cell by two spherical daughters can adopt very large (unphysiological) values leading to unphysiologically large cell displacements. This might partially be balanced by intermediately reducing the forces between the daughter cells. Secondly, cells in contact to the mother cell prior to replacing the mother cell by its two smaller daughter cells may loose contact as a consequence of the discontinuous local configuration change. For a cell at the border of a population this can lead to a detachment of the neighbor cell from the multicellular aggregate.

In another approach to cell division, cells entering the cell cycle are immediately represented as dumbbells which gradually grow by increasing their axis from zero length to twice the radius of the dumbbell spheres, where they separate into two spheres [84, 138]. This algorithm has been found to generate only minor differences on time scales much larger than the cell cycle time compared to the algorithm distinguishing between a spherical growth and a dumbbell deformation phase.

During the cell cycle, cells take several decisions depending on certain criteria. Commonly used criteria for cell cycle progression control by contact inhibition on the whole cell level are cell deformation, volume compression, a force threshold or a stress threshold, each of which have proven to qualitatively predict growth limitation in tumor models [144]. A machinery of molecular factors underly the decision processes (e.g., [145, 146]). Drasdo and co-workers considered in a number of papers [61, 82, 147] a mechanism in which progression of a cell in the cell cycle was possible only if the distance to any of its neighbor did not become smaller than a certain threshold value. For a cell in the interior of a multicellular aggregate, balance of forces usually leads to small distances of a cell to its neighboring cells, so a too large overlap is likely to correspond to a too large volume compression of that cell. Cells at the border can move and escape the compression which is why border cells were not limited by this criterion. Contact inhibition of growth, as well as other growth control mechanisms mark the difference between normal and abnormal cells. Aggressive cell lines pile markedly up in monolayer culture i.e., growing on a flat substrate, as their are able to grow and divide without contact to the underlying substrate [62, 63]. For example, in the absence of contact inhibition, cells localized in the interior of the monolayer may above a critical size not be able anymore to relax their mechanical stress by pushing other cells towards the monolayer border. Small stochastic differences in the cells’ distance from the underlying growth plane may then lead to normal force components, which for some of the cells are directed out-of-plane (Fig. 13). If those force components are not balanced by adhesion forces between substrate and these cells anymore, then these cells can be pushed out of the layer. CBMs capture this effect very naturally as it emerges out of the force-based dynamics, while in a cellular automaton model such an effect needs to be coded explicitly and is otherwise missed out.

The criterion that cell cycle progression is inhibited at any time in the cycle leaving cells in the interior of a multicellular aggregate arbitrarily long in the cell cycle does not correspond to biological observation. Rather a number of checkpoints are reported and generally accepted (e.g., $G_1$, $G_2$, M-checkpoint). Cells passing $G_1$-checkpoint are committed to pass the entire cell cycle. Hoehme et al. [33, 84] implemented this by assuming that cell cycle entrance was possible only after division in a certain time point or small window in $G_1$ phase mimicking the restriction point. Once a cell passes the restriction point it progresses until division with probability one. This condition gave a very good quantitative explanation of many observations in growing multicellular spheroids of EMT6/Ro cells [84]. A similar choice had been taken by Schaller and Meyer-Hermann [143]. Schaller and Meyer-Hermann [143] considered cell-cycle entrance if the local pressure was below a critical value. Drasdo and Hoehme [144] also compared the consequence of different cell cycle entrance criteria such as a compression threshold, a tension threshold (meaning cells enter cell cycle only if stretched), and different mechanical stress-related criteria, on monolayer growth (Fig. 14).

In Galle et al. [63] a volume-based criterion is used. The cell compression is computed using the target volume minus the volume loss from geometrical overlaps of neighboring cells, $V_i(t) = (4/3)\pi R_i(t)^3 - \sum _{j\;NN\;i} V_{ij}^\mathrm{cap}$ i.e., the caps lying in the neighbor cell. The overlap is overestimated if a cell overlaps simultaneously with two other cells that already overlap i.e., if three cells share a volume segment.

The different criteria influence the expansion speed and the cell density profile but seem to reproduce qualitatively the same dynamics.

There are several possibilities to compute stresses in CBM. For mechanical stress, one can define several pressure, or pressure-related quantities. Byrne and Drasdo [147] used the following definition to incorporate contact inhibition effects in tumor growth:

$$\begin{aligned} p_{i}=\sum _{j}{\frac{||\underline{F}^\mathrm{rep}_{ij}\underline{u}_{ij} ||}{A_{ij}}}, \end{aligned}$$

(36)

with $\underline{u}_{ij}=(\underline{r}_i -\underline{r}_j)/|\underline{r}_i-\underline{r}_j|$. $\underline{F}_{ij}^\mathrm{rep} =F_{ij}^\mathrm{rep}\underline{u}_ {ij}$ are the repulsive forces from cells j on cell i due to contacts with cells j and $A_{ij}$ is the contact area derived from the contact model. This measure associates every deformation or compression of the cell with a (positive) pressure, including cell deformation occurring as a consequence of adhesion. The measure was calculated as pressure equivalent in a growing monolayer, where cells are constantly under compression hence the cell-to-cell distance $\langle d_{ij}\rangle <R_i+R_j$.

Schaller and Meyer-Hermann [143] considered instead

$$\begin{aligned} p_{i}=\sum _{j}{\frac{||\underline{F}_{ij}^\mathrm{int} \underline{u}_{ij} ||}{A_{ij}}}, \end{aligned}$$

(37)

i.e., adhesive and repulsive forces in growing multi-cellular spheroids. In this case, the pressure would be zero if the total interaction force between a cell and its neighbors is zero, so the obtained values are always below those in the former definition in presence of cell cohesion. However, as the latter measure was considered in growing cell populations with all cells permanently under compression, it differs only by a pressure shift. Defining instead [144]:

$$\begin{aligned} p_{i}=\sum _{j}{\frac{\underline{F}_{ij}^\mathrm{int}\underline{u}_{ij}}{A_{ij}}} \end{aligned}$$

(38)

takes into account tensions that can occur if a cell stretches its cell–cell contacts during its attempt to actively migrate away from the cell aggregate it is attached to, which has been one of the cases studied in that reference in addition to the aforementioned case of cells permanently being under compression. This study was motivated by the experimental findings of Trepat et al. [18] suggesting that border cells migrate actively into their environment pulling interior cells behind. Reference [144] studied if the experimental observations could be explained if cell cycle entrance could occur only for cells under negative pressure, here defined such that the cell contacts are subject to stretch.

For growing monolayers with all cells under compression, all three measures differ only by moderate shifts of the pressure curves. In the first and second definition, $p_i \ge 0$. Using Eq. 37 has the disadvantage that cells for which are slightly compressed so that the interaction force is repulsive and cells for which the interaction force is negative due to adhesion can have the same pressure value. In the previous definition given by Eq. 36, the force is by definition always positive so this problem cannot occur. The disadvantages of all these measures is that they only give a scalar measure but do not allow to calculate the entire stress tensor. Moreover, they cannot easily quantitatively compared to the Cauchy stress tensor obtained in continuum theories.

An alternative estimation of stress can be obtained by introducing the concept of virial stress, a well-known concept in particle-based simulations. It was drawn from theories in molecular dynamics to compute macroscopic stress originating from the interaction of many particles [148]. Compared to the abovementioned definitions, more information can be extracted using the viral stress formula. The general formulation of virial stress defines the average stress^{Footnote 1} in a box of particles due to contact forces $F_{ij}^\mathrm{int}$ with other particles by:

$$\begin{aligned} \sigma _{i}= \frac{1}{V}\sum _{j}\left( \underline{F}_{ij}^\mathrm{int}\otimes \underline{r}_{ij}\right) \end{aligned}$$

(39)

where $\underline{r}_{ij}$ is the vector pointing from the center of cell i to the contact plane with cell j, and V is the sampling volume. In principal this volume needs to be large enough to obtain reliable averages over an ensemble of particles, yet one can still use it as an estimator for individual cell stresses, by taking the sampling volume as the volume of the cell. The hydrodynamic pressure $p_{i}$ on a cell can now be computed by:

$$\begin{aligned} p_{i} = \frac{1}{3}\text {tr}(\sigma _{i}). \end{aligned}$$

(40)

Equation 39 defines a full stress tensor for the cell, which will allow to extract more information than just a pressure. An example and the applicability of this method is further demonstrated in Appendix 1. From the stress tensor, we can also derive the principal stress directions $ \left\{ n_1, n_2, n_3 \right\} $ and eigenvalues $ \left\{ \sigma _1, \sigma _2, \sigma _3 \right\} $ by diagonalization of Eq. 39, yielding important information on how stresses on the cells are oriented. The algorithm for diagonalization of Eq. 39 is trivial in 2D and not demanding in 3D, using e.g., Cardano’s method.

3.1.2 Achievements, limitations and practical use

CBM now exist for several decades and despite their simplicity, they have been used with quite large success. Typical application fields are studies 2D monolayers (Fig. 14), (small) tumors (Fig. 14), intestinal crypts (Fig. 18), and liver (e.g., liver lobules, Fig. 15), with cell numbers starting from few and ranging up to $10^6$ [83, 150]. Refined image analysis techniques and available software increasingly permit a quantitative analysis of tissues in 3D at micro-meter resolution [38, 151]. The data infered can subsequently be transfered into modeling software such as CellSys [150], to simulate the tissue organization process of interest as demonstrated for the example of the regenerating liver (Fig. 15, [33, 152]). The experimental technology improvements regarding confocal microscopy and live microscopy now provides the necessary data on the histological level to calibrate agent-based models (see also [81]).

Although the mechanical information that can be extracted from these models is rather approximate, various biomechanical aspects in monolayers, tumor spheroids and organ tissues could be explained. For example, Drasdo and Hoehme [82] have shown that contact inhibition due to mechanical stress can explain the equal expansion speed at significantly different glucose medium concentrations in growing multicellular spheroids of EMT6/Ro cells. The same mechanism could explain different growth speeds in expanding one-cell thick dense monolayers observed by Bru et al. [90], and growth saturation of multicellular spheroids in agarose gel [62, 82, 144] observed by Helmlinger et al. [11], and Galle et al. [153]. Recently Delarue et al. [154] confirmed inhibition of proliferation in tumor spheroids by compressive stress.

Bassan et al. [138] simulated tissue rheology by measuring yield stresses, shear viscosities, complex viscosities as well as the loss tangents as a function of model parameters, concluding that cell division and apoptosis lead to a vanishing yield stress and fluid-like tissues. Macklin et al. [155] have addressed ductile carcinoma’s. A number of researcher has also been focused on intestinal crypts and epithelial tissues, studying e.g., cell organization and the role of basement membranes, and predicting the behavior of the tissue during steady state as well as after cell damage [35, 37, 52, 156–163].

An import aspect in cell migration is the interaction with ECM, as it is responsible for both the support and obstructing cell movement. For instance, in invasive cancers one observes that migrating cells breakdown the ECM, forming guidance tunnels for other cells. In addition, distant cells can feel one the other’s presence through the modified strain fields that they cause. In the basic CBM formulation the interaction with ECM is usually mimicked in an implicit fashion using friction terms and random motility forces. A 3D migration model was proposed by Zaman and co-workers [164, 165], based on internally generated forces transmitted into external traction forces (and considering a timescale during which multiple attachment and detachment events are integrated). They investigated the migration speed as a function of properties like ligand density of the adhesion receptors and ECM stiffness, and found configurations where the cells would optimize their migration speed. A similar approach was used later, by Rey and García-Aznar [166] to study cell migration in wound healing. Vermolen and Gefen [167] proposed a model with cell movement assuming a mechanical stimulus arises from the cells sensing a change in strain energy density in the ECM. An explicit interaction model of CBM with ECM fibers has been introduced by Schlueter et al. [66] in where the cells are able to interact and remodel the fiber orientation, showing the tendency of a cell to move to stiffer substrates. Harjanto et al. [168] used an approach where the cells can degrade, deposit, or pull on local fibers, depending on the fiber density around each cell. We note that the fibers themselves in principle are not inert stiff structures but have their specific dynamics which can be modeled as well using e.g., Brownian dynamics [169].

To address collective motion of cells, Bassan et al. [170] proposed a simple flocking-type mechanism, in which cells tend to align their motility forces with their velocity. Their model could explain the experimentally observed long-range alignment of motility forces in highly disordered patterns, as well as the buildup of tensile stress throughout the tissue. Also Sepulvada et al. [136] considered a CBM in which cells move in a stochastic manner and try to adapt their motion to that of their neighbors. Introducing leader cells, they found that fingers develop as observed in experiments when leader cells more actively invade free environment than following cells and regulate their motion according to their contacts with following cells.

Several extensions have been proposed to alleviate the limitations of CBM. To model more realistic cell shapes, Palsson et al. [171, 172] and Dallon et al. [173] used ellipsoidal shapes with axial degrees of freedom resulting in deformability in three principal directions. Herein the deformability was controlled by viscoelastic elements. The shape was adapted according to the forces that act on the cell. However, as the viscoelastic elements in these models are placed on the three axes of the ellipsoid, a cells’ deformationdepends on the relative orientation of the axes of the ellipsoids compared to the direction of the interaction force in the contact region of the two interacting ellipsoids. This can be seen as follows (see Fig. 16). Assume two cells that in the moment where they get into contact have a spheroidal shape, meaning all axes have the same length. If the cell–cell interaction force acts along the axes of two interacting ellipsoids, they deform easily. If on the other hand the cell axes are rotated by 45$^{\circ }$ compared to the direction of the interaction force, they would not or at most moderately deform (Fig. 16b).

Anisotropy and polarity in cells can be incorporated into the CBM in a rather simple fashion by assigning a polarity vector $\underline{P}$ to each cell. The latter indicates that a cell has a preferential direction, polar adhesion or biased motion [33, 66, 137]. Non-adherent neighboring cells can sense each others presence and move towards each other through mechano-chemical signals, originating for instance by the formation of filopodia or cell–cell communication through a medium. To incorporate such effects in the simulations, one can simply extend the range of cell–cell interaction force using a function that represents the intensity of this attraction as a function of the separation distance. Palsson et al. [172] proposed to characterize this attraction by a Gaussian:

$$\begin{aligned}&F_{ij}^\mathrm{int}\nonumber \\&\quad \!=\! {\left\{ \begin{array}{ll} -\pi \sigma \sqrt{\hat{R}} &{} \text {if}\; d_{ij} \ge R_i+R_j \\ -\pi \sigma \sqrt{\hat{R}}\exp \left( -\frac{(R_i\!+\!R_j\!-\!d_{ij})^2}{2b^2}\right) &{} \text {if}\; R_i\!+\!R_j \!<\! d_{ij}\! \le \! d_\mathrm{cutoff} \end{array}\right. } \nonumber \\ \end{aligned}$$

(41)

The parameter b represents the spreading of the filopodia length while the maximal interaction range is $d_\mathrm{cutoff}$. Another possibility are Morse potential functions, which were used by e.g., Rey et al. [166]. An artifact can occur if the interaction range mimicking the filopodia-mediate attraction is so large that two cells attract each other even if the local cell configuration does not leave free space for the filopodia to be stretched out from one to the other cell. In this case, explicit representation of filopodia even if computationally costly, may be necessary, as e.g., performed in Ref. [174] who modeled filopodia with a CBM by lines radially stretched from the cell border in such a way that they may touch but not intersect with neighbor cells (Fig. 17).

A drawback in CBM based upon pair-wise forces is that the contact forces and contact area become largely inaccurate when cells become densely packed. To understand this, consider a cell surrounded by many other cells in a dense packing as it may occur if neighbor cells grow. The volume of the central cell can be estimated by local Voronoi tessellation, and the calculation of the volume of the Voronoi polygon associated with the central cell (for a Voronoi tessellation, see Fig. 1). Consider the case of an incompressible homogeneous isotropic elastic cell. Incompressibility would correspond to a Poisson ratio of $\nu = 0.5$. Because the central cell is incompressible, even in the most dense arrangement of cells, its volume (let’s denote it by $V_0$) would remain the same all the time. This is a consequence of the volume control that generates a multibody interaction of the central cell with each of its neighbors. In contrast to this, in simulations the pairwise interaction forces (e.g., Hertz without adhesion, JKR in presence of cell–cell adhesion) do not include any notion of cell volume. Hence, even for the setting $\nu = 0.5$ in the Hertz (or JKR) force models, the cell–cell distances of the central cell with its neighbors could become so small, that the volume that can be associated with the central cell (e.g., by Voronoi tesselation) would become smaller than $V_0$. One way to account for the volume is adding a force term that takes into account a volume constraint. In this case the cell volume may be approximated by a Voronoi tessellation which, however, generates some other problems that are discussed in more detail below.

Another, related, point is that in a dense cell packing the Hertz contact area (or in presence of adhesion the JKR contact area) does not represent a good estimate of the contact area anymore, as in a dense packing the contact regions of the central cell with those of the neighbor cells, that are themselves neighbors of each other, overlap. In this case the interaction between the neighbors of the central cell leads to deformations of the neighbor cells, which impact on the contact area between the neighbors and the central cell. The contact area may then better be estimated by the geometrical area in the contact region, calculated again from a Voronoi tessellation. However, there is no evident smooth consistent transition between Hertz contact area at cell configurations with low densities, and geometrical area for cell configurations at high densities. To better see this, consider the case of no cell–cell adhesion i.e., a Hertz interaction force. The Hertz interaction area for two cells is $A_{ij}^\mathrm{H} = \pi (R_i+R_j-d_{ij})(R_iR_j)/(R_i+R_j)$. The geometrical contact area for two interacting cells, which is the area the two overlapping cells have in common, is $A_{ij}^\mathrm{g}=\pi R_i^2-(d_{ij}^2+R_i^2-R_j^2)/(2d_{ij}^2)$. In both cases, $d_{ij}$ is the distance of the centers of the two cells i, j. The Hertz interaction area is always smaller than the geometrical interaction area. For example in case $R_i=R_j\equiv R$ one finds $A_{ij}^\mathrm{H} = \pi (2R-d_{ij})R/2 < A_{ij}^\mathrm{g}=\pi (R^2-d_{ij}^2)$ (for $d_{ij} < 2R$). This has a number of shortcomings. The Hertz model is expected to be adequate only at sufficiently low cell deformations, which is not the case in cell configurations at high density. The latter cannot be avoided if no additional forces are added that take into account volume constraints. Below a certain distance, the geometrical contact area is more adequate but inconsistent with the Hertz force, as Hertz force and Hertz contact area are inherently linked. Moreover, the transition between the two areas is only smooth at $d_{ij}=R_i+R_j$, so at point contact. Therefore, the description of cell–cell interaction forces by Hertz (or JKR) force at large cell densities is only a crude approximation.

Finally, the Hertz and JKR-force emerge from linear elasticity assuming small deformations which in dense packings of cells can easily be violated, so that these force models can only be regarded as crude approximations in dense cell aggregates also with regard to this aspect.

The Voronoi tessellation approach originally proposed by Honda et al. [80, 175] have been considered also by other authors in center-based type modeling approaches, as it allows well approximating the surface and volume of a cell in dense aggregates, such as epithelial layers of intestinal crypts [34] (Fig. 18).

However, in some multi-cellular configurations such as multicellular spheroids or monolayers cells may detach from the main tumor body and enter the tumors’ environment [65, 178]. In a Voronoi tessellation the cell shape is constructed uniquely from the position of its neighbor cells hence the concept of an isolated cell within a Voronoi tessellation does not exist. In reference [83, 134] a possible solution to this problem using a modified Voronoi tessellation was sketched in which the biological cell shape was constructed as follows. A circle of radius R was drawn around each cell construction point (for a center-based cell, this would be the cell center). Then the Voronoi tessellation was constructed. If for a cell the Voronoi border was closer to the construction point than the circle, the biological cell border was associated with the Voronoi cell border, otherwise the biological cell border was associated with the circle (Fig. 18). Meyer-Hermann and co-workers [143, 179] explored this approach in more depth, and in three dimension for growing multi-cellular spheroids. As long as deformations of the cells remain moderate, the difference between the Voronoi approach and the presentation of the cells in isolation by spheres that deform due to Hertz or JKR force upon compression are moderate, as supported by very similar simulation results for EMT6/Ro multicellular spheroids with pairwise JKR-force [83] and Voronoi tessellation with a modified Hertz force [143].

However, volume estimates permit to consider explicitly compression terms. Using $K_\mathrm{V}=-V(\partial p/\partial V)$, a volume pairwise force term may be constructed by [180]

$$\begin{aligned} \underline{F}_{ij}^\mathrm{vol}\!= & {} \!\left( \frac{E_i}{3(1\!-\!2\nu _i)}log \frac{V_i}{V^0_i} \right. \nonumber \\&\left. \!+\, \frac{E_j}{3(1\!-\!2\nu _j)}log \frac{V_j}{V^0_j}\right) A_{ij}\frac{\underline{r}_i \!-\!\underline{r}_j}{|\underline{r}_i-\underline{r}_j|}. \end{aligned}$$

(42)

For $V_i\approx V^0_i$ and $V_j \approx V^0_j$, $log(V_i/V_i^0)$ can again be replaced by $(V_i-V_i^0)/V_i$. However, the approach has a shortcoming. The pairwise forces (Hertz and JKR) do already account for both, compression and deformation. The upper equation would represent compression hence the compression part would have to be eliminated from the Hertz (JRK) model which is not straightforward.

Pathmanaghan et al. [139] also used a Voronoi tessellation approach for a computational study of tissue mechanics. Despite the general higher realism for mechanics, these algorithms can be computationally expensive and also show some flaws for tissues under high mechanical stress [139].

In CBM a large part of the computational effort will be spend in the contact detection and the matrix solving of Eq. 27. Contact resolution is easily to implement while contact detection (i.e., finding list of neighbors) is relatively cheap $({\sim }NlogN)$ using grid-based methods or tree search algorithms. Once the forces are computed in Eq. 27 and convergence is found in the Conjugate Gradient algorithm, the velocities of the cells can be calculated to obtain the positions of the cells, by using e.g., a simple explicit Euler integration scheme. The stability of the simulation is then formally be determined by the timestep $\Delta t$.

$$\begin{aligned} \Delta t \le C\frac{\gamma }{k} \end{aligned}$$

(43)

where C is a prefactor, k is the stiffness of the system (e.g., Young modulus of the cells), and $\gamma $ is the friction (e.g., cell–cell friction). Thus, the higher the friction, the lower one can choose the timestep. In practice, one chooses a fixed timestep that results in a stable simulation. Alternatively, one may also opt for an adaptive time stepping scheme [143].

In some cases explicit matrix solving can be avoided if one considers only friction with a medium [143, 181], which speeds up the computation significantly. The maximal timestep is limited by the minimal ratio of the damping parameters to the contact stiffnesses. Time stepping may put significant restrictions to the solvability in a mechanically stiff system when the time scale of interest is long (e.g., several days or weeks in tumor simulations), because the number of time steps may become prohibitively high. To alleviate this, one can work with rescaled parameters such that the simulation time is shortened but the dynamics remains unaffected. This can be achieved, for example, if there are two timescales $\tau _1$ and $\tau _2$, with the first representing the characteristic mechanical relaxation of the cellular system and $\tau _2$ the division time of the cells. In case it holds that $\tau _2 \gg \tau _1$, then one can shorten $\tau _2$, such that this separation in timescales remains large yet the simulation time is decreased significantly.

Below we briefly summarize advantages and disadvantages of center-based models (CBMs) from our perspective. Evaluations on the computation times are done with regard to a standard PC.

3.2 Deformable cell models and vertex models

3.2.1 Basic concepts

CBM show several disadvantages. Basically their constitution does not allow to represent arbitrary cell shapes. Also they cannot give detailed information of the mechanical signals (tensile, compressive forces, shear forces) that are transmitted by the ECM via integrin receptors linked with the cytoplasm and CSK. On the other hand, detailed models that focus on the mechanics of the CSK as such (e.g., [169]) may be unable to capture the behavior of the whole cell. An intermediate scale model is thus desirable to bridge this gap. A number of authors introduced a series of models which can be categorized as Deformable Cell Models (DCM). Rejniak [182] proposed a 2D modeling framework based on deformable cells represented by connected line segments as units. The cells were capable of growing, dividing, and adhering to other cells, allowing for the formation of clusters of cells. The first 3D models did not address such complexity (only the cell surface is discretized) and were mainly used to simulate the large deformation mechanics in erythrocytes [183].

In DCM the cell body is discretized by a number of nodes which are connected by viscoelastic elements interacting via pairwise potential functions, creating a flexible scaffolding structure with arbitrary degrees of freedom per cell. The nodes at the boundary can be used to triangulate a cell surface (line segments in 2D), accounting for the mechanical response of the membrane and cortical CSK (see Fig. 19a). The forces in DCM originate from both cell–cell interactions and intracellular interactions due to the viscoelastic elements, bending resistance in the membrane, and global cell properties such as cell volume and surface area conservation.

A special class called vertex models (VM) are similar to DCM but the vertices form a polygonal tessellation (usually Voronoi) for the cells. These models are therefore rather suitable for tightly packed cell ensembles with negligible intercellular space. The terms that contribute to the mechanics in each cell are the line tension and adhesion along its common edges with other cells, the contractility of the cortical ring along the cell perimeter, and hydrostatic pressure related to cell area/volume. Also here one can define rules that govern cell proliferation, migration and apoptosis. Analyzing static equilibrium configurations such as the optimal packing in epithelial cells can be done using energy-based methods while for dynamical systems one can opt for equations of motion, similar as in CBM.

Often, the individual elements are represented by classical linear spring-damper systems. When the elastic and dissipative components are positioned in parallel, one has the well-known Kelvin–Voight element representing a solid like behavior. The vector force between two nodes i and j then reads (see Fig. 19b):

$$\begin{aligned} \underline{F}_{ij} = \underline{F^\mathrm{e}}_{ij} + \underline{F^\mathrm{v}}_{ij} = -k_\mathrm{s}(\underline{d}_{ij}-\underline{d}^0_{ij}) -\gamma \underline{v}_{ij} \end{aligned}$$

(44)

where $\underline{F^\mathrm{e}}_{ij}$ and $\underline{F^\mathrm{v}}_{ij}$ are the elastic and dissipative forces, $k_\mathrm{s}$ is the spring stiffness, $\gamma $ represents the level of dissipation, $\underline{d}^0_{ij}, \underline{d}_{ij}$ are the initial and actual distance vectors, and $\underline{v}_{ij}$ is the relative velocity between the nodes respectively. In some cases however, it is necessary to incorporate nonlinear elastic behavior which better describe the polymeric molecular structures [140, 184]. A well-known model for this is the Finitely Extensible Nonlinear Elastic (FENE) force which behaves soft at low stretch but become very stiff at high stretches (similar to uncoiling a rope):

$$\begin{aligned} F^\mathrm{e}_{ij}= -{k}_\mathrm{FENE} d_{ij} \left[ 1 - \left( \frac{d_{ij}}{d_{\max }}\right) ^{-2} \right] ^{-1}, \end{aligned}$$

(45)

where ${k}_\mathrm{FENE}$ is the stiffness, and $d^{\max }$ is the maximum stretch of the spring, at which the force will diverge. Van Liedekerke et al.[185] implemented a Neo-Hookean incompressible polymer model for a cell surface:

$$\begin{aligned} F^\mathrm{e}_{ij}={k}_\mathrm{poly}\left( \lambda - \lambda ^{-5} \right) , \end{aligned}$$

(46)

where $\lambda ={d}_{ij}/{d}^0_{ij}$ is the stretch ratio. This model takes the variation of the thickness of the material during stretching into account.

For a liquid like behavior the basic representation is the Maxwell model where a spring and damper are in series. The force equation then reads

$$\begin{aligned} \frac{\dot{F_{ij}}}{k} + \frac{F_{ij}}{\gamma } =v_{ij} \end{aligned}$$

(47)

This equation contains a derivative of the force, which can be approximated by $\dot{F} = (F(t)+F(t-\Delta t)) /\Delta t $ and involves the storage the force at the previous time step. Note that the spring-damper elements can be combined or rearranged [186] to obtain a more complex dynamics with multiple relaxation times. The surface bending resistance can be incorporated by computing the angle between two adjacent triangles $\alpha $ and $\beta $ (or line segments in 2D). This defines an out of plane energy $E_\mathrm{b}$:

$$\begin{aligned} E_\mathrm{b}= k_\mathrm{b}(1 - \cos (\mathbf n _{\alpha },\mathbf n _{\beta })) \end{aligned}$$

(48)

where $k_\mathrm{b}$ is the bending constant and $\mathbf n $ is the normal to each triangle. Using $\underline{F}_{i}=-\partial {E_\mathrm{b}}/\partial {\underline{r}_{i}}$, this can be transformed to an equivalent quadruplet force system acting on the nodes of the triangles $\{ \underline{F}_{1},\underline{F}_{12},\underline{F}_{21},\underline{F}_{2} \}$ (see Fig. 19b), for which momentum conservation demands that $\underline{F}_{1}+\underline{F}_{12}+\underline{F}_{21}+\underline{F}_{2} = 0$ [140].

The area and volume constraint forces arise from the limited compressibility of the membrane/cortex and cytoplasm of the cell. If we would consider a vesicle as a model (i.e., a spherical lipid bilayer encompassing a fluid), then additional forces $\underline{F}^\mathrm{area}$ will arise from the fact that the change of surface area of the membrane is limited. In addition, one should also try to keep the changes of the individual area of the triangles as low as possible, because similar to FE simulations, too strong deformed triangles can lead to mesh artefacts, instabilities and inaccurate results [184]. The magnitude of these forces can be tuned introducing an area conservation constant $K_\mathrm{A}$, with $K_\mathrm{A}=Eh$ where E and h are the apparent Young modulus and thickness of the membrane, and $\epsilon _\mathrm{A}=(A-A_{0})/A_{0}$ a strain measure (see Fig. 19). The forces $F^\mathrm{a}$ are then applied in the plane of each triangle in the direction from the barycenter of the triangle [140]. Note that the area conservation forces contribute to the membrane/cortex stiffness in a similar way as the spring forces in Eqs. 44 and 45. However, the latter ones are used to describe the elasticity cortical CSK, whereas the former one reflects the elasticity of the lipid bilayer of the cell. The linear spring constant for a sixfold symmetric triangulated lattice can be related approximatively to the cortex Young modulus $E^\mathrm{cor}$ with thickness h by [185, 187]

$$\begin{aligned} k_\mathrm{s}\approx \frac{2}{\sqrt{3}}E^\mathrm{cor}h^{} \end{aligned}$$

(49)

The bending stiffness of a (thin) cortex is in the same way approximated by

$$\begin{aligned} k_\mathrm{b}\approx \frac{E^\mathrm{cor}h^3}{12(1-\nu ^2)} \end{aligned}$$

(50)

where $\nu $ is the Poisson ratio (=1/3 for an equilateral 2D network of linear springs). For a more a more in depth analysis relating the spring force parameters to macroscopic constants, we refer to literature [140, 184, 185, 188].

The volume penalty forces are computed from the internal pressure using the volume change and compressibility $K_\mathrm{V}$. The internal pressure in the cell is given by

$$\begin{aligned} p = K_\mathrm{V}\epsilon _\mathrm{V} \end{aligned}$$

(51)

whereby $\epsilon _\mathrm{V}$ is the volume change strain (for small volume changes, we have $\epsilon _\mathrm{V}=(V-V_{0})/V_{0}$). The volume V of the cell can be computed summing up the volumes of the individual tetrahedra that build up the cell. The forces $F^\mathrm{vol}$ on the nodes are obtained by multiplying the pressure with the area of the triangle.

Finally, we arrive at the equation of motion for a node i of a cell in free suspension in its most simple form, taking all the former terms in account. This equation is formally the same as Eq. 27:

$$\begin{aligned} \sum _{j}\varGamma ^\mathrm{nn}_{ij} (\underline{v}_i-\underline{v}_j) +\varGamma ^\mathrm{ns}_{i} \underline{v}_i= & {} \sum _{j}\underline{F}^\mathrm{e}_{ij} + \sum _{\mathrm{quadruplets} \beta }\underline{F}^\mathrm{b}_{i, \beta } \nonumber \\&+ \sum _{\mathrm{triangles} \alpha }\underline{F}^\mathrm{a}_{i,\alpha } + \underline{F}^\mathrm{vol}_{i}.\nonumber \\ \end{aligned}$$

(52)

with $\varGamma ^\mathrm{nn}$ and $\varGamma ^\mathrm{ns}$ the matrices representing node-node friction and node substrate friction. The latter one can be obtained by considering the CBM substrate friction coefficient divided by the number of nodes at the surface of the cell.

The solving of this system is formally the same as Eq. 27 where the positions of the cells are replaced by the nodes. A Conjugate Gradient method will usually solve this system efficiently. The restrictions to the timestep are similar as in Eq. 43.

Having the forces at sub cellular scale at hand, one can now compute triangle-averaged stresses using the concept of viral stress. For example, the membrane/cortex tension in a node $\sigma _{i}$ can be calculated by [184]

$$\begin{aligned} \sigma _{i}\approx \frac{\sqrt{3}}{N_{j}}\sum _j{\frac{F^\mathrm{e}}{d_{ij}}} + K_\mathrm{A}\epsilon _\mathrm{A} \end{aligned}$$

(53)

where $F^\mathrm{e}$ is the elastic component of the force between two nodes, $N_{j}$ is the number of connections of node i (usually 6 -fold connectivity), and the second term is the contribution due to area conservation. Of course, these formulas need to be adapted if more detail such as an internal structure is added. In this regard it is also useful to mention that in some simulations of cell models without an internal structure (i.e., only fluid inside) unrealistic buckling instabilities can appear [189]. Adding an internal structure (CSK) or choosing a higher bending energy will likely solve this issue.

It is clear that the number of parameters and physical complexity in these model types is larger than in CBM, yet as the forces can be computed up to subcellular level, discerning a nucleus, organelles, membranes and internal structure, these models are candidates for investigation of mechanotransduction.

3.2.2 Achievements, limitations and practical use

The work of Rejniak and co-workers has treated multicellular spheroids, various cellular patterns in developing ductal carcinoma in situ, invasive tumors as well as normal development of epithelial ductal monolayers and their various mutants [190–193]. Deformable models have also been extensively used for erythrocytes to increase our insight in blood flow and related diseases such as malaria, with some authors using discretization up to the level of the spectrin CSK [194], and others using coarse graining models with realistic mechanical and dynamical behavior, see e.g., [184, 195–198]. It was also shown that the model parameters of the viscoelastic elements can be calibrated using optical tweezer or micro-pipetting experiments.

Models of cells with an internal structure or CSK are still relatively scarce, yet they are likely to become essential for addressing the more fundamental questions in cell migration and mechanotransduction. A basic, coarse representation of the CSK was proposed in Ingberg [199] who regarded a cell as a tensegrity structure making the connection with macroscopic structure models where forces are transmitted by flexible cables and stiff struts. In the so called subcellular element method (SEM), Sanderius et al. [186, 200] approached the cell as a 3D network using specific pairwise potentials and springs, to simulate the rheology of a single-cell and small cell clumps. This model can be regarded the off-lattice equivalent of the CPM. In both abovementioned methods however, the cell membrane and cortex was not discerned. An extended 2D deformable model was proposed by Jamali et al. [201], where in each cell a certain subcellular detail was incorporated, including the elastic plasma membrane, elements that perform the function of the CSK, and elements representing the cell nucleus. This model allows interaction and communication with its environment, changing its morphology, and performing basic cellular events such as growth, division, death, and polarization.

To analyze stress distribution and predict damage in compressed cells and impacted parenchyma tissue, Van Liedekerke et al. [185, 202] introduced a full Lagrangian particle model using the continuum smoothed particle hydrodynamics (SPH) technique, coupled with viscoelastic network model representing a solid hyperelastic membrane (Fig. 20), see Sect. 3.2.

Deformable cell models can incorporate models for specific adhesion (e.g., electrostatic force, depletion forces, van der Waals attractions) and non-specific adhesion (involving proteins and adhesion complexes). Often pairwise potential functions (Van der Waals, Morse) [186, 193, 201, 203] between nodes are used to mimic the effect of short range repulsion and long ranged attraction force. Odenthal et al. [140] on the other hand derived an interaction method based on the Maugis–Dugdale adhesion theory, using limited number of scalable parameters. Their model could successfully reproduce the experimental data of the dynamics of vesicles spreading to a substrate, and predicts a zone of inward traction stress on the substrate which was later observed experimentally [204], see Fig. 21c.

For cell migration, the same principals and methods originating from colloidal physics applied in CBM could be applied to DCM (see Eq. 29). Only recently attempts to model cell migration in a mechanistic way have been introduced. Kim et al. [205] integrated focal adhesion dynamics, cell membrane and nuclear remodeling, actin motor activity, and lamellipodia protrusion reflecting the 3D spatiotemporal dynamics of both cell spreading and migration. A comprehensive and detailed 2D model for cell migration in ECM was developed by Tozluoglu et al. [206], who analyzed the cell migration speed governed by the interplay between ECM structure, adhesion strength and CSK contractility. The sprouting mechanism in angiogenesis was investigated by Bentley et al. [207] who developed an agent-based model addressing both the internal actin-based cortical tension and filopodia driven migration in single cells, but also the higher, cell–cell level interactions of multiple adhered cells.

Using VM, a measure of how cell boundaries can be approximated by VM was thoroughly investigated by Honda and coworkers [80, 175, 208, 209] who studied how cells undergo rearrangement. Farhadifar and co-workers [210] used the approach to study the role of cell mechanics and cell division in determining network packing geometry in epithelial cells. Hilgenfeldt et al. [211] studied cell geometric order in an epithelial tissue. Manning et al. [212] deciphered the contributions of adhesion strength and cortical contractility in cells to tissue surface tension. Finally, Rudge and Haseloff [213] applied the vertex dynamics approach to study spatial mechanical factors and signal transduction in morphogenenis of plant tissue. An overview for the VM techniques can be found in [214].

Until now, deformable cell models have especially been proven to give accurate insight in cell mechanics on short timescale (seconds to minutes). The eventual dynamic that emerges will be dependent on the rheological model that is incorporated. Lim et al. [215, 216] examined several mechanical models that have been developed to characterize mechanical responses of living cells when subjected to both transient and dynamic loads, including cortical shell–liquid core (or liquid drop) models, the solid models, and power-law structural damping models. The latter seem to offer a more realistic description of the CSK dynamics [217]. It is important to mention that active phenomena involving actin contractility and CSK remodeling may require other internal models and parameters compared to those governing in short times experiments [218]. The choice of experiments for model calibration is therefore critical and can determine the type of parameters and mechanical properties.

The accuracy of the mechanics in DCM largely depends on the number of nodes assigned to each cell [184]. Depending on the application, a too small number of nodes per cell can create an unrealistic dynamics [140]. Computationally these models are quite expensive and the applications will limit themselves to relatively small cell numbers. The typical aspects of computational complexity (contact detection algorithms) in center-based models apply to deformable models as well.

Concluding, deformable cell models have already been shown to explain a large variety of phenomena in cell mechanics, especially related to short timescales. The computational effort and increased complexity make these models less widespread but they have the potential and might even become necessary to explain and predict certain phenomena related to cell mechanics and mechanobiology. Despite these useful developments in various directions, we believe that in order to investigate mechanotransduction in a quantitative way, more studies are needed to develop models which can capture the mechanical interactions between the cell interior and ECM accurately.

Below we briefly summarize advantages and disadvantages of deformable cell models (DCMs) from our perspective. Evaluations on the computation times are done with regard to a standard PC.

4 Hybrid discrete-continuum models

In case one faces large multicellular systems, discrete agent-based models can become computationally infeasible. Continuum models do not have this problem as they are not concerned with individual cells and they have no intrinsic scale. Continuum models basically solve the governing mass and momentum balance PDEs for the tissue dynamics, including source terms to for cell growth, birth and death, cell migration, as well as reaction-diffusion PDEs to compute the concentrations of chemical substances such as oxygen, glucose, chemokines, growth factors or cytokines. Continuum Models for tumors are often multi-phase mixture descriptions encompassing a tumor cells phase, intercellular fluid phase, and ECM phase, using Darcy’s law to determine the dynamics of each phase. It is not our intention to review these models here, instead we refer the reader to some basic key literature in this field [71, 219, 220]. Models that pay special attention to mechanical interaction with ECM and surrounding host tissue are generally more recent developments (see e.g., [221–225]). Continuum models generally have less parameters and they offer some advantages with regard to the information that needs to be analyzed. Nevertheless, the constitutive equations that characterize e.g., migration dynamics, cell adhesion, and stress-strain in the tissue and ECM will determine the results that will be obtained. The determination of those equations is a non-trivial task and in general they should be determined from experiments or might be inferred from agent-based models using averaging techniques [49, 83, 149, 226–229]. This branch is usually refered to as multiscale models within the mathematical and engineering community, and has to be distinguished from multi-level model which integrate components, mechanisms and information from different scales in one model [230].

It is likely however, that if one wants to capture the whole complexity of tissue dynamics, hybrid models will come into play. The formulation of hybrid models heavily draws on the conservation laws (mass, momentum, energy) and the key to success of the coupling to discrete models is an adequate (possibly on-the-fly) mapping between the constitutive behavior to the averaged discrete system with respect to growth kinetics, rheology, etc. [147, 224, 231–234]. A particularly appealing approach to bridge the gaps across the scales is based on the density functional theory (DFT) [229], which would allow calibration and validation of the cell-scale parameters and equations to be obtained directly from individual microscale measurements, and formulating rigorous upscaling techniques of the continuum equations at the larger scales.

Combining continuum and discrete (agent-based) models can be advantageous, e.g., if large parts in the tissue require less detail and can be treated with continuum approaches, while in other parts agent-based models can provide more detailed information. A coupling scheme was explained in Kim et al. [235] who proposed a domain decomposition for a tumor region whereby the active regions (proliferating cells) are represented by agents while the passive tissues (ECM, necrotic zone) were modeled by the Finite Element Method discretization of a PDE governing the behavior of an elastic material. The transfer of forces exerted by the cells to the mesh at the discrete-continuum interface are resolved by interpolation functions.

Hybrid approaches have been formulated, for example, to investigate wound healing [236–238], interaction of cells (CBM) and the basal membrane and ECM [52], angiogenesis [233, 239], and multicellular spheroids with invasion [233, 235].

In the approach where one regards tissue as deformable cells immersed in a continuum field (ECM, fluid) [190, 192], the dynamics of the cytoplasm and the medium (interstitial fluid, ECM) in which the cells are embedded is treated as a continuum and usually computed on an Eulerian mesh (see Fig. 21a). Interactions of the cell boundary with the cytoplasm as well as external medium are computed with the Immersed Boundary Method, a classic scheme for studying fluid–structure interactions. The interaction calculation between the Lagrangian and Eulerian points relies on the smearing out of vector or scalars fields computed in one single point over the neighboring computational points lying in a finite domain $\varOmega $, by considering the convolution integral:

$$\begin{aligned} \mathcal {F}(\underline{x})= \int _{\varOmega }\mathcal {F}(\underline{x'})\mathcal {K}(\underline{x'} -\underline{x}) \text {d}\underline{x'} \end{aligned}$$

(54)

in where $\mathcal {K}(\underline{x'}-\underline{x})$ is a symmetric kernel function with the property $\int _{\varOmega }\mathcal {K}(\underline{x'}) \text {d}\underline{x'} = 1$. In practice of course this integral will be replaced by a summation over all neighboring particles j:

$$\begin{aligned} \mathcal {F}(\underline{x}_{k})=\sum _{j}\mathcal {F}(\underline{x}_{j})\mathcal {K}(\underline{x}_{k}-\underline{x}_{j}) \mathcal {V}_{j} \end{aligned}$$

(55)

where $\mathcal {V}_{j}$ is the volume associated with that particle. One can also define $\mathcal {F}=\mathcal {M} / \mathcal {V}$, such that, for example if $\mathcal {F}$ represents the mass density $\rho $, then we have

$$\begin{aligned} \rho (\underline{x}_{k})= \sum _{j}m(\underline{x}_{j})\mathcal {K}(\underline{x}_{k} -\underline{x}_{j}) \end{aligned}$$

(56)

which is an ensemble averaging to derive a continuum field from discrete variables.

It might be noteworthy that although mesh-based methods like finite elements (FE) or finite differences (FD) are well established to solve the PDEs in continuum models, we notice that in the last decades alternative mesh free methods are gaining increasing attention. These methods use particles as computational nodes without underlying mesh or fixed connectivity. The vector fields and their gradients associated with a particle are computed using radial kernel functions rather than shape functions, see Eq. 56. SPH [240] is a mesh free method originally developed in astronomy and later fluid problems, while more recently it has also been applied in microscopic problems [241]. The main advantages are that no remeshing is required and media with a specific difficulties such as discontinuities or complex architectures can be dealt with more naturally. In addition these methods naturally allow a coupling with discrete particle models such as CBM or DCM. Disadvantages of SPH are the requirement a relatively high number of particles to attain the same accuracy as grid-based methods, and the non-trivial boundary treatment. Based on SPH, Van Liedekerke et al. [185, 202] proposed model to simulate dynamics of impact an damage in tissue (Fig. 21a). Tanaka et al. [242] and later Hosseini et al. [196] coupled the SPH Navier–Stokes solver to a 2D deformable cell model to simulate fluid flow associated with red blood cells. Van Liedekerke et al. [141] further extended the SPH method to solve the overdamped (Stokes) equations for a fluid, thereby drastically reducing the computational time.

Finally, we would like to emphasize that hybrid models do not necessarily restrict to the coupling of discrete and continuum models. In Fig. 21b we present a working simulation example where a CBM and DCM are coupled using the interaction model proposed in [140]. Moreover, one could even enrich the classical hybrid CBM—continuum model scheme with a higher detailed model such as DCM, as is schematically depicted in Fig. 21c.

5 Available tools and software development for modeling

While a lot of researchers develop and maintain their own code for agent-based simulations, several research groups have released tools for agent-based modeling in biology. These tools implement the CPM: CompuCell3D [131] Morpheus [130], Simmune [243] and Chaste [132, 244]; center-based models: CellSys [150], Chaste [132, 244], EPISIM [245], Timothy [246, 247] and Biocellion [248], or a deformable cell model: VirtualLeaf [249], LBIBCell [250] and SEM$++$ [181].

Most of these tools have one or more specific features, such as multi-scale modeling, usability, or performance. However, a few of them just provide the first, freely available, implementation of a specific model. The first of these is LBIBCell, which provides a framework for modeling cells as a viscous, Newtonian fluid encapsulated by an elastic polygon. The interactions between these two components are resolved with the immersed boundary method [190]. Furthermore, intracellular networks and extracellular concentrations can be modeled using reaction diffusion models. Another framework that provides a unique model implementation is VirtualLeaf, which implements a vertex-based model focused at plant development. Intracellular processes may be modeled using ODEs as well as transport between cells or between cells and their environment. SEM$++$ provides a framework for the subcellular element model [186] with several extensions. These extension include actin polymerization-based migration, cell–cell mediated phenotype determination and subcellular organelles. CellSys and Chaste are versatile tools for running center-based models. For CA models types (A, B, D) public codes seems to be less available for cell-based modeling although they exist for general applications [101].

All of the tools listed above are multiscale in a sense that they enable intracellular and/or extracellular signaling. Most commonly, ODEs are used to model intracellular signaling and reaction-diffusion equations are used to model extracellular concentrations. Simmune, a CPM framework, strongly focuses on linking spatial resolved signaling networks to cell-based models. For this, a network is defined for each lattice site and the networks are connected via their boundary conditions. In this manner, Simmune provides an efficient way of combining cell-based models with high resolution spatial signaling. In both CompuCell3D and EPISIM intracellular models can be added using SBML [251]. In CompuCell3D this is achieved via the BionetSolver plugin [252] which can handle any number of SBML models per cells. In EPISIM the SBML models are evaluated with Copasi [253]. Morpheus, one of the CPM frameworks, provides a variety of modeling formalisms for intracellular modeling, such as ODEs, CAs and fine state machines.

Several tools provide a basic graphical user interface (GUI) for parameter modification and/or simulation visualization. A few developers took a step further and developed much more elaborate GUIs that enable users to create and extend models with no or minimal programming. In the CPM framework CompuCell3D simulations are specified in XML and/or Python scripts. To ease setting up these scripts, CompuCell3D comes with a specialized editor: Twedit++, in which users can generate simulation scripts using a wizard and extend these scripts by adding existing code-snippets. Morpheus goes even a step further than this. Besides setting up, running, and visualizing simulations, Morpheus’ GUI also facilitates batch-processing for parameter sweeps and sensitivity analysis, and model analysis at run time. A tool for center-based models that focuses on usability is EPISIM. This tool employs a graphical modeling system (GMS) to set up the model and automatically converts that model to efficient code. To set up such a model EPISIM provides the Graphical Model Editor to specify an agent-based model using predefined components and functions from the Function Library. The model parameters and cell variables can edited with the Variable-sheet Editor.

Chaste, a tool that implements both the CPM and center-based models, aims to provide an extensive application program interface (API) for agent-based modeling. In contrast to most of the other tools discussed here, Chaste is a library that modelers can use to develop their own models. Combining either the CPM or center-based models with ODEs and/or PDEs. By utilizing automated testing, the developers ensure the correctness of the code while it is being developed. In this manner, the developers provide a tool for rapidly implementing models in a way that is reliable and reproducible.

Because biological systems can contain large numbers of cells, performance is an import issue with agent-based modeling tools. To improve performance, several tools are designed to multiple core: CompuCell3D, LBIBCell, SEM$++$, CellSys, Biocellion, and Timothy. Recently, two tools have been developed with a strong focus on performance: Timothy and Biocellion. Timothy is a tool for center-based modeling, which includes ODE and PDE models for intracellular and extracellular modeling. To be able to simulate cell populations of $10^9$ cells, Timothy was designed to take advantage of a specific supercomputer. Biocellion is a generic framework for center-based models. Users must supply the exact code body that governs the agent-based model as well as models for intracellular and extracellular signaling. Biocellion then distributes this code over the resources available on the specific platform. As result, Biocellion can run on a desktop PC with a few cores, or on a supercomputers with dozens or hundreds of cores.

6 Conclusion

Agent-based models of multicellular systems are becoming increasingly popular. An important reason for this is that they provide a natural direct description at the interface between the whole body of intracellular studies spanning genomics, transcriptomics, proteomics, metabolomics, lipidomics, and the study of intracellular signal transduction pathways on one hand, and organ and body physiology on the other hand. Intracellular mechanisms can readily been integrated into each individual cell agent, and the collective impact of cell responses on organ and body physiology be studied. In particular, computer simulations with agent-based models have been shown to be useful in analyzing development stages in cancer or transitions between them, cell-to-cell variability during cancer development, selection processes operating on differences between cell states, as well as transitions from aggregated population stages (as long as the tumor forms a single mass) to infiltrative invasion and growth, for example the epithelial to mesenchymal transition.

Recently, the role of mechanical stimuli in cell responses have increasingly moved into research focus by physics, complementing the more molecular mechanistic view of biologists. This has also raised the question of the pro’s and contra’s of the current agent-based approaches to biological cells. So far still insufficient quantitative comparisons between the model types and experiments exist to address this question. No optimal model exists unless one would classify as optimal a 1:1 in silico copy of the cell, which would be very computationally intense, and disagree to the often cited quote by Einstein, adapted: “Everything (a model) should be as simple as possible, but not simpler”. So the question is rather how much simplicity is adequate in a certain application. We present a simple comparative table overview of the models and their capabilities with regard to the feature they can address in the cellular system in Fig. 22.

Fundamentally, two classes can roughly be distinguished: lattice-based and lattice-free approaches. As reference throughout the paper, we have considered growing monolayers or growing multicellular spheroids as they have been modeled by almost each of these classes. With respect to this, all these models predict at least qualitatively similar dynamics.

We distinguished four lattice-based models refered to, in a wider sense, as cellular automaton models from which we enumerated types A, B, C according to their spatial resolution. Type A permits a certain number of cells $N_{\max } > 1$ on a lattice site, type B one cell on one lattice site, while in type C, known as Cellular Potts model, a cell can occupy many lattice sites simultaneously. In addition we considered LGCA as fourth class (type D). Different from types A-C, LGCA introduce the cell velocity as model parameter.

On regular lattices, the position of cells in the lattice can be described by integers. Operations, searching, etc.. is efficient so that the computation time for CA models is usually much faster than for their lattice-free counterparts addressing a comparable spatial resolution. The price to be paid for the high speed seems to be the possible artifacts in the emerging multicellular figures. Whether they occur or not depends on the parameter setting, but in some situations they might not be readily visible even though they could distort the value of observables.

Using an unstructured lattices (in our case a Voronoi diagram, representing the dual graph of a Delaunay triangulation) eliminates these artifacts, at the expense of higher computing time (Fig. 4).

A fundamental disadvantage of CA types A, B, D is that the representation of physical laws is rule-based, and that models on lattices have a lower length scale i.e., displacements smaller than the lattice spacings are not possible. This can lead to artifacts making the model being non-predictive. It is for example a huge challenge for the aforementioned models to properly simulate detachment of cells as a consequence of proliferation pressure in the interior of a monolayer (Fig. 13). Moreover, folding of tissue layers that are stabilized by polar adhesion forces as it is observed in intestinal crypts after irradiation (Fig. 18) or during the invagination of the blastula in early animal development [134] emerges naturally in a CBM (center-based model), which belongs to the class of lattice-free models, but it is very difficult to model properly with CA types A, B, D, and even for type C may represent a difficult task. The first reason is that CBMs represent a direct physical approach while in cellular automaton models physics has to be implemented in terms of rules such that on sufficiently large scales the correct mechanics is reproduced. The second reason is that CBMs, as they are off-lattice models, have now upper length scale and therefore permit gradual displacements.

LGCA (type D) for fluids, for example, can be shown to behave at large scales as the Navier–Stokes equation. Recent work by Deutsch and co-workers attempt to infer the large-scale behavior of LGCA for multicellular organization to better understand the large scale effect of the microscopic (cell scale) LGCA rules [99, 100, 102].

Gradual displacements cannot be modeled by type A, B, D cellular automaton models which makes it difficult to represent properly the cell mechanics, at least on spatial scales comparable to the lattice spacing. In model types A and B, a shifting length needs to be introduced to mimic pushing of cells that may occur if a growing and dividing cell exerts forces on its neighbor cells to generate the free space necessary for its division. By formal comparison with the expansion speed of monolayers, in type D (LGCA) models the velocity channel may take the function of the shifting length in model types A and B. Finally, the dynamics in type A, B, D models is stochastic and at least partially of the hopping type. For types A and B all neighbor configurations of a given configuration can be calculated and is usually not too large so that the master equation, that describes the underlying dynamics of the multicellular configuration in the model types A and B, can be numerically solved. The solution represents the dynamics of the probability distribution function for a certain multicellular configuration. For type C (CPM) and D (LGCA) CA models usually (Monte-Carlo) sampling methods are used. These have the disadvantage, that they can distort the time scale, making it difficult to associate an absolute, precise time scale to the number of Monte-Carlo steps. This represents one of the limitations of the CPM. Another limitation is that its physical parameters can couple in an unnatural way to cell migration. An incompressible elastic cell for example would not be able to move as moves are inherently linked to transient cell volume changes by flipping of lattice sites. This limitation might be removed either by separation between cell body and filopodia, or by flipping lattice sites only pairwise such that the volume of each cell remains constant. On the other hand, the modeling framework is very flexible, and complex cell shape changes can be captured. The latter has been proven successful in simulations of cell sorting in case of differential cell–cell adhesion [103] which turned out to provide an obstacle for those models that do not allow sufficiently large cell deformations. However, as all lattice-based models intrinsically base upon a stochastic dynamics, they are likely to fail if deterministic effects largely outrange stochastic contributions.

Off-lattice models have been considered based both on a Monte-Carlo dynamics and on equations of motion. For growing monolayers within the CBM approach both methods have been compared and no significant deviations in the simulation results have been found [62]. Definition of an absolute time scale is straightforward when using equations of motion, while they maintain their realism in absence of stochastic effects. They can also be easily extended with complex forces which is usually more intuitive than writing down the equation for the corresponding energy functions, although the latter may be a matter of the scientific background and community.

Within the original concept, in center-based models (CBMs) the interaction force of a cell with its neighbors is calculated based upon the pairwise force between the cell and each of its individual neighbors. CBMs can be completely parameterized by biophysical and cell-kinetic parameters, that in principle can all be measured. This is an important strength as it permits identification of a realistic range for each model parameter, thereby facilitating simulated sensitivity analyses. Moreover, this makes simulated sensitivity analyses feasible even in complex tissue organization models despite the otherwise longer simulation times compared to cellular automaton models [39]. CBM with a fixed intrinsic shape and interacting with pairwise forces only capture mechanical effects well if the cell shapes do not largely deviate from their intrinsic shape (usually spheres) and if the local cell density is not too high. The reason for the latter is that the precise cell shape in CBMs is not known so forces depending on shape or volume need approximations and force terms in addition to the pairwise force between the cell centers. This is prone to inconsistencies, for example, if on one hand pairwise Hertz forces are used to mimic the repulsive force upon compression and deformation, and an additional volume—force term is added at high cell densities, whereby the volume is approximated from Voronoi tessellation with a cut-off radius (Fig. 18 [143, 176]). The cut-off radius permits definition of cell shape and volume for the cell in isolation. Omitting the cut-off as in a pure Voronoi tessellation model [34], the shape and volume of a cell cannot be defined anymore as the shape of a cell within a Voronoi tessellation is solely determined by the position of the cells’ neighbor cells. The combined approach permits a reasonable estimate of the cell volume at any cell density which makes estimating of compression forces feasible at any density but generates an inconsistency as the Hertz contact force is inherently linked to a specific area (the Hertz contact area), that is different from the contact area between the same cells within a Voronoi tessellation. These differences can lead to significant problems in coherently defining forces and mechanical stresses which may partially be solved by sophisticated heuristic corrections but these render the models complicated and more computing time-intensive [180]. An additional drawback is that the Voronoi polygonal shape has to be updated at any moment in the simulation coming at the expense of longer simulation times.

Some authors assume the geometrical shape as approximation of the entire cell shape and control the cell volume by an additional dynamical variable [63], but the approach suffers from the same type of inconsistency. Hence these type of models are semi-quantitative in the sense that they do not permit precise force calculations simultaneously at small and large cell densities. A possibility to solve this issue might be by calibrating the interaction forces of center-based models by higher detailed models (e.g., deformable cell models, DCMs) which represent cell shape more realistically and self-consistently.

Deformable cell models represent cell shape explicitly. This is a fundamental advantage as now cell volume and contact area can be calculated more accurately. However, at high resolution, the choice of model mechanisms and the calibration of parameters is very difficult. The choice of models and parameters is usually determined by experiments probing cell mechanical properties at low scale, such as atomic force microscopy (AFM), pipette aspiration experiments, and optical stretchers. However, given the cells’ complexity and the difficulty to infer the individual stress contributions from membrane, cellular cortex, CSK, cytoplasm, nucleus and other cell organelles, the degrees of freedom within the highly parameterized DCMs can usually not be uniquely removed i.e., several models and parameter combinations may fit the same experiments. Due to the complexity of the model, sensitivity analyses are very time consuming. Moreover, in principle every cell type would have to be studied separately experimentally and then the experiments be used to calibrate a DCM for that cell type. Another issue is the computational effort. Because the mechanical timescales of subcellular properties are much smaller than those associated with observables as migration, growth and division, one can simulate only a small number of cells at each CPU. It seems that distributed computing can bring a solution here. Another solution may come for hybrid model approaches, so that in certain tissue parts, CBM may be replaced by a calibrated CBM. Nevertheless, given the development of experimental methods to collect information on cells individually and organized in tissues, DCMs are likely to play an increasingly important role in modeling and understanding mechanotransduction pathways.

What will be the role of lattice (CA) models in the future? This is hard to say, but with increasing experimental information on the physical cell properties, which are hard to properly represent in CA models, the future success of CAs will likely depend on in how far it will be possible to match the microscopic rules and CA parameters with physical observables. In granular matter, cellular automaton models, originally very successful, are almost not used anymore.

Notes

This formula is not complete as the general formula of viral stress includes a kinetic term related to the temperature of the molecular particle system. For a system of cells however, the latter term can be neglected and the stress corresponds to the Cauchy stress of the microscopic system [149].

References

Iskratsch T, Wolfenson H, Sheetz MP (2014) Appreciating force and shape the rise of mechanotransduction in cell biology. Nat Rev Mol Cell Biol 15:825–833
Article Google Scholar
Fletcher DA, Mullins RD (2010) Cell mechanics and the cytoskeleton. Nature 463:485–492
Article Google Scholar
Schwarz US, Safran SA (2013) Physics of adherent cells. Rev Mod Phys 85:1327–1381
Article Google Scholar
Sinha B, Köster D, Ruez R, Gonnord P, Bastiani M, Abankwa D, Stan RV, Butler-Browne G, Vedie B, Johannes L, Morone N, Parton RG, Raposo G, Sens P, Lamaze C, Nassoy P (2011) Cells respond to mechanical stress by rapid disassembly of caveolae. Cell 144(3):402–413
Article Google Scholar
Goldmann WH (2012) Mechanotransduction and focal adhesions. Cell Biol Int 36(7):649–652
Article Google Scholar
Cai Y, Sheetz MP (2009) Force propagation across cells: mechanical coherence of dynamic cytoskeletons. Curr Opin Cell Biol 21(1):47–50
Article Google Scholar
Ivanovska I, Swift J, Harada T, Pajerowski JD, Discher DE (2010) Physical plasticity of the nucleus and its manipulation. Methods Cell Biol 98:207–220
Article Google Scholar
Chen CS (2008) Mechanotransduction–a field pulling together? J Cell Sci 121(Pt 20):3285–3292
Article Google Scholar
Engler AJ, Sen S, Sweeney HL, Discher DE (2006) Matrix elasticity directs stem cell lineage specification. Cell 126:677–689
Article Google Scholar
Guilak F, Cohen DM, Estes BT, Gimble JM, Liedtke W, Chen CS (2009) Control of stem cell fate by physical interactions with the extracellular matrix. Cell Stem Cell 5(1):17–26
Article Google Scholar
Helmlinger G, Netti PA, Lichtenbeld HC, Melder RJ, Jain RK (1997) Solid stress inhibits the growth of multicellular tumor spheroids. Nat Biotechnol 15:778–783
Article Google Scholar
Cheng G, Tse J, Jain RK, Munn LL (2009) Micro-environmental mechanical stress controls tumor spheroid size and morphology by suppressing proliferation and inducing apoptosis in cancer cells. PLoS ONE 4:e4632
Article Google Scholar
Basan M, Risler T, Joanny JF, Garau XS, Prost J (2009) Homeostatic competition drives tumor growth and metastasis nucleation. HFSP J 3(4):265–272
Article Google Scholar
Alessandri K, Sarangi BR, Gurchenkov VV, Sinha B, Kießling TR, Fetler L, Rico F, Scheuring S, Lamaze C, Simon A, Geraldo S, Vignjevic D, Doméjean H, Rolland L, Funfak A, Bibette J, Bremond N, Nassoy P (2013) Cellular capsules as a tool for multicellular spheroid production and for investigating the mechanics of tumor progression in vitro. Proc Natl Acad Sci USA 110:14843–14848
Article Google Scholar
Menon S, Beningo KA (2011) Cancer cell invasion is enhanced by applied mechanical stimulation. PLoS ONE 6:e17277
Article Google Scholar
Radisky DC, Nelson CM (2013) Regulation of mechanical stress by mammary epithelial tissue structure controls breast cancer cell invasion. Oncotarget 4:498–499
Article Google Scholar
Tse JM, Cheng G, Tyrrell JA, Wilcox-Adelman SA, Boucher Y, Jain RK, Munn LL (2012) Mechanical compression drives cancer cells toward invasive phenotype. Proc Natl Acad Sci USA 109:911–916
Article Google Scholar
Trepat X, Wasserman MR, Angelini TE, Millet E, Weitz DA, Butler JP, Fredberg JJ (2009) Physical forces during collective cell migration. Nat Phys 5:426–430
Article Google Scholar
Tang X, Bajaj P, Bashir R, Saif TA (2011) How far cardiac cells can see each other mechanically. Soft Matter 7:6151
Article Google Scholar
Vernon RB, Angello JC, Iruela-Arispe ML, Lane TF, Sage EH (1992) Reorganization of basement membrane matrices by cellular traction promotes the formation of cellular networks in vitro. Lab Investig 66(5):536
Google Scholar
Manoussaki D, Lubkin SR, Vemon RB, Murray JD (1996) A mechanical model for the formation of vascular networks in vitro. Acta Biotheor 44(3):271–282
Article Google Scholar
Manoussaki D (2003) A mechanochemical model of angiogenesis and vasculogenesis. Esaim Math Model Numer Anal 37(4):581–599
Article MathSciNet MATH Google Scholar
Murray JD (2003) On the mechanochemical theory of biological pattern formation with application to vasculogenesis. C R Biol 326:239–252
Article Google Scholar
Namy P, Ohayon J, Tracqui P (2004) Critical conditions for pattern formation and in vitro tubulogenesis driven by cellular traction fields. J Theor Biol 227(1):103–120
Article MathSciNet Google Scholar
van Oers RFM, Rens EG, LaValley DJ, Reinhart-King CA, Merks RMH (2014) Mechanical cell-matrix feedback explains pairwise and collective endothelial cell behavior in vitro. PLoS Comput Biol 10:e1003774
Article Google Scholar
Pelham RJ, Wang Y-L (1997) Cell locomotion and focal adhesions are regulated by substrate flexibility. Proc Natl Acad Sci USA 94:13661–13665
Article Google Scholar
Cukierman E, Pankov R, Stevens DR, Yamada KM (2001) Taking cell-matrix adhesions to the third dimension. Science (New York, NY) 294:1708–1712
Article Google Scholar
Wu P, Giri A, Sun SX, Wirtz D (2014) Three-dimensional cell migration does not follow a random walk. Proc Natl Acad Sci USA 26:3949–3954
Article Google Scholar
Zaman MH, Trapani LM, Sieminski AL, Siemeski A, Mackellar D, Gong H, Kamm RD, Wells A, Lauffenburger DA, Matsudaira P (2006) Migration of tumor cells in 3D matrices is governed by matrix stiffness along with cell-matrix adhesion and proteolysis. Proc Natl Acad Sci USA 103:10889–10894
Article Google Scholar
Spencer SL, Gaudet S, Albeck JG, Burke JM, Sorger PK (2009) Non-genetic origins of cell-to-cell variability in TRAIL-induced apoptosis. Nature 459(7245):428–432
Article Google Scholar
Bertaux F, Stoma S, Drasdo D, Batt G (2014) Modeling dynamics of cell-to-cell variability in TRAIL-induced apoptosis explains fractional killing and predicts reversible resistance. PLoS Comput Biol 10:e1003893
Article Google Scholar
Ramis-Conde I, Chaplain MAJ, Anderson ARA, Drasdo D (2009) Multi-scale modelling of cancer cell intravasation: the role of cadherins in metastasis. Phys Biol 6(1):16008
Article Google Scholar
Hoehme S, Brulport M, Bauer A, Bedawy E, Schormann W, Hermes M, Puppe V, Gebhardt R, Zellmer S, Schwarz M, Bockamp E, Timmel T, Hengstler JG, Drasdo D (2010) Prediction and validation of cell alignment along microvessels as order principle to restore tissue architecture in liver regeneration. Proc Natl Acad Sci USA 107(23):10371–10376
Article Google Scholar
Meineke FA, Potten CS, Loeffler M (2001) Cell migration and organization in the intestinal crypt using a lattice-free model. Cell Prolif 34:253–266
Article Google Scholar
Buske P, Przybilla J, Loeffler M, Sachs N, Sato T, Clevers H, Galle J (2012) On the biomechanics of stem cell niche formation in the gut modelling growing organoids. FEBS J 279(18):3475–3487
Article Google Scholar
van Leeuwen IMM, Mirams GR, Walter A, Fletcher A, Murray P, Osborne J, Varma S, Young SJ, Cooper J, Doyle B, Pitt-Francis J, Momtahan L, Pathmanathan P, Whiteley JP, Chapman SJ, Gavaghan DJ, Jensen OE, King JR, Maini PK, Waters SL, Byrne HM (2009) An integrative computational model for intestinal tissue renewal. Cell Prolif 42:617–636
Article Google Scholar
Pin C, Parker A, Gunning AP, Ohta Y, Johnson IT, Carding SR, Sato T (2015) An individual based computational model of intestinal crypt fission and its application to predicting unrestrictive growth of the intestinal epithelium. Integr Biol 7:213–228
Article Google Scholar
Hammad S, Hoehme S, Friebel A, von Recklinghausen I, Othman A, Begher-Tibbe B, Reif R, Godoy P, Johann T, Vartak A, Golka K, Bucur PO, Vibert E, Marchan R, Christ B, Dooley S, Meyer C, Ilkavets I, Dahmen U, Dirsch O, Böttger J, Gebhardt R, Drasdo D, Hengstler JG (2014) Protocols for staining of bile canalicular and sinusoidal networks of human, mouse and pig livers, three-dimensional reconstruction and quantification of tissue microarchitecture by image processing and analysis. Arch Toxicol 88:1161–1183
Article Google Scholar
Drasdo D, Hoehme S, Hengstler JG (2014) How predictive quantitative modelling of tissue organisation can inform liver disease pathogenesis. J Hepatol 61(4):951–956
Article Google Scholar
Roeder I, Loeffler M (2002) A novel dynamic model of hematopoietic stem cell organization based on the concept of within-tissue plasticity. Exp Hematol 30:853–861
Article Google Scholar
Eden M (1961) A two-dimensional growth process. In: Proceedings of the fourth Berkeley symposium on mathematical statistics and probability, vol 4: contributions to biology and problems of medicine. The Regents of the University of California
Batchelor M, Henry B (1991) Limits to Eden growth in two and three dimensions. Phys Lett A 157:229–236
Article Google Scholar
Kansal AR, Torquato S, Harsh GR IV, Chiocca EA, Deisboeck TS (2000) Simulated brain tumor growth dynamics using a three-dimensional cellular automaton. J Theor Biol 203(4):367–382
Mallet DG, De Pillis LG (2006) A cellular automata model of tumor-immune system interactions. J Theor Biol 239(3):334–350
Article MathSciNet Google Scholar
Deisboeck TS, Wang Z, Macklin P, Cristini V (2011) Multiscale cancer modeling. Annu Rev Biomed Eng 13:127–155
Article Google Scholar
Deutsch A, Dormann S (2004) Cellular automaton modeling of biological pattern formation: characterization, applications, and analysis. In: Modeling and simulation in science, engineering and technology. Birkhaeuser Verlag AG, Boston
Richardson D (2008) Random growth in a tessellation. Math Proc Camb Philos Soc 74:515
Article MathSciNet MATH Google Scholar
Alber MS, Kiskowski MA, Glazier JA, Jiang Y (2003) On cellular automaton approaches to modeling biological cells. Math Syst Theory Biol Commun Comput Finance 134:1–39
Article MathSciNet MATH Google Scholar
Alber M, Chen N, Lushnikov PM, Newman SA (2007) Continuous macroscopic limit of a discrete stochastic model for interaction of living cells. Phys Rev Lett 99(16):168102
Article Google Scholar
Lushnikov PM, Chen N, Alber M (2008) Macroscopic dynamics of biological cells interacting via chemotaxis and direct contact. Phys Rev E Stat Nonlinear Soft Matter Phys 78(6):061904
Article Google Scholar
Scianna M, Preziosi L (2012) A hybrid model describing different morphologies of tumor invasion fronts. Math Model Nat Phenom 7(1):78–104
Article MathSciNet MATH Google Scholar
D’Antonio G, Macklin P, Preziosi L (2013) An agent-based model for elasto-plastic mechanical interactions between cells, basement membrane and extracellular matrix. Math Biosci Eng 10:75–101
Graner F, Glazier JA (1992) Simulation of biological cell sorting using a two-dimensional extended Potts model. Phys Rev Lett 69:2013–2016
Article Google Scholar
Merks RM, Glazier JA (2005) A cell-centered approach to developmental biology. Physica A 352:113–130
Article Google Scholar
Savill NJ, Hogeweg P (1996) Modelling morphogenesis: from single cells to crawling slugs. J Theor Biol 184:229–235
Article Google Scholar
Marée AF, Hogeweg P (2001) How amoeboids self-organize into a fruiting body: multicellular coordination in Dictyostelium discoideum. Proc Natl Acad Sci USA 98(7):3879–3883
Article Google Scholar
Alarcón T, Byrne HM, Maini PK (2003) A cellular automaton model for tumour growth in inhomogeneous environment. J Theor Biol 225(2):257–274
Article MathSciNet Google Scholar
Alarcón T, Byrne HM, Maini PK (2004) A mathematical model of the effects of hypoxia on the cell-cycle of normal and cancer cells. J Theor Biol 229(3):395–411
Article MathSciNet Google Scholar
Ribba B, Alarcon T, Marron K, Maini PK, Agur Z (2004) The use of hybrid cellular automaton models for improving cancer therapy. In: ACRI. LNCS. Springer, Heidelberg
Alarcon T, Byrne HM, Maini PK (2010) A multiple scale model for tumor growth. Multiscale Model Simul 3(2):440
Article MathSciNet MATH Google Scholar
Drasdo D, Kree R, McCaskill JS (1995) Monte Carlo approach to tissue-cell populations. Phys Rev E 52:6635–6657
Article Google Scholar
Drasdo D, Hoehme S, Block M (2007) On the role of physics in the growth and pattern formation of multi-cellular systems: what can we learn from individual-cell based models? J Stat Phys 128(1–2):287–345
Article MathSciNet MATH Google Scholar
Galle J, Loeffler M, Drasdo D (2005) Modeling the effect of deregulated proliferation and apoptosis on the growth dynamics of epithelial cell populations in vitro. Biophys J 88(1):62–75
Article Google Scholar
Montel F, Delarue M, Elgeti J, Vignjevic D, Cappello G, Prost J (2012) Isotropic stress reduces cell proliferation in tumor spheroids. New J Phys 14:055008
Article Google Scholar
Ramis-Conde I, Drasdo D, Anderson ARA, Chaplain MAJ (2008) Modeling the influence of the E-cadherin-$\beta $-catenin pathway in cancer cell invasion: a multiscale approach. Biophys J 95(1):155–165
Article Google Scholar
Schluter DK, Ramis-Conde I, Chaplain MAJ (2012) Computational modeling of single-cell migration: the leading role of extracellular matrix fibers. Biophys J 103:1141–1151
Article Google Scholar
Radszuweit M, Block M, Hengstler J, Schöll E, Drasdo D (2009) Comparing the growth kinetics of cell populations in two and three dimensions. Phys Rev E 79:051907
Article MathSciNet Google Scholar
Dormann S, Deutsch A, Lawniczak AT (2001) Fourier analysis of Turing-like pattern formation in cellular automaton models. Future Gen Comput Syst 17:901–909
Article MATH Google Scholar
Dormann S, Deutsch A (2002) Modeling of self-organized avascular tumor growth with a hybrid cellular automaton. In Silico Biol 2:393–406
Google Scholar
Jones GW, Chapman SJ (2012) Modeling growth in biological materials. SIAM Rev 54:52–118
Article MathSciNet MATH Google Scholar
Lowengrub JS, Frieboes HB, Jin F, Chuang Y-L, Li X, Macklin P, Wise SM, Cristini V (2010) Nonlinear modelling of cancer: bridging the gap between cells and tumours. Nonlinearity 23:R1–R9
Article MathSciNet MATH Google Scholar
Block M, Schöll E, Drasdo D (2007) Classifying the expansion kinetics and critical surface dynamics of growing cell populations. Phys Rev Lett 99(24):248101
Article Google Scholar
Anderson ARA, Chaplain MAJ, Rejniak KA (eds) (2007) Single-Cell-Based Models in Biology and Medicine. Springer, Berlin
MATH Google Scholar
Anderson ARA, Weaver AM, Cummings PT, Quaranta V (2006) Tumor morphology and phenotypic evolution driven by selective pressure from the microenvironment. Cell 17(5):905–915
Article Google Scholar
Lee D-S, Rieger H, Bartha K (2006) Flow correlated percolation during vascular remodeling in growing tumors. Phys Rev Lett 96:058104
Article Google Scholar
Enderling H, Hahnfeldt P (2011) Cancer stem cells in solid tumors: Is ’evading apoptosis’ a hallmark of cancer? Prog Biophys Mol Biol 106(2):391–399
Article Google Scholar
Alfonso JCL, Jagiella N, Núñez L, Herrero MA, Drasdo D (2014) Estimating dose painting effects in radiotherapy: a mathematical model. PLoS One 9:e89380
Article Google Scholar
Kansal A, Torquato S, Harsh G IV, Chiocca E, Deisboeck T (2000) Cellular automaton of idealized brain tumor growth dynamics. Biosystems 55(1–3):119–127
Article Google Scholar
Honda H (1978) Description of cellular patterns by Dirichlet domains: the two-dimensional case. J Theor Biol 72(3):523–543
Article MathSciNet Google Scholar
Honda H (1983) Geometrical models for cells in tissues. Int Rev Cytol 81:191–248
Article Google Scholar
Jagiella N, Mueller B, Mueller M, Vignon-Clementel I, Drasdo D (2015) Inferring growth control mechanisms in growing multi-cellular spheroids of NSCLC cells from spatial-temporal image data. PLoS Comput Biol (accepted)
Drasdo D, Hoehme S (2005) A single-cell-based model of tumor growth in vitro: monolayers and spheroids. Phys Biol 2:133–147
Article Google Scholar
Drasdo D (2005) Coarse graining in simulated cell populations. Adv Complex Syst 8(2–3):319–363
Article MathSciNet MATH Google Scholar
Hoehme S, Drasdo D (2010) Biomechanical and nutrient controls in the growth of mammalian Cell populations. Math Popul Stud 17(166–187):37–41
MathSciNet MATH Google Scholar
Bortz A, Kalos M, Lebowitz J (1975) A new algorithm for Monte Carlo simulation of Ising spin systems. J Comput Phys 17:10–18
Gillespie DT (1977) Exact stochastic simulation of coupled chemical reactions. J Phys Chem 81:2340–2361
Article Google Scholar
Beysens DA, Forgacs G, Glazier JA (2000) Cell sorting is analogous to phase ordering in fluids. Proc Natl Acad Sci USA 97:9467–9471
Article Google Scholar
Kardar M, Parisi G, Zhang YC (1986) Dynamic scaling of growing interfaces. Phys Rev Lett 56(9):889–892
Article MATH Google Scholar
Halpin-Healy T, Zhang Y-C (1995) Kinetic roughening phenomena, stochastic growth, directed polymers and all that. Aspects of multidisciplinary statistical mechanics. Phys Rep 254(4–6):215–414
Article Google Scholar
Brú A, Albertos S, Luis Subiza J, García-Asenjo JL, Brú I (2003) The universal dynamics of tumor growth. Biophys J 85(5):2948–2961
Article Google Scholar
Huergo MAC, Pasquale MA, González PH, Bolzán AE, Arvia AJ (2011) Dynamics and morphology characteristics of cell colonies with radially spreading growth fronts. Phys Rev E Stat Nonlinear Soft Matter Phys 84:021917
Article Google Scholar
Huergo MAC, Pasquale MA, González PH, Bolzán AE, Arvia AJ (2012) Growth dynamics of cancer cell colonies and their comparison with noncancerous cells. Phys Rev E Stat Nonlinear Soft Matter Phys 85(1):011918
Article Google Scholar
Yates CA, Baker RE (2013) Isotropic model for cluster growth on a regular lattice. Phys Rev E Stat Nonlinear Soft Matter Phys 88(2):023304
Article Google Scholar
Jagiella N (2012) Parameterization of lattice-based tumor models from data. PhD Thesis, INRIA
Frisch U, Hasslacher B, Pomeau Y (1986) Lattice-gas automata for the Navier-Stokes equation. Phys Rev Lett 56:1505–1508
Article Google Scholar
Rothman DH, Zaleski S (2004) Lattice-gas cellular automata: simple models of complex hydrodynamics. Cambridge University Press, Cambridge
MATH Google Scholar
Rivet J-P, Boon JP (2005) Lattice gas hydrodynamics. Cambridge University Press, Cambridge
MATH Google Scholar
Gershenfeld NA (1999) The nature of mathematical modeling. Cambridge University Press, Cambridge
MATH Google Scholar
Böttger K, Hatzikirou H, Chauviere A, Deutsch A (2012) Investigation of the migration/proliferation dichotomy and its impact on avascular glioma invasion. Math Model Nat Phenom 7(1):105–135
Article MathSciNet MATH Google Scholar
Tektonidis M, Hatzikirou H, Chauvière A, Simon M, Schaller K, Deutsch A (2011) Identification of intrinsic in vitro cellular mechanisms for glioma invasion. J Theor Biol 287(1):131–147
Hoekstra AG, Kroc J, Sloot PM (2010) Simulating complex systems by cellular automata. Springer, Berlin
MATH Google Scholar
Hatzikirou H, Deutsch A (2010) Lattice-gas cellular automaton modeling of emergent behavior in interacting cell populations. Underst Complex Syst 2010:301–331
Article MATH Google Scholar
Glazier JA, Graner F (1993) Simulation of the differential adhesion driven rearrangement of biological cells. Phys Rev E 47(3):2128–2154
Article Google Scholar
Turner S, Sherratt JA (2002) Intercellular adhesion and cancer invasion: a discrete simulation using the extended Potts model. J Theor Biol 216:85–100
Rubenstein BM, Kaufman LJ (2008) The role of extracellular matrix in glioma invasion: a cellular Potts model approach. Biophys J 95:5661–5680
Article Google Scholar
Shirinifard A, Gens JS, Zaitlen BL, Popawski NJ, Swat M, Glazier JA (2009) 3D multi-cell simulation of tumor growth and angiogenesis. PLoS ONE 4(10):e7190
Article Google Scholar
Boghaert E, Radisky DC, Nelson CM (2014) Lattice-based model of ductal carcinoma in situ suggests rules for breast cancer progression to an invasive state. PLoS Comput Biol 10:e1003997
Article Google Scholar
Szabó A, Merks RMH (2013) Cellular Potts modeling of tumor growth, tumor invasion and tumor evolution. Front Oncol 3:87
Article Google Scholar
Li JF, Lowengrub J (2014) The effects of cell compressibility, motility and contact inhibition on the growth of tumor cell clusters using the Cellular Potts Model. J Theor Biol 343:79–91
Article Google Scholar
Gao X, McDonald JT, Hlatky L, Enderling H (2013) Acute and fractionated irradiation differentially modulate glioma stem cell division kinetics. Cancer Res 73(5):1481–1490
Article Google Scholar
Merks RMH, Glazier JA (2006) Dynamic mechanisms of blood vessel growth. Nonlinearity 19(1):C1–C10
Article MathSciNet MATH Google Scholar
Merks RMH, Perryn ED, Shirinifard A, Glazier JA (2008) Contact-inhibited chemotaxis in de novo and sprouting blood-vessel growth. PLoS Comput Biol 4(9):e1000163
Article MathSciNet Google Scholar
Palm MM, Merks RMH (2013) Vascular networks due to dynamically arrested crystalline ordering of elongated cells. Phys Rev E 87(1):12725
Article Google Scholar
Boas SEM, Merks RMH (2014) Synergy of cellcell repulsion and vacuolation in a computational model of lumen formation. J R Soc Interface 11(92):20131049
Article Google Scholar
Bauer AL, Jackson TL, Jiang Y (2009) Topography of extracellular matrix mediates vascular morphogenesis and migration speeds in angiogenesis. PLoS Comput Biol 5:e1000445
Article MathSciNet Google Scholar
Bauer AL, Jackson TL, Jiang Y (2007) A cell-based model exhibiting branching and anastomosis during tumor-induced angiogenesis. Biophys J 92:3105–3121
Article Google Scholar
Daub JT, Merks RMH (2013) A cell-based model of extracellular-matrix-guided endothelial cell migration during angiogenesis. Bull Math Biol 75(8):1377–1399
Leoncini E (2010) Applied mathematics to biology and medicine. PhD Thesis, INRIA
Ouchi N, Glazier JA, Rieu J, Upadhyaya A, Sawada Y (2003) Improving the realism of the cellular Potts model in simulations of biological cells. Physica A 329(3–4):451–458
Article MathSciNet MATH Google Scholar
Magno R, Grieneisen VA, Marée AF (2015) The biophysical nature of cells: potential cell behaviours revealed by analytical and computational studies of cell surface mechanics. BMC Biophys 8(8). doi:10.1186/s13628-015-0022-x
Albert PJ, Schwarz US (2014) Dynamics of cell shape and forces on micropatterned substrates predicted by a cellular Potts model. Biophys J 106(11):2340–2352
Article Google Scholar
Swat MH, Belmonte J, Heiland RW, Zaitlen BL, Glazier JA, Shirinifard A (2014) CompuCell 3D Reference Manual 3.7.3
Merks RMH, Brodsky SV, Goligorksy MS, Newman SA, Glazier JA (2006) Cell elongation is key to in silico replication of in vitro vasculogenesis and subsequent remodeling. Dev Biol 289(1):44–54
Article Google Scholar
Scianna M, Preziosi L (2013) A cellular Potts model for the MMP-dependent and -independent cancer cell migration in matrix microtracks of different dimensions. Comput Mech 53:485–497
Scianna M, Preziosi L, Wolf K (2013) A cellular Potts model simulating cell migration on and in matrix. Math Biosci Eng 10(1):235–261
Article MathSciNet MATH Google Scholar
Szabó A, Varga K, Garay T, Hegedus B, Czirók A (2012) Invasion from a cell aggregate–the roles of active cell motion and mechanical equilibrium. Phys Biol 9:16010
Article Google Scholar
Lemmon CA, Romer LH (2010) A predictive model of cell traction forces based on cell geometry. Biophys J 99:L78–L80
Article Google Scholar
Boas SEM, Jimenez MIN, Merks RMH, Blom JG (2015) A global sensitivity analysis approach for morphogenesis models. arXiv preprint, pp 1–29
Palm MM, Merks RMH (2015) Large-scale parameter studies of cell-based models of tissue morphogenesis using CompuCell 3D or VirtualLeaf. In: Nelson CM (ed) Tissue morphogenesis, vol 1189, Methods in molecular biology. Springer, New York, pp 301–322
Starruß J, de Back W, Brusch L, Deutsch A (2014) Morpheus: a user-friendly modeling environment for multiscale and multicellular systems biology. Bioinformatics 30:1331–1332
Article Google Scholar
Swat MH, Thomas GL, Belmonte JM, Shirinifard A, Hmeljak D, Glazier JA (2012) Multi-scale modeling of tissues using CompuCell 3D. In: Asthagiri AR, Arkin AP (eds) Computer methods in cell biology. Academic Press, London, pp 325–366
Chapter Google Scholar
Mirams GR, Arthurs CJ, Bernabeu MO, Bordas R, Cooper J, Corrias A, Davit Y, Dunn S-J, Fletcher AG, Harvey DG, Marsh ME, Osborne JM, Pathmanathan P, Pitt-Francis J, Southern J, Zemzemi N, Gavaghan DJ (2013) Chaste: an open source C++ library for computational physiology and biology. PLoS Comput Biol 9:e1002970
Article MathSciNet Google Scholar
Marée AFM, Grieneisen VA, Hogeweg P, Maree AFM (2007) The Cellular Potts Model and Biophysical Properties of Cells, Tissues and Morphogenesis. In: Anderson ARA, Chaplain MAJ, Rejniak KA (eds) Single-cell-based models in biology and medicine. Springer, Basel, pp 107–136
Chapter Google Scholar
Drasdo D, Forgacs G (2000) Modeling the interplay of generic and genetic mechanisms in cleavage, blastulation, and gastrulation. Dev Dyn 219(2):182–191
Article Google Scholar
Odell GM, Oster G, Alberch P, Burnside B (1981) The mechanical basis of morphogenesis. I. Epithelial folding and invagination. Dev Biol 85(2):446–462
Article Google Scholar
Sepúlveda N, Petitjean L, Cochet O, Grasland-Mongrain E, Silberzan P, Hakim V (2013) Collective cell motion in an epithelial sheet can be quantitatively described by a stochastic interacting particle model. PLoS Comput Biol 9(3):e1002944
Article MathSciNet Google Scholar
Drasdo D, Loeffler M (2001) Individual-based models to growth and folding in one-layered tissues: intestinal crypts and early development. Nonlinear Anal Theory Methods Appl 47:245–256
Article MathSciNet MATH Google Scholar
Basan M, Prost J, Joanny J-F, Elgeti J (2011) Dissipative particle dynamics simulations for biological tissues: rheology and competition. Phys Biol 8:026014
Article Google Scholar
Pathmanathan P, Cooper J, Fletcher A, Mirams G, Montahan L, Murray P, Osborne J, Pitt-Francis J, Walter A, Chapman SJ (2009) A computational study of discrete mechanical tissue models. Phys Biol 6(3):36001
Article Google Scholar
Odenthal T, Smeets B, Van Liedekerke P, Tijskens E, Van Oosterwyck H, Ramon H (2013) Analysis of initial cell spreading using mechanistic contact formulations for a deformable cell model. PLoS Comput Biol 9:e1003267
Article Google Scholar
Van Liedekerke P, Smeets B, Odenthal T, Tijskens E, Ramon H (2013) Solving microscopic flow problems using Stokes equations in SPH. Comput Phys Commun 184:1686–1696
Article Google Scholar
Turlier H, Audoly B, Prost J, Joanny J-F (2014) Furrow constriction in animal cell cytokinesis. Biophys J 106:114–123
Article Google Scholar
Schaller G, Meyer-Hermann M (2005) Multicellular tumor spheroid in an off-lattice Voronoi-Delaunay cell model. Phys Rev E 71:51910
Article MathSciNet Google Scholar
Drasdo D, Hoehme S (2012) Modeling the impact of granular embedding media, and pulling versus pushing cells on growing cell clones. New J Phys 14(5):55025
Article Google Scholar
Aragona M, Panciera T, Manfrin A, Giulitti S, Michielin F, Elvassore N, Dupont S, Piccolo S (2013) A mechanical checkpoint controls multicellular growth through YAP/TAZ regulation by actin-processing factors. Cell 154:1047–1059
Article Google Scholar
Low BC, Pan CQ, Shivashankar GV, Bershadsky A, Sudol M, Sheetz M (2014) YAP/TAZ as mechanosensors and mechanotransducers in regulating organ size and tumor growth. FEBS Lett 588:2663–2670
Article Google Scholar
Byrne H, Drasdo D (2009) Individual-based and continuum models of growing cell populations: a comparison. Math Biol 58:657–680
Article MathSciNet MATH Google Scholar
Irving JH, Kirkwood JG (1950) The statistical mechanical theory of transport processes. IV. The equations of hydrodynamics. J Chem Phys 18:817
Article MathSciNet Google Scholar
Ghysels P, Samaey G, Tijskens B, Liedekerke PV, Ramon H, Roose D (2009) Multi-scale simulation of plant tissue deformation using a model for individual cell mechanics. Phys Biol 6(1):16009
Article Google Scholar
Hoehme S, Drasdo D (2010) A cell-based simulation software for multicellular systems. Bioinformatics 26:2641–2642
Article MATH Google Scholar
Friebel A, Neitsch J, Johann T, Hammad S, Hengstler JG, Drasdo D, Hoehme S (2015) TiQuant: software for tissue analysis, quantification and surface reconstruction. Bioinformatics. doi:10.1093/bioinformatics/btv346
Godoy P, Hewitt N, Albrecht U, Andersen M, Ansari N, Bhattacharya S, Bode J, Bolleyn J, Borner C, Böttger J, Braeuning A, Budinsky R, Burkhardt B, Cameron N, Camussi G, Cho C-S, Choi Y-J, Craig Rowlands J, Dahmen U, Damm G, Dirsch O, Donato M, Dong J, Dooley S, Drasdo D, Eakins R, Ferreira K, Fonsato V, Fraczek J, Gebhardt R, Gibson A, Glanemann M, Goldring C, Gómez-Lechón M, Groothuis G, Gustavsson L, Guyot C, Hallifax D, Hammad S, Hayward A, Häussinger D, Hellerbrand C, Hewitt P, Hoehme S, Holzhütter H-G, Houston J, Hrach J, Ito K, Jaeschke H, Keitel V, Kelm J, Kevin Park B, Kordes C, Kullak-Ublick G, LeCluyse E, Lu P, Luebke-Wheeler J, Lutz A, Maltman D, Matz-Soja M, McMullen P, Merfort I, Messner S, Meyer C, Mwinyi J, Naisbitt D, Nussler A, Olinga P, Pampaloni F, Pi J, Pluta L, Przyborski S, Ramachandran A, Rogiers V, Rowe C, Schelcher C, Schmich K, Schwarz M, Singh B, Stelzer E, Stieger B, Stöber R, Sugiyama Y, Tetta C, Thasler W, Vanhaecke T, Vinken M, Weiss T, Widera A, Woods C, Xu J, Yarborough K, Hengstler J (2013) Recent advances in 2D and 3D in vitro systems using primary hepatocytes, alternative hepatocyte sources and non-parenchymal liver cells and their use in investigating mechanisms of hepatotoxicity, cell signaling and ADME. Arch Toxicol 87(8):1315–1530
Galle J, Aust G, Schaller G, Beyer T, Drasdo D (2006) Individual cell-based models of the spatial-temporal organization of multicellular systems-achievements and limitations. Cytometry 69:704–710
Article Google Scholar
Delarue M, Montel F, Vignjevic D, Prost J, Joanny J-F, Cappello G (2014) Compressive stress inhibits proliferation in tumor spheroids through a volume limitation. Biophys J 107:1821–1828
Article Google Scholar
Macklin P, Edgerton ME, Thompson AM, Cristini V (2012) Patient-calibrated agent-based modelling of ductal carcinoma in situ (DCIS): from microscopic measurements to macroscopic predictions of clinical progression. J Theor Biol 301:122–140
Article MathSciNet Google Scholar
Buske P, Galle J, Barker N, Aust G, Clevers H, Loeffler M (2011) A comprehensive model of the spatio-temporal stem cell and tissue organisation in the intestinal crypt. PLoS Comput Biol 7(1):e1001045
Article Google Scholar
Dunn SJ, Näthke IS, Osborne JM (2013) Computational models reveal a passive mechanism for cell migration in the crypt. PLoS ONE 8(11):e80516
Article Google Scholar
van der Wath RC, Gardiner BS, Burgess AW, Smith DW (2013) Cell organisation in the colonic crypt: a theoretical comparison of the pedigree and niche concepts. PLoS ONE 8:e73204
Article Google Scholar
Fletcher AG, Breward CJW, Jonathan Chapman S (2012) Mathematical modeling of monoclonal conversion in the colonic crypt. J Theor Biol 300:118–133
Article MathSciNet Google Scholar
Dunn SJ, Fletcher AG, Chapman SJ, Gavaghan DJ, Osborne JM (2012) Modelling the role of the basement membrane beneath a growing epithelial monolayer. J Theor Biol 298:82–91
Article Google Scholar
Van Leeuwen IMM, Mirams GR, Walter A, Fletcher A, Murray P, Osborne J, Varma S, Young SJ, Cooper J, Doyle B, Pitt-Francis J, Momtahan L, Pathmanathan P, Whiteley JP, Chapman SJ, Gavaghan DJ, Jensen OE, King JR, Maini PK, Waters SL, Byrne HM (2009) An integrative computational model for intestinal tissue renewal. Cell Prolif 42(5):617–636
Article Google Scholar
Smallwood R (2009) Computational modeling of epithelial tissues. Wiley Interdiscip Rev Syst Biol Med 1(2):191–201
Article Google Scholar
van Leeuwen IMM, Edwards CM, Ilyas M, Byrne HM (2007) Towards a multiscale model of colorectal cancer. World J Gastroenterol 13(9):1399–1407
Article Google Scholar
Zaman MH, Kamm RD, Matsudaira P, Lauffenburger DA (2005) Computational model for cell migration in three-dimensional matrices. Biophys J 89:1389–1397
Article Google Scholar
Rangarajan R, Zaman MH (2008) Modeling cell migration in 3D. Cell Adhes Migr 2(2):106–109
Article Google Scholar
Rey R, García-Aznar JM (2013) A phenomenological approach to modelling collective cell movement in 2D. Biomech Model Mechanobiol 12:1089–1100
Article Google Scholar
Vermolen FJ, Gefen A (2012) A semi-stochastic cell-based formalism to model the dynamics of migration of cells in colonies. Biomech Model Mechanobiol 11:183–195
Article Google Scholar
Harjanto D, Zaman MH (2013) Modeling extracellular matrix reorganization in 3D environments. PLoS ONE 8:e52509
Article Google Scholar
Kim T, Hwang W, Lee H, Kamm RD (2009) Computational analysis of viscoelastic properties of crosslinked actin networks. PLoS Comput Biol 5:e1000439
Article Google Scholar
Basan M, Elgeti J, Hannezo E, Rappel W-J, Levine H (2013) Alignment of cellular motility forces with tissue flow as a mechanism for efficient wound healing. Proc Natl Acad Sci USA 110(7):2452–2459
Article Google Scholar
Palsson E, Othmer HG (2000) A model for individual and collective cell movement in Dictyosteliumdiscoideum. Proc Natl Acad Sci USA 97:10448–10453
Article Google Scholar
Palsson E (2001) A three-dimensional model of cell movement in multicellular systems. Future Gen Comput Syst 17:835–852
Article MATH Google Scholar
Dallon JC, Othmer HG (2004) How cellular movement determines the collective force generated by the Dictyostelium discoideum slug. J Theor Biol 231:203–222
Article MathSciNet Google Scholar
Krinner A (2010) Spherical individual cell-based models: limits and applications. PhD Thesis, University of Leipzig
Honda H, Yamanaka H, Dan-Sohkawa M (1984) A computer simulation of geometrical configurations during cell division. J Theor Biol 106:423–435
Article Google Scholar
Drasdo D, Höhme S (2003) Individual-based approaches to birth and death in avascu1ar tumors. Math Comput Model 37:1163–1175
Article MATH Google Scholar
Drasdo D (2000) Buckling instabilities of one-layered growing tissues. Phys Rev Lett 84:4244–4247
Article Google Scholar
Schlüter DK, Ramis-Conde I, Chaplain MAJ (2015) Multi-scale modelling of the dynamics of cell colonies: insights into cell-adhesion forces and cancer invasion from in silico simulations. J R Soc Interface 12:20141080
Article Google Scholar
Kempf H, Bleicher M, Meyer-Hermann M (2010) Spatio-temporal cell dynamics in tumour spheroid irradiation. Eur Phys J D 60(1):177–193
Article Google Scholar
Boulanger AC (2010) Agent-based model–continuum model in tumor growth (INRIA intership report). Tech. Rep., August 2009
Milde F, Tauriello G, Haberkern H, Koumoutsakos P (2014) SEM++: a particle model of cellular growth, signaling and migration. Comput Part Mech 1(2):211–227
Article Google Scholar
Rejniak KA (2005) A single-cell approach in modeling the dynamics of tumor microregions. Math Biosci Eng (MBE) 2:643–655
Article MathSciNet MATH Google Scholar
Discher DE, Boal DH, Boey SK (1998) Simulations of the erythrocyte cytoskeleton at large deformation. II. Micropipette aspiration. Biophys J 75:1584–1597
Article Google Scholar
Fedosov DA, Caswell B, Karniadakis GE (2010) Systematic coarse-graining of spectrin-level red blood cell models. Comput Methods Appl Mech Eng 199(2932):1937–1948
Article MathSciNet MATH Google Scholar
Van Liedekerke P, Tijskens E, Ramon H, Ghysels P, Samaey G, Roose D (2010) Particle-based model to simulate the micromechanics of biological cells. Phys Rev E 81:61906–61915
Article Google Scholar
Sandersius SA, Newman TJ (2008) Modeling cell rheology with the subcellular element model. Phys Biol 5:015002
Article Google Scholar
Boal D (2012) Mechanics of the cell, 2nd edn. Cambridge University Press, Cambridge
Book Google Scholar
Hansen JC, Skalak R, Chien S, Hoger A (1996) An elastic network model based on the structure of the red blood cell membrane skeleton. Biophys J 70:146–166
Article Google Scholar
Buenemann M, Lenz P (2008) Elastic properties and mechanical stability of chiral and filled viral capsids. Phys Rev E 78:051924
Article MathSciNet Google Scholar
Rejniak KA (2007) An immersed boundary framework for modelling the growth of individual cells: an application to the early tumour development. J Theor Biol 247:186–204
Article MathSciNet Google Scholar
Rejniak KA, Dillon RH (2007) A single cell-based model of the ductal tumour microarchitecture. Comput Math Methods Med 8(1):51–69
Article MathSciNet MATH Google Scholar
Dillon R, Owen M (2008) A single-cell-based model of multicellular growth using the immersed boundary method. AMS Contemp Math 466:1–15
Article MathSciNet MATH Google Scholar
Rejniak KA, Anderson ARA (2008) A computational study of the development of epithelial acini: I. Sufficient conditions for the formation of a hollow structure. Bull Math Biol 70:677–712
Li J, Dao M, Lim CT, Suresh S (2005) Spectrin-level modeling of the cytoskeleton and optical tweezers stretching of the erythrocyte. Biophys J 88(5):3707–3719
Article Google Scholar
Pivkin IV, Karniadakisa GE (2008) Accurate coarse-grained modeling of red blood cells. Phys Rev Lett 101:118105
Article Google Scholar
Hosseini M, Feng JJ (2009) A particle-based model for the transport of erythrocytes in capillaries. Chem Eng Sci 64:4488–4497
Article Google Scholar
Fedosov D, Caswell B, Suresh S, Karniadakis GE (2011) Quantifying the biophysical characteristics of Plasmodium-falciparum-parasitized red blood cells in microcirculation. Proc Natl Acad Sci USA 108:35–39
Article Google Scholar
Peng Z, Li X IV, Pivkin Dao M, Karniadakis GE, Suresh S (2013) Lipid bilayer and cytoskeletal interactions in a red blood cell. Proc Natl Acad Sci USA 110(33):13356–13361
Article Google Scholar
Ingber DE (2003) Tensegrity I. Cell structure and hierarchical systems biology. J Cell Sci 116:1157–1173
Article Google Scholar
Sandersius SA, Weijer CJ, Newman TJ (2011) Emergent cell and tissue dynamics from subcellular modeling of active biomechanical processes. Phys Biol 8:45007
Article Google Scholar
Jamali Y, Azimi M, Mofrad MRK (2010) A sub-cellular viscoelastic model for cell population mechanics. PLoS ONE 5(8):e12097
Article Google Scholar
Van Liedekerke P, Roose D, Ramon H, Ghysels P, Tijskens E, Samaey G (2011) Mechanisms of soft cellular tissue bruising. A particle based simulation approach. Soft Matter 7:3580
Article Google Scholar
Tamura K, Komura S, Kato T (2004) Adhesion induced buckling of spherical shells. J Phys Condens Matter 16:L421–L428
Article Google Scholar
Murrell MP, Voituriez R, Joanny J-F, Nassoy P, Sykes C, Gardel ML (2014) Liposome adhesion generates traction stress. Nat Phys 10:163–169
Article Google Scholar
Kim M-C, Neal DM, Kamm RD, Asada HH (2013) Dynamic modeling of cell migration and spreading behaviors on fibronectin coated planar substrates and micropatterned geometries. PLoS Comput Biol 9(2):e1002926
Article MathSciNet Google Scholar
Tozluolu M, Tournier AL, Jenkins RP, Hooper S, Bates PA, Sahai E (2013) Matrix geometry determines optimal cancer cell migration strategy and modulates response to interventions. Nat Cell Biol 15(7):751–762
Article Google Scholar
Bentley K, Mariggi G, Gerhardt H, Bates PA (2009) Tipping the balance: robustness of tip cell selection, migration and fusion in angiogenesis. PLoS Comput Biol 5:e1000549
Article Google Scholar
Nagai T, Honda H (2001) A dynamic cell model for the formation of epithelial tissues. Philos Mag Part B 81:699–719
Article Google Scholar
Nagai T, Honda H (2009) Computer simulation of wound closure in epithelial tissues: cell basal-lamina adhesion. Phys Rev E 80:061903
Article Google Scholar
Farhadifar R, Röper J-C, Aigouy B, Eaton S, Jülicher F (2007) The influence of cell mechanics, cell-cell interactions, and proliferation on epithelial packing. Curr Biol (CB) 17:2095–2104
Article Google Scholar
Hilgenfeldt S, Erisken S, Carthew RW (2008) Physical modeling of cell geometric order in an epithelial tissue. Proc Natl Acad Sci USA 105:907–911
Article Google Scholar
Manning ML, Foty RA, Steinberg MS, Schoetz E-M (2010) Coaction of intercellular adhesion and cortical tension specifies tissue surface tension. Proc Natl Acad Sci USA 107(28):12517–12522
Article Google Scholar
Rudge T, Haseloff J (2005) Advances in artificial life, vol 3630. Lecture notes in computer science. Springer, Berlin
Fletcher AG, Osborne JM, Maini PK, Gavaghan DJ (2013) Implementing vertex dynamics models of cell populations in biology within a consistent computational framework. Prog Biophys Mol Biol 113:299–326
Article Google Scholar
Lim CT, Zhou EH, Quek ST (2006) Mechanical models for living cells–a review. J Biomech 39:195–216
Article Google Scholar
Zhou EH, Xu F, Quek ST, Lim CT (2012) A power-law rheology-based finite element model for single cell deformation. Biomech Model Mechanobiol 11:1075–1084
Article Google Scholar
Trepat X, Lenormand G, Fredberg JJ (2008) Universality in cell mechanics. Soft Matter 4:1750
Article Google Scholar
Wottawah F, Schinkinger S, Lincoln B, Ananthakrishnan R, Romeyke M, Guck J, Kas J (2005) Optical rheology of biological cells. Phys Rev Lett 94:98103
Article Google Scholar
Roose T, Chapman SJ, Maini PK (2007) Mathematical models of avascular tumor growth. SIAM Rev 49:179–208
Article MathSciNet MATH Google Scholar
Byrne HM, Alarcon T, Owen MR, Webb SD, Maini PK (2006) Modelling aspects of cancer dynamics: a review. Philos Trans Ser A Math Phys Eng Sci 364:1563–1578
Article MathSciNet Google Scholar
Roose T, Netti PA, Munn LL, Boucher Y, Jain RK (2003) Solid stress generated by spheroid growth estimated using a linear poroelasticity model. Microvasc Res 66:204–212
Article Google Scholar
Ambrosi D, Preziosi L (2009) Cell adhesion mechanisms and stress relaxation in the mechanics of tumours. Biomech Model Mechanobiol 8:397–413
Article Google Scholar
Preziosi L, Tosin A (2009) Multiphase modelling of tumour growth and extracellular matrix interaction: mathematical tools and applications. J Math Biol 58:625–656
Article MathSciNet MATH Google Scholar
Preziosi L, Ambrosi D, Verdier C (2010) An elasto-visco-plastic model of cell aggregates. J Theor Biol 262(1):35–47
Article MathSciNet Google Scholar
Sciumè G, Santagiuliana R, Ferrari M, Decuzzi P, Schrefler BA (2014) A tumor growth model with deformable ECM. Phys Biol 11:065004
Article Google Scholar
Alber M, Chen N, Glimm T, Lushnikov PM (2006) Multiscale dynamics of biological cells with chemotactic interactions: From a discrete stochastic model to a continuous description. Phys Rev E Stat Nonlinear Soft Matter Phys 73(5):051901
Article MathSciNet Google Scholar
De Masi A, Luckhaus S, Presutti E (2007) Two scales hydrodynamic limit for a model of malignant tumor cells. Ann Henri Poincare (B) Probab. Stat 43(3):257–297
MATH Google Scholar
Stevens A (2000) The derivation of chemotaxis equations as limit dynamics of moderately interacting stochastic many-particle systems. SIAM J Appl Math 61(1):183–212
Article MathSciNet MATH Google Scholar
Chauviere A, Hatzikirou H, Kevrekidis IG, Lowengrub JS, Cristini V (2012) Dynamic density functional theory of solid tumor growth: preliminary models. AIP Adv 2:11210
Article Google Scholar
D’Alessandro LA, Hoehme S, Henney A, Drasdo D, Klingmüller U (2015) Unraveling liver complexity from molecular to organ level: challenges and perspectives. Prog Biophys Mol Biol 117:78–86
Article Google Scholar
Galle J, Preziosi L, Tosin A (2009) Contact inhibition of growth described using a multiphase model and an individual cell based model. Appl Math Lett 22:1483–1490
Article MathSciNet MATH Google Scholar
Murray P, Edwards C, Tindall M, Maini P (2009) From a discrete to a continuum model of cell dynamics in one dimension. Phys Rev E 80:031912
Article Google Scholar
Frieboes HB, Jin F, Chuang Y-L, Wise SM, Lowengrub JS, Cristini V (2010) Three-dimensional multispecies nonlinear tumor growth-II: tumor invasion and angiogenesis. J Theor Biol 264:1254–1278
Article MathSciNet Google Scholar
Byrne HM, Osborne JM, Walter A, Kershaw SK, Mirams GR, Fletcher AG, Pathmanathan P, Gavaghan D, Jensen OE, Maini PK (2010) A hybrid approach to multi-scale modelling of cancer. Philos Trans R Soc A 368:5013–5028
Article MathSciNet MATH Google Scholar
Kim Y, Stolarska MA, Othemer HG (2007) A hybrid model for tumor spheroid grwoth in vitro I: theoretical development and earliy results. Math Models Methods Appl Sci 17:1773–1798
Article MathSciNet MATH Google Scholar
Dallon JC, Sherratt JA, Maini PK (1999) Mathematical modelling of extracellular matrix dynamics using discrete cells: fiber orientation and tissue regeneration. J Theor Biol 199:449–471
Article Google Scholar
Cumming BD, McElwain DLS, Upton Z (2010) A mathematical model of wound healing and subsequent scarring. J R Soc Interface 7:19–34
Article Google Scholar
Yang L, Witten TM, Pidaparti RM (2013) A biomechanical model of wound contraction and scar formation. J Theor Biol 332:228–248
Article MathSciNet MATH Google Scholar
Milde F, Bergdorf M, Koumoutsakos P (2008) A hybrid model for three-dimensional simulations of sprouting angiogenesis. Biophys J 95:3146–3160
Article Google Scholar
Monaghan J (2012) Smoothed particle hydrodynamics and its diverse applications. Annu Rev Fluid Mech 44:323–346
Article MathSciNet MATH Google Scholar
Gholami B, Comerford A, Ellero M (2013) A multiscale SPH particle model of the near-wall dynamics of leukocytes in flow. Int J Numer Methods Biomed Eng 30(1):83–102
Article MathSciNet Google Scholar
Tanaka N, Takano T (2005) Microscoic scale simulation of blood flow using SPH method. Int J Comput Methods 02(04):555–568
Article MATH Google Scholar
Angermann BR, Klauschen F, Garcia AD, Prustel T, Zhang F, Germain RN, Meier-Schellersheim M (2012) Computational modeling of cellular signaling processes embedded into dynamic spatial contexts. Nat Methods 9(3):283–289
Article Google Scholar
Figueredo GP, Joshi TV, Osborne JM, Byrne HM, Owen MR (2013) On-lattice agent-based simulation of populations of cells within the open-source Chaste framework. Interface Focus 3:20120081
Article Google Scholar
Sütterlin T, Kolb C, Dickhaus H, Jäger D, Grabe N (2013) Bridging the scales: semantic integration of quantitative SBML in graphical multi-cellular models and simulations with EPISIM and COPASI. Bioinformatics (Oxford, England) 29:223–229
Cytowski M, Szymanska Z (2014) Large-scale parallel simulations of 3D cell colony dynamics. Comput Sci Eng 16:86–95
Article Google Scholar
Cytowski M, Szymanska Z (2015) Large scale parallel simulations of 3-D cell colony dynamics. II. Coupling with continuous description of cellular environment. Comput Sci Eng 99:1–6
Google Scholar
Kang S, Kahan S, McDermott J, Flann N, Shmulevich I (2014) Biocellion: accelerating computer simulation of multicellular biological system models. Bioinformatics (Oxford, England) 30:3101–3108
Merks RMH, Guravage M, Inzé D, Beemster GTS (2011) VirtualLeaf: an open-source framework for cell-based modeling of plant tissue growth and development. Plant Physiol 155(2):656–666
Article Google Scholar
Tanaka S, Sichau D, Iber D (2015) LBIBCell: a cell-based simulation environment for morphogenetic problems. Bioinformatics 31(14):2340–2347
Article Google Scholar
Hucka M, Finney A, Sauro HM, Bolouri H, Doyle JC, Kitano H, Arkin AP, Bornstein BJ, Bray D, Cornish-Bowden A, Cuellar AA, Dronov S, Gilles ED, Ginkel M, Gor V, Goryanin II, Hedley WJ, Hodgman TC, Hofmeyr JH, Hunter PJ, Juty NS, Kasberger JL, Kremling A, Kummer U, Le Novère N, Loew LM, Lucio D, Mendes P, Minch E, Mjolsness ED, Nakayama Y, Nelson MR, Nielsen PF, Sakurada T, Schaff JC, Shapiro BE, Shimizu TS, Spence HD, Stelling J, Takahashi K, Tomita M, Wagner J, Wang J (2003) The systems biology markup language (SBML): A medium for representation and exchange of biochemical network models. Bioinformatics 19(4):524–531
Article Google Scholar
Andasari V, Roper RT, Swat MH, Chaplain MAJ (2012) Integrating intracellular dynamics using CompuCell 3D and Bionetsolver: applications to multiscale modelling of cancer cell growth and invasion. PLoS ONE 7:e33726
Article Google Scholar
Hoops S, Gauges R, Lee C, Pahle J, Simus N, Singhal M, Xu L, Mendes P, Kummer U (2006) COPASI–a COmplex PAthway SImulator. Bioinformatics 22(24):3067–3074
Marsaglia G (1972) Choosing a point from the surface of a sphere. Ann Math Stat 43:645–646
Article MATH Google Scholar
Marmottant P, Mgharbel A, Käfer J, Audren B, Rieu J-P, Vial J-C, van der Sanden B, Marée AFM, Graner F, Delanoë-Ayari H (2009) The role of fluctuations and stress on the effective viscosity of cell aggregates. Proc Natl Acad Sci USA 106:17271–17275
Article Google Scholar
Sciumè G, Shelton S, Gray W, Miller C, Hussain F, Ferrari M, Decuzzi P, Schrefler B (2013) A multiphase model for three-dimensional tumor growth. New J Phys 15:015005
Article Google Scholar
Subramaniyan AK, Sun C (2008) Continuum interpretation of virial stress in molecular simulations. Int J Solids Struct 45:4340–4346
Article MATH Google Scholar
Harvey DG, Fletcher AG, Osborne JM, Pitt-Francis J (2015) A parallel implementation of an off-lattice individual-based model of multicellular populations. Comput Phys Commun 192:130–137
Article MathSciNet Google Scholar
Bittig T, Wartlick O, Kicheva A, González-Gaitárr M, Jülicher F (2008) Dynamics of anisotropic tissue growth. New J Phys 10:063001
Article Google Scholar
Landau LD, Pitaevskii LP, Lifshitz EM, Kosevich AM (1986) Theory of elasticity, 3rd edn, vol 7 (Theoretical physics). Butterworth-Heinemann, Oxford
Montel F, Delarue M, Elgeti J, Malaquin L, Basan M, Risler T, Cabane B, Vignjevic D, Prost J, Cappello G, Joanny J-FBC (2011) Stress clamp experiments on multicellular tumor spheroids. Phys Rev Lett 107:188102
Hatzikirou H, Brusch L, Deutsch A (2014) Form cellular automaton rules to a macroscopic mean-field description. Acta Phys Pol Ser B Suppl 3:399–416
Google Scholar
Scianna M, Preziosi L (2013) A cellular Potts model for the MMP-dependent and -independent cancer cell migration in matrix microtracks of different dimensions. Comput Mech 53(3):485–497
Article MATH Google Scholar
Savill NJ, Hogeweg P (1997) Modelling morphogenesis: from single cells to crawling slugs. J Theor Biol 184:229–235
Article Google Scholar
Yao M, Goult BT, Chen H, Cong P, Sheetz MP, Yan J (2014) Mechanical activation of vinculin binding to talin locks talin in an unfolded conformation. Sci Rep 4:4610
Google Scholar
Brú A, Pastor J, Fernaud I, Brú I, Melle S, Berenguer C (1998) Super-rough dynamics on tumor growth. Phys Rev Lett 81(18):4008–4011
Article Google Scholar
Eissing T, Kuepfer L, Becker C, Block M, Coboeken K, Gaub T, Goerlitz L, Jaeger J, Loosen R, Ludewig B, Meyer M, Niederalt C, Sevestre M, Siegmund H-U, Solodenko J, Thelen K, Telle U, Weiss W, Wendl T, Willmann S, Lippert J (2011) A computational systems biology software platform for multiscale modeling and simulation: integrating whole-body physiology, disease biology, and molecular reaction networks. Front Physiol 2:4
Article Google Scholar

Download references

Acknowledgments

The authors gratefully acknowledge fruitful discussions with Pierre Nassoy and Irene Vignon-Clementel, and funding by EU NOTOX, ANR INVADE, ANR PHYSCANCER, BMBF Virtual Liver, and BMBF Cancersys/Lungsys. The authors acknowledge the contributions of Emanuele Leoncini to Figs. 9 and 10, Axel Krinner to Fig. 11, and helpful discussions with Anne-Celine Boulanger, either within their master thesis at INRIA (EL, ACB) or PhD thesis at University of Leipzig (AK) with DD.

Author information

Authors and Affiliations

INRIA, Domaine Rocquencourt, B.P. 105, 78153, Le Chesnay Cedex, France
P. Van Liedekerke, M. M. Palm, N. Jagiella & D. Drasdo

Authors

P. Van Liedekerke
View author publications
You can also search for this author in PubMed Google Scholar
M. M. Palm
View author publications
You can also search for this author in PubMed Google Scholar
N. Jagiella
View author publications
You can also search for this author in PubMed Google Scholar
D. Drasdo
View author publications
You can also search for this author in PubMed Google Scholar

Corresponding author

Correspondence to D. Drasdo.

Additional information

P. Van Liedekerke and D. Drasdo contributed equally to this work.

Appendices

Appendix 1: Validation of different stress measures

To illustrate the possibilities in this approach, let us consider different cases where a cell is in contact with other cells (see Fig. 23). In the first case, a cell is squeezed between two other cells along the X axis. Applying Eq. 39, the stress tensor becomes

$$\begin{aligned} \sigma = \begin{pmatrix} \sigma _{xx} &{} \quad 0 &{} \quad 0 \\ 0 &{} \quad 0 &{} \quad 0 \\ 0 &{} \quad 0 &{} \quad 0 \end{pmatrix} \end{aligned}$$

(57)

where $\sigma _{xx}=\frac{3f}{2\pi r^2}$, from which the pressure p can be extracted:

$$\begin{aligned} p = \frac{1}{3}\text {tr}(\mathbf \sigma )= \frac{f}{2\pi r^2}. \end{aligned}$$

(58)

In addition we can define the “equivalent” stress $\sigma _\mathrm{eq}$ as

$$\begin{aligned}&\sigma _\mathrm{eq}\nonumber \\&\quad \!= \, \sqrt{0.5\left[ (\sigma _{xx}-\sigma _{yy})^2 \!+\! (\sigma _{yy}-\sigma _{zz})^2 \!+\! (\sigma _{xx}-\sigma _{zz})^2 \right] \!+\! 3\left[ \sigma ^2_{xy} \!+\!\sigma ^2_{yz} \!+\!\sigma ^2_{xz} \right] },\nonumber \\ \end{aligned}$$

(59)

which only counts the deformation. In case of the line contact this yields $\sigma _\mathrm{eq}=\frac{3f}{2\pi r^2}$. In the second case the cell has also contact along the Y axis, and then we find an increased pressure $p=\frac{1}{3}(\sigma _{xx}+\sigma _{yy})=\frac{f}{\pi r^2}$ while $\sigma _\mathrm{eq}$ remains the same. However, an additional contact along the Z axis will yield $p=\frac{2f}{\pi r^2}$ but $\sigma _\mathrm{eq} = 0 $ so the deformational stress vanishes because of symmetry. This is the equivalent situation to a cube compressed equally from all sides. Thus in agreement with one would expect, the pressure on the cell can increase with the number of contacts, while the deviatoric part may even decrease or vanish depending on how the contacts are distributed along the cell’s surface. This extra information is an important advantage over the previous methods. Note however, that this tensor has no off diagonal components, as no shear forces are taken into account. Note that if the cells would have isotropic contact area (infinite contact points), the sum in Eq. 39 can be replaced by an integral while the force should be replaced by a force density $\mathbf P (x,y,z) = P \mathbf e _\mathrm{r}(x,y,z)$. In this case we retrieve the pressure:

$$\begin{aligned} p= \frac{1}{3}\mathrm{tr}\left( \frac{\int {\mathbf{r}\otimes \mathbf{P}\mathrm{d}S}}{V}\right) = \frac{1}{3} \cdot \frac{r \cdot 4\pi r^2}{4/3 \pi r^3} P = P \end{aligned}$$

(60)

To test the virial stress implementation, we have performed a study simulating a spheroid with “inert” cells, exhibiting no growth nor division, which is gradually compressed by a capsule with decreasing radius. The cells have a Young modulus of $250\,\hbox {Pa}$ and Poisson ration of 0.49, and do not adhere to each other. We assume that the capsule wall is rigid. At the point of confluence, the cells start exerting a total force $F^\mathrm{cap}=\sum _{j}F^{\mathrm{cap, cell} j}_{j}$ on the capsule which will increase as the capsule shrinks further. This force is converted to an average pressure $P^\mathrm{cap}=F^\mathrm{cap}/A$ where A is the inner surface area of the capsule. Figure 24 indicates that the maximal directions in the cells at the boundary point are perpendicular to the capsule wall, while the minimal stress directions are parallel. In Fig. 24c, we depict a plot of the average pressure increase versus volume decrease of the capsule, giving an indication of the differences in stress measures that can be expected between the average viral stress $\langle p \rangle =\sum _{j}p_j/N$ per cell, and the capsule wall pressure $P^\mathrm{cap}$. It can be seen that the deviations are relatively small for moderate compression. At high compression a discrepancy arises because the cell overlaps become too high, a problem that is well known in CBM under compression (see Sect. 3.1). The stress measure Eq. 38 results in a much higher value because a “local” stresses are computed and summed.

Appendix 2: Numerical estimation of the stochastic force term in the Langevin equation

To obtain a force-based numerical representation of a Brownian motion for a particle or cell in a low Reynolds number environment in a isotropic medium, we start from the approximation for the square of the expectation value of the force vector norm in n dimensions:

$$\begin{aligned} \langle ||\underline{F}^\mathrm{R}||^2 \rangle \approx \frac{ 6(2 n)\pi K_\mathrm{B} T \eta R}{\Delta t}, \end{aligned}$$

(61)

where R is the cell radius, $\gamma $ is the friction coefficient, $\eta $ is the viscosity of the medium, and$\Delta {t}$ is the timestep in the simulation. Using Stokes’ formula, $\gamma = 6 \pi \eta R$, we get

$$\begin{aligned} \langle ||\underline{F}^\mathrm{R}||\rangle \approx \sqrt{\frac{ (2 n) K_\mathrm{B} T \gamma }{\Delta t}}, \end{aligned}$$

(62)

To get the actual force vector, at each time-step in the simulation, we need to pick a random force unit vector $\underline{n}^\mathrm{R}$, and scale it in the following way:

$$\begin{aligned} \underline{F}^\mathrm{R}= \sqrt{\frac{2 n\ K_\mathrm{B} T \gamma }{\Delta t}}\underline{n}^\mathrm{R}, \end{aligned}$$

(63)

Using $D=k_\mathrm{B}T/\gamma $ and introducing $S= 2 n D \gamma ^2$ we can also write

$$\begin{aligned} \underline{F}^\mathrm{R}= \sqrt{({S}/\Delta {t})}\underline{n}^\mathrm{R}, \end{aligned}$$

(64)

where D is the diffusion coefficient. The vector $\underline{n}^\mathrm{R}$ must be taken with care to obtain a true random direction (just normalizing three random numbers will introduce artefacts), see e.g., [254]. Note finally that for cells the micromotility, represented by the diffusion coefficient D, and the frictional resistance $\gamma $, cannot be expected to be related by the Einstein relation $(D=K_\mathrm{B}T/\gamma )$ as a strict $K_\mathrm{B}T$ equivalent in cells does not exist. In active cell movement, the equivalent of this parameter would be controlled by the cell itself as well as by other influences (e.g., the cells’ environment).

Rights and permissions

Reprints and permissions

About this article

Cite this article

Van Liedekerke, P., Palm, M.M., Jagiella, N. et al. Simulating tissue mechanics with agent-based models: concepts, perspectives and some novel results. Comp. Part. Mech. 2, 401–444 (2015). https://doi.org/10.1007/s40571-015-0082-3

Download citation

Received: 12 June 2015
Revised: 25 September 2015
Accepted: 20 October 2015
Published: 26 November 2015
Issue Date: December 2015
DOI: https://doi.org/10.1007/s40571-015-0082-3

Simulating tissue mechanics with agent-based models: concepts, perspectives and some novel results

Abstract

Similar content being viewed by others

Free and Interfacial Boundaries in Individual-Based Models of Multicellular Biological systems

Introduction to Computational Modeling of Multicellular Tissues

Dynamical Rheological Properties of In-Silico Epithelial Tissue by Vertex Models

1 Introduction

2 Lattice models

2.1 Cellular automaton with one cell per lattice site (type B)

2.2 Cellular automaton with many cells per lattice site (type A)

2.3 Lattice gas cellular automata (type D)

2.4 Cellular automaton with many lattice sites per cell (type C; cellular Potts models)

3 Off-lattice models

3.1 Center-based models

3.1.1 Basic concepts

3.1.2 Achievements, limitations and practical use

3.2 Deformable cell models and vertex models

3.2.1 Basic concepts

3.2.2 Achievements, limitations and practical use

4 Hybrid discrete-continuum models

5 Available tools and software development for modeling

6 Conclusion

Notes

References

Acknowledgments

Author information

Authors and Affiliations

Corresponding author

Additional information

Appendices

Appendix 1: Validation of different stress measures

Appendix 2: Numerical estimation of the stochastic force term in the Langevin equation

Rights and permissions

About this article

Cite this article

Keywords

Navigation

Simulating tissue mechanics with agent-based models: concepts, perspectives and some novel results

Abstract

Similar content being viewed by others

Free and Interfacial Boundaries in Individual-Based Models of Multicellular Biological systems

Introduction to Computational Modeling of Multicellular Tissues

Dynamical Rheological Properties of In-Silico Epithelial Tissue by Vertex Models

1 Introduction

2 Lattice models

2.1 Cellular automaton with one cell per lattice site (type B)

2.2 Cellular automaton with many cells per lattice site (type A)

2.3 Lattice gas cellular automata (type D)

2.4 Cellular automaton with many lattice sites per cell (type C; cellular Potts models)

3 Off-lattice models

3.1 Center-based models

3.1.1 Basic concepts

3.1.2 Achievements, limitations and practical use

3.2 Deformable cell models and vertex models

3.2.1 Basic concepts

3.2.2 Achievements, limitations and practical use

4 Hybrid discrete-continuum models

5 Available tools and software development for modeling

6 Conclusion

Notes

References

Acknowledgments

Author information

Authors and Affiliations

Corresponding author

Additional information

Appendices

Appendix 1: Validation of different stress measures

Appendix 2: Numerical estimation of the stochastic force term in the Langevin equation

Rights and permissions

About this article

Cite this article

Share this article

Keywords

Search

Navigation