Adhesion dynamics regulate cell intercalation behaviour in an active tissue

Cell intercalation is a key cell behaviour of morphogenesis and wound healing, where local cell neighbour exchanges can cause dramatic tissue deformations such as body axis extension. Here, we develop a mechanical model to understand active cell intercalation behaviours in the context of an epithelial tissue. Extending existing descriptions, such as vertex models, the junctional actomyosin cortex of every cell is modelled as a continuum morphoelastic rod, explicitly representing cortices facing each other at bicellular junctions. Cells are described directly in terms of the key subcellular constituents that drive dynamics, with localised stresses from the contractile actomyosin cortex and adhesion molecules coupling apposed cortices. This multi-scale apposed-cortex formulation reveals key behaviours that drive tissue dynamics, such as cell-cell shearing and flow of junctional material past cell vertices. We show that cell neighbour exchanges can be driven by purely junctional mechanisms. Active contractility and viscous turnover in a single bicellular junction are sufficient to shrink and remove a junction. Next, the 4-way vertex is resolved and a new, orthogonal junction extends passively. The adhesion timescale defines a frictional viscosity that is an important regulator of these dynamics, modulating tension transmission in the tissue as well as the speeds of junction shrinkage and growth. The model additionally predicts that rosettes, which form when a vertex becomes common to many cells, are likely to occur in active tissues with high adhesive friction. SIGNIFICANCE Cell intercalation, or neighbour exchange, is a crucial behaviour that can drive tissue deformations, dissipate stress and facilitate wound healing. Substantial experimental work has identified the key molecular players facilitating intercalation, but there remains a lack of consensus and understanding of their physical roles. Existing biophysical models that represent cell-cell contacts with single edges cannot study the continuous dynamics of intercalation, involving shear between coupled cell cortices. Deriving a continuum description of the cell cortex, explicitly coupling neighbouring cortices with adhesions, we define the biophysical conditions required for successful neighbour exchanges. Furthermore, we show how the turnover of adhesion molecules specifies a viscous friction that regulates active tissue dynamics.


