Robust statistical properties of T1 transitions in a multi-phase field model of cell monolayers

Large-scale tissue deformation which is fundamental to tissue development hinges on local cellular rearrangements, such as T1 transitions. In the realm of the multi-phase field model, we analyse the statistical and dynamical properties of T1 transitions in a confluent monolayer. We identify an energy profile that is robust to changes in several model parameters. It is characterized by an asymmetric profile with a fast increase in energy before the T1 transition and a sudden drop after the T1 transition, followed by a slow relaxation. The latter being a signature of the fluidity of the cell monolayer. We show that T1 transitions are sources of localised large deformation of the cells undergoing the neighbour exchange, and they induce other T1 transitions in the nearby cells leading to a chaining of events that propagate local cell deformation to large scale tissue flows.


Introduction
Collective motion of cells is essential to several processes including development of an embryo, tissue morphogenesis, wound healing, homeostasis and cancer metastasis [1][2][3] .These biological processes are highly complex and orchestrate mechanical, chemical and biochemical interactions across multiple scales [4][5][6][7] .Through the interplay between directed motion, neighbour alignment and mechanical interactions, cell tissues exhibit emergent structures and dynamics that are crucial for their biological function.A fundamental underlying process for emergent large-scale behavior is the topological rearrangement of neighbouring cells, also known as T1 transition.It is a local, dissipative event that leads to remodelling of the tissue architecture and influences the large-scale flow properties of cell tissues that affects tissue homeostasis and epithelial morphogenesis [8][9][10] .
In confluent tissues, the tissue architecture can change in several ways.To isolate the tissue dynamics driven by spontaneous T1 transitions, we consider an idealised situation where apoptosis and cell division are neglected, cells have a constant volume, identical mechanical properties, and their total number is fixed.During a T1 transition, typically two neighbouring cells move apart, while two of their neighbours come towards each other and make contact as illustrated in Figure 1.The average number of neighbours before and after the T1 transition is invariant.Through T1 transitions, cells undergo large deformations and shape changes, and encounter an energy barrier that they have to overcome through their activity 11,12 .Albeit, there are several competing scenarios of the mechanical-chemical-biological feedback involved in a T1 transition, our understanding of these coupled processes is still elusive 13 .
T1 transitions are common features also in granular matter and foams under external forcing [14][15][16] .The energy relaxation after a T1 transition has been studied in foams by measuring the length of T1 junctions 17 .This concept was adapted for active tissues 18 , where the length of the T1 junction before and after a T1 transition has been measured.During a T1 transition in dry foams 15 , the cells form a rosette where either four or more edges meet.It has been shown that a junction is energetically stable for three edges incident at 120 degrees.So, while undergoing a T1 transition, the cells pass from one metastable state to another via an unstable state comprising of a rosette.In confluent tissue extracellular spaces (gaps) change this process 19 .Rosettes and tri-junctions can no longer be defined by the number of edges meeting but are placed where gaps are formed.Also, various mathematical models have been used to study different facets of T1 transitions in foams and tissues 11,[19][20][21][22][23] .They are mainly based on vertex models, which approximate cells by polygonal shapes.However, cell shape plays a crucial role in T1 transitions and the ability to accurately describe complex cell shapes, beyond polygons, might be advantageous 19,24,25 .We consider a multi-phase field model that allows for spontaneous T1 transitions while capturing the cell shape at a high resolution and allowing for large shape deformations.Multi-phase field models have been used to probe several questions pertaining collective motion of cells [26][27][28][29][30][31][32][33][34] .These models consider cells as active incompressible droplets and unlike vertex models, T1 transitions emerge spontaneously as a result of shape deformations.In vertex models, also extracellular spaces (gaps) need explicit modelling with ad-hoc assumptions 19 , whereas in multi-phase field models they are emergent.
In this paper, we focus on characterising the energy profile preceding and succeeding T1 transitions.We show that this energy profile is statistically robust to changes in several model parameters.It is characterized by an asymmetric profile with a fast increase in energy before the T1 transition, a sudden drop after the T1 transition, followed by a slow relaxation.The relaxation profile provides insights into the flow properties of the tissue.Previously, the relaxation has been indirectly studied in tissues by examining the relaxation of an ellipsoid droplet immersed in a tissue 17,19,35 .The relaxation profile was attributed to yield stress due to limitations in the measurement timescales 35 , and also associated with the fluidization of the tissue 19,36 .We further consider the duration of T1 transitions and find that the average duration scales inversely to the maximum average energy attained during the T1 transition.Also, we show that T1 transitions may trigger the creation of other T1 transitions nearby and the chaining of T1 transitions leads to large-scale deformation and fluid like behaviour.
We introduce the multi-phase field model in Section 2 and discuss results on local statistical properties of T1 transitions in Section 3. We further analyse the dependency of these statistical properties on various model parameters.The effect of cell deformability and activity is considered in detail.We study the impact of chaining of T1 transitions on flow at larger scales.In Section 4 we relate these finding to mechanical and rheological properties of the tissue and postulate that they can be used to characterize fluidization.Details on the numerical methods, the initialization and characterization of T1 transitions are provided in Section 5.

