The biophysical nature of cells: potential cell behaviours revealed by analytical and computational studies of cell surface mechanics
 Ramiro Magno^{1, 2},
 Verônica A Grieneisen^{2} and
 Athanasius FM Marée^{2}Email author
DOI: 10.1186/s136280150022x
© Magno et al.; licensee BioMed Central. 2015
Received: 21 November 2014
Accepted: 24 April 2015
Published: 12 May 2015
Abstract
Background
The biophysical characteristics of cells determine their shape in isolation and when packed within tissues. Cells can form regular or irregular epithelial structures, round up and form clusters, or deform and attach to substrates. The acquired shape of cells and tissues is a consequence of (i) internal cytoskeletal processes, such as actin polymerisation and cortical myosin contraction, (ii) adhesion molecules within the cell membrane that interact with substrates and neighbouring cells, and (iii) processes that regulate cell volume. Although these processes seem relatively simple, when combined they unleash a rich variety of cellular behaviour that is not readily understandable outside a theoretical framework.
Methods
We perform a mathematical analysis of a commonly used class of model formalisms that describe cell surface mechanics using an energybased approach. Predictions are then confirmed through comparison with the computational outcomes of a Vertex model and 2D and 3D simulations of the Cellular Potts model.
Results
The analytical study reveals the complete possible spectrum of single cell behaviour and tissue packing in both 2D and 3D, by taking the typical core elements of cell surface mechanics into account: adhesion, cortical tension and volume conservation. We show that from an energybased description, forces and tensions can be derived, as well as the prediction of cell behaviour and tissue packing, providing an intuitive and biologically relevant mapping between modelling parameters and experiments.
Conclusions
The quantitative cellular behaviours and biological insights agree between the analytical study and the diverse computational model formalisms, including the Cellular Potts model. This illustrates the generality of energybased approaches for cell surface mechanics and highlights how meaningful and quantitative comparisons between models can be established. Moreover, the mathematical analysis reveals direct links between known biophysical properties and specific parameter settings within the Cellular Potts model.
Keywords
Cell surface mechanics Cellular Potts model Cell shape Cortical tension Membrane tension Adhesion Tissue packing Cellular vertexmodels Grid anisotropy 3D modelling Cell sorting Cellular dynamicsBackground
The assembly of a multicellular organism is achieved, alongside cell division, apoptosis and differentiation, by changing cell shape and tissue topology. Such cellular and tissue dynamics are highly determined by the cell’s physical properties. Therefore, to increase our understanding of morphogenesis, we need quantitative analyses and descriptions of the biomechanical properties of cells and tissues. However, despite the recent advances in this field [112], we still lack detailed knowledge as well as computational power to approach cell mechanics from a molecular perspective. Accordingly, there are several multicellular modelling approaches that waive the fine detail of specific molecular interactions by describing cell properties on a biophysical cellular level (see [13] for a review). Such approaches are analogous to the simplification of considering the surface tension of a liquid droplet, instead of explicitly specifying the molecular interactions involved. These cellular models maintain the desired mesoscopic cell features, while facilitating the comparison between simulations and experimental data, given its reduced number of variables and effective parameters [1416]. Although such mesoscopic approaches lack direct links between the molecular level and the macroscopic biological outcome, they remain powerful instruments for comprehending cellular behaviour, since they allow the investigation of novel regulatory elements and biochemical processes. Lecuit and Lenne [17] evaluated an overarching class of such mesoscopic biophysical models, the Cell Surface Mechanics (CSM) models. These describe how biologicallygenerated tensions – namely adhesionderived tension, cortical tension and tension due to the cell’s internal pressure – affect (i) the displacement of the cell’s membrane, (ii) modify the cell’s shape and (iii) influence its dynamics. Nevertheless, major issues prevent us from gaining optimal understanding from these models.
Firstly, it is not always clear how to quantitatively link experimental observations to model implementations. We will therefore discuss how biophysical properties are linked to subcellular components, such as the structure and mechanics of the cytoskeleton and the interactions between adhesion molecules among neighbouring cells. On a pragmatic level, this leads us to derive a protocol for CSM class simulations [17], to determine and modify expected cell behaviour in simulations and help establish a link with experimental observations.
A second source of confusion relates to the wide variety of computational implementations used to describe tissues. It is not yet evident how a modeltomodel comparison can be drawn, confusing the biology community when trying to compare common insights and features between different theoretical studies. Likewise, modellers are frequently not aware of the common features and artefacts of their implementations, making them susceptible to overlooking valuable information from parallel studies. We address this second concern in two ways: (i) we initially discuss in detail how a commonly used latticebased model formalism – the Cellular Potts model (CPM) – can (and should) be seen as being part of the CSM class, and how to determine its parameters and expected dynamics, whilst ensuring appropriate scaling; and (ii) we compare CPM simulations to analytical predictions as well as other CSM implementations, and highlight overlaps and interesting deviations.
We use the CPM as a reference model formalism from the CSM class. Originally, CPM was proposed to describe cellular dynamics based on differential adhesion, thereby defining tissue surface tension relations [18,19], to consolidate Steinberg’s Differential Adhesion Hypothesis (DAH) [20]. Classical CPM simulations, in which cells display an internal pressure and positive adhesion coupling constants, successfully captured cell sorting phenomena [18,19,21]. However, in the classical formulation the CPM did not include the effects of cortical tension. Ouchi et al. [22] have proposed an extension with a term representing the surface constraint tension (to modify the relationship between cell motility within aggregates and adhesion strength [23]), alongside using negative instead of positive coupling constants for the adhesion (discussed later). Recently, such a surface constraint term has been identified as the cortical tension, the third “ingredient” underlying cell surface mechanics [17,24]. Since then, cortical tension within the extended CPM has been used in several multicellular modelling studies, for example, in the description of corticaltension dependent cell shape alterations in dendritic cells [14]; in reverse engineering of the critical parameters determining Drosophila’s eye geometry [24]; and to understand germlayer organisation in zebrafish [15]. A number of other computational model formalisms have been proposed which, analogous to the CPM, use an energybased description for volume conservation, adhesion and cortical tension [2528]. They adopt, however, different ways of describing the cell shape as well as its dynamics. Such models have been used to explore the role of the interfacial tension, reproducing, through parameter changes of the CSM, different forms of epithelial cell packing [25]; plausible mechanisms underlying cell intercalation [27]; and again, the Drosophila’s eye geometry [26].
Surprisingly, despite the growing acceptance of the concept and importance of cortical tension from both modelling and cell biology communities, a clear explanation for the role of cortical tension in cell surface mechanics is still lacking. In fact, although several analytical studies regarding CSM have been published [29,30], it is still unclear what cell behaviour can be expected even when just a single cell is considered in terms of its cortical tension. For multicellular modelling, such a baseline understanding is vital to pinpoint the role of each biophysical parameter, allowing us to distinguish and predict beforehand individual and collectivecelldriven mechanisms. We used the CPM formalism to numerically verify our analytical predictions and to determine their implications in actual cell dynamics, highlighting the fact that the analytical results are independent of the specific simulation framework.
Methods
The cell surface mechanics model
where p and a are the perimeter and area of the cell (see Figure 1A). The function uses five parameters for the cellular properties: J, an energy per contact length due to adhesion to other cells or the surrounding medium; P, the membrane resting length; and A, the target cell area (resting area); while the constraints are modulated by the Lagrange multipliers λ _{ p } and λ _{ a } (comparable to elastic constants), which weigh the relative tension contributions of actinmyosin contraction and cell deformations, respectively. Although modifications of the above energy function could and have been proposed (see, e.g., [34]), almost all studies on CSM have been using this basic framework, sometimes further simplified (see, e.g., [19,25]), or extended with additional terms that, for example, capture chemotaxis, the microstructure of the extracellular matrix or fluid dynamics [3537]. These extensions, such as combining CSM with chemotaxis, can trigger highly intricate and sophisticated dynamics [38]. Nevertheless, understanding the dynamics of the core CSM model is an essential ground step to enable understanding of the full process and in interpreting the meaning and consequences of any subsequent model extension.
Note that the above equation is a simplification which assumes that the cell is completely surrounded by homogeneous contacts (which could be other cells or medium). In the case of an heterogeneous cell environment, the first term, in its most general form, should be written as \(\sum _{\sigma '}J_{\sigma '}p_{\sigma '}\), where σ ^{′} indicates surrounding medium and/or neighbouring cells; \(J_{\sigma ^{\prime }}\) the energy per contact length specifically for the interface between the cell and σ ^{′}; and \(p_{\sigma ^{\prime }}\) the length of that given interface. The sum over all of the cell interfaces’ lengths gives its perimeter \(p=\sum _{\sigma '}p_{\sigma '}\). For the initial analysis we will assume the environment to be homogeneous. There has been ample discussion whether it is more reasonable to use positive or negative Jvalues [22,24]. Part of the issue is that Jvalues may be used to encode several biological processes, making it a challenge to quantify or ascribe to them a biological meaning (see [29] for a detailed discussion). We therefore assume that its sign (as well as the sign of P, see discussion on mapping between J and P below) is undetermined. It is nonsensical, however, to consider negative values for the perimeter and area constraints, and it seems unreasonable to use a negative target area. Moreover, while in many modelling studies no perimeter constraint is being used (corresponding to λ _{ p }=0), it is nonsensical to have no area constraint. We therefore assume that λ _{ p } and A are always nonnegative and λ _{ a } is positive. We initially focus on a 2D cell, and later extend our analysis to 3D tissues.
Note that the formalism, besides discarding any intracellular detail, also describes cell surfaces without explicit “surface elements”, whose movement could be followed over time and would require energy to move closer/away from each other (when not affecting its perimeter or area). While clearly being a coarse simplification, this reduced level of membrane complexity is what allows CSM models to capture complex tissue dynamics involving many cells. (Note that while numerically CSM dynamics might be calculated through displacements of introduced surface elements, they are not relevant for the energy calculation of the configuration, and hence for the dynamics itself.)
where τ is defined as the lengthindependent component of the interfacial tension. The sign of τ is undetermined, while the lengthdependent component is always nonnegative.
The force \(\vec {F}\) represents a vector at a point \(\vec {x}\) on the membrane, which can be decomposed into an interfacial tension force \(\vec {F_{\gamma }}=\gamma \nabla p\) and a pressure driven force \(\vec {F_{\Pi }}=\Pi \nabla a\). An imbalance between these forces should lead to a dynamical change of the cell shape. It is important to note that biological cells are highly dissipative objects, for which viscous forces greatly exceed inertial forces (low Reynolds number, see [39]). Thus, the motion of cell deformations due to the force \(\vec {F}\) acting on the cell membrane is usually assumed to be overdamped, with inertial terms being negligible when compared to dissipative terms, leading to firstorder dynamics, i.e. \(\vec {F}\propto d\vec {x}/dt\) [13,40].
As one of our aims is to clarify how model parameters determine cell and tissue behaviour, we will first characterize the equilibrium states of single cells. To do so, we use the above definitions and derived variables.
Results
Single cell analysis
To understand the role of the biophysical parameters on the cell mechanics we started asking two basic questions. Firstly, what are all the possible equilibrium cell sizes, what is their stability, and how are these defined by the parameters; and secondly, what is the parameter regime for which, at equilibrium, the cell’s interfacial tension (γ) is positive. The answer to the former question tells us for which cell size(s) the forces acting on the cell are in balance, while the answer to the latter question enables us to determine whether the cell shape is stable. In 2D, when γ>0 the cell tends to minimise its perimeter relative to the enclosed area a. Under this condition, the circle is the cellular shape that best minimises the a/p ratio [41].
Positive interfacial tension at the equilibrium cell radii
Under positive interfacial tension, energy reduction leads to surface minimization. When the interfacial tension is negative, the tendency of reducing the surface for a determined area is lost, as there is no energy cost to additional surface – in fact, the tendency is to increase the length of the interface, causing cells which are constrained in area to adopt different noncircular and ruffled shapes. Thus, to predict the qualitative behaviour of the cell regarding shape, it is therefore important to determine the sign of the interfacial tension at the equilibrium radius.
The condition τ>−4π λ _{ p } R _{ a } thus determines an important parameter range in which the circular shape is stable. When the circular shape is unstable, the shape of the cell is hard to predict and will critically depend on the specific implementation of the cell surface and surface dynamics within the diverse CSM model formalisms. Even though the cell shape will be unpredictable, in this parameter regime there is no inherent conflict between fulfilling the area constraint (equivalent to solving Π(a ^{∗})=0) and fulfilling the effective perimeter constraint (equivalent to solving γ(p ^{∗})=0); hence, it is expected that for the shape that will be taken up, a ^{∗}=A and p ^{∗}=P−J/(2λ _{ p }). In contrast, when the circular shape is stable, we can analyse the cell radius equilibria for any CSM implementation, as follows next.
Cell equilibrium radii and their stability
when ε>0. (As discussed above, when ε<0, \(\frac {\partial E}{\partial r}\) is monotonically increasing, which means that there is no such transition, see Figure 2A).
The stability of the solutions is given by the second derivative (Eq. 17). When τ<0, the nontrivial solution ρ _{1} is stable while the trivial solution, r ^{∗}=ρ _{0}=0, is unstable; when τ>0, both the nontrivial solution ρ _{1} and the trivial solution, r ^{∗}=ρ _{0}=0 are stable, while the nontrivial solution r ^{∗}=ρ _{2} is unstable (Figure 2C).
Parameter conditions for the four possible dynamic behaviours of a 2D circular shaped cell
Region  Condition(s) 

I  τ<−4π λ _{ p } R _{ a } 
II  −4π λ _{ p } R _{ a }<τ<0 
III  \(0<\tau <4\pi \lambda _{a}\epsilon ^{\frac {3}{2}}\) 
IV  \(\tau >4\pi \lambda _{a}\epsilon ^{\frac {3}{2}}\) 
3D spherical cell
When an equivalent analysis is performed for the 3D case, we again find that there are only four qualitatively different behaviours possible for a single cell. Moreover, these behaviours can again be fully determined by only two aggregate parameters. Below, we briefly present the outline of the analysis, leaving the details for Appendix A.
Energy function, interfacial tension and pressure
Again, the interfacial tension has a surfaceareadependent and a surfaceareaindependent component, 8π λ _{ s } r ^{2} and \(J8\pi \lambda _{s}{R_{s}^{2}}\), respectively. As in 2D, we denote the surfaceareaindependent component of the interfacial tension by τ, i.e. \(\tau =J8\pi \lambda _{s}{R_{s}^{2}}\) and γ=τ+8π λ _{ s } r ^{2}. (Note that the expression for τ differs from the 2D case.)
Cell equilibrium radii, their interfacial tension and stability
which is very comparable to the condition expressed in Eq. 14b for the 2D case.
where τ is, as previously, the surfaceareaindependent part of the interfacial tension, here given by \(J8\pi \lambda _{s}{R_{s}^{2}}\), which can be either positive or negative; and the other lumped parameters are respectively \(a=\frac {4}{3}\pi \lambda _{v}\), b=8π λ _{ s }, and \(c=\frac {4}{3}\pi \lambda _{v}{R_{v}^{3}}\), with a strictly positive and b and c strictly nonnegative for any reasonable parameter choice.
The first observation is that \(\left.\frac {\partial E}{\partial r}\right _{r=0}=0\) for any parameter choice. This is an important difference with the 2D case. It implies that in 3D, a cell with a radius close to zero slows down its rate of shrinkage or expansion (depending on the stability of the zeroradius equilibrium), while in 2D, the rate of shrinkage or expansion typically remains nonzero for very small radii (given by \(\left.\frac {\partial E}{\partial r}\right _{r=0}=\tau \)). This difference has a clear impact on the cell dynamics, and plays a role in computational implementations of CSM formalisms.
The closedform solutions to \(\frac {\partial E}{\partial r}=0\) for the 3D case are not included here because their full expressions are hardly informative. Note, however, that its numerical evaluation can be easily performed if needed.
As in 2D, the number of equilibria and their stability depend on the sign of τ. If τ is negative, there is only one stable positive equilibrium radius, whereas if τ is positive two scenarios are possible: no equilibrium at all or two nontrivial equilibria, one lower and unstable, and another greater and stable, as summarized in Figure 3B, and derived in Appendix A.
where \(f(\mu)=\sinh \left (\frac {1}{3}\operatorname {arcsinh}\left (\mu \right)\right)\). Figure 3B shows the full 3D bifurcation diagram.
Thus, by calculating μ and ν, which are dependent on the parameters of the CSM, one can immediately establish (without the need for simulations) the expected behaviour of a single spherical cell. The distinguishing feature from the 2D case regards the dynamics (but not the stability) close to zero radius.
Cellular Potts model (CPM)
The CPM, a prominent member of the class of cell surface mechanics formalisms, distinguishes itself by being a latticebased formalism. This property requires special attention on how to relate the CPM numerical descriptors to the variables of the analytical description of CSM. Specifically, it is crucial to explain how the variables p and a of the analytical description relate to the latticebased descriptors of perimeter and area in a CPM simulation.
A cell within the CPM consists of a set of lattice (or grid) points that together define a region in a two dimensional (2D) or three dimensional (3D) lattice. In 2D, the set of lattice points (i,j) that define a specific cell are indicated by a unique identifier σ _{ i,j }. There are as many possible cell shapes as the degrees of freedom defined by the number of permutations of lattice points composing the cell σ, i.e. the bulk area \(n_{\sigma }=\sum _{\textit {ij}}\delta _{\sigma _{\textit {ij}},\sigma }\) (where δ is Kronecker’s delta function). This is an important feature of the CPM compared to other methods that treat cells as points or centreofmass based entities, as it allows an indefinite freedom for cell deformation, being only limited by the spatial resolution. Although many CPM studies explicitly relate the lattice spatial unit Δ x to a physical length (e.g. 1 μ m), relevant model parameters, such as target area or adhesion energies are still often expressed in terms of the lattice arbitrary spatial unit only. To link such latticepointbased parameters to the above analysis, as well as to understand how a change in spatial resolution affects the parameters, we here explicitly take the resolution of the lattice k=1/Δ x into account. We thus distinguish between the bulk area n _{ σ }, expressed in number of lattice points, and the actual area of a CPM cell a _{ σ }, expressed in Δ x ^{2}. These two areas only differ in that n _{ σ } is the number of lattice points composing the CPM cell σ, whereas a _{ σ }=n _{ σ }/k ^{2} is the physical area of a CPM cell, and k the factor describing the resolution of the lattice. The area a in the analytical description maps to the CPM area a _{ σ }. The explanation on how the CPM cell perimeter p _{ σ } relates to the perimeter p of the analytical description is not as simple as for the area, and is therefore deferred to the next section.
The energy function E from which the dynamics are computed is either Eq. 1 for 2D lattices or Eq. 22 for 3D lattices, as explained previously. In Eq. 28, the parameter Y represents a yield – the ability of the membrane to resist a force, in many studies, including this one, set to zero – and T a simulation temperature which determines the extent of fluctuations. The higher the temperature, the more probable it is to observe energetically unfavourable cell shape deformations.
Contact length perimeter and contact surface area in the CPM
Theoretical and measured values of the perimeter scaling factor ξ for small neighbourhood sizes, needed for calculating the effective J value in the CPM
2D  3D  

neigh.  neigh.  ξ , estimated  ξ , numer.  % error  neigh.  ξ , estimated  ξ , numer.  % error 
number  radius  using Eq. 29  determined  radius  using Eq. 33  determined  
1  1  0.67  1  50.00  1  0.79  1  27.32 
2  \(\sqrt {2}\)  1.89  3  59.10  \(\sqrt {2}\)  3.14  5  59.15 
3  2  5.33  5  6.25  \(\sqrt {3}\)  7.07  9  27.32 
4  \(\sqrt {5}\)  7.45  11  47.58  2  12.57  11  12.46 
5  \(2\sqrt {2}\)  15.08  15  0.56  \(\sqrt {5}\)  19.64  23  17.14 
6  3  18.00  18  0  \(\sqrt {6}\)  28.27  39  37.93 
7  \(\sqrt {10}\)  21.08  26  23.33  \(2\sqrt {2}\)  50.27  47  6.50 
8  \(\sqrt {13}\)  31.25  36  15.21  3  63.62  70  10.03 
12  \(2\sqrt {5}\)  59.63  68  14.04  \(\sqrt {13}\)  132.73  134  0.96 
20  \(2\sqrt {10}\)  168.66  173  2.58  \(\sqrt {22}\)  380.13  410  7.86 
100  \(3\sqrt {29}\)  2811.06  2850  1.39  \(\sqrt {118}\)  1.09·10^{4}  11127  1.75 
1000  \(\sqrt {3365}\)  130133  130181  0.037  \(3\sqrt {133}\)  1.13·10^{6}  1127539  0.20 
If one aims for a fair comparison between CPM simulations and either the analytical results shown throughout this paper, or the simulations performed with other CSMbased formalisms, the factor ξ has to be determined for the chosen neighbourhood and the parameters adapted accordingly. Moreover, without applying such a correction no direct quantitative link can be made between biophysical parameters used in CPM simulations and experimental measurements.
The analytical results given in Eq. 29 and Eq. 33 assume a perfectly circular neighbourhood. Small neighbourhoods can deviate substantially from this approximation, as shown in Figure 4B. We therefore measured the residual error in perimeter length after applying the correction factor. Figure 4F and H show the residual error for 2D and 3D, respectively, while Figure 4G and I show the same data magnified. For neighbourhoods with a radius equal or larger than 5, both in 2D and in 3D, the error remains within 6%. For small neighbourhoods, however, the errors can be large, with a 59% mismatch for a 2nd level neighbourhood (see Table 2). Good choices are the 5th level neighbourhood (0.56% error) and the 6th level neighbourhood (0% error). In 3D, the often used 3rd level neighbourhood presents a 27% error, while the 7th level neighbourhood presents a 6% error. In general it is best to use a numerically established correction factor ξ for small neighbourhoods, rather than directly applying Eq. 29 and Eq. 33 (see Table 2).
The choice of the neighbourhood is therefore very important, because it defines both the “localness” of the interfacial tension computation as well as its dependence on the grid geometry, which can be a source of angular bias in the cell’s shape [35]. In general, a radial neighbourhood spanning a wider region will take more lattice points into account and will be closer to the circle (or sphere in 3D), leading to a better approximation of the perimeter. However, the neighbourhood should be sufficiently small compared to the size of the cells themselves, as with increasing neighbourhood radius the information integration becomes less local. Hence, there is a tradeoff between being local, which is desired as it improves spatial resolution, and isotropy of the computation itself. Moreover, larger neighbourhoods are computationally more costly.
In addition, the same neighbourhood level/radius needs to be used in all computations involving perimeter, such as the energy change calculations stemming from adhesion as well as from perimeter conservation.
In order to prevent anisotropic bias stemming from the neighbourhood choice, one may opt to increase its radius while concurrently increasing the spatial resolution (e.g. a 20th level neighbourhood would already be very close to being perfectly isotropic for most CPM studies). This way, a good compromise between the localness and isotropy of the interfacial tension computation can be achieved.
CPM spatial resolution
Onetoone mapping between adhesion energies J and target perimeter P
where \(\overline {\Delta J_{w}}\,=\,\left (\frac {1}{2}\sum _{\sigma '}\Delta J_{\tau (\sigma),\tau (\sigma ')}p_{\sigma,\sigma '}+\Delta J_{\tau (\sigma),M}p_{\sigma,M}\right)/p_{\sigma }\) is the weighted mean adhesiondriven interfacial tension. This simple expression shows the balance that has to exist to keep the dynamics the same if the adhesion energies are to be changed, negative or positivewise. The equation shows that in general this compensation depends on the extent of the contact length of cell σ with other neighbouring cells, i.e. \(p_{\sigma,\sigma ^{\prime }}\) (this dependence is implicit in \(\overline {\Delta J_{w}}\)), and it is therefore not possible to find a unique mapping between a change in one specific adhesion energy \(\phantom {\dot {i}\!}J_{\tau (\sigma),\tau (\sigma ^{\prime })}\) and the perimeter constraints P _{ σ }. The reason is that if one changes a specific adhesion energy by adding \(\Delta J_{\tau (\sigma),\tau (\sigma ^{\prime })}\phantom {\dot {i}\!}\) (\(J_{\tau (\sigma),\tau (\sigma ^{\prime })}\to J_{\tau (\sigma),\tau (\sigma ')}+\Delta J_{\tau (\sigma),\tau (\sigma ')}\phantom {\dot {i}\!}\)), the equilibrium contact lengths may change for that cell. This changes the weights (of the weighted average \(\overline {\Delta J_{w}}\)), i.e. \(p_{\sigma,\sigma ^{\prime }}\phantom {\dot {i}\!}\), and probably also the cell shape, as the cell will now have a new equilibrium partitioning of its contact lengths with the neighbouring cells.
for the specific case of a concurrent change of all Jvalues with Δ J _{ c e l l,c e l l } (but cellmedium adhesion energies with half of this value).
When applying this method, it is important to realize that there is no such scaling possible for the parameter λ _{ p }. One consequence is that one cannot rescale a model which does not take a perimeter constraint into account (λ _{ p }=0) into one which does. That is, the rescaling only works once the perimeter constraint is already taken into account. Likewise, to be explicit, rescaling the model such that P _{ σ }=0 does not mean that perimeter constraint has been taken out of the equation. It rather implies a perimeter constraint with a rest length of zero.
Single cell simulations
In region I, a cell of initial radius 50Δ x grows rapidly in size (Figure 7C1). In this region, the cell is expected to have a negative interfacial tension at the circular shape. Therefore, the circular shape is unstable and the area and perimeter constraints can be both fulfilled independently, while the cell takes up a complex shape. In such a case, the stable equilibrium area a ^{∗} is directly defined by the cell’s target area (a ^{∗}=A), while the equilibrium perimeter p ^{∗} is determined by an interfacial tension of γ=0 (\(p^{*}=\frac {\tau }{2\lambda _{p}}=P\frac {J}{2\lambda _{p}}\)). Figure 7C2 shows the time dynamics of the perimeter and radius of four cells with different initial sizes, the blue line indicating the cell (Figure 7C1). The dashed lines indicate the expected perimeter and radius at the unstable equilibrium with circular cell shape, while the solid lines indicate the expected perimeter and radius when both constraints are fulfilled independently in a complex cell shape. After 500 Monte Carlo Step (MCS), the predicted p ^{∗}=2000 Δ x and a ^{∗}=31416Δ x ^{2} (or r ^{∗}=100Δ x) perfectly match the simulations, while the cell shapes are clearly deviating from being circular.
In region II, the simulated cell shows a more isotropic shape that approaches a circular shape (Figure 7D1). The simulated cells starting at different radii all converge to the same equilibrium perimeter and area (Figure 7D2). While the analytical equilibrium area a ^{∗} is approached very closely by the numerical simulations, the analytical equilibrium perimeter p ^{∗}=619Δ x deviates considerably from what is observed in the simulations (which is around 800Δ x, see Figure 7D2). This deviation can be explained by the membrane fluctuations that are due to the stochastic nature of the CPM update scheme. Close to region I the interfacial tension at equilibrium is only marginally larger than zero. Consequently, stochasticity can easily cause the formation of excess perimeter, causing its length to deviate from the analytical, deterministic prediction. This effect strongly depends on the simulation temperature parameter T. When this parameter is lowered, the difference between the predicted and numerical observed perimeter at equilibrium is reduced. Note, however, that low simulation temperatures can also introduce undesirable grid effects [35].
In region III, the cell dynamics exhibits bistability. In Figure 7E2, the dashed lines indicate the analytical unstable equilibrium (p ^{∗}=194Δ x, a ^{∗}=2995Δ x ^{2}) and the solid lines the stable nontrivial equilibrium (p ^{∗}=482Δ x, a ^{∗}=18468Δ x ^{2}), with the other stable equilibrium given by p ^{∗}=0Δ x, a ^{∗}=0Δ x ^{2}. Figure 7E1 shows the simulations of two different cells, one starting below the unstable steadystate (initial radius r=30Δ x), the other starting above the unstable equilibrium (initial radius r=36Δ x). While the larger cell increases in size, the smaller cell eventually disappears. Note that in this region there is hardly any discrepancy between the analytical and observed values, due to a much higher interfacial tension at equilibrium.
In region IV, the simulated cells exhibit a very smooth membrane, due to an even higher interfacial tension. The interfacial tension dominates over the cell’s pressure at any cell size, triggering a continuous size decrease until the cell vanishes (Figure 7F1). Both the perimeter and area concomitantly decrease until the trivial equilibrium p ^{∗}=0Δ x, a ^{∗}=0Δ x ^{2} is reached (Figure 7F2).
Secondly, in the analysis we find a sharp transition from region I (with a negative interfacial tension for an equilibrium circular cell, leading to a noncircular cell shape) and region II (with a positive interfacial tension, in which a circular shape is preferred). In the computer simulations we find that this transition is more blurred: when the interfacial tension becomes small, the forces restoring a round shape from thermallydriven shape deformations become small as well. Hence, close to the region I/II boundary the cell shape increasingly deviates from a round shape, resulting in a smooth transition. Still, overall, those simulations present a close match between the expected perimeter and area from the mathematical analysis. It illustrates that the CPM, albeit a latticebased, discrete and stochastic model formalism, still presents equilibrium dynamics as expected from the basics of cell surface mechanics.
Epithelial cell packing
In the previous section we have pointed out deviations from the analytical results when cell shapes differ too much from a perfect circle (in 2D) or perfect sphere (in 3D). Clearly, within a packed tissue context such an assumption of circularity is invalid. We therefore questioned whether the analysis could be extended to noncircular shapes.
For 2D we focus on epithelial tissue. Simple epithelium is normally a cell monolayer tissue made of cells strongly adhering to each other in a plane of cell junctions. The cellular organization of the epithelium derives from the mechanical properties, particularly at the level of cell junctions. Several modelling studies have looked whether the morphological properties of epithelial tissue can be derived from simple rules describing its mechanical properties [24,25]. Here, we try to approach the characterization of epithelial cell packing using the CSM description. We explore the assumption that the behaviour of an epithelial sheet is no more than the collective behaviour of its individual cells. Due to high cell density within epithelial tissues, cells are typically densely packed and exhibit hexagonal shapes. We therefore seek how a cell behaves when it is constrained to a hexagonal shape.
For example, a hexagonal cell can be described by k _{ p }=6 and \(k_{a}=\frac {3\sqrt {3}}{2}\), which can be derived when taking for l the polygon’s side length (but any choice for l will lead to the same result). Likewise, a circular cell can be described by k _{ p }=2π and k _{ a }=π, typically taking for l the radius r.
While the parameters k _{ p } and k _{ a } are constants that relate the perimeter and area with l, L _{ p } and L _{ a } are the reciprocal target lengths of the target perimeter P and area A.
We next determined whether this analytical result can indeed correctly describe cell behaviour of hexagonal cells in a multicellular context. We therefore compared the analytical predictions to published simulation results of a vertex model [25], as well as to a series of tissue simulations using the CPM formalism.
Vertex model cell packing simulations
To show that the results regarding the epithelial cell packing are a general consequence of the CSM model and not a quality particular to the CPM or any other numerical implementation, we took published data about the modelling of epithelial cell packing that uses a different mathematical formalism and compared it to our analytical predictions. For this, we used the work reported in Farhadifar et al. [25], where the authors used the Vertex model to investigate the mechanical properties of epithelial cell packing. Their approach describes the CSM tensions using the same energy description as presented here, compare Eq. 1 with Box 1 in [25].
To compare with the results from [25], we introduce the correspondence between their notation and ours: J=Λ, λ _{ p }=Γ, λ _{ a }=K, A=A ^{0} and P=0. The fact that they assumed P=0 in their approach means that our description of the CSM holds more generally, because with such an assumption the influence of the cortical tension on the interfacial tension is always positive. The authors also introduced a normalised tension \(\overline {\Lambda }\) and normalised contractility \(\overline {\Gamma }\) (see Eq. 84), which were used to depict their biophysical parameter space.
By representing the bifurcation diagram of Figure 10A in the equivalent (\(\overline {\Lambda },\overline {\Gamma }\))space, we show that the analytical results are perfectly matched by the Vertex model simulations of [25] (Figure 10B): Eq. 53a, which separates the positive from the negative interfacial tension regime, matches the boundary obtained by simulation (green dots in Figure 10B, extracted from [25]). Eq. 53b, which separates the regime with a single positive stable equilibrium from the bistable regime in our analysis, is not present in their simulation bifurcation diagram. Nevertheless, when we follow our analysis, the socalled Case II of their parameter space should be located within the bistable regime. Indeed, a snapshot of their simulation reveals a huge variation in cell size, strongly suggesting bistability (see Figure 2H in [25]). Finally, Eq. 53c, which separates region III and IV, lies right on top of their numericallyderived line that separates the “Hexagonal network” and “cell vanishing” regimes. Please note that a comparable analytical effort has been published to understand the boundaries between the different observed behaviours [30]. That study, however, did not note the distinction between region II and III, and they derived an approximation of the boundary between region III and IV, while the solution presented here is exact. Consequently, unlike [30], our analytically derived boundary between region III and IV matches perfectly the numerical simulations over their full domain.
These results show that the tissue properties which were defined in [25] as “Soft Network”, “Hexagonal Network” and “Cell Vanishing” perfectly correspond to region I, region II/III, and region IV, respectively, as presented in this paper. It highlights that such tissue properties are basically a manifestation of the mechanical properties of individual cells and not due to collective cell behaviour per se (except that collective cell behaviour is needed to give rise to the hexagonal cell shapes). In fact, knowledge on the expected cell shape combined with single cell analysis turned out to be sufficient to account for the whole range of observed tissue behaviour.
CPM cell packing simulations
In region I the cells in the cluster show a fringed border (Figure 11C1). In this region, the cells are expected to have a negative interfacial tension when taking up a hexagonal shape. Consequently, the hexagonal shape is unstable, and the tissue is formed by dynamically changing, nonhexagonal cells, with tissue topology alterations by neighbour changes, during which some cells loose their contact, while others come into contact with each other (socalled T1 transitions, see [19]). In this regime, both perimeter and area constraint are fulfilled independently, the equilibrium area a ^{∗} defined by the cell’s target area and the equilibrium perimeter p ^{∗} determined by a zero cell interfacial tension, as shown in Figure 11B, C2. The predicted p ^{∗}=2,000Δ x and a ^{∗}=31,416Δ x ^{2} (or l ^{∗}=110 Δ x) are perfectly matched by the simulations.
In region II, the simulated cell cluster takes up a honeycomb shape, with each individual cell approaching the expected hexagonal shape (Figure 11D1). After initiation of the simulation, all cells rapidly converge to the same equilibrium perimeter and area (Figure 11D2). As was also observed for the single cell simulations, the observed equilibrium perimeter is larger than the expected equilibrium perimeter p ^{∗}=648Δ x, due to the stochasticity of the CPM update scheme. For this particular parameter choice the difference is around 20%, but it increases when the equilibrium is followed towards region I and decreases when followed towards region IV, with lower simulation temperatures yielding a much closer match.
In region III, the cells exhibit bistability (Figure 11E1, F1). Figure 11E1 illustrates that when a simulation is initiated with cells that are larger than the unstable equilibrium (p ^{∗}=221Δ x, a ^{∗}=3,509Δ x ^{2}), the cells evolve towards the nontrivial stable equilibrium (p ^{∗}=490Δ x, a ^{∗}=17,336Δ x ^{2}, Figure 11E2). In contrast, in a simulation initiated with cells below the unstable equilibrium (Figure 11F1), most cells vanish (p ^{∗}=0Δ x, a ^{∗}=0Δ x ^{2}), but some, due to stochastic fluctuations, pass the unstable equilibrium to end up in the nontrivial stable one (Figure 11F2). Note that due to the increased interfacial tension compared to the previous simulations, fluctuations that deviate the cell shape from perfect hexagonality are much less pronounced, and hence the simulations match the analytical results more precisely.
In region IV, the cell edges are even smoother due to the high interfacial tension (Figure 11G1). However, all cells vanish rapidly, reaching the trivial equilibrium (p ^{∗}=0Δ x, a ^{∗}=0Δ x ^{2}, Figure 11G2).
Parameters used in the 2D analysis and their meaning
Parameter  Meaning  General  Circle  Hexagon 

E  Energy function or Hamiltonian  J p+λ _{ p }(p−P)^{2}+λ _{ a }(a−A)^{2}  
J  Adhesion energy (per contact length)  J  J  J 
p  Cell perimeter  k _{ p } l  2π r  6l 
a  Cell area  k _{ a } l ^{2}  π r ^{2}  \(\frac {3\sqrt {3}}{2}l^{2}\) 
P  Membrane rest length  k _{ p } L _{ p }  2π R _{ p }  6L _{ p } 
A  Target cell area  \(k_{a}{L_{a}^{2}}\)  \(\pi {R_{a}^{2}}\)  \(\frac {3\sqrt {3}}{2}{L_{a}^{2}}\) 
λ _{ p }  Perimeter constraint  λ _{ p }  λ _{ p }  λ _{ p } 
λ _{ a }  Area constraint  λ _{ a }  λ _{ a }  λ _{ a } 
l  Basic length scale  l  r (radius)  l 
k _{ p }  Perimeter scaling factor  \(\frac {p}{l}\)  2π  6 
k _{ a }  Area scaling factor  \(\frac {a}{l^{2}}\)  π  \(\frac {3\sqrt {3}}{2}\) 
L _{ p }  Membrane rest length, using basic length scale  \(\frac {P}{k_{p}}\)  \(R_{p}=\frac {P}{2\pi }\)  \(\frac {P}{6}\) 
L _{ a }  Target cell area,using basic length scale  \(\sqrt {\frac {A}{k_{a}}}\)  \(R_{a}=\sqrt {\frac {A}{\pi }}\)  \(\sqrt {\frac {2A}{3\sqrt {3}}}\) 
E  Energy function or Hamiltonian, using basic length scale  \({Jk}_{p}l+\lambda _{p}(k_{p}lk_{p}L_{p})^{2} +\lambda _{a}(k_{a}l^{2}k_{a}{L_{a}^{2}})^{2}\)  \(2\pi rJ+\lambda _{p}(2\pi r2\pi R_{p})^{2} +\lambda _{a}(\pi r^{2}\pi {R_{a}^{2}})^{2}\)  \(6lJ+\lambda _{p}(6l6L_{p})^{2} +\lambda _{a}(\frac {3\sqrt {3}}{2}l^{2}\frac {3\sqrt {3}}{2}{L_{a}^{2}})^{2}\) 
\(\frac {\partial E}{\partial l}\)  Energy variation per length change  \(k_{p}(\gamma 2\frac {k_{a}}{k_{p}}l\Pi)\)  2π(γ−r Π)  \(6(\gamma \frac {\sqrt {3}}{2}l\Pi)\) 
γ  Enterfacial tension  J+2k _{ p } λ _{ p }(l−L _{ p })  J+4π λ _{ p }(r−R _{ p })  J+12λ _{ p }(l−L _{ p }) 
Π  Pressure  \(2k_{a}\lambda _{a}(l^{2}{L_{a}^{2}})\)  \(2\pi \lambda _{a}(r^{2}{R_{a}^{2}})\)  \(3\sqrt {3}\lambda _{a}(l^{2}{L_{a}^{2}})\) 
τ  Lengthindependent component of interfacial tension  J−2k _{ p } λ _{ p } L _{ p }  J−4π λ _{ p } R _{ p }  J−12λ _{ p } L _{ p } 
ε  l ^{2} at which \(\frac {\partial ^{2}E}{\partial l^{2}}=0\)  \(\frac {{L_{a}^{2}}}{3}\frac {{k_{p}^{2}}\lambda _{p}}{6{k_{a}^{2}}\lambda _{a}}\)  \(\frac {{R_{a}^{2}}}{3}\frac {2\lambda _{p}}{3\lambda _{a}}\)  \(\frac {{L_{a}^{2}}}{3}\frac {8\lambda _{p}}{9\lambda _{a}}\) 
α  \(\left.\frac {\partial E}{\partial l}\right _{l=\sqrt {\epsilon },\tau =0}\)  \(8{k_{a}^{2}}\lambda _{a}\left (\frac {{L_{a}^{2}}}{3}\frac {{k_{p}^{2}}\lambda _{p}}{6{k_{a}^{2}}\lambda _{a}}\right)^{\frac {3}{2}}\)  \(8\pi ^{2}\lambda _{a}\left (\frac {{R_{a}^{2}}}{3}\frac {2\lambda _{p}}{3\lambda _{a}}\right)^{\frac {3}{2}}\)  \(54\lambda _{a}\left (\frac {{L_{a}^{2}}}{3}\frac {8\lambda _{p}}{9\lambda _{a}}\right)^{\frac {3}{2}}\) 
β  Aggregate parameter  \(4\frac {{k_{a}^{2}}\lambda _{a}}{k_{p}}\)  2π λ _{ a }  \(\frac {9\lambda _{a}}{2}\) 
ζ  Aggregate parameter  \(\frac {64{k_{a}^{4}}{\lambda _{a}^{2}}}{{k_{p}^{2}}}\left (\frac {{k_{p}^{2}}\lambda _{p}}{6{k_{a}^{2}}\lambda _{a}}\frac {{L_{a}^{2}}}{3}\right)^{3{\vphantom {\frac {1}{2}}}}\)  \(16\pi ^{2}{\lambda _{a}^{2}}\left (\frac {2\lambda _{p}}{3\lambda _{a}}\frac {{R_{a}^{2}}}{3}\right)^{3}\)  \(81{\lambda _{a}^{2}}\left (\frac {8\lambda _{p}}{9\lambda _{a}}\frac {{L_{a}^{2}}}{3}\right)^{3}\) 
Bifurcation 1 (γ(l ^{∗})=0)  Transition from negative to positive interfacial tension at equilibrium  τ=−2k _{ p } λ _{ p } L _{ a }  τ=−4π λ _{ p } R _{ a }  τ=−12λ _{ p } L _{ a } 
Bifurcation 2 (pseudotranscritical)  Transition of l ^{∗}=0 from unstable to stable  τ=0  τ=0  τ=0 
Bifurcation 3 (fold)  Transition from 2 to 0 nontrivial equilibria  \(\tau =\frac {8{k_{a}^{2}}\lambda _{a}}{k_{p}}\epsilon ^{\frac {3}{2}}\)  \(\tau =4\pi \lambda _{a}\epsilon ^{\frac {3}{2}}\)  \(\tau =9\lambda _{a}\epsilon ^{\frac {3}{2}}\) 
\(\overline {\Lambda }\)  Normalised tension,as used in [25]  \(\frac {J}{k_{a}^{\frac {3}{2}}\lambda _{a}{L_{a}^{3}}}\)  \(\frac {J}{\pi ^{\frac {3}{2}}\lambda _{a}{R_{a}^{3}}}\)  \(\frac {2\sqrt {2}J}{9\sqrt [4]{3}\lambda _{a}{L_{a}^{3}}}\) 
\(\overline {\Gamma }\)  Normalised contractility,as used in [25]  \(\frac {\lambda _{p}}{k_{a}\lambda _{a}{L_{a}^{2}}}\)  \(\frac {\lambda _{p}}{\pi \lambda _{a}{R_{a}^{2}}}\)  \(\frac {2\lambda _{p}}{3\sqrt {3}\lambda _{a}{L_{a}^{2}}}\) 
Bifurcation 1 (γ(l ^{∗})=0)  Transition from negative to positive interfacial tension at equilibrium  \(\overline {\Gamma }=\frac {\sqrt {k_{a}}}{2k_{p}}\overline {\Lambda }\)  \(\overline {\Gamma }=\frac {1}{4\sqrt {\pi }}\overline {\Lambda }\)  \(\overline {\Gamma }=\frac {1}{4\sqrt {2}\sqrt [4]{3}}\overline {\Lambda }\) 
Bifurcation 2 (pseudotranscritical)  Transition of l ^{∗}=0 from unstable to stable  \(\overline {\Lambda }=0\)  \(\overline {\Lambda }=0\)  \(\overline {\Lambda }=0\) 
Bifurcation 3 (fold)  Transition from 2 to 0nontrivial equilibria  \(\overline {\Gamma }=\frac {4k_{a}3\left (k_{a}k_{p}\overline {\Lambda }\right)^{\frac {2}{3}}}{2{k_{p}^{2}}}\)  \(\overline {\Gamma }=\frac {4\pi 3\left (2\pi ^{2}\overline {\Lambda }\right)^{\frac {2}{3}}}{8\pi ^{2}}\)  \(\overline {\Gamma }=\frac {23\sqrt [6]{3}\left (\overline {\Lambda }\right)^{\frac {2}{3}}}{8\sqrt {3}}\) 
3D cell packing
where \(\nu =\frac {\left (J2k_{s}\lambda _{s}{L_{s}^{2}}\right)9{k_{v}^{2}}\lambda _{v}}{4{k_{s}^{3}}{\lambda _{s}^{2}}}=\frac {\tau }{12a\phi ^{2}}\), \(\mu =\frac {27{k_{v}^{3}}\lambda _{v}^{\frac {3}{2}}{L_{v}^{3}}}{8{k_{s}^{3}}\lambda _{s}^{\frac {3}{2}}}=\frac {\psi }{\phi ^{\frac {3}{2}}}\), and \(f(\mu)=\sinh \left (\frac {1}{3}\operatorname {arcsinh}\left (\mu \right)\right)\).
Again, when the analytical results are compared to CPM simulations of a cell tissue, we find a very close correspondence between them. Figure 9B shows 3D CPM simulations of a small tissue for each of the four possible regions. It shows the individual cells within a cluster that is not in direct contact with the medium. As for the 2D case, in region I the cells take up an irregular shape, fulfilling both the surface area and the volume requirement. In region II, the cell shape is very regular, indeed very closely corresponding to a rhombic dodecahedron. In region III, the initial size of the individual cells determines whether the cells disappear or not, while in region IV all cells eventually disappear. The boundaries between those regions of different qualitative behaviour, as observed in the 3D CPM simulations, again closely match the analytical predictions (data not shown).
Discussion
Our analytical results, in combination with computer simulations, have elucidated how the biophysical parameters of CSM models can be translated into qualitative cell dynamics. Some observations and underlying assumptions made during this process merit further discussion. To start with, in region I of the cellular dynamics, the observed behaviour is driven by a negative interfacial tension. However, it is questionable whether such a parameter regime is biologically and physically reasonable; the main issue being that negative interfacial tension implies the intermingling of neighbouring cells. This phenomenon is not observed in low resolution cell shape CSM models, which often use, for example, a single line element for the whole interface between two cells (e.g. [25]), effectively inhibiting cells from intermingling. However, in the CPM the cell’s interface is resolved at much higher detail, clearly exposing the intermingling between cells (see Figure 11C1). From a biological perspective, such dynamics might be considered pathological and unrealistic, and one should be weary of using such a dynamical regime to capture normal biological processes.
On the other range of the spectrum, it is also worth discussing whether it is reasonable to assume that cells can completely disappear (i.e. shrink to zero) due to their cell surface mechanics. One could argue that the elastic description of area and perimeter conservation becomes unreasonable when deviations become large, and that while water can freely enter and leave the cell, ultimately the dry cellular mass should prevent full disappearance. Nevertheless, experimental biological studies have been showing pressuredriven apoptosis [48,49], which has been used as a possible interpretation for vanishing cells in the CPM [50]. Another explanation has been proposed by Marinari et al. [16], in a combined experimental and modelling study of the Drosophila notum. They showed that in a 2D CSM description of the notum an increased pressure resulted in the local disappearance of cells, which was then attributed to the process of delamination, when cells leave the epithelial plane and enter the (nondescribed) underlying tissue.
Given that the most reasonable dynamics are found for (a large range of) wellbalanced parameter choices, it will be important for most theoretical studies which are based upon CSM concepts — including the CPM — to apply the insights derived from our study: (i) to choose parameters sensibly, (ii) to ensure model behaviour matches experimental data, and (iii) to prevent potential modelling artefacts.
Although the analysis presented here is relevant for the large class of CSM models that use the same underlying building blocks, there is obviously a wide range of different model descriptions in use to capture cell dynamics. These alternative models cannot always be easily encapsulated within the same theoretical framework. For instance, other common strategies to describe the biophysics of cells is using the socalled subcellular elements method [51,52] or by modelling of cellular material as a continuous medium taking stress and strain tensors explicitly into account [53]. These modelling strategies grant a finer description of forces, not just at the membrane, as with CSM formalisms, but also within the cell. Although this level of detail is important when studying the cell rheology, quite often the simplifications defining CSM models make them better suited for morphodynamical simulations, where part of the individual cell detail may be exchanged for the feasibility of manycell simulations.
Finally, we argue that the proposal suggested by Ouchi et al. [22] — to use negative instead of positive Jvalues (adhesion energy constant) to solve the apparent mismatch (regarding the relationship between adhesion and cell motility) between models and experiments — fails to reconcile experiments and simulations. Using both theory and simulations, we have shown here that this limitation is not particular to the CPM, but rather present in any modelling framework that describes the strength of adhesion through an effective lower tension.
Conclusions
For the perplexed: how to determine the expected dynamics in any CSM model
Step 1.  Write the Hamiltonian or energy function in the standard form as used in this paper (i.e. E=J p+λ _{ p }(p−P)^{2}+λ _{ a }(a−A)^{2} for 2D; E=J s+λ _{ s }(s−S)^{2}+λ _{ v }(v−V)^{2} for 3D). 

Step 2.  When using CPM, first the effective values of J, and P and λ _{ p } in 2D, or S and λ _{ s } in 3D, have to be determined, or inversely the correct values of \(J_{_{\mathit {CPM}}}\) etc. have to be assigned. This is done by calculating the perimeter scaling factor ξ, which depends on the neighbourhood radius, in its turn depending on the neighbourhood level. The value of ξ can be obtained either by using the theoretical formulae Eq. 29, Eq. 33, for 2D and 3D respectively, in the case that the neighbourhood is sufficiently large, or by using the numerical estimates given in Table 2, for small neighbourhoods that can present significant deviations from those theoretical values. Transformations between the effective and CPM parameter values are then given by Eq. 31, Eq. 32 for 2D, and Eq. 34 for 3D. Note that because of the way interfacial tension is implemented in the CPM, \(J=\xi J_{_{\mathit {CPM}}}\) for cellmedium interfaces and \(\frac {\xi }{2}J_{_{\mathit {CPM}}}\) for cellcell interfaces. 
Step 3.  Calculate k _{ p }, k _{ a } (for 2D) cq. k _{ s }, k _{ v } (for 3D) for the cell shape of interest. The values for typical shapes in 2D and 3D, namely circle and hexagon, cq. sphere and rhombic dodecahedron, are given in Table 3 and Table 6. 
Step 4.  Calculate τ and ε (for 2D), ν and μ (for 3D when λ _{ s }≠0), or ν ^{′} and μ ^{′} (for 3D when λ _{ s }=0), using the formulae given in Table 3 (for 2D) and Table 6 (for 3D). 
Step 5.  The expected behaviour of a single cell or tissue is given by Figure 3; the bifurcation lines which form the transitions between the different zones are given in Table 3 (for 2D) and Table 6 (for 3D). 
For the perplexed: how to correctly rescale a CPM model
Rescale spatial resolution  When the spatial resolution of a CPM model is kfold increased, the required changes in the standard set of kinetic parameters are given in Eq. 38 and Eq. 39. 

Resizing neighbourhood  When the neighbourhood used in a simulation is changed, the value of \(J_{_{\mathit {CPM}}}\phantom {\dot {i}\!}\), and \(P_{_{\mathit {CPM}}}\phantom {\dot {i}\!}\) and \(\lambda _{p\,{}_{\mathit {CPM}}}\phantom {\dot {i}\!}\) in 2D, or \(S_{_{\mathit {CPM}}}\phantom {\dot {i}\!}\) and \(\lambda _{s\,{}_{\mathit {CPM}}}\phantom {\dot {i}\!}\) in 3D, have to be modified, such that the effective values remain the same (Eq. 31, Eq. 32, Eq. 34). This can be achieved by setting \(J_{_{\mathit {CPM}}}'=\frac {\xi _{\textit {old}}}{\xi _{\textit {new}}}J_{_{\mathit {CPM}}}\phantom {\dot {i}\!}\), where ξ _{ old } and ξ _{ new } are the perimeter scaling factor before and after resizing the neighbourhood, respectively. Likewise, \(P_{_{\mathit {CPM}}}'=\frac {\xi _{\textit {new}}}{\xi _{\textit {old}}}P_{_{\mathit {CPM}}}\phantom {\dot {i}\!}\), \(\lambda _{p\,{}_{\mathit {CPM}}}'=\frac {\xi _{\textit {old}}^{2}}{\xi _{\textit {new}}^{2}}\lambda _{p\,{}_{\mathit {CPM}}}\phantom {\dot {i}\!}\), \(S_{_{\mathit {CPM}}}'=\frac {\xi _{\textit {new}}}{\xi _{\textit {old}}}S_{_{\mathit {CPM}}}\phantom {\dot {i}\!}\), and \(\lambda _{s\,{}_{\mathit {CPM}}}'=\frac {\xi _{\textit {old}}^{2}}{\xi _{\textit {new}}^{2}}\lambda _{s\,{}_{\mathit {CPM}}}\phantom {\dot {i}\!}\). Details on calculating ξ are in Step 2 of Table 4. 
Concurrent rescaling of J and P  It is possible to concurrently change all J values (for example from all being positive to all being negative) in a CPM, or in fact any CSM simulation, without causing any change in the dynamics of the model (Figure 6), by means of wellchosen shifts in the membrane rest lengths of all cells (P _{ σ } and S _{ σ } in 2D and 3D, respectively). The required shifts in the rest lengths are given in Eq. 47 (for 2D) and Eq. 49 (for 3D). In contrast, it is not possible to change only a subset of the J values without causing changes in dynamics. Specifically, in such a case, it is still possible to keep for a specific configuration the weighted mean adhesiondriven interfacial tension constant (using Eq. 46 for 2D and Eq. 48 for 3D), but those weighted means are expected to change over time (for example due to cell sorting), generating an imbalance and hence a change of dynamics on the long run. 
Appendix
A 3D Spherical Cell
Positive interfacial tension at the equilibrium radii
which is very comparable to the condition expressed in Eq. 14b of the 2D case.
Equilibrium radii and their stability
where \(\tau =J8\pi \lambda _{s}{R_{s}^{2}}\), which can be both positive or negative; \(a=\frac {4}{3}\pi \lambda _{v}\), strictly positive; and b=8π λ _{ s } and \(c=\frac {4}{3}\pi \lambda _{v}{R_{v}^{3}}\), both strictly nonnegative.
We can use Descartes’ rule of signs to determine the highest number of equilibria larger than zero (remember that r ^{∗}<0 is nonsensical). When τ<0, there is only one sign change, limiting the number of equilibria larger than zero to at most one, while when τ>0, two sign changes limit it to at most two.
The slope of the function \(\frac {\partial E}{\partial r}\) at the equilibrium r ^{∗}=0 is 8π τ, implying that this equilibrium is always stable for τ>0 (positive slope), and always unstable for τ<0 (negative slope). (Recall that the force on the cell membrane is positive in the direction of the energy gradient steepest descent.)
Moreover, because \(\left.\frac {\partial E}{\partial r}\right _{r=0}=0\) and \(\left.\frac {\partial E}{\partial r}\right _{r\rightarrow \infty }=\infty \), a negative slope implies an odd number of equilibria larger than zero, while a positive slope implies an even number. Combined with the insights derived from Descartes’ rule of signs, we can therefore conclude that when τ<0, there is at least and at most one equilibrium for r>0, i.e. there is always one single equilibrium. Also, this equilibrium has to be stable, since the slope at the equilibrium has to be positive. When τ>0, there can be either no equilibria, or two positive equilibria. When there are two equilibria, the lower one will be the unstable and the higher one the stable.
Parameters used in the 3D analysis and their meaning
Parameter  Meaning  General  Sphere  Rhombic dodecahedron 

E  Energy function or Hamiltonian  J s+λ _{ s }(s−S)^{2}+λ _{ v }(v−V)^{2}  
J  Edhesion energy (per contact length)  J  J  J 
s  Cell surface  k _{ s } l ^{2}  4π r ^{2}  \(8\sqrt {2}l^{2}\) 
v  Cell volume  k _{ v } l ^{3}  \(\frac {4}{3}\pi r^{3}\)  \(\frac {16}{3\sqrt {3}}l^{3}\) 
S  Rest surface area  \(k_{s}{L_{s}^{2}}\)  \(4\pi {R_{s}^{2}}\)  \(8\sqrt {2}{L_{s}^{2}}\) 
V  Ttarget cell volume  \(k_{v}{L_{v}^{3}}\)  \(\frac {4}{3}\pi {R_{v}^{3}}\)  \(\frac {16}{3\sqrt {3}}{L_{v}^{3}}\) 
λ _{ s }  Surface constraint  λ _{ s }  λ _{ s }  λ _{ s } 
λ _{ v }  Volume constraint  λ _{ v }  λ _{ v }  λ _{ v } 
l  Basic length scale  l  r (radius)  l 
k _{ s }  Surface scaling factor  \(\frac {s}{l^{2}}\)  4π  \(8\sqrt {2}\) 
k _{ v }  Volume scaling factor  \(\frac {v}{l^{3}}\)  \(\frac {4}{3}\pi \)  \(\frac {16}{3\sqrt {3}}\) 
L _{ s }  Rest surface area, using basic length scale  \(\sqrt {\frac {S}{k_{s}}}\)  \(R_{s}=\frac {1}{2}\sqrt {\frac {S}{\pi }}\)  \(\sqrt {\frac {S}{8\sqrt {2}}}\) 
L _{ v }  Target cell volume,using basic length scale  \(\sqrt [3]{\frac {V}{k_{v}}}\)  \(R_{v}=\sqrt [3]{\frac {3V}{4\pi }}\)  \(\sqrt [3]{\frac {3\sqrt {3}V}{16}}\) 
E  Energy function or Hamiltonian, using basic length scale  \({Jk}_{s}l^{2}+\lambda _{s}(k_{s}l^{2}k_{s}{L_{s}^{2}})^{2} +\lambda _{v}(k_{v}l^{3}k_{v}{L_{v}^{3}})^{2}\)  \(4\pi Jr^{2}+\lambda _{s}(4\pi r^{2}4\pi {R_{s}^{2}})^{2} +\lambda _{v}(\frac {4}{3}\pi r^{3}\frac {4}{3}\pi {R_{v}^{3}})^{2}\)  \(8\sqrt {2}Jl^{2}\) \(+\lambda _{s}(8\sqrt {2}l^{2}8\sqrt {2}{L_{s}^{2}})^{2} +\lambda _{v}(\frac {16}{3\sqrt {3}}l^{3}\frac {16}{3\sqrt {3}}{L_{v}^{3}})^{2}\) 
\(\frac {\partial E}{\partial l}\)  Energy variation per length change  \(2k_{s}l\left (\gamma \frac {3k_{v}}{2k_{s}}l\Pi \right)\)  4π r(2γ−r Π)  \(16\sqrt {2}l\left (\gamma \frac {l}{\sqrt {6}}\Pi \right)\) 
γ  Interfacial tension  \(J+2k_{s}\lambda _{s}(l^{2}{L_{s}^{2}})\)  \(J+8\pi \lambda _{s}(r^{2}{R_{s}^{2}})\)  \(J+16\sqrt {2}\lambda _{s}(l^{2}{L_{s}^{2}})\) 
Π  Pressure  \(2k_{v}\lambda _{v}(l^{3}{L_{v}^{3}})\)  \(\frac {8}{3}\pi \lambda _{v}(r^{3}{R_{v}^{3}})\)  \(\frac {32}{3\sqrt {3}}\lambda _{v}(l^{3}{L_{v}^{3}})\) 
\(\frac {\partial E}{\partial l}\)  Energy variation per length change, full expansion  2k _{ s }(a l ^{5}+b l ^{3}−c l ^{2}+τ l)  8π(a r ^{5}+b r ^{3}−c r ^{2}+τ r)  \(16\sqrt {2}\left (al^{5}+bl^{3}cl^{2}+\tau l\right)\) 
a  Aggregate parameterin \(\frac {\partial E}{\partial l}\) equation  \(\frac {3{k_{v}^{2}}\lambda _{v}}{k_{s}}\)  \(\frac {4}{3}\pi \lambda _{v}\)  \(\frac {16}{9}\sqrt {2}\lambda _{v}\) 
b  Aggregate parameterin \(\frac {\partial E}{\partial l}\) equation  2k _{ s } λ _{ s }  8π λ _{ s }  \(16\sqrt {2}\lambda _{s}\) 
c  Aggregate parameterin \(\frac {\partial E}{\partial l}\) equation  \(\frac {3{k_{v}^{2}}\lambda _{v}{L_{v}^{3}}}{k_{s}}\)  \(\frac {4}{3}\pi \lambda _{v}{R_{v}^{3}}\)  \(\frac {16}{9}\sqrt {2}\lambda _{v}{L_{v}^{3}}\) 
τ  Lengthindependent component of interfacial tension  \(J2k_{s}\lambda _{s}{L_{s}^{2}}\)  \(J8\pi \lambda _{s}{R_{s}^{2}}\)  \(J16\sqrt {2}\lambda _{s}{L_{s}^{2}}\) 
ϕ  \(\frac {b}{6a}\)  \(\frac {{k_{s}^{2}}\lambda _{s}}{9{k_{v}^{2}}\lambda _{v}}\)  \(\frac {\lambda _{s}}{\lambda _{v}}\)  \(\frac {3\lambda _{s}}{2\lambda _{v}}\) 
ψ  \(\frac {c}{8a}\)  \(\frac {{L_{v}^{3}}}{8}\)  \(\frac {{R_{v}^{3}}}{8}\)  \(\frac {{L_{v}^{3}}}{8}\) 
ν (when ϕ>0)  \(\frac {\tau }{12a\phi ^{2}}\)  \(\frac {\left (J2k_{s}\lambda _{s}{L_{s}^{2}}\right)9{k_{v}^{2}}\lambda _{v}}{4{k_{s}^{3}}{\lambda _{s}^{2}}}\)  \(\frac {\left (J8\pi \lambda _{s}{R_{s}^{2}}\right)\lambda _{v}}{16{\pi \lambda _{s}^{2}}}\)  \(\frac {\left (J16\sqrt {2}\lambda _{s}{L_{s}^{2}}\right)\lambda _{v}}{48\sqrt {2}{\lambda _{s}^{2}}}\) 
μ(when ϕ>0)  \(\frac {\psi }{\phi ^{\frac {3}{2}}}\)  \(\frac {27{k_{v}^{3}}\lambda _{v}^{\frac {3}{2}}{L_{v}^{3}}}{8{k_{s}^{3}}\lambda _{s}^{\frac {3}{2}}}\)  \(\frac {{R_{v}^{3}}\lambda _{v}^{\frac {3}{2}}}{8\lambda _{s}^{\frac {3}{2}}}\)  \(\frac {{L_{v}^{3}}\lambda _{v}^{\frac {3}{2}}}{6\sqrt {6}\lambda _{s}^{\frac {3}{2}}}\) 
ν ^{′}(when ϕ=0)  \(\frac {\tau }{12a}\)  \(\frac {{Jk}_{s}}{36{k_{v}^{2}}\lambda _{v}{\vphantom {\frac {1}{2}}}}\)  \(\frac {J}{16\pi \lambda _{v}}\)  \(\frac {3J}{64\sqrt {2}\lambda _{v}}\) 
μ ^{′}(when ϕ=0)  ψ  \(\frac {{L_{v}^{3}}}{8}\)  \(\frac {{R_{v}^{3}}}{8}\)  \(\frac {{L_{v}^{3}}}{8}\) 
Bifurcation 1 (γ(l ^{∗})=0)  Transition from negative to positive interfacial tension at equilibrium  \(\nu =\frac {9{k_{v}^{2}}\lambda _{v}{L_{v}^{2}}}{2{k_{s}^{2}}\lambda _{s}}\)  \(\nu =\frac {\lambda _{v}{R_{v}^{2}}}{2\lambda _{s}}\)  \(\nu =\frac {\lambda _{v}{L_{v}^{2}}}{3\lambda _{s}}\) 
ν ^{′}=0  ν ^{′}=0  ν ^{′}=0  
Bifurcation 2 (pseudotranscritical)  Transition of l ^{∗}=0 from unstable to stable  ν=0  ν=0  ν=0 
ν ^{′}=0  ν ^{′}=0  ν ^{′}=0  
Bifurcation 3 (fold)  Transition from 2 to 0 nontrivial equilibria  ν=f(μ)(μ−f(μ)), where \(f(\mu)=\sinh \left (\frac {1}{3}\text {arcsinh} \left (\mu \right)\right)\)  ν=f(μ)(μ−f(μ))  ν=f(μ)(μ−f(μ)) 
\(\nu '=\frac {\mu '^{\frac {4}{3}}}{2^{\frac {2}{3}}}\)  \(\nu '=\frac {\mu '^{\frac {4}{3}}}{2^{\frac {2}{3}}}\)  \(\nu '=\frac {\mu '^{\frac {4}{3}}}{2^{\frac {2}{3}}}\) 
where \(\nu '=\frac {\tau }{12a}=\frac {J}{16\pi \lambda _{v}}\) (taking into account that the specific case of ϕ=0 only occurs when λ _{ s }=0, and hence τ is equal to J) and \(\mu '=\psi =\frac {{R_{v}^{3}}}{8}\). When ϕ=0, due to the absence of a surface constraint, the γ(r ^{∗})=0 bifurcation takes place at ν ^{′}=0. This is the same condition as for the pseudotranscritical bifurcation, while the bifurcation line \(\nu '=\frac {\mu '^{\frac {4}{3}}}{2^{\frac {2}{3}}}\) separates the region with two positive equilibria from the region without equilibria, again located in the first quadrant of the (μ ^{′},ν ^{′})plane.
B Cell packing: derivation of the equations for the hexagonal cell
Using a same reasoning as for the circular shape (see Eq. 11b, Eq. 12b and text below it), it follows that the interfacial tension at any nontrivial equilibrium is always positive when τ>−2λ _{ p } k _{ p } L _{ a }, and negative otherwise.
In conclusion, two equilibria with positive, real values requires that \(0<\tau <2\beta \epsilon ^{\frac {3}{2}}\,\).
For the hexagon, the perimeter p and area a can be parametrised as p=6l and \(a=\frac {3\sqrt {3}}{2}l^{2}\). Thus, k _{ p }=6 and \(k_{a}=\frac {3\sqrt {3}}{2}\). The first bifurcation then becomes τ=−12λ _{ p } L _{ a }. The second bifurcation is given by τ=0 and the third by \(\tau =2\beta \epsilon ^{\frac {3}{2}}\), with k _{ p }=6 and \(k_{a}=\frac {3\sqrt {3}}{2}\) substituted into τ, β and ε. Table 3 gives an overview of all parameters used and their meaning, for the general cell shape as well as specifically for the circle and hexagon.
C Additional parameter information regarding the CPM simulations
Figure 5
In the single cell CPM simulations, a 6th level neighbourhood was used (and hence a correction factor ξ=18 for scaling corrections, see Table 2). For the first simulation, k=1; J=80,000 (\(J_{_{\mathit {CPM}}}=4,444\phantom {\dot {i}\!}\)); A=7,854 (R _{ a }=1,250); P=0 (\(P_{_{\mathit {CPM}}}=0\phantom {\dot {i}\!}\)); λ _{ a }=2; and λ _{ p }=25 (\(\lambda _{p\,{}_{\mathit {CPM}}}=0.07716\phantom {\dot {i}\!}\)) (for switching between actual and CPM parameters, see Eq. 31, Eq. 32). For all other simulations the parameters were rescaled using Eq. 38.
In the cell sorting CPM simulations, a Moore neighbourhood was used. For the first simulation, k=1; J _{ G,G }=J _{ Y,Y }=J _{ B,B }=400, J _{ G,M }=600, J _{ Y,M }=1,200, J _{ B,M }=1,800, J _{ G,Y }=J _{ Y,B }=800, J _{ G,B }=1,400; A=30; P=100; λ _{ a }=1,000; λ _{ p }=20; and T=600. (Here, only the CPM parameters are given, which can be translated into CSM values using Eq. 31.) In the k≠1 simulations the parameters were rescaled using Eq. 38.
Figure 7
The initial cell radii used to illustrate the dynamics in the different regions were as follows: For Region I (C2), r={10, 50, 70, 150}; for Region II (D2), r={10, 50, 70, 150}; for Region III (E2), r={10, 20, 30, 36, 50, 70, 115, 150}; and for Region IV (F2) r={50, 70, 100, 150}. All simulations used a 6th level neighbourhood and A=31,416 (R _{ a }=100); P=2,000 (\(P_{_{\mathit {CPM}}}=36,000\phantom {\dot {i}\!}\)); λ _{ a }=0.0625; λ _{ p }=25 (\(\lambda _{p\,{}_{\mathit {CPM}}}=0.07716\phantom {\dot {i}\!}\)); and T=20,000. The Jvalue was modified to set τ; in regions I–IV, J=0 (\(J_{_{\mathit {CPM}}}=0\phantom {\dot {i}\!}\)), J=80,000 (\(J_{_{\mathit {CPM}}}=4,444\phantom {\dot {i}\!}\)), J=200,000 (\(J_{_{\mathit {CPM}}}=11,111\phantom {\dot {i}\!}\)), and J=300,000 (\(J_{_{\mathit {CPM}}}=16,667\phantom {\dot {i}\!}\)), respectively.
Figure 11
As the initial condition, all cells were positioned within a region of 300×300 lattice points, except for (F1), which utilized 100×100 lattice points. All simulations used a 6th level neighbourhood and A=31,416 (L _{ a }=110, R _{ a }=100); P=2,000 (\(P_{_{\mathit {CPM}}}=36,000\phantom {\dot {i}\!}\)); λ _{ a }=0.0625; λ _{ p }=25 (\(\lambda _{p\,{}_{\mathit {CPM}}}=0.07716\phantom {\dot {i}\!}\)); and T=20,000. The Jvalues were modified to set τ, in which \(J_{_{\mathit {CPM}}}=J_{C,M}=2J_{C,C}\phantom {\dot {i}\!}\). In regions I–IV, J=0 (\(J_{_{\mathit {CPM}}}=0\phantom {\dot {i}\!}\)), J=80,000 (\(J_{_{\mathit {CPM}}}=4,444\phantom {\dot {i}\!}\)), J=200,000 (\(J_{_{\mathit {CPM}}}=11,111\phantom {\dot {i}\!}\)), and J=300,000 (\(J_{_{\mathit {CPM}}}=16,667\phantom {\dot {i}\!}\)), respectively.
D Notation correspondence to Farhadifar et al. (2007)
The bifurcation line defined by Eq. 89 has the \(\overline {\Gamma }\)intercept at \(\overline {\Gamma }=\frac {2}{8\sqrt {3}}\) and \(\overline {\Lambda }\)intercept at \(\overline {\Lambda }=\frac {2\sqrt {2}}{3\cdot 3^{\frac {3}{4}}}\).
E 3D cell packing: derivation of the equations for the rhombic dodecahedronal cell
When ϕ=0, the bifurcation is simply at ν ^{′}=0. Having derived the general equations, the specific equations for a rhombic dodecahedronal cell can be straightforwardly derived using \(k_{s}=8\sqrt {2}\) and \(k_{v}=\frac {16}{3\sqrt {3}}\), and are given in Table 6.
Abbreviations
 2D:

2dimensional Euclidean space
 3D:

3dimensional Euclidean space
 CPM:

Cellular Potts model
 CSM:

Cell surface mechanics
 DAH:

Differential adhesion hypothesis
 MCS:

Monte Carlo step
Declarations
Acknowledgements
RM received support from the Portuguese Fundação para a Ciência e
Tecnologia — fellowship SFRH/BD/32966/2006 — and grant PTDC/BIMMED/1063/2012. VAG gratefully acknowledges support from the Royal Society Dorothy Hodgkin fellowship. The work was supported by the UK Biological and Biotechnology Research Council (BBSRC) via grant BB/J004553/1 to the John Innes Centre.
Authors’ Affiliations
References
 Bausch AR, Ziemann F, Boulbitch AA, Jacobson K, Sackmann E. Local measurements of viscoelastic parameters of adherent cell surfaces by magnetic bead microrheometry. Biophys J. 1998; 75(4):2038–049.View ArticleGoogle Scholar
 Bausch AR, Moller W, Sackmann E. Measurement of local viscoelasticity and forces in living cells by magnetic tweezers. Biophys J. 1999; 76(1):573–9.View ArticleGoogle Scholar
 Lau AW, Hoffman BD, Davies A, Crocker JC, Lubensky TC. Microrheology, stress fluctuations, and active behavior of living cells. Phys Rev Lett. 2003; 91(19):198101.View ArticleADSGoogle Scholar
 Heidemann SR, Wirtz D. Towards a regional approach to cell mechanics. Trends Cell Biol. 2004; 14(4):160–6. doi:10.1016/j.tcb.2004.02.003.View ArticleGoogle Scholar
 Discher DE, Janmey P, Wang YL. Tissue cells feel and respond to the stiffness of their substrate. Science. 2005; 310(5751):1139–43. doi:10.1126/science.1116995.View ArticleADSGoogle Scholar
 Brangwynne CP, MacKintosh FC, Kumar S, Geisse NA, Talbot J, Mahadevan L, et al. Microtubules can bear enhanced compressive loads in living cells because of lateral reinforcement. J Cell Biol. 2006; 173(5):733–41. doi:10.1083/jcb.200601060.View ArticleGoogle Scholar
 Kumar S, Maxwell IZ, Heisterkamp A, Polte TR, Lele TP, Salanga M, et al. Viscoelastic retraction of single living stress fibers and its impact on cell shape, cytoskeletal organization, and extracellular matrix mechanics. Biophys J. 2006; 90(10):3762–773. doi:10.1529/biophysj.105.071506.View ArticleGoogle Scholar
 Kasza KE, Rowat AC, Liu J, Angelini TE, Brangwynne CP, Koenderink GH, et al. The cell as a material. Curr Opin Cell Biol. 2007; 19(1):101–7. doi:10.1016/j.ceb.2006.12.002.View ArticleGoogle Scholar
 Loesberg WA, te Riet J, van Delft FC, Schon P, Figdor CG, Speller S, et al. The threshold at which substrate nanogroove dimensions may influence fibroblast alignment and adhesion. Biomaterials. 2007; 28(27):3944–951. doi:10.1016/j.biomaterials.2007.05.030.View ArticleGoogle Scholar
 Massiera G, Van Citters KM, Biancaniello PL, Crocker JC. Mechanics of single cells: rheology, time dependence, and fluctuations. Biophys J. 2007; 93(10):3703–13. doi:10.1529/biophysj.107.111641.View ArticleGoogle Scholar
 Mizuno D, Tardin C, Schmidt CF, Mackintosh FC. Nonequilibrium mechanics of active cytoskeletal networks. Science. 2007; 315(5810):370–3. doi:10.1126/science.1134404.View ArticleADSGoogle Scholar
 Saez A, Ghibaudo M, Buguin A, Silberzan P, Ladoux B. Rigiditydriven growth and migration of epithelial cells on microstructured anisotropic substrates. Proc Nat Acad Sci USA. 2007; 104(20):8281–6. doi:10.1073/pnas.0702259104.View ArticleADSGoogle Scholar
 In: (Anderson ARA, Chaplain MAJ, Rejniak KA, editors.)SingleCellBased Models in Biology and Medicine. Basel: Birkhäuser Verlag; 2007.MATHGoogle Scholar
 Beltman JB, Marée AFM, Lynch JN, Miller MJ, de Boer RJ. Lymph node topology dictates T cell migration behavior. J Exp Med. 2007; 204(4):771–80. doi:10.1084/jem.20061278.View ArticleGoogle Scholar
 Krieg M, ArboledaEstudillo Y, Puech PH, Käfer J, Graner F, Muller DJ, et al. Tensile forces govern germlayer organization in zebrafish. Nat Cell Biol. 2008; 10(4):429–36. doi:10.1038/ncb1705.View ArticleGoogle Scholar
 Marinari E, Mehonic A, Curran S, Gale J, Duke T, Baum B. Livecell delamination counterbalances epithelial growth to limit tissue overcrowding. Nature. 2012; 484(7395):542–5. doi:10.1038/nature10984.View ArticleADSGoogle Scholar
 Lecuit T, Lenne PF. Cell surface mechanics and the control of cell shape, tissue patterns and morphogenesis. Nat Rev Mol Cell Biol. 2007; 8(8):633–44. doi:10.1038/nrm2222.View ArticleGoogle Scholar
 Graner F, Glazier JA. Simulation of biological cell sorting using a twodimensional extended Potts model. Phys Rev Lett. 1992; 69(13):2013–6.View ArticleADSGoogle Scholar
 Glazier JA, Graner F. Simulation of the differential adhesion driven rearrangement of biological cells. Phys Rev E. 1993; 47(3):2128–54.View ArticleADSGoogle Scholar
 Steinberg MS. Reconstruction of tissues by dissociated cells: some morphogenetic tissue movements and the sorting out of embryonic cells may have a common explanation. Science. 1963; 141:401–8.View ArticleADSGoogle Scholar
 Käfer J, Hogeweg P, Marée AFM. Moving forward moving backward: directional sorting of chemotactic cells due to size and adhesion differences. PLoS Comput Biol. 2006; 2(6):56. doi:10.1371/journal.pcbi.0020056.View ArticleADSGoogle Scholar
 Ouchi NB, Glazier JA, Rieu JP, Upadhyaya A, Sawada Y. Improving the realism of the cellular Potts model in simulations of biological cells. Physica A. 2003; 329(3–4):451–8.View ArticleADSMATHMathSciNetGoogle Scholar
 Rieu JP, Upadhyaya A, Glazier JA, Ouchi NB, Sawada Y. Diffusion and deformations of single Hydra cells in cellular aggregates. Biophys J. 2000; 79(4):1903–14.View ArticleGoogle Scholar
 Käfer J, Hayashi T, Marée AFM, Carthew RW, Graner F. Cell adhesion and cortex contractility determine cell patterning in the Drosophila retina. Proc Nat Acad Sci USA. 2007; 104(47):18549–54. doi:10.1073/pnas.0704235104.View ArticleADSGoogle Scholar
 Farhadifar R, Roper JC, Aigouy B, Eaton S, Julicher F. The influence of cell mechanics, cellcell interactions, and proliferation on epithelial packing. Curr Biol. 2007; 17(24):2095–104. doi:10.1016/j.cub.2007.11.049.View ArticleGoogle Scholar
 Hilgenfeldt S, Erisken S, Carthew RW. Physical modeling of cell geometric order in an epithelial tissue. Proc Nat Acad Sci USA. 2008; 105(3):907–11. doi:10.1073/pnas.0711077105.View ArticleADSGoogle Scholar
 Honda H, Nagai T, Tanemura M. Two different mechanisms of planar cell intercalation leading to tissue elongation. Dev Dyn. 2008; 237(7):1826–36. doi:10.1002/dvdy.21609.View ArticleGoogle Scholar
 Du X, Osterfield M, Shvartsman SY. Computational analysis of threedimensional epithelial morphogenesis using vertex models. Phys Biol. 2014; 11(6):066007. doi:10.1088/14783975/11/6/066007.View ArticleADSGoogle Scholar
 Manning ML, Foty RA, Steinberg MS, Schoetz EM. Coaction of intercellular adhesion and cortical tension specifies tissue surface tension. Proc Nat Acad Sci USA. 2010; 107(28):12517–22. doi:10.1073/pnas.1003743107.View ArticleADSGoogle Scholar
 Staple DB, Farhadifar R, Roper JC, Aigouy B, Eaton S, Julicher F. Mechanics and remodelling of cell packings in epithelia. Eur Phys J E. 2010; 33(2):117–27. doi:10.1140/epje/i2010106770.View ArticleGoogle Scholar
 Iwamoto M, Sugino K, Allen RD, Naitoh Y. Cell volume control in Paramecium: factors that activate the control mechanisms. J Exp Biol. 2005; 208(3):523–37. doi:10.1242/jeb.01417.View ArticleGoogle Scholar
 de Vries WN, Evsikov AV, Haac BE, Fancher KS, Holbrook AE, Kemler R, et al. Maternal βcatenin and Ecadherin in mouse development. Development. 2004; 131(18):4435–45. doi:10.1242/dev.01316.View ArticleGoogle Scholar
 Forgacs G, Foty RA, Shafrir Y, Steinberg MS. Viscoelastic properties of living embryonic tissues: a quantitative study. Biophys J. 1998; 74(5):2227–34.View ArticleGoogle Scholar
 Grieneisen VA. Estudo do estabelecimento de configurações em estruturas celulares. Porto Alegre: Master’s thesis, Universidade Federal do Rio Grande do Sul. November 2004.
 Marée AFM, Grieneisen VA, Hogeweg P. The Cellular Potts Model and biophysical properties of cells, tissues and morphogenesis In: Anderson ARA, Chaplain MAJ, Rejniak KA, editors. SingleCellBased Models in Biology and Medicine. Basel: Birkhäuser Verlag: 2007. p. 107–36.Google Scholar
 Rubenstein BM, Kaufman LJ. The role of extracellular matrix in glioma invasion: a cellular Potts model approach. Biophys J. 2008; 95(12):5661–80. doi:10.1529/biophysj.108.140624.View ArticleADSGoogle Scholar
 Xu Z, Lioi J, Mu J, Kamocka MM, Liu X, Chen DZ, et al. A multiscale model of venous thrombus formation with surfacemediated control of blood coagulation cascade. Biophys J. 2010; 98(9):1723–32. doi:10.1016/j.bpj.2009.12.4331.View ArticleGoogle Scholar
 Marée AFM, Hogeweg P. How amoeboids selforganize into a fruiting body: multicellular coordination in Dictyostelium discoideum. Proc Nat Acad Sci USA. 2001; 98(7):3879–83. doi:10.1073/pnas.061535198.View ArticleADSGoogle Scholar
 Purcell EM. Life at low Reynolds number. Am J Phys. 1977; 45(1):3–11.View ArticleADSMathSciNetGoogle Scholar
 Nordlund TM. Quantitative Understanding of Biosystems: an Introduction to Biophysics. Boca Raton: CRC Press; 2011.Google Scholar
 Osserman R. The isoperimetric inequality. Bull Amer Math Soc. 1978; 84(6):1182–238.View ArticleMATHMathSciNetGoogle Scholar
 Glazier JA. Dynamics of cellular patterns: PhD thesis, University of Chicago; 1989.
 Deutsch A, Dormann S. Cellular Automaton Modeling of Biological Pattern Formation. Boston: Birkhäuser; 2004.Google Scholar
 Weisstein EW. Sum of squares function. Visited on 14/05/2015. http://mathworld.wolfram.com/SumofSquaresFunction.html.
 Zajac M, Jones GL, Glazier JA. Simulating convergent extension by way of anisotropic differential adhesion. J Theor Biol. 2003; 222(2):247–59.View ArticleGoogle Scholar
 Marée AFM, Jilkine A, Dawes A, Grieneisen VA, EdelsteinKeshet L. Polarization and movement of keratocytes: a multiscale modelling approach. Bull Math Biol. 2006; 68(5):1169–211. doi:10.1007/s1153800691317.View ArticleGoogle Scholar
 Grieneisen VA, Xu J, Marée AFM, Hogeweg P, Scheres B. Auxin transport is sufficient to generate a maximum and gradient guiding root growth. Nature. 2007; 449(7165):1008–13. doi:10.1038/nature06215.View ArticleADSGoogle Scholar
 Chen CS, Mrksich M, Huang S, Whitesides GM, Ingber DE. Geometric control of cell life and death. Science. 1997; 276(5317):1425–8.View ArticleGoogle Scholar
 Ruoslahti E. Stretching is good for a cell. Science. 1997; 276(5317):1345–6.View ArticleGoogle Scholar
 Hogeweg P. Evolving mechanisms of morphogenesis: on the interplay between differential adhesion and cell differentiation. J Theor Biol. 2000; 203(4):317–33. doi:10.1006/jtbi.2000.1087.View ArticleGoogle Scholar
 Newman TJ. Modeling multicellular systems using subcellular elements. Math Biosci Eng. 2005; 2(3):613–24.View ArticleMATHMathSciNetGoogle Scholar
 Harvey CW, Morcos F, Sweet CR, Kaiser D, Chatterjee S, Liu X, et al. Study of elastic collisions of Myxococcus xanthus in swarms. Phys Biol. 2011; 8(2):026016. doi:10.1088/14783975/8/2/026016.View ArticleADSGoogle Scholar
 Oakes PW, Banerjee S, Marchetti MC, Gardel ML. Geometry regulates traction stresses in adherent cells. Biophys J. 2014; 107(4):825–33. doi:10.1016/j.bpj.2014.06.045.View ArticleADSGoogle Scholar
 Namias V. Simple derivation of the roots of a cubic equation. Am J Phys. 1985; 53:775. doi:10.1119/1.14311.View ArticleADSGoogle Scholar
 Murray JD. Mathematical Biology I: An Introduction, 3rd edn. New York: Springer; 2002.MATHGoogle Scholar
Copyright
This is an Open Access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly credited. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated.