Interacting active surfaces: A model for three-dimensional cell aggregates

We introduce a modelling and simulation framework for cell aggregates in three dimensions based on interacting active surfaces. Cell mechanics is captured by a physical description of the acto-myosin cortex that includes cortical flows, viscous forces, active tensions, and bending moments. Cells interact with each other via short-range forces capturing the effect of adhesion molecules. We discretise the model equations using a finite element method, and provide a parallel implementation in C++. We discuss examples of application of this framework to small and medium-sized aggregates: we consider the shape and dynamics of a cell doublet, a planar cell sheet, and a growing cell aggregate. This framework opens the door to the systematic exploration of the cell to tissue-scale mechanics of cell aggregates, which plays a key role in the morphogenesis of embryos and organoids.


Notations of differential geometry
Here we summarise notations of differential geometry used in this study, which follow Refs. [1] and [2]. The surface S is described by the parametrisation X(s 1 , s 2 ), with X a 3D vector and s 1 , s 2 two surface coordinates. The cartesian basis of the 3D space in which the surface is embedded is denotedẽ α , with α = x, y, z. In the following we use greek indices to label 3D coordinates, and latin indices to label coordinates on the surface. The tangent vectors e i , normal vector n, metric and curvature tensors g ij and C ij and Christoffel symbols Γ k ij are respectively defined by: The determinant of the metric tensor is denoted g = detg ij , and the surface element is dS = √ gds 1 ds 2 . The vectors e 1 , e 2 of the dual basis are defined by: The inverse of the metric tensor g ij is denoted g ij , such that g ik g kj = δ i j . We also introduce the Levi-Civita tensor on a curved surface: which verifies ij j k = −g ik . A 3D vector a = a αẽ α has tangential and normal components to the surface, according to: a = a i e i + a n n = a i e i + a n n .
The metric tensor g ij and its inverse g ij can be used to raise and lower indices of vector and tensor components of the surface; for instance a i = g ij a j and a i = g ij a j . We also introduce covariant derivatives of a tangent vector a and tangent second-rank tensor b on the surface as follows: The covariant derivatives of the metric and Levi-Civita tensor vanish, ∇ i g jk = 0 and ∇ i jk = 0. In the following we also use the notation ∂ i = g ij ∂ j and ∇ i = g ij ∇ j .
2 Variations of the metric, curvature and Christoffel symbols Given a variation of the surface parametrisation δX, we can compute the variation of the different geometrical quantities of the surface. For the tangent vectors, we have The variation of the metric tensor can then be computed as: δg ij = δ(e i · e j ) = δe i · e j + e i · δe j = ∂ i δX · e j + e i · ∂ j δX .
For the inverse of the metric, we can use that g ij g jk = δ i k , and then δ(g ij g jk ) = 0 =⇒ δg ij = −g ik g jl δg kl .
The variation of the metric determinant g = det (g ij ) is given by Jacobi's formula: As a result, we have the useful relation for the infinitesimal variation of the square root of the metric determinant: Eq. (12) can be used then to compute variations of the basis vectors of the dual basis δe i = δ g ij e j = δg ij e j + g ij ∂ j δX .
To compute variations of the normal, we note that n · e i = 0 and n · n = 1 so δn · e i = −δe i · n, and δn · n = 0, leading to: The variations of the curvature tensor are then, Finally, variations of the Christoffel symbols have the form Expressions (11), (17), and (18) can be rewritten in terms of infinitesimal tangential displacement δX i and normal displacement δX n by substituting δX = δX i e i + δX n n. This leads to the expressions [1] which can be used to obtain the virtual work expression Eq. (34). When the surface changes according to a flow field v, and coordinates follow the flow in a Lagrangian way, the surface shape change can be written δX = vdt. Therefore, using Eq. (19), the rate of change of the components of the metric tensor is: where v ij is the strain-rate tensor, defined in Eq. (4) in the main text. Eq. (22) is the definition for the components of the Lie derivative of g along the flow v [3].

Balance of linear and angular momentum and the principle of virtual work
Here we show that the condition of vanishing infinitesimal virtual work, Eq. (2) in the main text, follows from the force balance and torque balance equations on the surface S, at low Reynolds number where inertial terms can be neglected. We consider a more general case than in the main text, where the surface subjected to an external force density f and an external torque density Γ. Internal forces and torques are described by the tension tensor t i = t ij e j + t i n n and the moment tensor m i = m ij e j + m i n n, such that: are respectively the force and torque acting on a cut on the surface with length dl, with ν the unit vector tangent to the surface and normal to the cut. The statement of conservation of linear momentum (or force balance) for a thin sheet, at low Reynolds number where inertial terms are neglected, is [1] where the first equation stems from balance of forces tangent to S and the second from balance of forces normal to S. On the other hand, balance of angular momentum (or torque balance) tangent and normal to S reads In the following, for convenience we introduce the symmetric and antisymmetric part of the tension tensor t ij : where t ij S = t ji S and t ij A = −t ji A . Eqs. (26) and (27) can be rewritten as: We now multiply Eq. (24) by an infinitesimal tangential displacement field −δX j and integrate the result over S, to obtain: where we have used Eqs. (29) and (30) and integrated by parts terms taking into account that S is a closed surface. Multiplying Eq. (25) by an infinitesimal normal displacement field −δX n and again integrating over S, we obtain: Adding Eqs. (31) and (32), we get Finally, using the analytical expressions for the variations of the metric, curvature, and Christoffel symbols δg ij , δC ij , and δΓ k ij in terms of tangential and normal variations (Eqs. (19)-(21)), we get where we have defined ∇ × δX = 2 ij (∂ j δX n − C jk δX k )e i + (∇ i δX j ) ij n (see Refs. [1,2]),m ij = −m ik j k , and the tension tensort which accounts for the work generated by the in-plane tension tensor t ij as well as the work generated by the bending moments m ij upon an in-plane deformation characterised by δg ij . The expression for the differential virtual work, Eq. (2) in the main text, corresponds to Eq. (34) with vanishing external torque density (Γ = 0) and vanishing normal part of the moment tensor (m i n = 0).

Equilibrium tension and bending moments in the Helfrich model
Here we derive the tensiont ij e and bending momentm ij e arising from a Helfrich energy of the form To identifyt ij e andm ij e , we take variations of Eq. (36), where we have used that dS = √ gds 1 ds 2 with g the determinant of the metric. The differential of the trace of the curvature tensor is: where we have used that δg ij = −g ik g jl δg kl . Noting that the variation of the square root of the metric determinant is given by Eq. (14), we can then write: which, by comparing with Eq. (34) and using that δW = δF at equilibrium, allows us to identifŷ which contribute to the total tension and bending moments in the constitutive equations, Eq. (3) in the main text.

Microscopic motivation for the cell-cell adhesion potential
Here we discuss a microscopic motivation for the characterisation of cell-cell interactions with a potential. Following Eq. 17 in the main text, we consider the free energy of an ensemble of linkers, described by a twopoint concentration field c IJ (X I , X J ) and concentration c I denotes the concentration of unbound, free linkers in cell I: where parameter definitions are discussed in the main text. Variation of the free energy (41) with respect to the concentrations gives: To obtain chemical equilibrium, we note that concentrations can change due to the reaction of binding of two linkers on a pair of surfaces I, J with reaction coordinate r IJ = r JI , and can change due to exchange of free linkers with the bulk, with reaction coordinates r 0I : such that minimisation with respect to δr IJ for each pair of surfaces I, J gives: and minimisation with respect to δr 0I gives c I = c 0 . One then obtains: is the binding affinity of linkers for two element of surface separated by the distance Taking into account that a surface deformation dilutes locally the density of linkers as δc I = − 1 2 c I g ij I δg I,ij and δc IJ = −c IJ 1 2 g ij I δg I,ij + 1 2 g ij J δg J,ij , the variation of the free energy (41) with respect to a shape change leads to (48) We now assume that linkers can bind and unbind the surfaces quickly so that c I and c IJ are equal to their equilibrium values. The contribution to the first sum in Eq. (47) is then equivalent to a homogeneous surface tension, which we assume can be absorbed in the surface tension term in Eq. 3 in the main text. The contribution in the second sum of Eq. (47) can be rewritten using the equilibrium concentration (46): Thus, we can write an effective free energy with In the main text we consider instead an effective Morse potential of the form which like Eq. (51), has a minimum −D at the position |X I − X J | = r min . The effective potential defined in Eq. (52) also has a strong repulsive component at short separation distance for l r min .
Variations of Eq. (50) lead to Eq. (24) in the main text, using Eq. (14). In addition, we also show below that variation of the effective free energy F IJ with respect to the surface shapes I, J, only depends on normal infinitesimal displacement of the surfaces, and so only depends on the surface shapes: where we have used that: and we generalise notations of differential geometry (S1 Appendix 1) to surfaces I, J, ... by indicating with a I, J, ... subscript the surface considered. The variation δF IJ in Eq. (53) only depends on δX In and δX Jn , which shows that the net tangential force (accounting for both the internal tension and the external force density) from the interaction potential is zero.

Discretised finite element equations
We discuss here the form of the finite element equations that we solve numerically. We consider a parameterisation of the shape S given in Eq. (8) in the main text. For the discretised surface, where notations are defined as in Eq.(8) in the main text. Using this decomposition and Eqs. (11), (17), and (18) for the variations of the metric and curvature tensors and Christoffel symbols, we have where we have defined and here

Active viscous layer with bending rigidity
We start by looking at the discretisation of an active viscous layer with bending rigidity. To obtain the equations governing the dynamics of the nodes of the mesh, we substitute Eq. (12) into Eq. (2) from the main text. Using the expressions in Eqs. (56) and (57), we get Note that because of the compact nature of the basis functions B a , the force F a only depend on nodes that have non-zero basis functions in the neighbouring elements of node a, which are those formed by the first, second, and third nearest neighbours in the mesh (those connected to node a by a 3-edge, 2-edge, or 1-edge path in the mesh) and we denote by a . As it will be apparent in the next section, these equations also depend on P (n+1) , the intracellular pressure difference. We split the different terms leading to F a individually for clarity: These terms are We note here that we integrate the dissipative terms Eqs. (65) and (66) on S (n) , whereas Eqs. (67) and (68) are integrated on S (n+1) ; this guarantees that, when γ is homogeneous and time-independent, the effective energy F = S (γ + κ(C k k − C 0 ) 2 /2)dS decreases over time, which endows the discretisation with stability regardless of the value of ∆t (n) , see [3]. The term F a,volume is discussed in the next section. To calculate the integrals in Eqs. (65)-(68), and since our surface is parametrised locally by finite element parametrisations of the form of Eq. (9) in the main text, we split integrals into a sum over integrals in each element. Integrals over elements are then transformed into integrals in the parametric domain of barycentric coordinates, i.e. given an integrand f (X): where s 1 e , s 2 e denote the barycentric coordinates of element e and S e the surface spanned by the element e. Using Gaussian integration, we approximate this integral by where (s 1 ek , s 2 ek ) is the parametric position of Gauss point k and w k its associated weight, on element e. In practice, we loop over elements e and over their Gauss points k, and calculate the partial sums to the expressions Eq. (65)-(68), which are then assembled following a typical finite element routine. To solve Eq. (63), which is a non-linear equation of the degrees of freedom X (n+1) a , we use a Newton-Raphson method, which requires to solve Eq. (15) of the main text. For this method, we need to compute the matrix ∂F a /∂X . To obtain this matrix, we split the expressions as in Eq. (65)-(68), ∂F a,active ∂F a,Helfrich where I is the identity matrix in three dimensions, and where k e kα n β + e iγ e jα n β + e iα e jβ n γ = (∂ i B a )(∂ j B b ) e iα e jβ n γ + e jα n β e iγ − g ij n α n γ n β . (77)

Adding cell pressure as a Lagrange multiplier
To add the cell pressure P imposing the constraint in volume, we first note that the volume enclosed by S can be computed as a surface integral by applying the divergence theorem, where V denotes the volume enclosed by S. We note that the variation of volume can be written, by application of Reynold's transport theorem, as Thus, by plugging the definition of the pressure force, Eq. (6) in the main text, into the principle of virtual work, Eq. (2) in the main text, we can write The right hand side can then be interpreted as the second term in the variation of the Lagrangian −δ(P (V − V 0 )) = −δP (V −V 0 )−P δV , with V 0 the cell target volume. This allows us to enforce the constraint V (n+1) = V 0 through P = P (n+1) . After discretisation, the right hand side of Eq. (80) leads to the following contribution to Eq. (64): The constraint of prescribed volume can then be written as which adds an extra equation to be solved. For the Newton-Raphson method, we need to compute the derivatives ∂F a,volume (84)