Multi-phase field model
We represent a two-dimensional confluent cell tissue within a multi-phase field model following formulations 28,[32][33][34] .We consider a system of N cells of equal area occupying a square domain of size [0, L] × [0, L] and use periodic boundary conditions.Each cell is represented by a scalar phase field φ i (x,t) as an indicator function of the domain occupied by each cell labeled by i = 1, 2, • • • , N. Namely, the bulk phase values φ i ≈ 1 and φ i ≈ −1 indicate the interior and exterior of the cell, respectively.The cell boundary is defined by the localised transition region between the two bulk values.The time evolution of the i-th phase field follows a conservative dynamics which preserves the cell areas and is given by where ∆ is the two-dimensional Laplacian applied to the variational derivative of a free energy functional F with respect to the phase field φ i .The free energy F = F CH + F INT contains the Cahn-Hilliard energy and the interaction energy 28, 34 The capillary number Ca and interaction number In are tuning parameters for the cell deformability and the strength of mutual repulsion/attraction interactions, respectively.In equation 2, the Cahn-Hilliard energy has a local free energy density given by the double well potential with the minima corresponding to the two bulk values and a gradient energy.The parameter ε controls the width of the diffuse interface.The Cahn-Hilliard energy ensures phase separation into two bulk regions which are separated by a thin, diffusive interface.This energy alone is minimised by cells with circular shapes.In equation 3, each cell's interior and interface (B(φ i ) = (φ i + 1)/2) is coupled with every other cell through a local interaction potential, where the parameter a = 1 models repulsion, while a > 1 models attraction and repulsion (see 34 for a detailed analysis of role of a).
Cell activity is introduced through the advection velocity v i (x,t) in equation 1 and is given by where v 0 is a constant parameter that controls the magnitude of the activity, e i = [cos θ i (t), sin θ i (t)] where θ i is the orientation of the self-propulsion which evolves as The first term on the right side of equation 5 is a rotational diffusion term with a Wiener process W i .The second term is a relaxation to the orientation of the cell's shape elongation.The cell elongation is identified by the principal eigenvector of the shape deformation tensor 29,32 which is symmetric and traceless and has the two components Its corresponding eigenvalues are λ ± i = ± S 2 i,0 + S 2 i,1 and eigenvectors are η ± i = ( . The vector η + i is parallel to the elongation axis of the cell and determines the preferred self-propulsion direction as Therefore, the second term on the right hand side of equation ( 5) aligns θ i (t) with β i (t).The parameter α controls the time scale of this alignment of the self-propulsion direction with the elongation axis of the cell.There are different possibilities to define the advection velocity v i (x,t) (see Ref. 32 for an overview and comparison).The current form includes approaches of Ref. 29 and, as the elongation is a result of the interaction with neighbouring cells, it accounts for contact inhibition of locomotion 37,38 .The model leads to properties appropriate to describe, e.g., Madin-Darby canine kidney (MDCK) cells 32,39 .The junctions that shrink and expand are called T1 junctions.We refer to Section 5 for the procedure to detect T1 transitions and their durations.A T1 transition not only leads to topological rearrangements of the four neighbouring cells, it also involves deformation of the cells.While details, such as the specific shape of the cells and their deformation, the duration of the T1 transition and the relaxation process differ between T1 transitions, we will demonstrate that robust statistical features of T1 transitions exist.We define the epicenter of a T1 transition as the point with the minimum total distance from the centers of the involved cells in the neighbour exchange process midway through the T1 transition.We define the immediate region around the epicenter as the core of a T1 transition, which is of essence because it is the region where T1 junctions shrink and expand, and the gap appears and disappears.
Figure 2a shows the total free energy density midway through a T1 transition.The epicentre is shown by the white dot and the estimated core is highlighted by the green circle.It has a radius r avg = 0.02L, where L is the side length of the computational domain.We compute a coarse-grained energy whose value at any point in the domain is the average of the energy density in a circular region centered at that point with radius r avg .Figure 2b shows this coarse grained energy field f r avg , which we will call 'energy' field in the following.The signature of triple-junctions and T1 transitions already becomes appealing due to their higher energy.The difference between both is enhanced by using a log scale, see Figure 2c.Considering this energy field in the epicenter over time provides a spatial-temporal description of T1 transitions.For discussions on the sensitivity of this procedure on r avg we refer to Section 5.
Figure 3a shows the time evolution of this energy averaged over 158 T1 transitions.The time is negative before a T1 transition and is positive after a T1 transition, and is denoted by t T 1 .The energy during the T1 transitions is excluded, which leads to a discontinuity at t T 1 = 0.The two values at t T 1 = 0 correspond to the averaged energies at the start and the end of the T1 transitions.As the duration of T1 transitions differs, an averaged energy as a function of time during the T1 transition does not provide any meaningful information.Details on the energy during the T1 transition are shown in Figure 3b using a normalized time.The energy profile in Figure 3a, 3b has a peak at the T1 transition.The profile is asymmetric with a strong increase in energy before the T1 transition and a sudden decrease after the T1 transition followed by a slow relaxation.The asymmetry can be quantified by considering the 75% of the maximum value, which is marked in Figure 3a.Figures 3c-3e illustrate the evolution for one T1 transition, the one depicted in Figure 1.These figures contain overlays of several snapshots as per the time marked in the figures.The darkest of these snapshots pertains to the latest time.The yellow region marks the estimated core of the T1 transition.The asymmetry before and after the T1 transition, Figure 3c, 3e, respectively, is clearly visible.The T1 junctions are longer at t T 1 = 5 compared with t T 1 = −5.During the T1 transition, Figure 3d, the asymmetry is less pronounced.Most of the deformations are concentrated in the core.These deformations arise as a result of the formation of the gap, and subsequently its disappearance.The shrinking and formation of T1 junctions and the deformations within the core are a signature of the T1 transition.However, they also influence the deformation of the four cells outside of the core, and their neighbours, which can be perceived by the overlayed cell shapes.Interestingly in the depicted T1 transition, the deformations of each of the four cells seems to be persistent before, during and after the T1 transition (see the arrows indicating the direction of deformations).We will elaborate on this and other coarse grained effects in Section 3.4.The energy profile indicates an accumulation of energy to reach the energy barrier at the T1 transition.This is due to probing several possibilities in local movement and cell shape deformation, which are coupled by the definition of activity, taking into account cell elongation and contact inhibition of locomotion.After the energy barrier has been overcome the fast relaxation of the energy can be associated with a steep gradient in the energy landscape in one direction.
The asymmetric shape of the energy profile is robust to changes in most model parameters, as demonstrated in Figure 4 where α, Ca, a, D and v 0 are varied and the energy profile associated with passive sheared foams is included for comparison.Figure 4b shows the energy rescaled by the maximum energy as changes in Ca directly affect the free energy, see equation (2).Within the range of parameters explored, the changes in the values of alignment parameter α, interaction coefficient a, and diffusivity D have minimal effects on the energy profile.We see that the profile is robust even in absence of noise (D = 0) (Figure 4d).On the other hand, the profile deviates from Figure 3a for low values of v 0 and Ca. Figure 4e shows that the cell activity v 0 affects the rate at which the cells approach a T1 transition which is indicated by the slower accumulation of energy for low v 0 .However, change in v 0 has a minor effect on the energy relaxation immediately after a T1 transition.The slow relaxation afterwards is largest for large values of v 0 .This can be associated with the definition of activity, which is related to cell elongation and at least on average cells elongate in the direction of movement after the T1 transition.The characteristic profile of the accumulation of energy before the T1 transition and the fast relaxation of energy after the T1 transition is also present for low values of Ca, see Figure 4b.However, as Figure 4b considers a rescaled energy the actual rates depend on Ca.The slow relaxation after the sudden decrease only slightly depends on Ca.We would like to point out that the results for low values of v 0 and Ca should be considered with care, as the number of T1 transitions considered in these cases is much lower.While the system is still in the fluid phase, the extreme values for v 0 = 0.1 and Ca = 0.05 already approach the transition to the solid phase.
In passive foams T1 transitions can be induced by applying shear.This is considered by an advection velocity field v i (x,t) = 0.5|x 1 − L/2| and the resulting energy profile is compared with the profile from Figure 3a, see Figure 4f.The profiles differ before the T1 transition and within the slow relaxation, but are similar in the sudden drop of energy right after the T1 transition.The latter reiterates that the energy relaxation right after a T1 transition is independent on activity.The differences in the accumulation of the energy can be associated with the persistent orientation of advection velocity due to shear, which results in collective deformation and a more deterministic approach of the T1 transition.Also the termination of the decay in the passive case results from the restricted possibilities of relaxation due to the applied shear.| while the active case corresponds to parameters in Table 1.

