Single and two-cells shape analysis from energy functionals for three-dimensional vertex models

Vertex models have been extensively used for simulating the evolution of multicellular systems, and have given rise to important global properties concerning their macroscopic rheology or jamming transitions. These models are based on the definition of an energy functional, which fully determines the cellular response and conclusions. While two-dimensional vertex models have been widely employed, three-dimensional models are far more scarce, mainly due to the large amount of configurations that they may adopt and the complex geometrical transitions they undergo. We here investigate the shape of single and two-cells configurations as a function of the energy terms, and we study the dependence of the final shape on the model parameters: namely the exponent of the term penalising cell-cell adhesion and surface contractility. In single cell analysis, we deduce analytically the radius and limit values of the contractility for linear and quadratic surface energy terms, in 2D and 3D. In two-cells systems, symmetrical and asymmetrical, we deduce the evolution of the aspect ratio and the relative radius. While in functionals with linear surface terms yield the same aspect ratio in 2D and 3D, the configurations when using quadratic surface terms are distinct. We relate our results with well-known solutions from capillarity theory, and verify our analytical findings with a three-dimensional vertex model.

At the cellular scale, due to the explicit representation of cell boundaries, vertex models have been exploited to simulate tissue rheology and evolution from cellular mechanics. 10Initially, the vertex models were developed to understand the groups of joint bubbles in foams. 14,15Vertex models discretise epithelial cells by a set of vertices 10,16 and impose static or dynamic equilibrium jointly by volumetric constraints.As such, they are equally valid to describe cell packings or capillarity problems composed of bubble-like bodies.Although our aim is to analyse cellular geometries, our result can be also applied to non-biological problems driven by surface area minimisation and volumetric constraints.Two-dimensional vertex models have been widely exploited to describe fluidity of tissues, 17 topological transitions [18][19][20][21] or jamming. 22Three-dimensional extensions have been developed for epithelial rheology, 23,24 folding, 25 blastocytes analysis, 26 wound healing, 27 cyst formation, 28 or cell division, 29 to cite a few.
Vertex models require the definition of a functional that depends on the vertex coordinates y ¼ y 1 , …, y N vert È É , which is minimised with respect to each vertex position y I , I ¼ 1, …, N vert .In fact, the form of this functional determines the forces at equilibrium, and the global predictions that we can derive from the model.The general form of the energy functional depends explicitly on geometrical measures of the cell, and may be expressed in two dimensions (2D) as where A i and P i are respectively the area and perimeter of the cell i, A 0 a reference preferential area, λ and Γ the area and perimeter modulus, l mn is the length between vertices m and n and Λ mn the tension of segment y m À y n .The equivalent form of the functional in three dimensions (3D) reads, with V 0 a preferential reference volume.Parameters λ c and λ e represent penalisations of different surface areas (cell-cell vs. cell-environment, lateral vs. basal, etc.).Depending on the biological problem, some authors have considered additional terms for bending, 30 or hybrid energy definitions with 3D vertex coordinates but with cells located on a two dimensional manifold, and with additional normal pressure forces and constraints. 25,31he functionals W 2D y ð Þ and W 3D y ð Þ additionally depend on the exponents n V , n P , n L , n A r and n A c .In the literature, there is a general agreement in using quadratic terms for area and volume penalisation in 2D and 3D, respectively, that is, in using n V ¼ 2, which is tantamount to reduce volume variations.However, the exponents of other terms may vary depending on the authors and applications.Tables 1 and 2 show some of the values of the exponents for surface energy T A B L E 1 Order of exponent in different terms in 2D vertex models.See Equation (1).

Reference
Line (n L ) Perimeter (n P ) Area (n V ) Duclut et al., 32 Barton et al., 33 Osterfield et al. 34 1 2 2 Fletcher et al., 35 Farhadifar et al. 19 Bi et al., 22,36 -1,2 2 Spencer et al., 37   term (perimeter in 2D).The use of quadratic terms, that is, n P ¼ 2, may be interpreted as an elastic response, where the resulting forces g I ¼ ∂ y I W are proportional to the area (or perimeter in 2D).Instead, in linear functionals with n A ¼ n P ¼ 1, the forces do not vary as area or perimeter increases, and are thus more pertinent for fluid-like behaviour, where the constants Γ, λ, λ c , λ e represent surface tensions.The objective of the present work is to analyse the effects on the resulting shape when using linear or quadratic functionals, for simple cell configurations, and study their differences in 2D and 3D models.We will derive in Section 2 some analytical results for single and two-cells configurations, and in the latter case, focus on the aspect ratio or relative area of the cell-cell surface.We will assume axial symmetry, and thus our system is fully defined by the radius r of a single cell, or the radius and parameter h, corresponding to half distance of the maximum overlapping, as indicated in Figure 1.In Section 3 we will verify our theoretical results with a numerical 3D vertex model.