Adding cell-cell interactions
We now consider cell-cell interactions. Here, one needs to consider a larger system of equations, taking into account vertex forces calculated in the previous sections, but also force stemming from cell-cell interaction. For cell I and node a, the condition δW = 0 defined in Eq. (23) in the main text now gives Here I, a identifies the set of vertices (identified as the pair of labels J for the cell considered and b for the vertex considered) that are neighbours of vertex a in cell I. These contain all vertices that are neighbours to a in cell I and also all vertices that interact with node a through the interaction potential. The function F I,a then has the following decomposition: F I,a = F I,a,friction + F I,a,viscous + F I,a,active + F I,a,Helfrich + F I,a,volume + F I,a,adhesion .
The 5 first terms in the right-hand side follow definitions obtained in the previous subsections. The remaining term, the contributions from adhesion, reads, using Eq. (24) in the main text: .
We note that, together with Eq. (85), for each cell I we solve F I,P = V 0 with the expression in Eq. (82) to impose the volume contraint through P a )/∆t [n] where [n] identifies the n-th step of the mesh relaxation dynamics and solve Eqs. (30) and (31) of the main text. We let the time-step in this process to increase as ∆t [n] ← g∆t [n] with g > 1 a factor (g = 1.1 in our simulations) if (1) the number of iterations in the Newton-Raphson method employed to solve Eqs. (30) and (31) of the main text is smaller than a maximum number (5 in our simulations) and (2) the error in the normal motion, evaluated as S (v · n) 2 dS/A is smaller than a tolerance (10 −3 in our simulations). If conditions (1) or (2) are not satisfied after a step in the dynamics, we repeat that step by reducing the time-step ∆t [n] ← ∆t [n] /g until conditions (1) and (2) are satisfied. We stop this mesh improvement dynamics when |∂F mesh /∂X α a | ∞ = max a (|∂F mesh /∂X α a |) is smaller than a given tolerance (10 −3 in our simulations). In S1 video we show a prototypic dynamics of mesh relaxation driven by this method (left), where we show that the distance between the original and the reparametrised smooth surfaces is below 5 · 10 −3 (to be compared with the radius of the surface 1); this error can be controlled by changing the tolerance for maximum normal motion (right, where tolerance for normal motion is 10 −4 ), which however leads to a larger number of steps in the reparametrisation. We note that in simulations with many cells, this reparametrisation can be applied in parallel to improve the mesh of each cell individually and, as the shape is mantained to a good degree of precision, meshes do not self-intersect in the course of remeshing.