Duration and other properties of T1 transitions
As mentioned earlier, the duration of T1 transitions strongly depends on the specific cell arrangements.We now discuss the statistical properties of the duration.Figure 5a shows the probability distributions of the duration of T1 transitions.The distributions peak at smaller values and have a long tail for larger values.The profiles corresponds to repulsive and adhesive (a > 1), and only repulsive interactions (a = 1), and are fitted by Gamma distributions.The average duration of T1 transitions for repulsive interactions (3.418 measured for 539 T1 transitions across 3 simulations) is smaller compared to that for repulsive and adhesive interactions (3.826 for 631 T1 transitions across 4 simulations).Keeping other parameters fixed, the average number of T1 transitions in the repulsive and adhesive case was 157.75 while for the repulsive case was 179.66, respectively.Therefore, in the repulsive case, cells undergo neighbour exchanges faster and more often.transitions as a function of the maximum energy reached during a T1 transition.While the data is scattered, it qualitatively shows that high energy T1 transitions are faster.This qualitative result holds for both cases and can be explained by a larger accumulation of energy in the core, which increases the spatial energy gradients and in turn speeds up the relaxation of the energy which leads to the shorter duration.Figure 5c shows the averaged shape index (perimeter/ √ area) of the four cells involved in a T1 transition as function of time relative to a T1 transition.The asymmetry found for the energy profile and the discontinuity at t T 1 = 0 is also present for this quantity.The cells deform and elongate as they approach a T1 transition and relax afterwards.This increases and decreases their shape index, respectively.The faster relaxation leads to the asymmetry in the evolution of the shape index.The asymmetry around a T1 transition is also seen in the average velocity of the center of mass of the cells involved in a T1 transition as shown in Figure 5d.While the velocity is almost constant before the T1 transition, the velocity peaks at the T1 transition and slows down afterwards until it reaches the average value before the T1 transition.The peak in the average velocity of the center of mass is due to the large deformations of the portions of cells within the core and their fast relaxation after the T1 transition.Both quantities, the shape index and the cell velocity of the four cells involved in a T1 transition are also experimentally accessible.These quantities can be related to the energy considered above.