28
In both developing and adult animal tissues, cell rearrangements are a common mechanism by which cells actively drive 29 tissue deformation and passively dissipate stress [1][2][3][4][5]. In epithelia, directed neighbour exchanges between four cells The close ups highlight a bicellular junction (1), a vertex (2) and a rosette center (3). (C) Representation of bicellular junctions (i) and vertices (ii) in the apposed-cortex model. Adhesions can bind to neighbouring cortices within δ max (orange shading). For the given point on cell 3, blue arrows indicate the inverse lengths of possible adhesions within δ max . The red arrow represents the mean-field force resulting from the binding kinetics of adhesion when τ ad < τ cortex (equation (22) in the Supplementary Document). If τ ad ≥ τ cortex unbound cortex nodes bind to the nearest available connection and persist for an average duration τ ad . (D) The three configurations of a cell cortex in the morphoelastic framework: the undeformed configuration (parameterised by S 0 ) is taken to the virtual/reference configuration (stress-free; parameterised by S) by an active contraction, γ (the right half of the cortex has γ < 1, marked by a dark blue line). The imposed boundary conditions and body loads bring the cortex to the physically realised configuration, called the current configuration (parameterised by s), with an elastic stretch, α.
cortex in a bicellular junction separately (Figure 1Bi). Furthermore, direct access to adhesion dynamics requires 80 modelling each cell cortex explicitly with adhesion bonds coupling the cortices of neighbouring cells ( Figure 1C). 81 We begin by considering the passive mechanical properties of a cell and assume that the cortex is the leading mechanical 82 driver of cell behaviour [40]. Each cell cortex is modelled as an extensible, unshearable and torsion-free rod [41], whose 83 centreline is parameterised by a Lagrangian coordinate, S, in the reference (stress free) configuration ( Figure 1D). The 84 cortex is assumed to have resistance to bending and extension, with dimensionless elastic mechanical energy, U, given where ε = α − 1 is the strain with elastic stretch, α, and c/α is the Frenet curvature; κ = B/E/δ 0 is a dimensionless 87 mechanical parameter representing the ratio of bending, B, to extensional, E, modulus, relative to a reference length, 88 δ 0 . At rest, an isolated segment of cortex would lie as a straight rod. 89 We allow the elastic forces generated within the cortex to be dissipated viscously. Cortex relaxation is mainly driven by 90 actin turnover time, on timescales of ∼ 50 s [10,42]. Active cell rearrangement occurs on tens-of-minutes timescales 91 [17]. Given that the cortex relaxation time is significantly shorter than the timescale over which rearrangement occurs, 92 we assume that the cortex behaves as a purely fluid material and update the rest length to the current length, between 93 mechanical equilibria, at every simulation step. 94 In addition to the passive constitutive material properties of the cortex, we consider forces acting on the cortex due to 95 adhesions. Adhesion is modelled as a single agent that accounts for the composite effect of all molecules associated 96 with the adhesion complex, such as E-cadherin, αand β-catenin and vinculin [43]. Adhesion bonds acting over a unit 97 length of cortex, coupling it to another cortex, are each modelled as a simple spring, with energy ( Figure S1A) where ω is the stiffness of an adhesion bond, times the lineic density of bonds, and scaled by the extensional modulus, 99 E, of the cortex; δ is the dimensionless length of the adhesion bond (the distance between the two cortices that it is 100 attached to; Figure 1C). Lengths have been scaled relative to the rest length of an adhesion bond, δ 0 . Adhesion bonds where f ad is the external force due to adhesions. The cortical force, n, is given by (see the Supplementary Document) 2.3 Active mechanics in the apposed-cortex model 134 We use the theory of morphoelasticity to incorporate generalised active stresses in the cell cortex. We must consider 135 three configurations for the cortex: first, the initial configuration in which the cortex is undeformed and has no active 136 stress; second, the reference configuration in which the cortex is unstressed, though subjected to active (pre)stress that 137 deforms it. This is a virtual configuration that cannot, in general, be achieved in Euclidian space; third, the current 138 configuration, which is the physically realised configuration that the cortex adopts in response to active stresses and 139 external forces. These configurations are visualised in Figure 1D. The cortex is taken from its undeformed configuration, 140 parameterised by S 0 , to the virtual configuration, parameterised by S, by an active stretch: The active stretch represents local changes in material length, with γ > 1 for active growth and γ < 1 representing active 142 contraction. Though we use traditional terminology, referring to γ as a stretch, we deal only with active contraction in a 143 tissue. Setting γ < 1 could represent loading regions of the cortex with active myosin, the main driver of contractility 144 in many tissues. The effect of loading contractile machinery on the cortex can hence be conceptualised as substituting 145 cortex segments with segments that have a smaller rest length. The elastic stretch is then the local change of length 146 between the virtual configuration and the current configuration, s, due to external physical constraints: from which elastic stresses arise. We next perform simulations with active contractility. The parameters used for each 148 simulation are given in Table 1. The geometry of cell vertices, set by κ 2 /ω, is perturbed when cells are mechanically active. Simulating isotropic 151 whole-cortex contractility (setting γ < 1) in a minimal 3-cell tissue, we find that active stresses perturb vertex stability.

152
For sufficiently strong active contractility, the vertex opening is forced beyond the maximum adhesion binding length, 153 δ max , and the tissue fractures ( Figure 2C and Movie 1A). Tracking the size of the vertex opening as active contractility 154 is increased, we see that larger values of the adhesion strength, ω, produce more stable vertices that can sustain more 155 active stress before fracturing ( Figure 2D). Cell vertices are thus potential sites of weakness in a tissue and there is a 156 critical contractility threshold, beyond which the tissue will fracture.

