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

Interacting active surfaces: a model for three-dimensional cell aggregates | bioRxiv /* */ /* */ <!-- <!-- /*! * yepnope1.5.4 * (c) WTFPL, GPLv2 */ (function(a,b,c){function d(a){return"[object Function]"==o.call(a)}function e(a){return"string"==typeof a}function f(){}function g(a){return!a||"loaded"==a||"complete"==a||"uninitialized"==a}function h(){var a=p.shift();q=1,a?a.t?m(function(){("c"==a.t?B.injectCss:B.injectJs)(a.s,0,a.a,a.x,a.e,1)},0):(a(),h()):q=0}function i(a,c,d,e,f,i,j){function k(b){if(!o&&g(l.readyState)&&(u.r=o=1,!q&&h(),l.onload=l.onreadystatechange=null,b)){"img"!=a&&m(function(){t.removeChild(l)},50);for(var d in y[c])y[c].hasOwnProperty(d)&&y[c][d].onload()}}var j=j||B.errorTimeout,l=b.createElement(a),o=0,r=0,u={t:d,s:c,e:f,a:i,x:j};1===y[c]&&(r=1,y[c]=[]),"object"==a?l.data=c:(l.src=c,l.type=a),l.width=l.height="0",l.onerror=l.onload=l.onreadystatechange=function(){k.call(this,r)},p.splice(e,0,u),"img"!=a&&(r||2===y[c]?(t.insertBefore(l,s?null:n),m(k,j)):y[c].push(l))}function j(a,b,c,d,f){return q=0,b=b||"j",e(a)?i("c"==b?v:u,a,b,this.i++,c,d,f):(p.splice(this.i++,0,a),1==p.length&&h()),this}function k(){var a=B;return a.loader={load:j,i:0},a}var l=b.documentElement,m=a.setTimeout,n=b.getElementsByTagName("script")[0],o={}.toString,p=[],q=0,r="MozAppearance"in l.style,s=r&&!!b.createRange().compareNode,t=s?l:n.parentNode,l=a.opera&&"[object Opera]"==o.call(a.opera),l=!!b.attachEvent&&!l,u=r?"object":l?"script":"img",v=l?"script":u,w=Array.isArray||function(a){return"[object Array]"==o.call(a)},x=[],y={},z={timeout:function(a,b){return b.length&&(a.timeout=b[0]),a}},A,B;B=function(a){function b(a){var a=a.split("!"),b=x.length,c=a.pop(),d=a.length,c={url:c,origUrl:c,prefixes:a},e,f,g;for(f=0;f<d;f++)g=a[f].split("="),(e=z[g.shift()])&&(c=e(c,g));for(f=0;f<b;f++)c=x[f](c);return c}function g(a,e,f,g,h){var i=b(a),j=i.autoCallback;i.url.split(".").pop().split("?").shift(),i.bypass||(e&&(e=d(e)?e:e[a]||e[g]||e[a.split("/").pop().split("?")[0]]),i.instead?i.instead(a,e,f,g,h):(y[i.url]?i.noexec=!0:y[i.url]=1,f.load(i.url,i.forceCSS||!i.forceJS&&"css"==i.url.split(".").pop().split("?").shift()?"c":c,i.noexec,i.attrs,i.timeout),(d(e)||d(j))&&f.load(function(){k(),e&&e(i.origUrl,h,g),j&&j(i.origUrl,h,g),y[i.url]=2})))}function h(a,b){function c(a,c){if(a){if(e(a))c||(j=function(){var a=[].slice.call(arguments);k.apply(this,a),l()}),g(a,j,b,0,h);else if(Object(a)===a)for(n in m=function(){var b=0,c;for(c in a)a.hasOwnProperty(c)&&b++;return b}(),a)a.hasOwnProperty(n)&&(!c&&!--m&&(d(j)?j=function(){var a=[].slice.call(arguments);k.apply(this,a),l()}:j[n]=function(a){return function(){var b=[].slice.call(arguments);a&&a.apply(this,b),l()}}(k[n])),g(a[n],j,b,n,h))}else!c&&l()}var h=!!a.test,i=a.load||a.both,j=a.callback||f,k=j,l=a.complete||f,m,n;c(h?a.yep:a.nope,!!i),i&&c(i)}var i,j,l=this.yepnope.loader;if(e(a))g(a,0,l,0);else if(w(a))for(i=0;i (function(w,d,s,l,i){w[l]=w[l]||[];w[l].push({'gtm.start':new Date().getTime(),event:'gtm.js'});var f=d.getElementsByTagName(s)[0];var j=d.createElement(s);var dl=l!='dataLayer'?'&l='+l:'';j.src='//www.googletagmanager.com/gtm.js?id='+i+dl;j.type='text/javascript';j.async=true;f.parentNode.insertBefore(j,f);})(window,document,'script','dataLayer','GTM-M677548'); Skip to main content Home About Submit ALERTS / RSS Search for this keyword Advanced Search New Results Interacting active surfaces: a model for three-dimensional cell aggregates View ORCID Profile Alejandro Torres-Sánchez , View ORCID Profile Max Kerr Winter , View ORCID Profile Guillaume Salbreux doi: https://doi.org/10.1101/2022.03.21.484343 Alejandro Torres-Sánchez 1 Theoretical Physics of Biology laboratory, The Francis Crick Institute , 1 Midland Road, NW1 1AT, London, United Kingdom Find this author on Google Scholar Find this author on PubMed Search for this author on this site ORCID record for Alejandro Torres-Sánchez Max Kerr Winter 1 Theoretical Physics of Biology laboratory, The Francis Crick Institute , 1 Midland Road, NW1 1AT, London, United Kingdom Find this author on Google Scholar Find this author on PubMed Search for this author on this site ORCID record for Max Kerr Winter Guillaume Salbreux 1 Theoretical Physics of Biology laboratory, The Francis Crick Institute , 1 Midland Road, NW1 1AT, London, United Kingdom 2 University of Geneva , Quai Ernest Ansermet 30, 1205 Genève, Switzerland Find this author on Google Scholar Find this author on PubMed Search for this author on this site ORCID record for Guillaume Salbreux For correspondence: guillaume.salbreux{at}unige.ch Abstract Full Text Info/History Metrics Supplementary material Preview PDF Abstract 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 constitutive equations using a finite element method, and provide a parallel implementation in C++. We discuss examples of application of this framework to simulations involving 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. Author summary Understanding how tissue-scale morphogenesis arises from cell mechanics and cell-cell interactions is a fundamental question in developmental biology. Here we propose a mathematical and numerical framework to address this question. In this framework, each cell is described as an active surface representing the cell acto-myosin cortex, subjected to flows and shape changes according to active tensions, and to interaction with neighbouring cells in the tissue. Our method describes cellular processes such as cortical flows, cell adhesion, and cell shape changes in a deforming three-dimensional cellular aggregate. To solve the equations numerically, we employ a finite element discretisation, which allows us to solve for flows and cell shape changes with arbitrary resolution. We discuss applications of our framework to describe cell-cell adhesion in doublets, three-dimensional cell shape in a simple epithelium, and three-dimensional growth of a cell aggregate. 1 Introduction Tissue morphogenesis relies on the controlled generation of the cellular forces that collectively drive tissue-scale flows and deformation [ 1 , 2 ]. The interplay between cell-cell adhesion, cellular mechanics and the cytoskeleton plays a key role in determining how biological tissues self-organise [ 3 ]. These ingredients are also crucial for the growth of in vitro organoids, organ-like structures derived from stem cells which can self-organise into complex structures reminiscent of actual organs [ 4 – 6 ]. Several classes of models have been proposed to describe the mechanics of multicellular aggregates, such as cellular Potts models [ 7 – 13 ], phase field models [ 14 – 18 ] and vertex models [ 19 ]. Vertex models, and the closely related Voronoi vertex models [ 20 , 21 ], describe cells in a tissue as polyhedra that share faces, edges and vertices forming a three-dimensional junctional network [ 19 , 22 – 24 ]. Cell deformations are encoded by the displacement of the vertices X a of the network. These displacements are dictated by vertex forces F a stemming from cell pressures P c , surface tensions t f , and line tensions Γ e that are coupled to virtual changes in cell volume δV c , face area δA f , and edge length δl e respectively in a work function where to get to the last expression one needs to express δV c , δA f and δl e in terms of a virtual displacement of the vertices δ X a . Two and three-dimensional versions of the vertex model have been employed for many applications, for instance to study cell packing [ 23 , 25 ], cell sorting [ 26 ], wound closure [ 27 ], cyst formation [ 28 ], tumourigenesis in tubular epithelia [ 29 ] among many other [ 19 ]. However because of their definition, vertex models do not explicitly resolve cortical flows on the cell surface. The effect of cell-cell adhesion is also implicitly introduced in the surface tension t f , which mix together physical processes arising from molecular bonds between cells and surface forces arising in the cell membrane and in the actomyosin cytoskeleton. In vertex models with vertices positions as degrees of freedom, topological transitions leading to cell-neighbour exchange are encoded explicitly by formulating rules to change edges in the network. At the single cell scale, a number of studies have shown the relevance of coarse-grained, continuum models to describe the mechanics of the cell surface. In this approach, an active fluid theory taking into account cellular cortical flows, gradients of active cytoskeletal tension and their regulation, and orientation and filament alignment in the actin cortex, has proven successful to describe the mechanics of cell polarisation, cell motility or cell division [ 30 – 39 ]. From a computational perspective, there has been a growing attention to the simulation of the dynamics of fluid interfaces both with prescribed [ 40 – 43 ] and with time-evolving shape [ 38 , 39 , 44 – 50 ]. However, to our knowledge no computational framework has attempted to provide a physical description of three-dimensional cellular aggregates taking into account explicitly the mechanics of a single cell surface described as an active surface, as well as cell-cell adhesions. Here we bridge this gap and introduce a new modelling and simulation framework, and a freely available code [ 51 ], for the mechanics of cell aggregates in three dimensions ( Fig. 1 ). We describe cells as interacting active surfaces [ 52 ]. The governing equations for the cell surface mechanics are discretised using a finite element method. In this method, each cell is represented by a three-dimensional mesh with vertices positions X a . In analogy with Eq. (1) in vertex models, we start from the virtual work theorem for interfaces ( Eq. (2) ) and find the net forces at the vertices F a ( Eq. (11) ), which vanish in the absence of inertia. This condition imposes cortical flows and cell shape changes. In this framework, topological transitions appear as a natural output of the remodelling of cell-cell interactions and are not treated explicitely. Download figure Open in new tab Fig 1. Schematic of the interacting active surface framework. A tissue is described by a collection of triangular meshes representing each cell. The dynamics of the tissue is described by the dynamics of the vertices making the cellular meshes, similar to how the movement of the vertices of a vertex model describe the deformation of a tissue. The motion of the cell mesh is obtained by coarse-graining continuum mechanics equations of a theory of active surfaces via the finite element method. In this theory, cortical flows, cortical tensions, intracellular pressures, bending moments and forces arising from cell-cell interactions are taken into account. The main benefits of our method are that it (1) can incorporate complex descriptions of the physics of the cell surface, including sources of tension, in-plane and normal moments [ 52 ], (2) accounts for cell-cell adhesion explicitly through constitutive laws that can be adapted to represent different biological scenarios, (3) resolves cell shape with arbitrary resolution given by the mesh size of the finite element discretisation. The nonlinear equations of the model are treated computationally with the use of nonlinear solvers. Since we consider only the discretisation of the cell surface, the number of degrees of freedom is considerably smaller than in 3D Cellular Potts or phase-field models, both of which require a 3D discretisation. On the other hand, to capture cell shape, cortical flows and cell-cell interactions accurately, we need to use many more vertices per cell than in a vertex model, which leads to a greater computational cost. To limit computational time, the code is parallelised to allow each cell to be stored on a different partition, each possibly using a pool of cores (hybrid MPI-OpenMP method). As such, it is possible to simulate several tens of cells on a computing cluster. We now turn to the description of our framework. In Section 2 we describe the mathematical formulation and the discretisation of our method. We show some examples of application in Section 3 and end in Section 4 with conclusions, summary and ideas for future work. 2 Materials and methods In this section we describe the main elements of our framework. We start by introducing the governing equations for a single, isolated cell described as a fluid active surface. We use a virtual work formulation for the mechanics of a surface, together with a finite element discretisation, to obtain the equations that dictate the movement of the vertices of the mesh. We introduce cell-cell interactions, represented by an interaction potential between pairs of surfaces, that result in a force density and a tension acting on each cell. Finally, we describe a mesh reparametrisation method that allows simulations to handle large tangential deformations of the surface which can arise for continuously flowing fluid interfaces such as the cell cortex. 2.1 Governing equations and discretisation for a single surface 2.1.1 Virtual work for interfaces Our starting point is the statement of virtual work for a closed interface representing a single cell. We describe with a parametrisation X ( s 1 , s 2 ) where s 1 , s 2 are two surface coordinates and X a point in the 3D space in which the surface is embedded. We denote the coordinates of the 3D cartesian basis by greek indices α, β ,…, and the coordinates on the surface by latin indices i,j …. Here and elsewhere in the manuscript we use Einstein summation convention for repeated indices. Given the tangent vectors e i = ∂ i X one can compute the metric tensor g ij = e i · e j and the curvature tensor C ij = – n · ∂ i e j , where n = ( e 1 × e 2 )/| e 1 × e 2 | is the outer normal to the surface. Other notations of differential geometry are given in S1 Appendix 1. Here and in the following, we consider the limit of low Reynolds number where inertial terms are negligible, a limit relevant to the physics at the cell scale which is of interest here [ 53 ]. The mechanics of a single surface can then be described by the following statement of the principle of virtual work for [ 52 ]: where δ X is an infinitesimal displacement of the surface, and δg ij and δC ij the associated infinitesimal variation of the metric and curvature tensors. Expressions for δg ij and δC ij in terms of δX α are given in S1 Appendix 2. Eq. (2) relates infinitesimal variations of geometric quantities of the interface to their work conjugates: the external force density f α coupled to the surface displacement, the tension tensor related to metric variations, and the bending moment tensor coupled to variations of the curvature tensor. Eq. (2) is equivalent to the statement of balance of linear and angular momentum at low Reynolds number (S1 Appendix 3). We have neglected external and internal normal moments for simplicity, which lead to extra terms in Eq. (2) (S1 Appendix 3). 2.1.2 Constitutive laws To describe the mechanics of a single surface, we need to specify constitutive equations for the mechanical tensors and . We denote by v the velocity field on the cell surface. Here, we assume that the cell surface can be represented as an active viscous layer, with a bending rigidity and a spontaneous curvature. As a result, we use the following constitutive equations: where v ij is the strain-rate tensor on the surface: with ∇ i the covariant derivative operator defined in S1 Appendix 1. Here, η is the surface viscosity; for simplicity we do not explicitly distinguish between the shear and bulk viscosity of the surface. γ is the surface tension of the cell, which is not necessarily homogeneous on the cell surface. As we expect this contribution to largely arise from active forces generated in the actomyosin cortex [ 54 ], we later refer to it as the active tension. The contributions and arise from an effective Helfrich free energy penalising the membrane curvature with bending modulus κ and spontaneous curvature C 0 (S1 Appendix 4). The external force density acting on a single cell is split into two contributions arising respectively from the action of the pressure difference across the surface and from an effective external friction force density where ξ is a friction coefficient. In the following, the pressure P is adjusted to impose a value of the cell volume. Eq. (2) , together with the constitutive equations (3) - (7) , leads to a complete set of equations to determine the velocity field v on the cell surface. The constitutive laws Eqs. (3) - (7) are a simple choice of physical description of the surface. Other terms, for instance additional viscous or active bending moments, could be introduced: linear irreversible thermodynamics provides with a set of additional possible terms that could play a role in the dynamics of isotropic, polar or nematic active surfaces and could be added to the framework described here [ 52 , 55 ]. 2.1.3 Finite element discretisation We now discretise the geometry of the surface using finite elements based on a triangular control mesh of N e triangles or elements and N v vertices ( Fig. 2 ). The parametrisation of is of the form where the sum is taken over vertices of the control mesh, labelled by a , X a is the position of vertex a on the control mesh, and is its basis function, defined in terms of coordinates ( s 1 , s 2 ) on the control mesh. Given that may have bending or in-plane moments coupled to variations of the curvature tensor and the Christoffel symbols in the differential virtual work (see Eq. (2) and Eq. (34) in S1 Appendix 3), the basis functions B a must be chosen to have second order derivatives that are square-integrable; here we follow [ 38 , 56 , 57 ] and use Loop subdivision surfaces, which lead to smooth surfaces that satisfy this condition by construction. We note that, in practice, a consistent parametrisation ( s 1 , s 2 ) of the entire control mesh is not needed: indeed one can consider the surface as a union of surface elements associated to each triangle e of the control mesh, which are given in terms of barycentric coordinates , of the element e by see Fig. 2 . Here we have denoted by 〈 e 〉 the set of vertices whose basis functions are nonzero in element e , and corresponds to the contribution of the basis function B a , within element e , parametrised by the barycentric coordinates of e . For Loop subdivision surfaces, 〈 e 〉 is formed by the vertices in the triangular element e , and all first neighbours of these vertices. Download figure Open in new tab Fig 2. A smooth surface representing a cell is obtained based on a triangular control mesh with vertex positions X a (A) and a set of basis functions per vertex a (B). (A) To define the mapping between the control mesh and the cell surface, the barycentric coordinates of points in a triangular element e in the control mesh , which span a reference triangle, are used to define a point on the cell surface ( Eq. (9) ). Points are obtained by summing basis functions , weighted by X a , over vertices a whose basis functions have a non-zero contribution to this element, an ensemble denoted 〈 e 〉. (B) Example of the basis function associated to a vertex in the mesh. For Loop subdivision surfaces basis functions, the basis function spans the first and second rows of elements surrounding the vertex (thicker white line). The vertices that interact with vertex a in the same cell, represented by the set 〈〈 a 〉〉 (green) are formed by the first, second, and third nearest neighbours in the mesh. The discretisation of the surface ( Eq. (8) ) transforms the virtual work principle ( Eq. (2) ) into an expression of the form Here F a , which can be interpreted as the net force on vertex a , can be obtained by substituting, in the differential virtual work, the analytical expressions for the variations δg ij and δC ij in terms of δ X = ∑ a B a ( s 1 , s 2 ) δ X a (S1 Appendix 6). We have denoted by 〈〈 a 〉〉 the set of vertices interacting with vertex a , which for Loop subdivision surfaces is formed by its first, second and third ring of neighbours, see Fig. 2B . Since Eq. (10) has to be satisfied for any δ X a , the net forces on the vertices need to vanish Through this discretisation, we transform the original continuum problem into a set of coupled ordinary differential equations (ODEs). These ODEs need to be discretised in time to be resolved computationally; in the following we denote by ( n ) the n –th time step of the time evolution. Here we employ a semi-implicit Euler discretisation, where terms arising from the effective bending energy and active tension are discretised in a fully implicit manner, whereas viscous and frictional terms are treated explicitly. This particular choice leads to a variational time-integrator that preserves the dissipative structure of the dynamics for a homogeneous and time-independent active tension [ 38 ]. This implies that the tension and bending moment tensors at step n are evaluated as: which can be written in terms of and through the relation between the surface definition and the position of vertices, Eq. (8) . To discretise the strain rate tensor, we have used its relation with the rate of change of components of the metric tensor (see S1 Appendix 2). Plugging these expressions in the differential virtual work Eq. (2) allows us to transform Eq. (11) into a set of (nonlinear) algebraic equations (S1 Appendix 6) which can be solved using a Newton-Raphson method together with the discretisation of the volume constraint, which is imposed through the nonlinear constraint Here P ( n +1) is the intracellular pressure, playing the role of a Lagrange multiplier imposing conservation of volume, see S1 Appendix 6.2. In this method, the solution is updated by , P ( n +1) ← P ( n +1) + Δ P ( n +1) where Δ X a and Δ P ( n +1) satisfy the linearised equations until the norm of F a and F P – V 0 is below a given tolerance. Here , and form the tangent matrix. The tangent matrix is symmetric, in particular . Because only vertices in 〈〈 a 〉〉 interact with a, this matrix is sparse and the linear system can be solved efficiently with an iterative solver. Our computational framework makes use of Trilinos [ 58 ] to handle all linear algebra objects, including sparse matrices, in parallel. 2.2 Forces arising from cell-cell interactions Cell-cell adhesion modifies the previous equations by introducing interactions between the vertices of interacting surfaces. We now consider a set of surfaces , I =1… N describing N interacting cells. We assume that a pair of distinct cells I, J interact via an effective energy: where φ is an effective adhesion potential. We discuss below a microscopic motivation for this potential. The virtual work differential of the cell aggregate can then be written as: where the contribution δW I is given by Eq. (2) with infinitesimal displacement ( δ X I ), metric variation ( δg I ij ) and curvature variation ( δC I ij ) of cell I , and with the tension tensor , bending moment tensor , and external force density excluding cell-cell interaction forces f , given by the constitutive equations (3) - (7) . Here and in the following, we denote distinct cell pairs 〈 I, J 〉, such that in sums taken over 〈 I, J 〉, each pair is counted only once. Variations of the interaction free energy lead to Comparing this expression with Eq. (2) , each cell interaction with a cell J is contributing an additional external force density on cell I , f I J , and an additional isotropic tension to cell I , : In addition, although this point is not directly apparent from Eq. (18) , the variation can be written only in terms of normal displacements, provided that φ only depends on | X I – X J | (Eq. (53) in S1 Appendix 3). This shows that the net driving force from the interaction potential, including the effect of both the force density f I J and tension , has a vanishing tangential component. Intuitively, the interaction energy does not change if cells do not change shape, so it cannot generate a driving force for tangential motion. Following the same procedure as in the previous section, the condition δW = 0, with δW defined in Eq. (17) , gives rise to a set of (nonlinear) algebraic equations of the form 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 interact with vertex a in cell I . The form of the effective energy ( Eq. (16) ) can be motivated microscopically by considering an ensemble of stretchable linkers connecting pairs of surfaces, which quickly equilibrate by binding and unbinding to cell surfaces, and whose free concentration is set by contact with a reservoir. We characterise such an ensemble by a two-point concentration field c IJ ( X I , X J ), which quantifies the number of bound linkers joining the points X I and X J per unit area of the first and second surfaces (thus, it has units of the inverse of an area squared). The concentration c I denotes the concentration of unbound, free linkers in cell I , with units of the inverse of an area. In the dilute limit, the free energy of this ensemble can be written as: where k B is the Boltzmann constant and T the temperature. The first sum corresponds to the free energy of free linkers, described as an ideal solution in contact with a reservoir imposing a chemical potential. This chemical potential determines the value of c 0 . The second sum corresponds to the free energy of bound linkers, described as an ideal solution, with an energy per linker dependent on the linker elongation, quantified by the potential ϕ . The potential is defined such that the force sustained by a linker of length r is – ϕ ′( r ). One can then show that if linkers are at equilibrium with respect to surfaces I, J with a fixed shape, the free energy of interaction of two surfaces I, J is (S1 Appendix 5) which gives a relation between the microscopic behaviour of the linkers and the potential φ introduced in Eq. (16) . For the particular choice of ϕ ( r ) = ϕ 0 + k ( r – r min ) 2 /2 with k the bond stiffness and r min its reference length, one obtains: which corresponds to an inverted Gaussian with centre at the equilibrium length r min , width and depth . This description however still does not take into account short range repulsion between two cells surfaces. This could be taken into account by introducing a second repulsive interaction potential between surfaces. Here we choose instead to introduce a convenient effective potential of interaction, the Morse potential which like the interaction potential in Eq. (23) , vanishes for r → ∞, has a minimum at r = r min with minimum value – D , and is also peaked around its minimum with characteristic length l . In addition, for l ≪ r min it exhibits a sharp short-range repulsion. Although φ Morse ( r ) decays with r rapidly, it is convenient to have a strict cut-off on its range, to limit interacting vertices of the meshes. Therefore, we further multiply φ Morse ( r ) by a smooth step function: where This ensures that the potential φ ( r ) goes to zero exactly at a distance r min + 3 l with first and second order continuous derivatives; a convenient property to solve numerically the non-linear equations ( Eq. (21) ). 2.3 Surface reparametrisation The method discussed in the previous sections is based on a Lagrangian scheme, such that a node of the mesh flows with the material particles of the interface. This, however, can lead to large in-plane distortions, notably if surface tension gradients generate in-plane flows [ 38 , 47 ]. To compensate for the resulting mesh distortion, we describe here a reparametrisation method, in the spirit of Refs. [ 56 , 57 ]. In this method, vertices of the control mesh move tangentially to the surface to minimise an effective mesh quality energy. We stress that this step does not bear any physical meaning. Here we discuss again a single cell. Given the mesh of a cell, we define the energy where the sum is performed over the triangles of the control mesh, and where A e is the area of the triangle e and l e 1 , l e 2 , l e 3 its side lengths. I e and J e represent the invariants (trace and square root of the determinant) of the Cauchy-Green deformation tensor [ 59 ] assuming a reference equilateral triangle of size 〈 A 〉 ∑ e A e /N e , where the sum is taken over the N e triangles of a meshed surface. To specify the free energy f , we use the Neohookean energy f ( I e , J e ) = μI e + λ( J e – 1) 2 [ 59 ]. Note that we define this energy on the mesh rather than on . We want to minimise Eq. (28) , but with the restriction that the cell shape does not change, so that this operation corresponds to a surface reparametrisation without surface deformation. For this, we evolve the position of the vertices of the mesh according to velocities . These velocities are obtained by introducing a continuous velocity field on the surface , and by solving the following equations for the vertices velocities: where p ( s 1 , s 2 ) = ∑ a p a B a ( s 1 , s 2 ) plays the role of a normal pressure, here a field enforcing the condition that the normal flow vanishes in a weak sense, Eq. (31) . This leads to the relaxation of the free energy , with vertices constrained to the shape of , where the constraint is enforced weakly, i.e. in a finite element sense. In practice, this reparametrisation step is performed after a number of steps of the physical evolution of the cell surfaces, and we stop this mesh improvement dynamics when is smaller than a given tolerance. 3 Results We now discuss applications of the interacting active surface framework. We first consider flows in a single spherical cell driven by gradient of cortical tension, a set-up which allows to compare simulation results to an analytical solution. We then examine the shape of the simplest multicellular aggregate, a doublet formed by two cells, when the two participating cells have equal or different tensions. Next, we consider an aggregate of cells assembled in a planar configuration, recapitulating the organisation of a small epithelial island. Finally, we introduce cell divisions in our framework and simulate the growth of a three-dimensional cell aggregate from a single cell. In the following we introduce a reference length scale which corresponds to cell radius of a spherical cell, a reference surface tension , and a reference time scale , which corresponds to the characteristic time scale of cortical flows. We use these reference quantities for normalisation of other quantities. 3.1 Single cell: convergence with mesh size We first consider flows driven by gradients of active tension in a single cell. This allows us to test the convergence of the numerical method for the dynamics of a single spherical cell, since we can compare the velocity field resulting from the method discretisation to an analytical solution using spherical harmonics. We consider a pattern of surface tension on a spherical surface of radius ℓ , given by where Y Aa is the spherical harmonic of degree A and order a . The volume enclosed by the surface is assumed to be subjected to a uniform pressure difference P . The resulting velocity field can then be written as (S1 Appendix 7) where the functions ϕ and v n can also be expanded in spherical harmonics, , and the coefficients read: In Fig. 3 , we consider flows resulting from a surface tension profile given by the coefficients if 1 4, and from the imposed inner pressure ( Fig. 3A ). We obtain the velocity field analytically ( v * , Fig. 3B ) and numerically ( v , Fig. 3C ), for different mesh sizes. We then compute the L 2 norm of the error v – v * , i.e. , with respect to the L 2 norm of v * , . We find an excellent agreement between the exact and numerically obtained velocity field ( Fig. 3B-D ). The corresponding error scales with ( h/ℓ ) 2 , where h is the average mesh size, in line with the reported convergence rate of subdivision surfaces for other systems of partial differential equations [ 60 ] ( Fig. 3D ). Download figure Open in new tab Fig 3. Single cell dynamics and convergence of the numerical framework. (A) An inhomogeneous pattern of active surface tension γ is imposed on a spherical surface. (B) Analytically computed velocity field generated by the active tension profile in A. The velocity field is decomposed into its normal (colormap) and tangential (arrows) components. (C) Numerical solution for the velocity profile generated by the active tension profile in A. (D) Discretisation error, evaluated here in terms of the L 2 norm of the difference between the analytical and numerical solutions for the velocity, as a function of the average mesh size h and for (blue), (orange) and (green). Other parameters are ξℓ 2 / η = 4, C 0 ℓ = 1. 3.2 Shape and dynamics of an adhering cell doublet We now discuss the equilibrium shape of an adhering cell doublet. In this and the following sections, the cell pressure difference P I is imposed as a Lagrange multiplier enforcing the condition V I = V 0 , with V I the cell volume and a reference volume, and we assume that C 0 = 0. With these choices, 5 normalised, non-dimensional parameters have to be specified for each simulation: , and . In the following we set ξ = 10 -3 so that the effect of friction is small. A cadherin bond has a typical length ~15-30 nm [ 61 ] and a typical actomyosin cortex thickness is 200nm [ 62 ], both much smaller than the typical radius of a cell, ~ from a few to tens of μm. Therefore, in the interaction of potential of cell surface we take . For simplicity, in the following we constrain . We typically choose values which for a cell radius of 5 μ m, correspond to r min = 300nm and l = 100nm. We note that if the reference cell volume V 0 is constant, the system can be viewed as a generalised gradient flow minimising the net free energy , subjected to the constraint V = V 0 , where the first term represents an effective energy for the active tension γ and is defined in Eq. (16) . We thus expect the system to eventually reach an equilibrium state with vanishing cortical flows. We first analyse the behaviour of a doublet of identical cells. We initialise the simulation by putting two spherical cells close to each other, such that they are within the interaction range of the potential φ ( r ) ( Eq. (26) ) without touching. As expected, after an initial transient and contact growth, the doublet reaches an equilibrium shape ( Fig. 4A-B ). Increasing the relative adhesion strength leads to an increasing adhesion patch and a lower cell pressure ( Fig. 4C-D and S2 Fig. B). The value of the normalised distance modulates the distance between the two cells ( Fig. 4E ). For beyond a threshold which depends on and , the adhesion patch develops a buckling instability. We found that this instability eventually leads to self-intersection of the computational mesh, as there is no energy contribution in our framework preventing such self-intersection (S2 Fig. A). Varying the bending modulus , but maintaining relatively small values , we find slightly smoother shapes at the boundary of the adhesion patch as well as a slight modulation of the threshold for buckling instability of the contact zone, with larger leading to higher stability ( Figs. 4C,E ). Download figure Open in new tab Fig 4. Dynamics and steady-state shape of a cell doublet. (A) Normalised area of contact as a function of time . (B) Snapshots of deforming doublet at times indicated in (A), with the cell surface velocity field superimposed. (C) Coloured lines: ratio of contact surface area A contact to cell surface area A , for simulations with different values of and . Black dotted line: theoretical approximation valid in the limit of . (D) Snapshots of doublet equilibrium shape for increasing adhesion strength; different parameters in (D1)-(D3) correspond to points labelled in C. (E) Comparison of a slice for the different simulations, with values of and indicate in C, and . The value of affects the shape smoothness of the edge of the adhesion patch. (F) Convergence of the method evaluated by computing the inner cell pressure P for different average mesh sizes h , and comparing the results with a simulation with h/∓ ≈ 2 · 10 -2 (finer). For each h , we compute a box plot using different values of and fixed . (G) Schematic of adhering doublet, with different active tensions γ 1 and γ 2 for each cell. (H) Position of the cell centre of mass X 1 and X 2 , as a function of active tension asymmetry between the two adhering cells, α = ( γ 1 – γ 2 )/( γ 1 + γ 2 ), where γ 1 and γ 2 are the surface tensions of the two cells. Beyond α = 0.69, the cell with lowest tension completely engulfs the one with highest tension. (I) Snapshots of doublet equilibrium shape, clipped by a plane passing by the line joining the cell centres, for increasing difference of active surface tension; corresponding to points labelled in H. In H, I: . To interpret these results, we turn to an approximate analysis of the shape of a doublet. Assuming a small bending rigidity , we approximate the cell doublet by two spherical caps of height h c and base of radius r c , forming an adhesion patch where the cells are separated by a distance d ( Fig. 4E ). The effective free energy of such a doublet configuration can then be written as where is the volume of one cell, P is the cell pressure, and ζ ( d ) an effective surface tension at the contact arising from cell-cell adhesion. For a sufficiently large patch compared to the interaction distances r min , and l , the effective surface tension can be approximated as (S1 Appendix 8): which can be evaluated numerically for a given potential φ . Minimising the effective free energy ( 36 ), one obtains equilibrium values , d * , P * (S1 Appendix 8), which depend on the effective surface tension of a cell at the contact: where the surface tension of the whole interface is 2 γ l . Here, we have introduced a numerical function , whose functional form depends on the potential φ ( r ). For the value of r min / l chosen here, β ≃ 10.7. When the net tension at the contact 2 γ l becomes negative, for we expect the system to develop a buckling instability. Indeed, the corresponding threshold for buckling is well predicted by the simulation with smallest l/ℓ and bending modulus κ ( Fig. 4C ). Before the buckling instability, the ratio between the contact area and the cell surface area as well as the cell pressure P * are well predicted by the approximate analysis, which become more accurate as ( Fig. 4C and S2 Fig. B). To further check the numerical method, we analyse the convergence of the intracellular pressure P as a function of the mesh size h , by comparing the pressure value P obtained at different values of h and (for a fixed ) with a simulation with a fine mesh h/ℓ ≈ 2 · 10 -2 ( Fig. 4F ). We observe that, on average, errors converge as ~ h 3 . Finally, we examine an asymmetric doublet system where cells have different tensions γ 1 and γ 2 ( Fig. 4G-I ). The corresponding equilibrium state has been considered previously (Ref. [ 63 ] and references therein). Fixing a relatively low value of , we change the ratio α = ( γ 1 – γ 2 )/( γ 1 + γ 2 ). As expected from previous studies, we observe that the cell with lowest tension progressively engulfs the cell with highest tension as the ratio α is increased. As a result, their centre of masses approach each other for increasing values of α ( Fig. 4H-I ). In our numerical simulations, beyond α ≃ 0.69, the cell with lowest tension self-intersects before completely engulfing the one with highest tension (S2 Fig. C). 3.3 Epithelial monolayer We now discuss simulations of aggregates containing a larger number of cells. We start by considering a sheet of cells with free boundary conditions, intended to resemble the organisation of a simple epithelial island. To obtain an initial condition for such a sheet, spherical cells with homogeneous active tension are positioned with their centres on a 6×6 hexagonal lattice with side s = 2 ℓ + r min + l . We then let the system relax to its equilibrium state, varying the adhesion parameter , keeping its value below the instability threshold identified in the doublet analysis. As the adhesion parameter is increased, the cellular shapes progressively deviate from loosely adhering spheres to packed columnar cells ( Fig. 5A-B ). The reduced volume with V the cell volume and A its surface area, a measure of the deviation of the cell shape from a sphere with v = 1 being a sphere, is close to 1 for small and then decreases ( Fig. 5B ). Download figure Open in new tab Fig 5. (A) Simulation results (equilibrium shapes) for a planar sheet for different values of the adhesion parameter . (B) Reduced volume v as a function of . (C) Schematic for the measured apical and basal cell surface area A ab , lateral surface area A l , side length a and cell thickness h c . (D) Results are compared to a simple 3D vertex model with a lateral surface tension γ l and apical and basal surface tension γ ab . (E) Box plots: ratio of apico-basal to lateral surface area 2 A ab /A l for the center cells, as a function of the adhesion parameter . Dashed black, blue and green lines: prediction of simplified theories describing the cell shape as an hexagonal prism, a cylinder with two spherical caps, and the union of a hexagonal prism with two spherical caps. Inset: cellular aspect ratio a/h c , with a measuring the side of the hexagonal face and h c the thickness of the sheet, as a function of . Here and h = A l /(6 a ). (F) Schematic for the polar vector P characterising the asymmetry of the cell shape. (G) Box plot for P z (blue) and | P z | (red). As adhesion increases, cells deform asymmetrically in the direction orthogonal to the planar sheet. We compare the simulation results with a theoretical prediction from a 3D vertex model on a perfect hexagonal lattice, i.e. formed by uniform hexagonal prisms with height h c and side length a ( Figs. 5C-D ). We consider that cells are subjected to an apico-basal tension , lateral tension , and an inner pressure P enforcing the cell volume to be equal to the reference volume V 0 = 4 πℓ 3 /3. We then write the corresponding effective free energy for a single cell in the tissue (S1 Appendix 9): For a hexagonal prism, the apical and basal surface area, A l = 6 ah c the lateral surface area, and the cell volume is . Minimising the effective free energy, we obtain the equilibrium ratio A ab /A l which can be compared to simulation results for the 16 inner cells ( Fig. 5E ). This area ratio is related to the cell aspect ratio, with smaller values corresponding to more columnar cells and larger values corresponding to more squamous cells. Although this simplified model captures the qualitative trend of the cell shape dependency on cell adhesion, it underestimates the simulated area ratio. We now consider alternative simplified descriptions where (i) each cell is considered as a cylinder connected to two spherical caps, or (ii) use an approximation where the surface area and volume of the cell is defined through the union of a hexagonal prism and two spherical caps (S1 Appendix 9). These choices better capture the actual simulated shapes ( Fig. 5E ), showing that taking the apical and basal curvatures play a significant role in the equilibrium cell shape. These refined models however still underestimate the area ratio A ab /A l . This can arise from the fact that these simplified descriptions do not take into account the surface bending modulus, the tissue-scale deformation due to the system finite size and free boundary conditions, and consider approximate cell shapes. In addition, we noticed that for higher values of , the cellular shapes appear more heterogeneous ( Fig. 5A ). This heterogeneity appears to be linked to an asymmetry in the cell apical and basal surface areas (see zoomed area comparing how lateral faces look for small and for larger in Fig. 5A ). To verify this, we introduce a polar order parameter for the shape of cell I ( Fig. 5F ): with the volumetric domain enclosed by , and S I and V I the surface area and volume of cell I . We calculate the order parameter for all cells in the simulation and consider its projection on the direction orthogonal to the plane containing the initial cell centers, z . We observe that, with increasing adhesion strength , the average projected cell polarity does not clearly deviate from zero, , but the average of the absolute value of the projected cell polarity, 〈| P z |〉, strongly increases ( Fig. 5G ). This suggests that at high enough adhesion, cells adopt polarised apico-basal shapes orthogonal to the plane of the tissue, with no consistent overall shape polarisation orthogonal to the tissue ( Fig. 5A ). Possibly, such a spatial arrangement favours larger contact areas, which is beneficial at large adhesion. 3.4 Adding cell divisions: growth of an organoid We now discuss simulations modelling the growth of a cell aggregate from a single cell, for which we introduce cell divisions in our framework. When a cell divides, the mother cell is replaced by two daughter cells as follows: a randomly oriented plane passing through the mother cell centre is selected, splitting the mother cell in two parts. The two daughter cells are then separated by this plane and simply fill the original shape of the mother, except for a small region that separates the daughter cells by a distance d *, perpendicular to the division plane ( Fig. 6A ). To determine timepoints of cell division, each cell is assigned a cell cycle time t D , which we take equal for all cells. We assume that in between divisions, cell volume follows a linear growth law , where V 0 is the volume of the cell at its birth. This effectively leads to cells doubling their volume during their lifetime - we note however that a small volume loss occurs at division due to the initial separation of the daughter cells by a distance d * ( Fig. 6A ). At each time point, cell volume is imposed through the Lagrange multiplier P I . Download figure Open in new tab Fig 6. Growth of small cell aggregates, driven by synchronous cell divisions and cell volume increase . (A) Each cell doubles its volume between birth and division, over a cell cycle time t D . Cell division is introduced by splitting the mother cell with a plane passing through the cell centre, and generating two daughter cells separated by a distance d * ( d */2 from the division plane). (B-D) Simulation results for two values of the ratio τ/t D . (B): Largest centre-to-centre cell distance max〈 I,J 〉 R IJ , as a function of time. Jumps correspond to cell division events. (C) Average reduced volume 〈 v 〉 as a function of time. (D) Snapshots of simulations of two growing aggregates. We then simulate the growth of an aggregate starting from a single cell ( Fig. 6B-D and S3-4 videos). The dynamics of the growing aggregate strongly depends on τ/t D , which measures the ratio between a characteristic time scale of cell shape relaxation, and the cell cycle time. For smaller values of τ/t D , the growing aggregate is more compact ( Fig. 6B ), cells have a smaller reduced volume and have therefore shapes further away from spheres ( Fig. 6C ). Here, the aggregate compactness is measured by calculating the maximum distance between cell centres ( Fig. 6B ). These observations indicate that the shape of a cell aggregate can strongly depend on a competition between its growth rate and internal mechanical relaxation times. 4 Discussion The framework of interacting active surfaces introduced here is a novel method to study the mechanics of cell aggregates such as early developing embryos or organoids, and opens the door to their systematic modelling and simulation. We have demonstrated here that it can be used to study in detail the shape of adhering cell doublets, simple epithelia, as well as growing cellular aggregates. Our method is well-suited to capture the mechanics of tissues and organoids connecting it to cell level processes such as cortical flows, cortical tension and cellular adhesion in the organisation of a cellular aggregate. In this study we have restricted ourselves to relatively simple constitutive equations for the tension and bending moment tensors ( Eqs. (3) ), and we have considered situations with a uniform and constant surface tension within each cell. Our method is based on using the virtual work principle ( Eq. (2) ), a very general statement of force and torque balance for a surface, to obtain a set of algebraic equations for the cell surface described with finite elements. As such it is versatile and we expect that more complex constitutive equations, corresponding to more detailed physical descriptions of the cell surface, can be easily introduced in our description. We now discuss some of these possible extensions of our framework. First, we have not included here apico-basal polarity, an axis of cell organisation which results from a spatially segregated protein distribution and inhomogeneous cytoskeletal structures [ 64 ]. To take this into account, one could introduce a polarity field in each cell and consider an active tension γ on the cell surface whose value at each point depends on the polarity field orientation. This could be used to introduce, for instance, differences in apical, basal and lateral surface tension which are taken into account in 3D vertex models [ 19 ]. It would be natural to introduce a concentration field on the cell surface, describing a regulator of the cortical tension, such as myosin concentration. At the level of a single surface, such coupling between cortical flows and its regulator can give rise to pattern formation, spontaneous symmetry breaking and shape oscillation [ 37 , 65 , 66 ]. The dynamics of the concentration per unit area, c , of such a regulator can be obtained from the balance equation on the surface: where D t c is the material derivative of c ( ∂ t c in a Lagrangian description), j is the flux of c relative to the centre of mass, and r is a reaction rate. A natural choice for the flux would be j = – D ∇ c to represent diffusion according to Fick’s law. A natural choice for the reaction rate would be r = k on – k off c for turnover dynamics, with target concentration c 0 = k on / k off and typical turnover time . The discretisation of such fields can be easily introduced in our framework, following the methods detailed in [ 38 , 67 ]. Importantly, the model introduced here is based on an interaction potential between cells, which can be motivated microscopically from a description of cell-cell linkers that equilibrate quickly to their Boltzmann distribution, with free linkers on the surface in contact with a reservoir imposing a constant concentration. From a computational perspective, the exact integration of the interaction potential requires the computation of double integrals, which have a large computational cost. Alternatively, one could approximate the double integrals further in the limit by considering the interaction of each point on surface only with its closest point projection on , following classical numerical approaches for the adhesion between interfaces (Ref. [ 68 ] and references therein). On the other hand, the interplay of adhesive cell-cell linkers such as E-cadherin with cortical dynamics also plays an important role in orchestrating cell adhesion [ 69 ]. Unlike the specific adhesion of solid interfaces, cell-cell adhesion dynamics involves a complex interplay between the diffusion, advection and binding dynamics of linkers [ 70 ]. Notably, E-cadherin junctions have been shown to be mechanosensitive [ 71 ] and to act to regulate the actomyosin levels at junctions [ 72 ]. To take these effects into account, an explicit description of E-cadherin concentration on the cell surface might be required. Thus, an extension of our model could introduce explicitly two-point density fields c IJ ( X 1 , X 2 ) representing the concentration of bound linkers between cells I and J , as well as a concentration field of free linkers on each cell c I . Alternatively, one could introduce cell-cell adhesion by considering a finite number of explicitely described individual linkers [ 73 ]. Our model does not account for the friction generated by relative surface flows between cells that adhere to each other, which is likely to play an important role during cell rearrangements. The effective friction stemming from an ensemble of transiently binding and unbinding linkers can be modelled effectively with a friction coefficient motivated by microscopic models such as a Lacker-Peskin model [ 74 ], which lead to predictions of force-velocity relations which depend on whether linkers are force-sensitive, e.g. slip or catch bonds [ 75 , 76 ]. One could include these terms systematically in our finite element discretisation following the ideas in [ 77 , 78 ]. In its current version and with these additions, we hope that the interacting active surface framework will be a useful tool to investigate the mechanics and self-organisation of cellular aggregates. 5 Author contributions Conceptualisation, Methodology, Writing -Original draft preparation, Writing -Review and editing: ATS, GS. Software: ATS. Validation: MKW. Funding acquisition, Supervision: GS. 6 Acknowledgements ATS, MKW and GS acknowledge support from the Francis Crick Institute, which receives its core funding from Cancer Research UK (FC001317), the UK Medical Research Council (FC001317), and the Wellcome Trust (FC001317). GS and ATS acknowledge support from a grant from the Engineering and Physical Sciences Research Council (EP/T003103/1). We thank Quentin Vagne and Guillermo Vilanova for comments on the manuscript. References 1. ↵ Heer NC , Martin AC. Tension, contraction and tissue morphogenesis . Development . 2017 ; 144 ( 23 ): 4249 – 4260 . OpenUrl Abstract / FREE Full Text 2. ↵ Stooke-Vaughan GA , Campàs O. Physical control of tissue morphogenesis across scales . Current opinion in genetics & development . 2018 ; 51 : 111 – 119 . OpenUrl 3. ↵ Lecuit T , Lenne PF. Cell surface mechanics and the control of cell shape, tissue patterns and morphogenesis . Nature reviews Molecular cell biology . 2007 ; 8 ( 8 ): 633 – 644 . OpenUrl CrossRef PubMed Web of Science 4. ↵ Huch M , Knoblich JA , Lutolf MP , Martinez-Arias A. The hope and the hype of organoid research . Development . 2017 ; 144 ( 6 ): 938 – 941 . OpenUrl Abstract / FREE Full Text 5. Kamm RD , Bashir R , Arora N , Dar RD , Gillette MU , Griffith LG , et al. Perspective: The promise of multi-cellular engineered living systems . APL bioengineering . 2018 ; 2 ( 4 ): 040901 . OpenUrl 6. ↵ Rossi G , Manfrin A , Lutolf MP. Progress and potential in organoid research . Nature Reviews Genetics . 2018 ; 19 ( 11 ): 671 – 687 . OpenUrl CrossRef PubMed 7. ↵ Graner F , Glazier JA. Simulation of biological cell sorting using a two-dimensional extended Potts model . Phys Rev Lett . 1992 ; 69 ( 13 ): 2013 – 2016 . OpenUrl CrossRef PubMed Web of Science 8. Glazier JA , Graner F. Simulation of the differential adhesion driven rearrangement of biological cells . Phys Rev E Stat Phys Plasmas Fluids Relat Interdiscip Topics . 1993 ; 47 ( 3 ): 2128 – 2154 . OpenUrl CrossRef PubMed 9. Marée AFM , Jilkine A , Dawes A , Grieneisen VA , Edelstein-Keshet L. Polarization and movement of keratocytes: a multiscale modelling approach . Bull Math Biol . 2006 ; 68 ( 5 ): 1169 – 1211 . OpenUrl CrossRef PubMed Web of Science 10. Anderson ARA , Chaplain MAJ , Rejniak KA 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. Single-Cell-Based Models in Biology and Medicine . Basel : Birkhäuser Basel ; 2007 . p. 107 – 136 . 11. Scianna M , Preziosi L , Wolf K. A Cellular Potts Model simulating cell migration on and in matrix environments . Math Biosci Eng . 2013 ; 10 ( 1 ): 235 – 261 . OpenUrl CrossRef PubMed 12. Vroomans RMA , Hogeweg P , ten Tusscher KHWJ. Segment-specific adhesion as a driver of convergent extension . PLoS Comput Biol . 2015 ; 11 ( 2 ): e1004092 . OpenUrl CrossRef 13. ↵ Louis PY , Nardi FR Boas SEM , Jiang Y , Merks RMH , Prokopiou SA , Rens EG. Cellular Potts Model: Applications to Vasculogenesis and Angiogenesis . In: Louis PY , Nardi FR , editors. Probabilistic Cellular Automata: Theory, Applications and Future Perspectives . Cham : Springer International Publishing ; 2018 . p. 279 – 310 . 14. ↵ Nonomura M. Study on multicellular systems using a phase field model . PLoS One . 2012 ; 7 ( 4 ): e33501 . OpenUrl CrossRef PubMed 15. Löber J , Ziebert F , Aranson IS. Collisions of deformable cells lead to collective migration . Sci Rep . 2015 ; 5 : 9172 . OpenUrl CrossRef PubMed 16. Najem S , Grant M. Phase-field model for collective cell migration . Phys Rev E . 2016 ; 93 ( 5 ): 052405 . OpenUrl 17. Moure A , Gomez H. Phase-Field Modeling of Individual and Collective Cell Migration . Arch Comput Methods Eng . 2021 ; 28 : 311 – 344 . OpenUrl 18. ↵ Kuang X , Guan G , Wong MK , Chan LY , Zhao Z , Tang C , et al. Computable early Caenorhabditis elegans embryo with a phase field model . PLoS Comput Biol . 2022 ; 18 ( 1 ): e1009755 . OpenUrl 19. ↵ Alt S , Ganguly P , Salbreux G. Vertex models: from cell mechanics to tissue morphogenesis . Philos Trans R Soc Lond B Biol Sci . 2017 ; 372 ( 1720 ): 20150520 . OpenUrl CrossRef PubMed 20. ↵ Bi D , Yang X , Marchetti MC , Manning ML. Motility-Driven Glass and Jamming Transitions in Biological Tissues . Phys Rev X . 2016 ; 6 ( 2 ): 021011 . OpenUrl CrossRef 21. ↵ Merkel M , Lisa Manning M. A geometrically controlled rigidity transition in a model for confluent 3D tissues . New J Phys . 2018 ; 20 ( 2 ): 022002 . OpenUrl 22. ↵ Onuki A , Kawasaki K Nagai T , Kawasaki K , Ohta S. Vertex Model of Cellular Pattern Growth in Two and Three Dimensions . In: Onuki A , Kawasaki K , editors. Dynamics and Patterns in Complex Fluids . Berlin, Heidelberg : Springer Berlin Heidelberg ; 1990 . p. 179 – 181 . 23. ↵ Nagai T , Honda H. A dynamic cell model for the formation of epithelial tissues . Philos Mag B . 2001 ; 81 ( 7 ): 699 – 719 . OpenUrl CrossRef 24. ↵ Fletcher AG , Osterfield M , Baker RE , Shvartsman SY. Vertex models of epithelial morphogenesis . Biophys J . 2014 ; 106 ( 11 ): 2291 – 2304 . OpenUrl CrossRef PubMed Web of Science 25. ↵ Staple DB , Farhadifar R , Röper JC , Aigouy B , Eaton S , Jülicher F. Mechanics and remodelling of cell packings in epithelia . Eur Phys J E Soft Matter . 2010 ; 33 ( 2 ): 117 – 127 . OpenUrl CrossRef PubMed 26. ↵ Landsberg KP , Farhadifar R , Ranft J , Umetsu D , Widmann TJ , Bittig T , et al. Increased Cell Bond Tension Governs Cell Sorting at the Drosophila Anteroposterior Compartment Boundary . Curr Biol . 2009 ; 19 ( 22 ): 1950 – 1955 . OpenUrl CrossRef PubMed Web of Science 27. ↵ Nagai T , Honda H. Computer simulation of wound closure in epithelial tissues: cell-basal-lamina adhesion . Phys Rev E Stat Nonlin Soft Matter Phys . 2009 ; 80 ( 6 ): 061903 . OpenUrl PubMed 28. ↵ Bielmeier C , Alt S , Weichselberger V , La Fortezza M , Harz H , Jülicher F , et al. Interface Contractility between Differently Fated Cells Drives Cell Elimination and Cyst Formation . Curr Biol . 2016 ; 26 ( 5 ): 563 – 574 . OpenUrl CrossRef PubMed 29. ↵ Messal HA , Alt S , Ferreira RMM , Gribben C , Wang VMY , Cotoi CG , et al. Tissue curvature and apicobasal mechanical tension imbalance instruct cancer morphogenesis . Nature . 2019 ; 566 ( 7742 ): 126 – 130 . OpenUrl CrossRef PubMed 30. ↵ Mayer M , Depken M , Bois JS , Jülicher F , Grill SW. Anisotropies in cortical tension reveal the physical basis of polarizing cortical flows . Nature . 2010 ; 467 ( 7315 ): 617 – 621 . OpenUrl CrossRef PubMed Web of Science 31. Callan-Jones AC , Voituriez R. Active gel model of amoeboid cell motility . New J Phys . 2013 ; 15 ( 2 ): 025022 . OpenUrl CrossRef 32. Turlier H , Audoly B , Prost J , Joanny JF. Furrow constriction in animal cell cytokinesis . Biophys J . 2014 ; 106 ( 1 ): 114 – 123 . OpenUrl CrossRef PubMed 33. Bergert M , Erzberger A , Desai RA , Aspalter IM , Oates AC , Charras G , et al. Force transmission during adhesion-independent migration . Nature cell biology . 2015 ; 17 ( 4 ): 524 – 529 . OpenUrl CrossRef PubMed 34. Reymann AC , Staniscia F , Erzberger A , Salbreux G , Grill SW. Cortical flow aligns actin filaments to form a furrow . Elife . 2016 ; 5 : e17807 . OpenUrl CrossRef 35. Callan-Jones AC , Ruprecht V , Wieser S , Heisenberg CP , Voituriez R. Cortical Flow-Driven Shapes of Nonadherent Cells . Phys Rev Lett . 2016 ; 116 ( 2 ): 028102 . OpenUrl 36. Callan-Jones AC , Voituriez R. Actin flows in cell migration: from locomotion and polarity to trajectories . Curr Opin Cell Biol . 2016 ; 38 : 12 – 17 . OpenUrl CrossRef PubMed 37. ↵ Mietke A , Jülicher F , Sbalzarini IF. Self-organized shape dynamics of active surfaces . Proc Natl Acad Sci U S A . 2019 ; 116 ( 1 ): 29 – 34 . OpenUrl Abstract / FREE Full Text 38. ↵ Torres-Sánchez A , Millán D , Arroyo M. Modelling fluid deformable surfaces with an emphasis on biological interfaces . J Fluid Mech . 2019 ; 872 : 218 – 271 . OpenUrl 39. ↵ da Rocha HB , Bleyer J , Turlier H. A viscous active shell theory of the cell cortex . arXiv . 2021 ;. 40. ↵ Reuther S , Voigt A. The Interplay of Curvature and Vortices in Flow on Curved Surfaces . Multiscale Model Simul . 2015 ; 13 ( 2 ): 632 – 643 . OpenUrl 41. Gross BJ , Atzberger PJ. Hydrodynamic flows on curved surfaces: Spectral numerical methods for radial manifold shapes . J Comput Phys . 2018 ; 371 : 663 – 689 . OpenUrl 42. Gross B , Atzberger PJ. Spectral Numerical Exterior Calculus Methods for Differential Equations on Radial Manifolds . J Sci Comput . 2018 ; 76 ( 1 ): 145 – 165 . OpenUrl 43. ↵ Nestler M , Nitschke I , Voigt A. A finite element approach for vector- and tensor-valued surface PDEs . J Comput Phys . 2019 ; 389 : 48 – 61 . OpenUrl 44. ↵ Dziuk G , Elliott CM. Finite element methods for surface PDEs . Acta Numer . 2013 ; 22 : 289 – 396 . OpenUrl CrossRef 45. Rodrigues DS , Ausas RF , Mut F , Buscaglia GC. A semi-implicit finite element method for viscous lipid membranes . J Comput Phys . 2015 ; 298 : 565 – 584 . OpenUrl 46. Barrett JW , Garcke H , Nürnberg R. A stable numerical method for the dynamics of fluidic membranes . Numer Math . 2016 ; 134 ( 4 ): 783 – 822 . OpenUrl 47. ↵ Sahu A , Omar YAD , Sauer RA , Mandadapu KK. Arbitrary Lagrangian–Eulerian finite element method for curved and deforming surfaces: I. General theory and application to fluid interfaces . J Comput Phys . 2020 ; 407 : 109253 . OpenUrl 48. Reuther S , Nitschke I , Voigt A. A numerical approach for fluid deformable surfaces . Journal of Fluid Mechanics . 2020 ; 900 : R8 . doi: 10.1017/jfm.2020.564 . OpenUrl CrossRef 49. Omar YAD , Sahu A , Sauer RA , Mandadapu KK. Nonaxisymmetric Shapes of Biological Membranes from Locally Induced Curvature . Biophys J . 2020 ; 119 ( 6 ): 1065 – 1077 . OpenUrl 50. ↵ Bächer C , Khoromskaia D , Salbreux G , Gekle S. A Three-Dimensional Numerical Model of an Active Cell Cortex in the Viscous Limit . Frontiers in Physics . 2021 ; 9 . doi: 10.3389/fphy.2021.753230 . OpenUrl CrossRef 51. ↵ Torres-Sánchez A , Winter MK , Salbreux G. ias ; 2021 . https://github.com/torressancheza/ias . 52. ↵ Salbreux G , Jülicher F. Mechanics of active surfaces . Phys Rev E . 2017 ; 96 ( 3-1 ): 032404 . OpenUrl 53. ↵ Purcell EM. Life at low Reynolds number . American journal of physics . 1977 ; 45 ( 1 ): 3 – 11 . OpenUrl CrossRef Web of Science 54. ↵ Salbreux G , Charras G , Paluch E. Actin cortex mechanics and cellular morphogenesis . Trends in cell biology . 2012 ; 22 ( 10 ): 536 – 545 . OpenUrl CrossRef PubMed Web of Science 55. ↵ Salbreux G , Jülicher F , Prost J , Callan-Jones A. Theory of nematic and polar active fluid surfaces . arXiv preprint arXiv:220109251. 2022 ;. 56. ↵ Feng F , Klug WS. Finite element modeling of lipid bilayer membranes . J Comput Phys . 2006 ; 220 ( 1 ): 394 – 408 . OpenUrl 57. ↵ Ma L , Klug WS. Viscous regularization and r-adaptive remeshing for finite element analysis of lipid membrane mechanics . J Comput Phys . 2008 ; 227 ( 11 ): 5816 – 5835 . OpenUrl 58. ↵ Trilinos ;. https://trilinos.github.io/index.html . 59. ↵ Marsden JE , Hughes TJR. Mathematical Foundations of Elasticity . New York, New York, USA : Courier Corporation ; 1994 . 60. ↵ Barendrecht PJ , Bartoň M , Kosinka J. Efficient quadrature rules for subdivision surfaces in isogeometric analysis . Comput Methods Appl Mech Eng . 2018 ; 340 : 1 – 23 . OpenUrl 61. ↵ Brasch J , Harrison OJ , Honig B , Shapiro L. Thinking outside the cell: how cadherins drive adhesion . Trends in cell biology . 2012 ; 22 ( 6 ): 299 – 310 . OpenUrl CrossRef PubMed Web of Science 62. ↵ Clark AG , Dierkes K , Paluch EK. Monitoring actin cortex thickness in live cells . Biophysical journal . 2013 ; 105 ( 3 ): 570 – 580 . OpenUrl CrossRef PubMed Web of Science 63. ↵ Maître JL , Turlier H , Illukkumbura R , Eismann B , Niwayama R , Nédélec F , et al. Asymmetric division of contractile domains couples cell positioning and fate specification . Nature . 2016 ; 536 ( 7616 ): 344 – 348 . OpenUrl CrossRef PubMed 64. ↵ Motegi F , Plachta N , Viasnoff V. Novel approaches to link apicobasal polarity to cell fate specification . Current opinion in cell biology . 2020 ; 62 : 78 – 85 . OpenUrl 65. ↵ Bois JS , Jülicher F , Grill SW. Pattern formation in active fluids . Physical review letters . 2011 ; 106 ( 2 ): 028103 . OpenUrl CrossRef PubMed 66. ↵ Mietke A , Jemseena V , Kumar KV , Sbalzarini IF , Jülicher F. Minimal model of cellular symmetry breaking . Physical review letters . 2019 ; 123 ( 18 ): 188101 . OpenUrl 67. ↵ Torres-Sánchez A , Santos-Oliván D , Arroyo M. Approximation of tensor fields on surfaces of arbitrary topology based on local Monge parametrizations . J Comput Phys . 2020 ; 405 : 109168 . OpenUrl 68. ↵ Sauer RA , De Lorenzis L. A computational contact formulation based on surface potentials . Comput Methods Appl Mech Eng . 2013 ; 253 : 369 – 395 . OpenUrl 69. ↵ Lecuit T , Yap AS. E-cadherin junctions as active mechanical integrators in tissue dynamics . Nature cell biology . 2015 ; 17 ( 5 ): 533 – 539 . OpenUrl CrossRef PubMed 70. ↵ Kaurin D , Bal PK , Arroyo M. Peeling dynamics of fluid membranes bridged by molecular bonds: moving or breaking . arXiv . 2022 ;. 71. ↵ Iyer KV , Piscitello-Gómez R , Paijmans J , Jülicher F , Eaton S. Epithelial viscoelasticity is regulated by mechanosensitive E-cadherin turnover . Current Biology . 2019 ; 29 ( 4 ): 578 – 591 . OpenUrl CrossRef 72. ↵ Sumi A , Hayes P , D’Angelo A , Colombelli J , Salbreux G , Dierkes K , et al. Adherens junction length during tissue contraction is controlled by the mechanosensitive activity of actomyosin and junctional recycling . Developmental cell . 2018 ; 47 ( 4 ): 453 – 463 . OpenUrl CrossRef PubMed 73. ↵ Nestor-Bergmann A , Blanchard GB , Hervieux N , Fletcher AG , Etienne J , Sanson B. Adhesion-regulated junction slippage controls cell intercalation dynamics in an Apposed-Cortex Adhesion Model . PLoS Comput Biol . 2022 ; 18 ( 1 ): e1009812 . OpenUrl 74. ↵ Lacker HM , Peskin CS. A mathematical method for the unique determination of cross-bridge properties from steady-state mechanical and energetic experiments on macroscopic muscle . In: Lecture on Mathematics in the Life sciences . vol. 16 . AMS ; 1986 . p. 121 – 153 . OpenUrl 75. ↵ Srinivasan M , Walcott S. Binding site models of friction due to the formation and rupture of bonds: state-function formalism, force-velocity relations, response to slip velocity transients, and slip stability . Phys Rev E Stat Nonlin Soft Matter Phys . 2009 ; 80 ( 4 Pt 2 ): 046124 . OpenUrl CrossRef PubMed 76. ↵ Sabass B , Schwarz US. Modeling cytoskeletal flow over adhesion sites: competition between stochastic bond dynamics and intracellular relaxation . J Phys Condens Matter . 2010 ; 22 ( 19 ): 194112 . OpenUrl CrossRef PubMed 77. ↵ Mergel JC , Sahli R , Scheibert J , Sauer RA. Continuum contact models for coupled adhesion and friction . J Adhes . 2019 ; 95 ( 12 ): 1101 – 1133 . OpenUrl 78. ↵ Mergel JC , Scheibert J , Sauer RA. Contact with coupled adhesion and friction: Computational framework, applications, and new insights . J Mech Phys Solids . 2021 ; 146 : 104194 . OpenUrl Back to top Previous Next Posted March 22, 2022. Download PDF Supplementary Material Email Thank you for your interest in spreading the word about bioRxiv. NOTE: Your email address is requested solely to identify you as the sender of this article. Your Email * Your Name * Send To * Enter multiple addresses on separate lines or separate them with commas. You are going to email the following Interacting active surfaces: a model for three-dimensional cell aggregates Message Subject (Your Name) has forwarded a page to you from bioRxiv Message Body (Your Name) thought you would like to see this page from the bioRxiv website. Your Personal Message CAPTCHA This question is for testing whether or not you are a human visitor and to prevent automated spam submissions. Share Interacting active surfaces: a model for three-dimensional cell aggregates Alejandro Torres-Sánchez , Max Kerr Winter , Guillaume Salbreux bioRxiv 2022.03.21.484343; doi: https://doi.org/10.1101/2022.03.21.484343 Share This Article: Copy Citation Tools Interacting active surfaces: a model for three-dimensional cell aggregates Alejandro Torres-Sánchez , Max Kerr Winter , Guillaume Salbreux bioRxiv 2022.03.21.484343; doi: https://doi.org/10.1101/2022.03.21.484343 Citation Manager Formats BibTeX Bookends EasyBib EndNote (tagged) EndNote 8 (xml) Medlars Mendeley Papers RefWorks Tagged Ref Manager RIS Zotero Tweet Widget Facebook Like Google Plus One Subject Area Biophysics Subject Areas All Articles Animal Behavior and Cognition (7995) Biochemistry (18704) Bioengineering (14815) Bioinformatics (44306) Biophysics (22537) Cancer Biology (19659) Cell Biology (26828) Clinical Trials (138) Developmental Biology (13922) Ecology (20948) Epidemiology (2067) Evolutionary Biology (25404) Genetics (16137) Genomics (23451) Immunology (18656) Microbiology (42369) Molecular Biology (18014) Neuroscience (93230) Paleontology (696) Pathology (2975) Pharmacology and Toxicology (5079) Physiology (8090) Plant Biology (15961) Scientific Communication and Education (2094) Synthetic Biology (4546) Systems Biology (10212) Zoology (2381) window.__CF$cv$params={r:'a3b5406b7d4093cc',t:'MTc4OTQ1MDM1NA==',u:'01a0a38d8a947b61971ef3cf8be8340a',ut:'nmkj677S3wCfp8voxpmcV7rXskqC9bB0A7OwzX2TGTg-1789450357-1.2.1.1-iy2hlxTAYW0Ieby4NYTV.elwzCSXhE2eKrVCo2McmuinSdUGXyI3ikKnAQcDaljr0l7CZOYh4eNFZLClh6yzDLLHyqgGil7ZWiusY2LsD0I',i:60};(function(){if(!document.body)return;var s=document.createElement('script');s.src='/cdn-cgi/challenge-platform/scripts/precursor/main.js';document.head.appendChild(s);})();