Effect of cell deformability, activity and gaps on T1 transitions
The asymmetric energy profile in Figure 3a is robust to tuning of most of the model parameters.Significant variations only occur for low values of Ca and v 0 , see Figure 4b, 4e.We now analyse the effect of cell deformability and activity on T1 transitions in more detail.This requires a detailed analysis of the influence of gaps.The gap fraction is related to the confluency as confluency = 100(1 − gap fraction).It essentially is a fixed quantity set by the initial data.We fix all parameters as per table 1 and compare two different initial cell sizes, denoted by 'low gap' with gap fraction 0.00048 and 'high gap' with gap fraction 0.00212.Both can be considered as confluent.The number of T1 transitions within the considered time frame is not influenced by this variation.The total numbers of T1 transitions are 162 and 158 for low and high gap cases, respectively.However, the  (g) -(l)).Total T1 considers the total number of T1 transitions within the considered time frame, T1 duration is the averaged time from start to end of all T1 transitions, Gap fraction is the extracellular space, considered as ∑ i B(φ i ) below a fixed threshold, again averaged over time, Shape index considers the averaged shape index of the four cells involved in the T1 transitions.Time between T1 is the average time a cell spends between successive T1 transitions, Max energy is the maximum energy reached at a T1 transition and v avg is the average velocity of center of mass of all cells.average duration of T1 transitions is reduced by reducing the gap fraction.The values are 2.559 and 3.794 for low and high gap cases, respectively.We measure the gap fraction as the fraction of domain where ∑ i B(φ i ) is less than a fixed threshold which is set to 0.2.This essentially excludes possible partial overlap of the diffuse interface region of cells and only accounts for gaps at tri-junctions and rosettes.This makes the measured gap fraction to depend on deformability and activity.For the considered cases low Ca leads to rounder cells with stronger overlap of the diffuse interfaces of the cells, which are in contact.This leads to an increase in the measured gap fraction, see Figure 6c.A similar dependency, but smaller in magnitude, is found for activity.Larger v 0 lead to stronger interactions between cells and thus more overlap of the diffuse interface region of cells in contact which again leads to an increase in measured gap fraction, see 6i.The gap fraction in both figures is the average quantity over the considered time frame.Both results and the dependencies discussed below are considered for the 'high gap' setting.
As shown in Figure 6a, the number of T1 transitions increases with increasing cell deformability parameter Ca.Cells that are more deformable can more easily acquire the shape deformations associated with T1 transitions.When Ca is low, these deformations are energetically more expensive resulting in fewer T1 transitions.Also the duration of T1 transitions depends on Ca, as shown in Figure 6b.T1 transitions are shorter when cells are more deformable.We suspect that this might be due to the presence of smaller gaps at T1 transitions, as this requires less shape deformation.Figure 6d shows the average cell shape index of the four involved cells in a T1 transition as function of cell deformability Ca.The shape index increases as deformability increases.The shape index of Ca = 0.05 is less than that of a regular pentagon.The shape index of regular pentagon (3.813) was attributed as the critical shape index for jamming transition in classical vertex models 40 without gaps.It has been argued that gaps influence the mechanical properties and solid-liquid transition 17 , which might explain this discrepancy, as our system is still within the fluid phase.Further details, which are related to the previous dependencies are shown in Figures 6e and 6f. Figure 6e shows the average time a cell spends between successive T1 transitions as function of Ca.This quantity is large for low Ca but decreases and plateaus to low values upon increasing Ca. Figure 6f shows the maximum energy reached during a T1 transition against Ca.We see from the dotted curve that the maximum energy is proportional to 1/Ca.Recall that 1/Ca scales the Cahn-Hilliard energy as per equation ( 2).This means that F r avg is primarily affected by the Cahn-Hilliard energy, which explains the correspondence of our results with the length of T1 junctions discussed earlier and considered in 17,18 .
The dependency on v 0 shows qualitatively similar behaviour for the number of T1 transitions, the duration of T1 transitions, the shape index of the cells involved in T1 transitions and the time a cell spends between successive T1 transitions, see Figures 6g, 6h, 6j, 6k, respectively.The increase in T1 transitions and decrease in the time between T1 transitions with activity is a property of active systems, which are driven out of equilibrium.T1 transitions are topological defects and thus an indication of out of equilibrium.The decrease in duration with increasing activity can again be associated with the decrease in measured gap fraction, see 6i, and also the increasing shape index with activity is a direct consequence of the form of active forcing considered.Figure 6l shows the average velocity of center of mass of all cells as a function of v 0 .As expected, activity is primarily converted into motion with an almost linear dependency.

