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

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.

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.

Fig. 1

Construction scheme of unstructured lattice for cellular automaton simulations with type B automaton [72]. From left to right (i) a square lattice is generated with one point in each square. (ii) the nearest neighbors are connected by lines generating a Delaunay triangulation; (iii) the Voronoi tessellation representing the dual graph of the Delaunay triangulation is generated from it, by drawing the perpendicular bisection on every edge of the Delaunay triangulation; (iv) the final Voronoi tessellation with the construction points. Notice that each point in space closer to one construction point (each of the black points is a construction point) than to any other construction point belongs to the Voronoi cell of that construction point. Points with equal distance to more than one construction point form the borders between neighboring Voronoi cells (black lines). If a biological cell changes site, the volume (area in Fig. 2d) changes too, as the volume (area in Fig. 2d) is an intrinsic property of the lattice site, not of the cell

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}$.

Fig. 2

a Determination of free lattice site with circle of radius $\Delta R$ [72]. b Fundamental processes during growth of multicellular spheroids [81]. c Typical growth scenario in absence of migration and death [72]

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).

Fig. 3

Left Simulations of a growing monolayer with a type B cellular automaton (points experiment, red line simulations with type B cellular automaton, blue line with squares off-lattice CBM (see later section). The parameters for the type B cellular automaton are $\Delta L=9$ (in lattice units, here cell diameters), $m=60$ (Erlang parameter), $\phi = 0$ (migration rate). The experimental growth curve can be explained also by other parameter sets, and other mechanisms [72]. Additional experimental information is needed to uniquely fix model mechanisms and parameters. Middle Cell cycle duration distribution for type B cellular automaton and CBM for the simulations depicted in the middle. Right Simulation of a xenograft in mouse with a type A cellular automaton (points experiment, line simulation) [67]. Cells can be shifted over $\Delta L = 2.85$ lattice units, here $b=10$ cell diameters. NIH3T3 xenografts show only small necrotic lesions, hence nutrient limitation can be neglected even inside the tumor. This means the speed of tumor expansion is controlled by the thickness of the proliferating rim of size $\Delta L = 2.85 b = 28.5$ cell diameters. Here, $N_{\max }=1000$ cells and zero micro-motility (no migration) have been assumed. (Color figure online)

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].

Fig. 4

Growth figures in absence of migration and in case of $m=60$ states that have to be passed until division. a Square lattice with division along the x and y axes. b Square lattice with division along x, y axes and in diagonals (Moore neighborhood), c hexagonal lattice, d Voronoi lattice

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.

Fig. 5

Matching type A and B cellular automaton models. a Lin–log plot of growth kinetics for $N^{1/3}(t)$ for $m=1$, $\phi =0$ ($\phi =$phi in the figure legends). For short times growth is exponential. b The resulting growth kinetics $N^{1/3}(t)$ for $m=1$, $\phi =0$ in lin–lin plot shows that $N^{1/d}\propto R \propto t$ for long times. c Scaling of the maximum number of cells in a compartment of a type-A cellular automaton versus $k=\Delta L-1$ in a type B automaton. d Growth kinetics for $m=100$, $\phi = 0$. e Growth kinetics as in (b) but for $\phi = 0.1/h$. f Growth kinetics for varying migration rate (In e and f, the number of cells per compartment)

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.

Fig. 6

Integration of type A and type B cellular automaton in a hybrid model. In the proliferating rim, individual cell positions a spatially resolved by a type B automaton, otherwise, a type A automaton is used. Inspired by octrees, we divide the whole space into large compartments which themselves can be further sub-divided

Fig. 7

LGCA propagation scheme for a cubic lattice. Shown are nine lattice squares (called nodes), each with four velocity channels oriented towards the $+x, -x, +y, -y$-direction, and a rest channel in the center of node. The left field of nine lattice sites show the state before, the right after one propagation step where the cell speed is $\tilde{m}=1$. Filled circles denote occupied channels, empty circles empty channels (redrawn from Hatzikirou and Deutsch in [101]). For example, the central node on the left side has two occupied channels, i.e., is occupied by two particles. The maximum number of particles at one node is $4+1=5$ i.e., several particles can occupy the same node in space (as for type A cellular automata introduced previously)

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.

Fig. 8

Simulation of a growth process in the LGCA (courtesy of Hatzikirou and Deutsch [102]). Shown are the lines of each cell density at times $t=50, 100, 150$, and the growth kinetics. The center is filled with cells

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].

Fig. 9

Simulation of a growth process in a cellular Potts model defined on an irregular (a Voronoi) lattice in $d=2$ spatial dimensions [118]. The found growth kinetics is the same as for the CA type A and B models and the LGCA. Note, that $R\propto N^{1/d}$

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].

Fig. 10

Simulation of a growth process in a cellular Potts model defined on an irregular (a Voronoi) lattice in $d=2$ spatial dimensions [118]. A smaller Erlang parameter k and volume threshold $\varTheta$ accelerates growth (The Erlang parameter is here denoted by k, not by m)

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].

Fig. 11

Cellular Potts simulation of a row of three cells with a pushing force on the rightmost cell ($\mu = 10$, $V = 2500$, $\lambda = 10$, $J(\text {cell,cell}) = 10$, $J(\text {cell,medium}) = 5$). a Evolution of the three cells in a single simulation. b Evolution if the length of each cell averaged over 50 simulations (filled areas around the curves reflect the standard deviation). The cell length is estimated from the cell’s moment of inertia [123]

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.

Fig. 12