157
If the forces acting at a vertex are not isotropic, in addition to displacing the vertex, cortical material can flow between 158 cell junctions, across the vertex (Movie 1B and Figure S2B). Indeed, from the view point of a given cell, a vertex can 159 slide along the cortex. Adhesions disconnect along the junction that is receding, while new adhesions form with the 160 cortex of the cell that is advancing. A pair of neighbouring bicellular junctions can thus commensurately increase 161 and decrease their lengths by passing material between each other, across the vertex ( Figure S2C). Importantly, this 162 behaviour cannot be manifested by models where junctions behave as spring-dashpot elements connected at vertices 163 ( Figure S3A). There, shrinkage can happen from contraction and viscous dissipation only. In this model, the cell cortex 164 behaves as an extensible, viscous rope that is anchored by cell-cell adhesions. We can conceptualise this as a rope 165 surrounding pulleys at vertices, with a friction that resists material flow between junctions ( Figure S3B). Now, shrinkage 166 can be achieved by both contraction and material flow between junctions. In the model presented here, the friction 167 associated with this shear flow depends on the adhesion timescale and is studied in further detail below. to genetically-specified asymmetric localisation of cell surface receptors between the apposed cortices of a bicellular 178 junction. We mimic this mechanism by making the local level of active contractility in a cortex dependent on the 179 identity of the cells with whom it can form such a coupling. In a similar manner to adhesion, the identity sensing 180 couplings have a maximum binding range δ γ . Unless otherwise stated, we set δ γ = δ max for simplicity. In this case, we 181 simulate the active contraction of the bicellular junction between cells 1 and 2 by setting γ < 1 in the portions of their 182 cortices where they can share adhesion bonds.
183 Figure 3 shows an example simulation of an actively driven neighbour exchange (see also Movie 2A). We find that 184 active contractility in a single junction is sufficient to shrink the junction to zero length. Importantly, however, the 185 neighbour exchange stalls at the 4-way configuration if the active contractility is lost before new adhesion connections 186 are made between the cells coming into contact (cells 3 and 4 in Figure 3). The model thus requires that δ γ is greater 187 than the spacing between the cells losing a junction (cells 1 and 2; Figure 3Bii, top), until the cells coming into contact 188 are within δ max and can form adhesion bonds (cells 3 and 4; Figure 3Bii, bottom). We highlight that using δ γ to specify 189 active contractility is a modelling choice. We could also imagine that cells 3 and 4 could be brought into contact by 190 other means e.g. active contractility lingering in the cortex after the specification is lost, or external cell-and tissue-level 191 forces helping to drive through the 4-way configuration.   apposed cortices in bicellular junctions ( Figure 4B). By contrast, for slow adhesion turnover, apposed-cortex shear will 211 gradually increase the angles between the cortex and adhesion bonds ( Figure 4D). This, in turn, causes adhesion forces 212 to oppose further shear, behaving as an effective frictional viscosity µ ad = ωτ ad . Both the area loss and apposed-cortex  The larger shear rate for slower adhesion timescales is indicative of more relative material flow between junctions. As 218 with the 3-cell tissue in Figure S2, we quantify this by tracking the locations of material points (Lagrangian tracers) in 219 the cells. Figure 4D shows kymographs of sampled material points in the whole cortex of cell 1 over the course of 220 a neighbour exchange. For both fast and slow adhesion turnover, we see that material points in the shrinking {1, 2} junction collapse as the junction contracts (see also blue circles on the cortex of cell 1 in Movies 3B&C). However, 222 this effect is exacerbated when the adhesion timescale is fast relative to the timescale of the cortex, τ ad < τ cortex , as 223 material is dragged from the rest of the cortex into the contracting {1, 2} junction. This is an extreme case where 224 vertices are acting as almost-frictionless pulleys, with significant cortical material flowing freely between junctions 225 ( Figure S3). In contrast, slower adhesion timescales increase the friction between apposed cortices and prevent material 226 transfer (notice that black lines mostly stay within the same shaded grey region in Figure 4D, lower). defines how quickly the 4-way is resolved. Note that resolution is slowest when τ ad = τ cortex because, although cells 1 252 and 2 uncouple quickly, there is less active stress to drive through the 4-way since the cortices are further apart (see 253 Figure S4A). This is an artefact of our implementation prescribing contractility using δ γ . After resolution of the 4-way 254 configuration, tension drops dramatically. The shortest timescale now present in the system is the adhesion timescale, 255 which therefore sets the rate of extension.