| ANALYTICAL RESULTS
In the following derivations we initially assume an energy functional defined for a three-dimensional cellular system, which we will adapt for comparing our results with two-dimensional configurations.A quadratic term for the volume penalisation is assumed n V ¼ 2 ð Þ, while we consider two possible values of the surface exponent, n A ¼ 1 (linear functional) and n A ¼ 2 (quadratic functional).More precisely, the following general expression of the energy functional is studied, W r, h ð Þ¼ 1 2 Þ the interface area between cells i and j, and A e i r, h ð Þ the external area of cell i.We write our functional in terms of the radius r of each cell and overlapping h between them, instead of distances in order to simplify our subsequent derivations, and we restrain our analysis to one cell and two cells, that is, N ¼ 1, 2 f g.For N ¼ 1, we set A c ¼ 0 and assume radial symmetry, while for N ¼ 2 we consider cylindrical symmetry with respect to the axis joining the two cells centres.The common area A c is counted twice, one for each cell sharing the surface.This assumption would not be realistic in the analysis of droplets, but is consistent with the presence of two cortical surfaces at cell-cell contacts.The values V 0 and A 0 are the reference volume and areas, respectively, and r 0 ¼ ffiffiffiffiffiffiffiffiffi ffi q the reference radius.Note that in our analysis, we use the same exponent n A for all surface area terms.
We note that the resulting configuration of the cell system, which corresponds to the minimisation of the functional W with respect to the independent variables (radius r and overlapping h), remains unchanged when the functionals is scaled by a constant factor.For this reason, we have normalised our functional and omitted the volume parameter λ v , and have left only the cell-cell and cell-environment factors λ c ij ð Þ and λ e i .These factors represent in fact relative penalisation of surface areas with respect to bulk (volumetric) modulus.Furthermore, the parameter n A in the denominator allows us to interpret the factors λ c ij ð Þ and λ e 1 as the conjugate forces (surface tension) in the linear case, or the surface elasticity constant in the solid case.That is, , where forces are proportional to the area.
2.1 | Single cell analysis

| Linear functional
Due to the absence of interface surface, we ignore the second term (3), and deduce the radius of equilibrium as the value of r that minimises the following functional, The optimality condition dW 1 l r ð Þ=dr ¼ 0 may be expressed as, Real roots of the equation above in the interval of interest r 0, that is, when the surface tension λ e is lower than a critical value, Above this value, the functional W 1 L is a monotonically increasing function for r ≥ 0, and thus the minimiser corresponds to the degenerated configuration r ¼ 0. For λ e < λ crit , there are two positive real roots.We consider the largest one, which converges towards r 0 as λ e decreases towards 0.
In this linear case, and after introducing the parameter λ e ¼ λ v λ e V 2 0 =A 0 as the (not normalised) surface traction, with λ v the bulk modulus, it can be verified that our functional in ( 4) is a normalised version of the energy In this case, the equilibrium condition W 0 ¼ 0 turns into which after identifying the term Þas the cell internal pressure P, this equation corresponds to relation P ¼ 2λ e =r, that is, the Young-Laplace relation. 45

| Quadratic functional
The energy functional for n A ¼ 2 reads which implies the following optimality condition, In view of this equation, and as expected, for λ e ¼ 0, the optimal radius is r ¼ r 0 , while for λ e > 0, there is always at least one real positive root.The roots of ( 5) and ( 8) as a function of λ e are shown in Figure 2. It can be checked that the optimal radius is reduced in a more pronounced manner with respect to changes in λ e in the quadratic case, while in the linear case the cell collapses for λ e ≥ λ crit .
Similar analysis in two dimensions yields an equivalent critical value of λ e for linear functionals with A 0 and P 0 reference area and perimeter.In contrast to the three-dimensional case though, a quadratic functional Þin 2D exhibits also a critical value Beyond this value, no solution exists other than the degenerated value r ¼ 0. These results show that reduced 2D models have strong qualitative differences, and care must be taken when calibrating surface tension.
Figure 2 shows the value of r that minimises the energy functional for two and three dimensions.The critical values of λ are indicated with thick dots.It can be observed that the sole case that admits unbounded values of λ e is the quadratic functional in 3D.
Also, a similar equation to the Young-Laplace relation may be written for the quadratic case by setting the internal pressure as , and the relation between pressure and stiffness at equilibrium turns into P ¼ 2λ e A=r.This is similar to the Young-Laplace relation but with an additional area factor A r ð Þ due to the surface elasticity.