Chaining of T1 transitions
So far, we have analysed robust statistical properties of T1 transitions within their cores.However, we have also seen that these local features influence the position and shape of the four cells involved in a T1 transition, and their neighbours.This can induce new T1 transitions and lead to the formation of chains of T1 transitions as illustrated in Figure 7.Each of these images consists of 10 tissue states captured at equally-spaced time instants and overlaid on top of each other.The cell shapes outlined in the darkest colors correspond to the latest time.The yellow circles mark the cores of the T1 transitions at those time instants.The chaining of T1 transitions is a result of the assumptions on constant cell area and a confluent tissue.Any cell deformation associated with a T1 transition induces deformation of the neighbouring cells and thereby increases the possibility of new T1 transitions.This is further enhanced by activity and the considered propulsion mechanism which favours the direction of elongation.
This chaining of T1 transitions is also observed experimentally in sheared foams 41 and in our simulations of passive foams which are sheared with a constant shear velocity profile.For v 0 = 0, typically one or two T1 transitions occur due to the initial non-equilibrium configuration of the tissue.As cells relax toward an equilibrium state, their motility is reduced which prevents any further T1 transitions.The situation for small v 0 is similar.The tissue becomes jammed by cells being caged amongst their neighbours and no T1 transitions occur 32 .Furthermore, when cell deformability (Ca) is low, the energetic cost for cell deformations that are necessary to undergo T1 transitions is high, which prevents or at least reduces T1 transitions and the tissue also becomes jammed 32 .This corresponds to the low number of T1 transitions in Figure 4b, 4e for low Ca and low v 0 , respectively.
However, in the considered case in Figure 7 we are far away from jamming and the chaining of T1 transitions leads to cell deformation propagating to larger scales.This is highlighted in Figure 8a, which shows the evolution of the cell tissue in the whole time window considered in Figure 7 together with the trajectory of the center of mass of the colored cells, which highlights the movement on larger spatial scales.The chaining of T1 transitions is also a source of large-scale flows as evidenced in Figure 8b.We consider the velocities of the centers of mass of all cells, average this quantity with the neighboring cells and construct a continuous velocity field by interpolating in space.The velocity field is shown together with the cell boundaries at t = 52.The mean direction corresponds with the direction of the black path shown in Figure 8a.However, as the variations in magnitude and direction of the flow field in Figure 8b indicate, T1 transitions can also induce fluctuations and could play an important role in sustaining chaotic flows (active turbulence) in cell tissues [42][43][44]    Large-scale tissue deformation requires cellular rearrangements.The simplest rearrangement in confluent cell tissue is a T1 transition.We have analysed these neighbour exchanges among cells in detail using a multi-phase field model and identified a characteristic asymmetric energy profile, see Figure 3.The energy profile has a peak at the T1 transition.The profile is asymmetric with a strong increase in energy before the T1 transition and a sudden decrease after the T1 transition which is followed by a slow relaxation.Detailed studies on the dependency of this profile on model parameters show robustness to variations in most parameters.They also allowed to associate the strong energy increase before the T1 transition with the strength in activity.This region is characterized by an accumulation of energy to reach the energy barrier at the T1 transition.This is achieved by probing several possibilities of direction of movement and shape deformation.This process is enhanced by activity, which is quantified by Figure 4e.In contrast to this the sudden relaxation after the T1 transition can clearly be associated with energy relaxation.It is almost independent of activity, see Figure 4e, and cell deformability, see Figure 4b, and also present in sheared foams, see Figure 4f.We would like to remark that the behaviour is independent but the actual slope and duration of this regime depends on deformability, as the energy is scaled in Figure 4b.The sudden decrease is associated with a steep gradient in the energy landscape in one direction set by the deformation of the cells in the core of the T1 transition.The third characteristic region, the slow relaxation, depends on activity and cell deformability.This relaxation profile provides insight in the mechanical properties of the tissue.Similar energy profiles have been obtained by actuation and relaxation of magnetic microdroplets which are injected into the tissue 17,19,35 .In these experiments a slow relaxation is associated with the fluidization of the tissue 19,35 , while stagnation of the relaxation indicates more solid-like behaviour 17 and is associated with irreversible (plastic) tissue rearrangements.We postulate that these mechanical characterizations can also be obtained from the energy decay of the T1 transitions.
In the considered confluent tissue the type of interaction between the cells, if repulsive or repulsive and attractive, seems to play a minor role on the characteristic energy profile of a T1 transition, see Figure 4c.However, the degree of confluency is known to influence the solid-fluid phase transition 35 .Increasing the extracellular space enhances fluidization.While we only consider the fluid phase, we observe an increased duration of T1 transitions for larger extracellular space.A finite duration of T1 transitions in cell tissues has been associated with molecular processes and is considered in an adhoc manner in vertex models 22 .Within the multi-phase field model a finite duration is a result of the mechanical properties of the cells and the their interactions.An increased duration of T1 transitions is observed for low deformability and low activity, see Figures 6b and 6h, respectively.Both indicating more solid-like behaviour, which is consistent with 22 , where increased duration of T1 transitions leads to decreasing the overall number of T1 transitions and a possible stiffening of the global tissue mechanics.However, these results don't take extracellular space into account.
Even if characterized locally, due to the confluent cell tissue, large enough deformations induced by T1 transitions lead to permanent cell deformations in the neighbourhood, which can trigger other T1 transitions, leading to a chaining effect.This behaviour is associated with the foam-like architecture and consistent with previously reported nonlinear tissue mechanics 35 .It is this chaining of T1 transitions which allows for large-scale tissue deformations and flow patterns which can be associated with sustaining chaotic flows, see Figure 8b.
We believe these results also to hold in more general situations, e.g. for varying cell sizes and varying mechanical cell properties.

