Effects of local curvature on epithelia tissue -- coordinated rotational movement and other spatiotemporal arrangements

Coordinated movements of epithelia tissue are linked with active matter processes. We here consider the influence of curvature on the spatiotemporal arrangements and the shapes of the cells. The cells are represented by a multiphase field model which is defined on the surface of a sphere. Besides the classical solid and liquid phases, which depend on the curvature of the sphere, on mechanical properties of the cells and the strength of activity, we identify a phase of global rotation. This rotation provides a coordinated cellular movement which can be linked to tissue morphogenesis. This investigation on a sphere is a first step to investigate the delicate interplay between topological constraints, geometric properties and collective motion. Besides the rotational state we also analyse positional defects, identify global nematic order and study the associated orientational defects.

Introduction. -Epithelia tissues display a large variety of forms, from flat to highly curved. They are involved in morphogenetic events leading to complex multicellular organisms and are intensively studied in-vivo, in-vitro and in-silico, see [1][2][3][4][5] for reviews. While the mechanics of epithelia tissues are well understood in flat space and also various properties have been identified which trigger outof-plane deformations to start the development towards highly curved structures [6,7], less is known about the influence of curvature on the spatiotemporal arrangements and shapes of the cells. A physical framework which accounts for the curvature of epithelia would therefore be highly desirable to better understand the development of multicellular organisms.
The influence of curvature on the motility of epithelial cells has been studied recently using Madin-Darby canine kidney (MDCK) cells [8][9][10]. The considered geometries are tubes, rods or wavy patterns. They all have in common a non-zero mean but zero Gaussian curvature. Furthermore topological constraints resulting from a closed surface, as e.g. an egg chamber in early embryogenesis, are not present in these investigations. These additional mathematical effects are accounted for in epithelial acini [11,12], which are spherical like structure with the epithe-lium surrounding a lumen. These experiments demonstrate not only an adjustment of cell shapes to curvature but also coordinated rotational movement. Similar persistent and continuous azimuthal rotation are observed during Drosophila oogenesis [13,14].
We will consider an idealized spherical surface to acquire deeper insights into these complex phenomena of active matter, tissue mechanics, and geometrical and topological constraints. Various modeling approaches have been considered to account for the multiscale aspects involved in the mechanics of epithelia tissue.
On a coarse-grained continuous description active gel models have shown to reproduce significant mechanical features, see e.g. [15] for a review. First attempts to formulate and simulate these models on surfaces also exist [16][17][18][19][20][21][22]. However, most of these approaches are restricted to axisymmetric surfaces or intrinsic coupling terms with the geometry. For the effect of additional extrinsic coupling terms we refer to [22][23][24][25]. Furthermore, these coarsegrained models do not allow to analyse the spatiotemporal arrangements and shapes of the cells.
Another class are vertex models. For flat geometries they are a versatile tool to reproduce e.g. solid-liquid phase transitions, number of neighbor distributions and p-1 arXiv:2202.03409v1 [cond-mat.soft] 7 Feb 2022 stress fields in epithelia tissue [26][27][28][29][30]. Various 3D approaches have been developed and also extensions towards formulations on curved surfaces [31], where the effect of curvature on the solid-liquid transition is analysed on spherical surfaces.
We here consider a multiphase field approach, see [32][33][34][35][36][37][38][39], which allows for cell deformations and detailed cell-cell interactions, as well as subcellular details to resolve the mechanochemical interactions underlying cell migration. Also topological changes, such as T1 transitions, follow naturally in this framework. Multiphase field models, together with efficient numerics and appropriate computing power, allow to analyse the connection of single cell behaviour with collective behaviour, at least for moderate numbers of cells. They already led to quantitative predictions of many generic features in multicellular systems [39][40][41]. We here formulate a multiphase field model on the surface of a sphere. We analyse the solid-liquid transition and identify a rotational mode between the solid and the liquid phase. We further consider the cell shapes, their neighbor relations and emerging nematic liquid crystal order.
Modeling. -We consider two-dimensional phase field variables φ i (x, t), one for each cell. Values of φ i = 1 and φ i = −1 denote the interior and the exterior of a cell, respectively. The cell boundary is defined implicitly by the region −1 < φ i < 1. The dynamics for each φ i is considered as for i = 1, . . . , N , where N denotes the number of cells. F is a free energy and v i is a vector field used to incorporate active components, with a self-propulsion strength v 0 . The operators ∇ S and ∆ S denote the surface divergence and Laplace-Beltrami operator on the sphere S, respectively. All quantities are non-dimensional quantities. As in previous studies [37,39,[42][43][44][45] we consider conserved dynamics. The free energy F = F CH + F IN T contains passive contributions, where with non-dimensional capillary and interaction number, Ca and In, respectively. The first is a Cahn-Hilliard energy, with W (φ i ) = 1 4 (φ 2 i − 1) 2 a double-well potential and a small parameter determining the width of the diffuse interface. Due to this energy non-interacting cells tend to become circular. The second is an interaction energy with B(φ i ) = φi+1 2 , a simple shift of φ i to values in [0, 1] and a cell-cell interaction potential approximating a short range potential, which is only active within the interior of the cell and its diffuse boundary. The approach offers the possibility to consider repulsive as well as attractive interactions which can be tuned by parameter a, see [45]. We here consider a = 2 which models repulsive and attractive interactions. For the definition of v i in eq.
(1), we follow the simplest possible approach, which can be considered as a generalization of active Brownian particles [46][47][48] to deformable objects [38]. In this approach the specified propulsion speed v 0 is the same for each cell, but the specified direction of motion, determined by the angle θ i is controlled by rotational noise dθ i (t) = √ 2D r dW i (t), with diffusivity D r and a Wiener process W i , which results in v i = (cos θ i , sin θ i ), within a local coordinate system for the tangent plane at the centre of mass of cell i. For consistency we define one of the orthonormal vectors with respect to the projected velocity vector of the previous time step.
Numerical issues. -The resulting system is solved by surface finite elements [49,50] within the toolbox AMDiS [51,52]. A piecewise flat surface triangulation is used which is adaptively refined within the diffuse interface. The other approximations in space and time follow the validated strategies for the flat case considered in [39]. Following [53] we include a degenerate term in the Cahn-Hilliard energy. This modification does not affect the asymptotic behaviour as → 0 but helps to ensure φ i ∈ [−1, 1] and to increase accuracy, see [53,54]. To efficiently solve the resulting system of equations each cell is considered on a different core. Due to the short-range interaction communication can be reduced and the implementation essentially scales with the number of cells [55]. In order to achieve a large area fraction the initial condition is constructed by regularly arranging cells with an equilibrium tanh-profile. All cells are of equal size. All simulations correspond to an area fraction of 94%. The considered parameters are comparable to those considered in the plane. We essentially use = 0.006, In = 6.0 and vary Ca ∈ [150, 450] and v 0 ∈ [0.0, 3.0]. We consider N = 6, 12, 32, 92 and adapt the radius of the sphere to allow for a constant cell size. For N = 32 the radius of the sphere is 1.
Results. -In order to understand the active system we first consider the passive case with v 0 = 0 and ask for the equilibrium configuration. In contrast to the flat situation, where this is given by a hexagonal ordering the topology of the sphere introduces defects in the ground state. The optimal arrangement of the cells is related to the Thomson problem [56] and various numerical methods have been used to find the minimal energy configuration for different N . While this remains a difficult task for large N , for the considered cases the configurations are known and have also been found using related phase field crystal models [57,58]. Figure 1 shows the corresponding cell configurations, with an appropriate color coding indicating the number of neighbors. These results correspond The color coding shows the number of neighbors (blue -four, green -five, yellow -six). The radius of the sphere is scaled to be equal. The fine details show the adaptive refined mesh for each cell.
to Ca = 150. Varying Ca does not influence the arrangements, it only modifies the shape of the cells. The equilibrium configuration can be used to define a standard deviation of the number of neighbors std eq . For N = 6, which corresponds to the arrangement of a cube and N = 12, corresponding to a dodecahedron, we have std eq = 0.0 as all cells have 4 or 5 neighbors, respectively. For N = 32 and 92, there are 12 cells with 5 neighbors, and 20 and 80 cells with 6 neighbors, corresponding to std eq = 0.47 and 0.34, respectively. A threshold on the deviation from these equilibrium values d = std − std eq ≤ 0.06 is used to identify the solid and the liquid state. This criteria is related to the coordination number in flat space, which is frequently used as an easily accessible structural property to identify solid-liquid transitions. We here consider a second criteria the mean square displacement (msd) of the centres of mass of all cells as used in [38]. Before the calculation of the msd the radius of the sphere is scaled to be equal for the different numbers of cells N . We again identify a threshold msd ≤ 0.02 to distinguish between solid and liquid. As known from flat space both activity v 0 and cell deformability, here considered through Ca, can alter the solid-liquid transition. On the sphere we have some additional degrees of freedom which can be identified by combining the two criteria. We end up with four cases: If both msd and d are below their thresholds, the cells remain in the stable configuration and the system is not moving as a whole. The system is in the solid state. If both quantities are above their thresholds, the cells are moving and are out of the stable configuration. This in-dicates a liquid state. If msd is above, but d below the threshold, the cells remain in the stable configuration but the system moves. This indicates coordinated rotational movement. We call this rotation state. The fourth possibility was never observed in our studies. Figure 2 shows the resulting phase diagram for the considered parameters. To further confirm the identified global rotation we As these angles are characteristic for the equilibrium configuration the constant values in the kymographs confirm a rotation if msd is above the considered threshold. At least for low numbers of cells, N = 6 and 12 this is the case and the rotation state is characteristic for the phase diagram. Due to the random selection of the direction of movement in our model, the rotation is not persistent, but randomly changes its axis. This randomness is also responsible for the disappearance of a rotation state for larger numbers of cells. For N = 32 the angles remain only constant over some time period, which confirms the presence of a rotation state also in this case. However, it is p-3 less stable. For N = 92 the angles change, which is consistent with the identified liquid state in the phase diagram. Other modeling approaches with a more deterministic relation of the direction to the cell shape or internal degrees of freedom, see [36,39] for realisations in flat space, will probably cure this issue and lead to global persistent rotation also for larger numbers of cells, similar to the rotation observed for active crystals on a sphere using a related active phase field crystal model [59]. Another property common for many active systems is also shown in Fig. 3. We plot the mean absolute velocity as a function of v 0 . Again this is only shown for Ca = 150. The average velocity is computed from the position of the centres of mass at each time step. The shaded region shows the variance. The deviation from the expected increase with v 0 corresponds to the transition from rotation to less ordered movement in the liquid phase.
We are next interested in the distribution of the number of neighbors as a function of v 0 . We compute this quantity in each time step and average over all times. Fig. 4 shows the resulting distributions. While the deviation for N = 6 and 12 are minor, for 32 and 92 cells we observe a clear broadening of the distribution and a shift of the maximal value towards 5. The increased heterogeneity of the numbers of neighbors for increased v 0 is also visible by looking at snapshots. Fig. 5 shows configurations with rosettes and also a sequence for a T1 transition. As in flat space these are natural features of the multiphase field model and their occurrence increases with v 0 . While in the equilibrium state the cells are isotropic, with no preferred orientation, this changes in the liquid phase. The snapshots in Fig. 5 already indicate a deformation of the cells. We quantify this by computing the aspect ratio AR, as the ratio of the major and the minor axis of an ellipse fitted to the cell. As the exact calculation on a surface is tedious we only use an approximation assuming the surface to be locally flat in the centre of mass of each cell. This is valid as long as the cell size is small compared to the surface of the sphere, which is the case for N = 92. We compute x = (AR − 1.0)/(AR − 1.0) with AR the mean of the aspect ratio of the cells. In [60] it is shown that x is universal and can be described by a k-Gamma distribution pdf (x, k) = k k x k−1 e −kx /Γ(k), with Legendre-Gamma function Γ. Fig. 6 shows the obtained results for 92 cells with Ca = 150 and various values for v 0 . According to [60] diverse epithelia systems, includ- ing Madin-Darby canine kidney (MDCK) cells, Human bronchial epithelia cells (HBECs) and the Drosophila embryo during ventral furrow formation, show a distribution with k in a narrow range between 2 and 3. The computed values for k in Fig. 6 are within this range.
This agreement further indicates the possibility to define a nematic order induced by the cell shapes. Taking the major axis in the centre of mass of each cell as a discrete orientated object allows to define a coarse-grained tangential director field. This induces nematic order and leads the presence of topological defects, see Fig. 7. This coarse-graining allows to establish a connection to surface active gel models, similar as considered for the flat case, where epithelia tissue are successively described by these theories [15]. We compute an averaged director field by averaging the major axis in each cell with their neighboring cells. This is done within the tangent planes of the centres of mass of each cell. We interpolate between the resulting directors using a Gaussian kernel where the standard deviation corresponds to the radius of one cell. This defines a global tangential nematic order on the surface. From this director field (and the corresponding orthogonal vector-field in the tangent plane) we compute a global p-4 tangential Q-Tensor field q. Defects in this field are identified by maxima of ||∇q|| 2 . To determine the topological charge we project q in the tangent plane around the defect and transform this into the xy-plane. This allows to use established methods in flat space for this issue. Following a physical approach developed in [61], which has been successfully used for multiphase field models in [39,44], we use the sign of δ = ∂q11 ∂x ∂q12 ∂y − ∂q11 ∂y ∂q12 ∂x to distinguish between +1/2 and −1/2 defects where q 11 and q 12 denote the two independent entries of a 2-dimensional Q-tensor. We track the trajectories of the defects using the software trackpy 0.5.0 [62] which implements the Crocker-Grier linking algorithm [63]. The resulting data can now be compared with results of surface active gel models on a sphere, see [64,65], and other modeling approaches and experimental data.
First, due to topological reasons the sum of the topological charges of all defects has to be equal to the Euler characteristic of the surface, which is 2 for the considered sphere. This is at least on average fulfilled, see Fig. 8(a). It shows the number of defects averaged over all time steps and for various v 0 . On average we have roughly six +1/2 defects and two −1/2 defects, the numbers slightly increase with v 0 , which gives the desired topological charge. The deviation from this average values also indicate a frequent creation and annihilation of defect pairs. The presence of −1/2 defects and the fluctuation in the overall number qualitatively differs from results for surface active gels models and also known experimental results. A prominent example is [66], considering microtubles and molecular motors encapsulated within a spherical lipid vescicle. In this setting the topological constraint suppresses the creation of additional defects and the systems remains close to the equilibrium configuration with four +1/2 defects. In the active setting they oscillate. The behaviour is reproduced by various simulation approaches [64][65][66][67]. At least the continuous approaches [64,65] also demonstrate that regular defect oscillation are only present for low activities. If the activity is increased they become irregular and in principle also the creating and annihilation of defects pairs becomes possible. Experimental data for such defects in spherical epithelia tissue are not know. Fig. 8(b) analyses the life time of the defects, which indicates a stronger persistence of +1/2 defects and Fig. 8(c) compares the average velocity between +1/2 and −1/2 defects for different v 0 . In agreement with results in flat space [39] no significant difference can be seen for the velocity. As expected the distributions shift towards larger velocities for increasing v 0 . However, also this effect is present for both +1/2 and −1/2 defects.
Conclusion. -We have considered a minimal multiphase field model for epithelia tissue on the surface of a sphere. The topological constraint induces positional defects in the equilibrium configuration. Deviations from this configuration and cell movements have been used to identify a global rotating state between the solid and the liquid state. Even for the considered random propulsion mechanism this rotation state is persistent, at least for small cell numbers. Further investigations for larger cell numbers in the liquid state indicate changes in the number of neighbor distribution and shape deformations. The elongation of cells is used to define a global nematic order. This allows to analyse the emerging orientational defects. They again follow the required topological constraint but p-5 in contrast with other systems and coarse-grained computational results obtained with surface active gel models the topological constraint does not suppress the creation and annihilation of additional defects. This might indicate fundamental differences between these systems and asks for experimental investigation on spherical epithelia tissue. * * * A.V. acknowledges support by the German Research Foundation (DFG) under Grant FOR3013 and Vo899/19-2. We further acknowledge computing resources provided at Jülich Supercomputing Center under Grant No. PFAMDIS and at ZIH at TU Dresden.