| Symmetric two-cells system analysis
After assuming an axis-symmetric configuration, where the volume and areas are defined as a function of spherical radius r and overlapping distance h (half of the maximum overlapping length, see Figure 1), the volume, the external area and the common interface area of the two cells can be expressed as, r as a function of λ e for single cell, linear and quadratic functionals.Left: results for 3D (A 0 ¼ 4π, V 0 ¼ 4 3 π, which implies λ crit ¼ 0:709, according to Equation ( 6)).Right: results for 2D A 0 ¼ π, P 0 ¼ 2π, which yields In order to analyse the resulting shape, we will use the shape factor h=r. Due to the symmetry assumptions, the energy due to the contributions of the two cells reads, In the next paragraphs, we deduce the values of h, r ð Þ that minimise this functional in the linear

| Linear functional
The optimal values of r and h are those that minimise the gradient, that is, rW These conditions can be written as, This is a system of two non-linear equations.By isolating the term V =V 0 À 1 ð Þ =V 0 from one of them, we obtain, which after using the expressions in (9) yields, Two important conclusions can be derived: the shape factor h=r is independent of A 0 and V 0 , and this factor solely depends on the ratio of the penalty terms λ e =λ c , but not on their actual values.Also, note that we are interested in optimal values r Ã , h Ã ð Þdifferent from the trivial solution r Ã , h Ã ð Þ¼ 0, 0 ð Þ, which is also a solution of (11).

| Quadratic functional
After setting n A ¼ 2 in (10), the equilibrium equations rW ¼ 0 now become, By replacing the factor V =V 0 À 1 ð Þ =V 0 from one of the equations, we find that It is interesting to compare the results in ( 13)-( 15) with those obtained in two dimensions.This simplified model is reasonable for elongated columnar cells, but is commonly used in 2D vertex models, which motivates our analysis.By replacing the volume by the area A ð Þ and the surface area by the perimeter P ð Þ in Equations 11 and 14, and using the following expressions of these measures, Þis half of the opening angle of the contact interface, (see Figure 1), the following relations between the ratio λ c =λ e and the shape factor are obtained, While the linear functional gives the same relation between λ c =λ e and the shape factor h=r for two and three dimensions, the quadratic functional differs, despite giving similar trends.Figure 3 shows the evolution of the shape factor h=r as a function of the ratio λ e , λ c ð Þfor all the studied cases.We remark that: while in the quadratic case, the cells split h=r ¼ 0 ð Þonly when λ e ¼ 0, in the linear case this happens whenever λ c ≥ λ e .Moreover, the shape factor h=r is proportional to λ c in the linear case, while in the quadratic case is non-linear, and converges towards h=r ¼ 0 as λ c =λ e increases.In both cases, the single spherical shape h ¼ r ð Þ is recovered when λ c ¼ 0 and λ e > 0. Figure 4A shows the value of h=r on the parameter plane λ c , λ e ð Þ.

| Contact angles
The contact angle of a contractile adhesive sphere is a usual measure for describing the mechanics of cells in contact with a substrate. 45,46We here compare our findings with these results, which assume a linear energy dependence, and are thus only valid for a linear functional n A ¼ 1 ð Þ.The relation of contact angle θ indicated in Figure 1 with the shape factor h=r can be deduced using simple geometrical arguments, F I G U R E 3 Evolution of shape factor h=r (A) and contact angle θ (B) in two cell-system as a function of ratio λ c =λ e for the linear and quadratic functional in 2D and 3D.
Furthermore, from capillarity theory, and in the absence of an adhesive solid substrate, we have from tension equilibrium that, 45 which together with (19) yields our results in ( 13) and ( 17) for a linear functional, where the penalty parameters play the role of surface tension.If a quadratic functional is used instead, relations in ( 17) and ( 19) allow us to write Therefore, when using a quadratic term in the energy functional, the factor λ c =λ e also fully determines the shape factor and the contact angle θ, but in a different manner to the results in capillary theory.Figure 4B shows the trend of the angle as a function of the factor λ c =λ e .