Model Parameters
Unless otherwise specified, we use the model parameters as per Table 1 τ τ save T L ε v 0 a Ca In D r α 0.005 0.5 150 100 0.15 0.5 1.5 0.2 0.1 0.1 0.1 Table 1.Default values of the model parameters.

Finite element simulations
The simulations are run for time interval [0, T ] discretised into N t units with a uniform timestep size τ, i.e.T = N t τ.We employ a semi-implicit discretization in time.Discretization in space follows the finite element method.We adaptively refine the diffuse interface and employ a parallelization approach which scales with the number of cells.For details we refer to 28, 32-34, 45, 46 .The algorithm is implemented in the open-source library AMDiS 47,48 .

Figure 1 .
Figure 1.Successive time snapshots of tissue section undergoing a T1 transition, a finite-time neighbour exchange process between cells A, B, C and D. The transition starts when cells B and D lose contact and is completed when cells A and C make contact.During the T1 transition an extracellular space (gap) is formed between cells A, B, C and D. Also see Supplementary Movie 1

Figure 2 .
Figure 2. (a).Free energy density, in a region surrounding a T1 transition.(b) and (c) Coarse-grained energy density in a linear and log-scale, respectively.The white dot represents the epicenter of the T1 transition while the green dotted circle represents the coarse graining radius r avg , the estimated core of the T1 transition.