Left Center-based model with cells represented as spheres, defining the basic contact variables. Right Division of two cells in a CBM

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.

Fig. 13

Growth scenario in expanding monolayers if cell–substrate adhesion is insufficient to keep the cells in contact with the substrate (left), hence cells in the interior are pushed out of the basal layer by the forces exerted by their neighbor cells during growth and division [63] [62]. Blue cells denote cells with contact to the underlying substrate, red cells do not contact to the substrate anymore. Right Cells trying to grow by increasing their size exert forces on their neighbors. Those can lead to force components normal to the substrate (thick red arrow, overemphasized). Such a process represents a difficult obstacle to cellular automaton models of types A, B, D, as gradual displacements in combination with an increasing force component pointing perpendicular outside the layer are responsible for the macroscopic detachment of cells. (Color figure online)

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.

Fig. 14

Left Growth of a monolayer for pure pushing forces and pulling forces emerging if border cells migrate actively into the environment. Right Growth of a multicellular spheroid in an environment of elastic particles to mimic the mechanical resistance of agarose gel [144]. As in the cellular automaton models, growth is first exponential turning into linear growth

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¹ 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]).

Fig. 15

Simulation of a regenerating liver lobule, the smallest functional and anatomical repetitive unit of liver after drug induced damage in the center of the lobule. This scheme illustrates the general strategy of how liver architectural parameters obtained by image analysis of confocal micrographs (left) are used together with a quantification of dynamic processes in the liver (top row regeneration after intoxication with the drug carbontetrachloride; brown hepatocytes, blue central necrosis) to construct a dynamic model of the in vivo situation (bottom row regeneration after intoxication in the model; light rose quiescent hepatocytes, dark rose proliferating hepatocytes, brown hepatocytes carrying the enzyme Glutamine sythetase, red sinusoids, central and portal vein, blue portal artery). Left picture confocal micrographs after image processing; blue hepatocytes nuclei, white sinusoids [33, 152]. (Color figure online)

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).

Fig. 16

Cells represented as ellipsoids with fixed major axes can introduce orientation artefacts when in contact with other cells (compare situation a and b)

Fig. 17

Left Sketch of a center-based cell with filopodia [174]. Right Comparison of the random migration dynamics with the stochastic equation of motion using a white noise term (red symbols) or filopodia. In the simulation, a cell had eight filopodia. Their length was chosen from a Gaussian distribution, and their orientation chosen randomly and uncorrelated in time. After each move, the orientation was again chosen randomly. The force was assumed to be proportional to the filopodia length but this parameter had no strong influence. A difference in movement to the standard equation of motion with a noise term to mimic active motion could be found in presence of other cells, as filopodia could not stretch if the path was blocked by a neighbor cell. (Color figure online)

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).

Fig. 18

Left model of an intestinal crypt using a Voronoi tessellation. The model assumes periodic boundary conditions in x-direction, and a hard border at the bottom. Cells dividing in the bottom push other cells towards the top from where there are shed in the intestinal lumen [34]. Middle modified Voronoi-tessellation to account for isolated and border cells [61, 176]. The configuration emerged from sorting of two cell types A (triangles) and B with (differential) adhesion energies strongest between A–A, second strongest between A–B and weakest between B–B cells. Right Cross section of a crypt in a model [137]. Here, the one-cell-thick layer of the intestinal crypt has been stabilized by a bending force. If either the proliferation pressure of the cells dividing in the bottom of the crypts become too high because of a too high proliferation rate, or the bending stiffness is reduced by reduction of the cell layers’ elastic modulus, the layer folds [177]. Such a folding or bending is very difficult to represent in cellular automaton models while it emerges naturally out of polar adhesion forces in CBMs

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.

Fig. 19

a DCM representation of two adhering cells. Each cell has a nucleus, a coarse central oriented cytoskeleton (CSK), and a membrane/cortex. b Detail of the elements in a cell membrane triangulation

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.

Fig. 20

Simulation of a deformable cell model (erythrocyte) adhering to an adhesive substrate [140]. As the cell adheres, tensile stresses develop in the cortex at the free surface of the cell (red color). At the cell–substrate interface, compressive stresses develop (blue color), which can in turn deform the substrate. (Color figure online)

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.

Fig. 21

a left Deformable cell model as proposed by Rejniak [182]. The cell surface is represented by connected viscoelastic elements between the freely moving particles, while the continuum fields are computed on an Eulerian grid. The interaction domain $(\varOmega )$ between Lagrangian and Eulerian points is indicated by the grey disc. Right Cell model as described in Van Liedekerke et al. [202]. Both the cell surface and the medium are represented by particles. b A snapshot of a simulation coupling a CBM and DCM. The adhesion model used here is JKR for the CBM, and the Maugis–Dugdale model for the DCM and DCM–CBM interactions [140]. The stresses in the cell cortex (Eq. 53) are represented a colormap. c Schematic view of a proposed three level hybrid modeling, using DCM, CBM, and a continuum model. (Color figure online)

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.

Fig. 22

Table overview of the pro’s and cons of the agent-based models with regard to the cellular property they address. We assume the models are constructed in 3D

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.

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)

Fig. 23

Evolution of virial stress and pressure as cells are in contact along the X axis, (X, Y) axes and (X, Y, Z) axes

Fig. 24

Top Simulation of a spheroid compressed by a capsule. The arrows indicate the maximal (a) and minimal (b) stress directions. c Pressure measured during compression of a spheroid for pressure formulations Eqs. 40 and 38

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).