| Asymmetrical two-cells system
We study here the stability of the two-cell system when the two contractility parameters λ e of the two cells are unequal.Let us call λ e 1 the value for one of the cells, and λ e 2 ¼ kλ e 1 the value for the neighbouring cell.Figure 5 depict the three independent geometrical variables, r 1 , r 2 and h 1 .The overlapping h 2 , indicated in Figure 5, can be deduced from geometrical considerations as In the subsequent derivations, we will employ the shape factors r 2 =r 1 and h 1 =r 1 , which fully describe the relative size and overlapping of the two-cells system, and restrict our study to the 3D case.
After assuming that reference volume and area V 0 and A 0 , the energy functional reads, and the equilibrium conditions rW ¼ 0 as, The gradients have now the derivatives with respect to the three independent variables, that is, The system of three equations above allows us to isolate the factors containing V 0 and A 0 , and deduce the following relation, Expressions in ( 21) and ( 22) should be completed with the following relations, which can be derived from (9).
An important preliminary result can be concluded from the expression in (22): contrarily to the symmetric case, the shape, which is defined by factors r 2 =r 1 and h 1 =r 1 , is not uniquely defined by the ratios λ c =λ e 1 and k ¼ λ e 2 =λ e 1 .Indeed, for arbitrary values of these ratios, Equation ( 22) is satisfied for multiple values of r 2 =r 1 and h 1 =r 1 .Said otherwise, the shape factors, which can be computed by solving (21), also depend on the reference volume V 0 and area A 0 .This is confirmed in Figure 6, where after solving Equation ( 21) numerically, it is shown that the evolution of the shape factors h 1 =r 1 and r 2 =r 1 are independent of V 0 in the symmetric case k ¼ 1 ð Þ, while for k ≠ 1, the shape factors are indeed different for different values of V 0 .
In Figure 7 we additionally show the evolution of the same shape factors as a function of λ c =λ e 1 , and for two values of k.For the symmetric case, k ¼ 1, we have the same curves shown in 3 for the 3D case, while the relative increase of λ e 2 with respect to λ e 1 produces in turn an increase of the overlapping h 1 =r 1 , but very minor effects on the relative sizes F I G U R E 5 Scheme of asymmetrical two-cells system with radius r 1 and r 2 .
r 2 =r 1 .The latter is only noticeable on the figure for the linear case (continuous line).In the quadratic functional, with λ c > 0, the factor r 2 =r 1 is just smaller than one, with differences close to 1e À 4.

| Vertex Model
We aim at verifying the previous results with a vertex model, where cell surface triangularised with a set of vertices y i .A discrete version of the energy functional can be then built from the set of vertices y ¼ y 1 , …, y N f gand by defining the volume discrete measures of cell area and volume:  The discretised energy functional mimics the definitions in ( 3), but with the previous measures of area and volume.In addition, and in order to ease the numerical convergence, the functional is minimised in consecutive pseudo-time steps Δt ¼ t n À t nÀ1 , which is tantamount to regularise the minimisation process.In summary, the following discrete energy functional is employed: with Δy ¼ y À y nÀ1 and η a regularisation parameter.For simplicity we denote by y all the vertex positions y 1 , …, y N vert at current time t n .In the results presented next, we let the geometry to evolve until Δy < 1e À 6 y .At each time step t n , the functional W h y ð Þ is minimised with respect to the new positions y n , which is equivalent to solving the following system of non-linear equations: These equations are solved in an implicit manner resorting to Newton-Raphson process.This requires the computation of the Jacobain.In Section (A) we detail the system of equations and the solution process.
The initial vertex positions y 0 and triangulation is constructed by defining a cloud of points, so called nodes x i , around the two estimated positions of the cell-centres.From these nodes, a Delaunay triangulation (tetrahedra in 3D) is built (shown as thin red lines in Figure 8A), and the vertices are placed at each barycentre of each tetrahedra.An additional vertex at each segment joining the cell centre y c and the connected nodes x i is added in order to have the external and common surfaces fully triangularised.The nodal network x 1 , …x N x also allows to track cell-cell connectivity in case of topological transitions are allowed.The study of these transitions is under investigation and out of scope of the present paper.
This model is an extension of the previous work presented in, 27 where the geometry was though limited to a monolayer, with two parallel apical and basal surfaces.Here, no restrictions on the geometry are imposed, other than the size of the discretisation, which is regulated by the density of the nodal cloud.Models with finer density can be found for instance in, 13,29 which may be considered as an extension of the coarser vertex models.