256
Our experimental data show that intercalating cells in the Drosophila germ-band resolve the 4-way configuration 257 with almost no delay, which is remarkable ( Figure 5B). We also find that junction extension is slower than shrinkage 258 ( Figure 5B). However, junction extension is proportionally faster, relative to shrinkage, than it is in our simulations.

259
It is likely that there are multiple mechanisms facilitating extension, including anisotropic forces from within the 260 bulk of the cells [10,52]. We asked whether this model could reproduce the experimental shrinkage and growth rates 261 without adding these additional forces. Indeed, we show that this can be done using a purely junctional mechanism by   control in reproducing different physiological conditions than models with a higher level of phenomenology. 287 We show the importance of considering whether the actomyosin cortex can flow past a vertex, from one junction to 288 the next, and that this can be regulated by the resistance to shear between apposed cell cortices. In this model, these 289 phenomena emerge as a consequence of -and are modulated by -adhesion turnover, whose timescale τ ad defines a 290 viscous friction, µ ad = ωτ ad , that regulates the penalty to tangential shearing between apposed cortices. Relative to the 291 cortical viscosity, µ ad defines regimes where we can consider cortices behaving as springs pinned together at vertices 292 versus extensible ropes passing through pulleys at vertices ( Figure S3). In the latter case, junctions can change length by 293 both contraction/expansion and by passing material past vertices. Interestingly resistance to cell area change emerges 294 from this mechanism, rather than as a phenomenological input. These results suggest that measurements of τ ad in vivo, 295 though difficult to extract, would provide valuable insights to the emergence of different tissue behaviours. 296 We provide relationships between the geometric configuration of cell vertices and the mechanical properties of the cortex 297 (κ) and adhesions (ω) that can be used to fit the model. adhesions on the biophysical properties of a tissue will be straightforward to explore in this framework.

304
Many existing models of epithelial mechanics, such as vertex-based models, treat cell intercalations as discrete, 305 instantaneous events and thus cannot study the dynamics of their behaviour. We use the continuum apposed-cortex

348
Xyz stacks were taken and a single plane is shown in Figure 1B at the position of adherens junctions. The source code implementing the apposed-cortex model is publicly available for use at github.com/

363
Alexander-Nestor-Bergmann/appcom. Documentation and quick-start tutorials can be be found at appcom.

517
Movie 1A (associated with Figure 2C). A simulation of a tissue with three cells, where active contractility is 518 increasing in the whole cortex of all cells (γ decreasing from 1) until the tissue fractures. Cell shading represents the 519 magnitude of isotropic cell stress, P cell , and cortex shading represents the strain, ε = α − 1.

520
Movie 1B (associated with Figure S2B). A simulation of a tissue with three cells, where junction shrinkage is driven 521 by asymmetric active contractility. There is active contractility (γ = 1 − 0.005) in only cells 1 and 2. Cell shading 522 represents the magnitude of isotropic cell stress, P cell , and cortex shading represents the strain, ε = α − 1.