1
Energy profile of T1 transitionsWithin our multi-phase field approach, T1 transitions are neighbour exchange processes with a finite duration.A prototypical time sequence of a T1 transition is illustrated in Figure1.Four cells A, B, C and D are involved.Before the T1 transition, the cell junction shared by cells B and D shrink.The T1 transition starts when the cells B and D break contact and move apart.This results in the formation of an extracellular space which we call 'gap'.Cells A and C move towards each other, close the gap, and form a new contact concluding the T1 transition.After the T1 transition, this new junction between cells A and D expands.

Figure 3 .
Figure 3. (a) Evolution of energy (averaged for 158 T1 transitions) at the epicenter of the T1 transitions.Negative time corresponds to time before a T1 transition and positive time corresponds to time after a T1 transition.The shaded region denotes a width of 1 standard deviation.The gray dashed line is the average energy across the whole domain.(b) Average energy profile during a T1 transition as function of percentage of T1 duration.The standard deviation is also indicated.(c), (d) and (e) Montages of deformed cells involved in a T1 transition.Each montage is made up of 5 images, that capture the cells at equidistant times, stacked over each other.The darkest colored overlay represents the latest time.(c) Cell shapes before the start of the T1 transition, (d) during the T1 transition, and (e) after the end of the T1 transition.Also see Supplementary Movie 2 for corresponding simulation