| Shape factor on discretised model
We have run the 3D vertex model using the linear and quadratic functionals.We have used a discretisation that is sufficient to estimate the shape factor h=r by fitting two spheres on the resulting set of vertices y i .The fitting was achieved by solving the following least squares problem: where P r,h y i ð Þ is the projection of vertex y i on the surface of two spheres of radius r and half-overlapping equal to h.The axis that joints the two centres of the spheres and the mid plane is deduced from the whole set of vertices y i .Alternatively, we note that from the expressions in ( 9), the ratio of A c and A e is given by, Therefore, the shape factor h=r can be also estimated from the value of 2A c =A e .Figure 9 shows the latter estimate, which coincides with a relative error smaller than 5% with respect to the first estimate.This figure also shows the theoretical results obtained in Section 2.2.As it can be observed, the shape factors of the vertex model confirm the theoretical results.

| CONCLUSIONS AND DISCUSSION
The present paper has focused in two main contributions: (i) derivation of analytical functional minimisers: radius in single cells, and shape factors in two-cells systems, and (ii) development and validation of numerical model for analysing cell systems, regarded as those driven by volume constraints and contact area minimisations.Motivated by the different energy functionals employed in the literature, we have deduced the influence of linear and quadratic surface area terms in the functional.
In single cell analysis, we have deduced the presence of limit values of the surface tension in linear and quadratic functionals for 2D, but only for the linear case in 3D.In this last case, there is a minimum radius beyond which the cell collapses.In the quadratic case instead, the resulting radius can be reduced by further increasing λ e .The presence of these limits is important prior to implementing numerical models for simulating cell responses.
In linear functionals, the contact angle is a direct experimental measure of the surface tension.We have shown that this angle reduces to 0 when λ c approaches λ e , while in the quadratic case such angle only vanishes in the limit case λ c =λ e !∞. when considering asymmetrical two-cells configurations, our analytical results indicate that the difference in external adhesion has effects on the overlapping between the cells, but very minor on their radius.The Equation in (A1) can be expressed as g y ð Þ ¼ 0, which is solved through a Newton-Raphson iterative process, which at each iteration k solves the following system of linear equations: with K ¼ ∂g ∂y y¼y k the Jacobian matrix.The vertex coordinates are iteratively updated as y kþ1 ¼ y k þ δy.Matrix K can be derived from the expressions in (A2).Convergence was accepted whenever max g y k À Á , δy < 1e À 10.We have used the initial values r 0 ¼ 1 and η ¼ 1:0.The latter parameter has been incrementally decreased until η ¼ 1e À 3. At point, the difference between time-steps satisfied Δy < 1e À 3.

X e C det y e 1
À y e c , y e 2 À y e c , y e 3 À y e c À Á ð23Þ where e denotes each triangle of the cell surface formed by three vertices y e 1 , y e 2 , y e 3 È É , and y e c is an arbitrary point, which we choose to be the geometrical centre y e c ¼ P N c i¼1 y e i =N c with N c the number of vertices in cell c. Figure 8 illustrates one of such triangles and tetrahedral.

2 F 1 F
I G U R E 6 Evolution of shape factors h 1 =r 1 and r 2 =r 1 as a function of material asymmetry k ¼ λ e2 =λ e1 for two values of reference volume V 0 .Continuous (À) and discontinuous (À À) line styles correspond to linear and quadratic cases, n A ¼ 1 and n A ¼ 2, respectively.All cases run with λ c ¼ 0:025 and λ e1 ¼ 0:05.I G U R E 7 Evolution of shape factors h 1 =r 1 and r 2 =r 1 as a function of ratio λ c =λ e1 for symmetric k ¼ 1 and unsymmetric (k ¼ 1:5) cases.Continuous (À) and dicontinuous (À À) line styles correspond to linear and quadratic cases, n A ¼ 1 and n A ¼ 2, respectively.

8
The geometry of 3D vertex model for two cells presented by a vertex network (dark lines) and nodal network (red thinner), shown in (A).Areas and volume are computed from the discretised triangularised surface of each cell, as depicted in (B).

1 F
I G U R E 9 Comparison of theoretical results and those obtained from vertex 3D model.
Order of exponent in different terms in 3D vertex models.See Equation(2).
V ) Honda et al.