523
Movie 2A (associated with Figure 3). A simulation where junction shrinkage is driven by active contractility in a 524 single bicellular junction. There is active contractility (γ = 1 − 0.04) in the junction shared between cells 1 and 2 525 (contractility is prescribed where the two cortices are within range δ γ = 4) and τ ad = 10τ cortex . Cell shading represents 526 the magnitude of isotropic cell stress, P cell , and cortex shading represents the strain, ε = α − 1. to cells 3 and 4; cortices next to darker (or lighter, resp.) blue lines have γ = 0.99 (or 0.9975, resp.). Cell shading 553 represents the magnitude of isotropic cell stress, P cell , and cortex shading represents the strain, ε = α − 1.

554
Movie 5A (associated with Figure 6). Multiple neighbour exchanges in a tissue where the adhesion timescale is  Movie 5B (associated with Figure 6). Rosette formation in a tissue with large adhesion viscosity, τ ad = 100τ cortex .

SUPPLEMENTARY DOCUMENT
• The reference (virtual) configuration, V, parameterised by S ∈ [0, L]. This is the unstressed configuration 9 adopted by the cortex at time t. It is a conceptual configuration that the cortex adopts in response to active 10 (pre-) stresses, but in the absence of external forces and boundary conditions. Due to physical constraints, the 11 reference configuration may not be physically realisable in Euclidean space and can be defined only locally.

12
For example, active stresses in the cortex may lead to a self-intersecting geometry, but this would be prevented 13 in physical space. 14 • The current configuration, C, parameterised by s ∈ [0, l]. This is the actual configuration that is adopted by 15 the cortex in physical space, balancing pre-stresses, body loads and boundary conditions. 16 These configurations are related through a multiplicative decomposition. Material elements in the cortex are taken from 17 the initial, undeformed configuration, C 0 , to the reference (virtual) configuration, V, by an active stretch: The active stretch represents local changes in material length, with γ > 1 for active growth and γ < 1 representing (2) The total stretch, λ, between the initial and current configuration then satisfies This decomposition is illustrated in Figure 1D.

25
In modelling the cortex as a morphoelastic rod, we follow the notation [1] and Chapter 6 in particular. Starting from the full 3D description, we provide simplified 2D equations describing the mechanics of the cortex in the apical plane. Under the rod representation, the cell cortex is defined geometrically by the position of the centre line, r(S, t), relative to the reference frame, mapping S to the fixed Cartesian basis {e 1 , e 2 , e 3 } via r = xe 1 + ye 2 + ze 3 . In the general frame, the rod is characterised by an orthonormal basis, {d 1 , d 2 , d 3 }. We set d 3 to be aligned with the tangent of the centre line, normal to the rod cross section, and d 1 and d 2 are chosen to lie in the principal directions of the cross section, with d 3 = d 1 × d 2 . The full kinematic description of the cortex is given by where v is a stretch vector (|v| > 0) and u is the Darboux vector. We detail these quantities under our specific 26 assumptions below. 27 We simplify the generalised three-dimensional formulation with some key assumptions. Since epithelial dynamics are often driven in the apical plane, we model only the apical cortex and assume that it remains planar (torsion free). Furthermore, we assume that the cortex is unshearable, has a circular cross section and is naturally straight (no reference curvature). Importantly, however, the cortex is extensible. Under these conditions, d 3 = cos θe 1 + sin θe 2 represents the rod tangent, d 1 is the rod normal and d 2 = e z . The angle, θ(S), is the deflection of the cortex relative to the e 1 axis. From these conditions, we have where u 2 = α ∂θ ∂s = ∂θ ∂S is related to the Frenet curvature, ∂θ ∂s , of the cortex. Thus, the full kinematic description (4) is simplified to where derivatives w.r.t. S (virtual) are denoted (·) .

Cortex mechanics 29
With respect to the reference (virtual) configuration the balance of linear and angular momentum in the cortex give where n(S) is the force in a cross-section of the cortex, n ext is the external force, m(S) is the moment and m ext is the 30 external moment. Following our assumptions above, the cortex behaves as a planar rod, extensible in the tangential d 3