Figure 4 .
Figure 4. Evolution of energy for different parameter values.The pink and cyan shaded region are used to denote time before and after the T1 transitions, respectively.The number of T1 transition used to obtain these results in indicated.(a) The aligning parameter α is varied.(b) The parameter to control cell deformability, Ca is varied.As Ca is a parameter that influences the overall total energy, for better comparison the energy is rescaled by division with the maximum energy.(c) Adhesion and repulsion corresponds to a = 1.5 and repulsion corresponds to a = 1.(d) The diffusivity D is varied.(e) The magnitude of the activity v 0 is varied.(f) The passive shear corresponds to advection field v i (x,t) = 0.5|x 1 − L 2 | while the active case corresponds to parameters in Table1.

Figure 5 .
Figure 5. (a) Probability distributions of the duration of T1 transitions for only repulsion interactions (magenta dots) and for both repulsion and adhesion (cyan dots).Both data sets are fitted by Gamma distributions highlighting the exponential tails.(b) Scatter plot of duration of T1 transition as function of the maximum energy reached during a T1 transition.(c) Evolution of average shape index and (d) Evolution of the average velocity of center of mass of the cells involved in the T1 transitions as function of time relative to a T1 transition.The shaded regions mark the standard deviations of both quantities.

Figure 6 .
Figure 6.Dependency of various properties on deformability Ca ((a) -(f)) and activity v 0 ((g) -(l)).Total T1 considers the total number of T1 transitions within the considered time frame, T1 duration is the averaged time from start to end of all T1 transitions, Gap fraction is the extracellular space, considered as ∑ i B(φ i ) below a fixed threshold, again averaged over time, Shape index considers the averaged shape index of the four cells involved in the T1 transitions.Time between T1 is the average time a cell spends between successive T1 transitions, Max energy is the maximum energy reached at a T1 transition and v avg is the average velocity of center of mass of all cells. .

Figure 7 .
Figure 7. Chaining of T1 transitions.Each panel is a montage of 10 snapshots of tissue configurations taken successively at constant times intervals.Latest time is represented by the cell shapes marked in the darkest color shades.The cores of the T1 transitions are highlighted in yellow.Also see Supplementary Movie 3.

Figure 8 .
Figure 8.(a) Montage of tissue snapshots from time t = 25 to t = 79 (see figure 7).The black path is the trajectory of the center of mass of the 11 coloured cells.(b) LIC visualization of streamlines, magnitude (color) and direction (black arrows) of the flow velocity.The velocity and the cell boundaries correspond to time t = 52.