Analytical solution for the flow on a spherical active surface
Here we derive an analytical expression for the flow on a spherical viscous surface resulting from gradient of active tensions, in the absence external torques. Our analysis is similar to the derivation in Ref. [4].
We separate here the velocity into its in-plane v i and normal v n components. We further consider the Hodge decomposition of the tangential velocity [5,6], where φ and ψ are the irrotational and solenoidal potentials, i.e. ∂ i φ is an irrotational field whereas j i ∂ j ψ is a solenoidal field. In this decomposition, we assume the surface is simply-connected as otherwise there would be another, harmonic component of v i [5,6]. Using the constitutive equations (12) in the main text, the tension tensor components are: where we have considered that the bending modulus κ and preferred curvature C 0 are uniform on the surface, while γ can be inhomogeneous. The external force density is given by Eq. (5) in the main text, and arises from an external pressure and an effective external friction force density. The in-plane and normal force balance equations (24) and (25) can then be written as: where ∆ = ∇ i ∇ i . We now introduce three arbitrary scalar fields w φ , w ψ and w n . We multiply Eq. (93) by ∂ i w φ and by j i ∂ j w ψ , and Eq. (94) by w n and integrate the resulting equations on S. We obtain: where we have integrated by parts some of the terms. Taking into account that ∇ i ∂ j A = ∇ j ∂ i A with A a scalar and that ∇ i ∇ j B k = ∇ j ∇ i B k − R l kij B l with B a vector and the Riemann tensor on a surface of Gaussian curvature K = det(C j i ), and assuming that S is a sphere of radius (C ij = g ij / , K = 1/ 2 ), we obtain: We now expand φ = A a Y Aa φ Aa , ψ = A a Y Aa ψ Aa , v n = A a Y Aa v Aa n , and γ = A a Y Aa γ Aa with Y Aa the spherical harmonic of degree A and order a, with −A ≤ a ≤ A. Defining spherical coordinates θ, φ, such that a point of the sphere of radius is given by X(θ, φ) = (sin θ cos φẽ x + sin θ sin φẽ y + cos θẽ z ) , the spherical harmonics are then given by with P Aa the associated Legendre polynomial of degree A and order a, and Z a normalising constant such that Imposing w φ = w ψ = w n = Y Bb , and taking into account the orthonormality of spherical harmonics, i.e. S Y Aa Y Bb * dS = δ AB δ ab with * indicating the complex conjugate, we obtain the following algebraic equations: where Eq. (107) holds for all A ≥ 0. We note that the zero-th mode A = a = 0 for φ and ψ is arbitrary since uniform values of φ and ψ do not change the velocity field defined by Eq. (90); we choose here that they vanish for simplicity. Given γ Aa and P , the solution is 9 Approximate description of a cell doublet Here we compute an approximate solution for the cell doublet valid in the limit ofr min ,l 1. We assume that the doublet is formed by two mirror spherical caps of height h c and base radius r c separated by a distance d, see Fig 1. Sincel 1, only the two circular discs at the base of each cell interact. We label a point X 1 on the contact disc of cell 1 by polar coordinates r 1 , θ 1 , and we define a point on the contact disc of the second cell by: where the z axis goes along the line joining the cell centers. The distance between two points on the two surfaces is then |X 2 − X 1 | = √ d 2 + r 2 . The interaction energy can be written as: where the second integral is taken over the domain S 2 (r 1 , θ 1 ) corresponding to points on cell 2 which are labelled by r, θ and are within the circular disc. Here we approximate this integral as an integral on the infinite plane; as the interaction potential ϕ is short-ranged, this is valid sufficiently far away from the boundaries of the contact discs. Therefore, the approximation applies when the contact disc is much larger than l, r min . One then obtains: whose minimum is at d * = r min − l ln 2. In the limit l r min , one then has ζ * Morse = −βDlr min with β = 4π. For the potential with a cutoff that we employ in our numerical simulations (Eq. (21) in the main text), the value of ζ * can be computed numerically. Taking r 1 = r min and r 2 = r min + 3l and r min = 3l, one finds ζ * = −βDr min l with β 10.7. The effective tension at the interface at equilibrium is then which becomes negative, and thus we expect to get a buckling instability, when 10 Approximate description of a planar, regular cell packing Here we discuss approximate models for the shape of cells organised within a planar sheet, as in Fig 6A in the main text. We consider that cells have a fixed volume V 0 , a surface tension for free interfaces γ ab , and a surface tension for interfaces in contact with other cells, γ l < γ ab . The effective energy of a single cell can then be written: where A ab is the apical surface area, equal to the basal surface area, A l is the lateral or contact surface area, V the cell volume and P the cellular pressure difference. Treating each cell as an hexagonal prism with side length a and height h c , one has A ab = 3 √ 3a 2 /2 the apical and basal surface area, A l = 6ah c the lateral surface area, and the cell volume is V = 3 √ 3a 2 h c /2. In that case, minimization of F leads to: Alternatively, one can approximate the shape of a cell as a cylinder of radius r c and height h c , connected to two spherical caps with height h s . Then A ab = π(r 2 c + h 2 s ), A l = 2πr c h c and V = πr 2 c h c + 2πh s (3r 2 c + h 2 s )/6. Minimisation of F leads to: One can also envision a "hybrid" model, where the cell is considered as a union of an hexagonal prism of height h c and side length r c , and two spherical caps or height h s and radius r c . We consider an approximate case, in that the extra surface area that arises from joining the prism hexagonal surfaces to the spherical caps circular discs, is neglected. In that case, one has A ab = π(r 2 c + h 2 s ), A l = 6r c h c and the cell volume is given by V = 3 √ 3r 2 c h c /2 + 2πh s (3r 2 c + h 2 s )/6. Minimisation of F leads to: The solution of these 3 models, together with results from simulations, are plotted in Fig 6E. In all cases, one takes γ ab =γ and which is the effective surface tension of an adhering interface, for one of the two participating cells, in the limit where the interface is large compared to l, r min (Eq. (122)).

Details about computational load
Here we provide information regarding computational resources used for the simulations of a planar sheet. Each cell is represented with a mesh with 806 vertices and 1608 elements, corresponding to a value of h/ ≈ 1 · 10 −1 .
The full non-linear system of 36 cells has 173,700 degrees of freedom (173,664 node velocities and 36 pressure Lagrange multipliers). Each cell is stored in a different partition using 4 cores for OpenMP parallelisation (for a total of 144 cores), so each partition handles the 2418 degrees of freedom of its assigned cell. On average, a time-step is solved in 2.5 minutes (on Intel Xeon E5-2670 processors), corresponding to 5 steps of a Newton-Raphson iteration; around 75% of this time is spent in filling the linear system of equations, most of which (95%) is spent in filling the terms corresponding to cell-cell interactions. After the solution for velocities and pressures is found, meshes are updated using the reparametrisation scheme described in Section 2.3, which, on average, requires another 2.5 minutes. To reach the equilibrium state from the initial distribution of spheres, we employ ∼200 time-steps, which correspond to a total of ∼17 h of simulation.