31
direction only, such that the force can be written where E is the extensional modulus of the cortex and n 1 is unknown. The cortex moment is reduced to where B is the bending modulus of the cortex. Using (8) and (9) with (7b), assuming no external moment, we find the Thus we have These equations describe the position of the cortex in the deformed configuration relative to the reference (virtual) 37 configuration. However, for our purposes, it is more convenient to work with reference to the initial, undeformed 38 configuration, upon which we can impose the active pre-stresses. We can make the change of variables, writing and where (·) * = ∂(·)/∂S 0 . With respect to the undeformed configuration, we obtain a sixth-order system with four unknowns (θ, α, x, y): Find (θ, α, x, y) : s ∈ [0, l] → R 4 for which where (15c) and (15d) provide the position of the cortex in the current configuration. The system (15) is closed with the 41 following periodic end-point conditions: (θ * * , θ * , θ, α, x, y)| s=l = (θ * * , θ * , θ, α, x, y)| s=0 + (0, 0, 2π, 0, 0, 0).

42
Cell-level stress 43 We define the stress tensor, σ, over the region, A , of the cell enclosed by the cell boundary, C . In the current configuration, the stress tensor is symmetric and divergence free in equilibrium, satisfying σ = ∇ · (r ⊗ σ). Taking an integral over the cell area and applying the divergence theorem we have [2, 3] such that the cell-level stress tensor, averaged over the cell area, A, can be evaluated as The symmetry of stress is assured by the moment balance. The trace of the stress tensor dictates the sign of its larger 45 eigenvalue. We can therefore use the total cell pressure, P cell = −tr(σ cell )/2, to characterise the state of stress within 46 the cell: cells with P cell > 0 (P cell < 0) are under net tension (compression), with the principal axis of cell stress being tensile (compressive).
The cortex relaxation timescale, τ cortex ∼ 50 s, is largely driven by the timescale of actin turnover [4,5]. This timescale 50 is much shorter than the timescale over which rearrangement occurs (∼ 20 min [6]). We therefore assume that the 51 cortex behaves as a purely fluid material and update the rest length to the current length at every time step. Dynamics in 52 the system can then be driven by changes in the mechanical properties of cells (e.g. pre-strain, rest length and stiffness) 53 or the boundary conditions.

54
Many discrete models, such as vertex-based models, impose friction directly on cell vertices or other discrete material 55 points [7-9]. This is traditionally imposed as a dissipative drag with the substrate, though attempts have been made to 56 link dissipation to changes in cell lengths and areas [10,11]. We extend this within the apposed-cortex framework to 57 account for both dissipation in the cell cortex, on timescales τ cortex , as well as friction/drag between apposed cortices 58 due to adhesion connections, on timescales of τ ad (discussed below). We neglect inertia and dissipative interactions 59 with neighbouring fluids and structures further apically and basally"/"in the direction d 2 .

60
External forces: cell-cell adhesions 61 We consider the nature of the distributed external body forces, f ext , acting on the cortex whose work must balance the 62 variation of mechanical energy in the cortex at equilibrium.

63
Adhesions are a key mechanical feature that regulate tissue stability and are a primary source of external forces acting which removes force contributions from adhesion connections at a distance greater than δ max . We see that the adhesion 77 bonds strongly couple the equilibrium equations of neighbouring cells (12). such that unbound cortex locations form new bonds instantaneously whenever there is a neighbouring cortex within 82 δ max . We allow for the possibility of connection to a node that already has another adhesion.

83
In addition to recovery time, there is an additional -less considered -adhesion timescale: the average lifetime of 84 an adhesion complex, τ ad . This timescale represents the average time taken from binding to unbinding. It is notably 85 distinct from, and more difficult to measure than, the adhesion recovery time. To our knowledge, there have been no 86 explicit measurements of this additional timescale. However, the aggregated force from adhesions acting on the cortex 87 will be strongly affected by it. We can consider it relative to the timescale of the cortex, τ cortex , and model separately 88 the cases for τ ad τ cortex and τ ad ≥ τ cortex below.
where the force is averaged over the segment lengths, ds i and ds j . The total force acting at s i on cortex i is then the 91 sum of adhesion force to all points s j , A(s i ), connected to that location : In the case where the adhesion timescale is slower than the cortical timescale, adhesion bonds persist and maintain 93 connections between neighbouring cortices between successive iterations of simulated cortex equilibria. In this case, all 94 existing adhesion bonds in the tissue, and their cortex-cortex pairs (s i , s j ), must therefore be tracked. For simplicity, 95 we assume that unbound cortex locations form new bonds to the nearest available neighbouring cortex.

96
Fast adhesions: τ ad τ cortex . When the adhesion timescale is significantly shorter than the cortical timescale, 97 adhesion forces are averaged to work in the mean-field. In this regime, individual adhesion bonds need not be tracked.
The sum is over all other cells in the tissue, but the Heaviside in (18)

113
Nondimensional governing equations 114 We nondimensionalise lengths on the adhesion rest length (equivalent to the bicellular spacing), δ 0 , and use using tildes to denote dimensionless variables. The dimensionless force balance in the normal and tangential directions (per unit S) is then whereκ = B/E/δ 0 is a nondimensional parameter corresponding to the relative length scale over which bending Solving with respect to the deformed configuration, we can define a dimensionless state vector of unknowns 118 Ψ(s) = (θ * * , θ * , θ, α,x,ỹ). From (15), we are required to solve:

123
Parameter selection

124
When initialising a cell, we must set a dimensionless reference length, relative to the spacing between apposed cortices

132
It is difficult to measure the cortex extensional and bending moduli and the stiffness of adhesion bonds in vivo. Figure 2B 133 of the main text demonstrates that the ratio of the dimensionless parameters, κ 2 /ω, can be fitted to the size of openings 134 around cell vertices. To our knowledge, due to difficulties in achieving the required resolution, these measurements 135 have not yet been reported. Using STEAD microscopy, we do not observe significant openings around cell vertices. 136 We therefore choose order-of-magnitude estimates κ = 1 × 10 −2 , ω = 1 × 10 −2 , which give δ vert ∼ 1.33, for this 137 representative study.

138
Deriving the mechanical balance from the energetic formulation 139 In the main text, we present the model in terms of the mechanical energy (returning to the dimensional model): where c = θ . We demonstrate that this is equivalent to the rod formulation described above by deriving the balance of We derive these terms by taking the first variation of (29) to get where we have used the torsion-free kinematic identities (6) and the fact that αδd 3 = δr − (d 3 · δr )d 3 . For the 148 extensional term, we note that ε = √ r · r − 1 such that 149 δε = r · δr α = d 3 · δr .
which is equivalent to (11). The terms evaluated at the boundaries impose moment and force continuity at the end-points. The initialisation of a new tissue can be done as described in the main text, by duplicating a single cell fitted within a 157 hexagon. Alternatively, it is also possible to initialise multiple circles at randomised locations within a global stencil 158 and relax them all simultaneously. However, this requires smaller changes in the adhesion strength since neighbouring 159 cells move simultaneously and must be prevented from overlapping, thus it requires increased computational time and 160 can be less numerically stable.  This process is repeated, updating the rest length of the cortices at every time step, until a global equilibrium is reached. 166 Adaptive mesh refinement was used at every relaxation step, to ensure that node spacing remained numerically stable 167 and fast. Additional nodes were added in regions where the mesh spacing was greater thanδ 0 /4 and nodes were 168 removed where the spacing was less thanδ 0 /20. The model parameters used for the simulations are given in Table 1.

169
Source code for running the published simulations can be found at github.com/Alexander-Nestor-Bergmann/ 170 appcom. Documentation and quick-start tutorials can be be found at appcom.readthedocs.io.