Effects of curvature on epithelial tissue —Coordinated rotational movement and other spatiotemporal arrangements

Epithelial 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–5] for reviews. While the mechanics of epithelial tissues are well understood in flat space and also various properties have been identified which trigger out-of-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–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 structures with the epithelium 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 epithelial 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–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–25]. Furthermore, these coarse-grained 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 stress fields in epithelial tissue [26–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. Also active Brownian particles have been considered to model the interaction of tissue cells, fluid-like properties and velocity correlations in dense systems, see, e.g., [32–40]. Formulations of these models on a sphere are considered in [41–43].

We here consider a multiphase field approach, see [44–51] for realizations in flat space. This approach 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 [51–53]. 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. The phase transition and the neighbor relations are analysed with respect to curvature and cell shapes and emerging nematic liquid crystal order are considered for a fixed curvature with respect to activity.

We consider two-dimensional phase field variables $\phi_i(\mathbf{x},t)$, one for each cell. Values of $\phi_i = 1$ and $\phi_i = -1$ denote the interior and the exterior of a cell, respectively. The cell boundary is defined implicitly by the region $-1 < \phi_i < 1$. The dynamics for each $\phi_i$ is considered as

$$\begin{equation} \partial_t \phi_i + v_0 (\mathbf{v}_i \cdot \nabla_{\cal{S}} \phi_i) = \Delta_{\cal{S}} \frac{\delta \mathcal{F} }{\delta \phi_i}, \end{equation} \tag{ 1 }$$

for $i = 1, \ldots, N$, where N denotes the number of cells. $\mathcal{F}$ is a free energy and $\mathbf{v}_i$ is a vector field used to incorporate active components, with a self-propulsion strength v₀. The operators $\nabla_{\cal{S}}$ and $\Delta_{\cal{S}}$ denote the surface divergence and Laplace-Beltrami operator on the sphere $\cal{S}$, respectively. All quantities are non-dimensional quantities. As in previous studies [49,51,54–57] we consider conserved dynamics. The free energy $\mathcal{F} = \mathcal{F}_{CH} + \mathcal{F}_{INT}$ contains passive contributions, where

$$\begin{eqnarray} \mathcal{F}_{CH} &= \sum_{i=1}^N \frac{1}{Ca}\int_{\cal{S}} \frac{\epsilon}{2}\|\nabla_{\cal{S}} \phi_i\|^2 + \frac{1}{\epsilon}W(\phi_i)\,\text{d}{\cal{S}},\\ \end{eqnarray} \tag{ 2 }$$

$$\begin{eqnarray} \mathcal{F}_{INT} &= \sum_{i=1}^N \frac{1}{In}\int_{\cal{S}} B(\phi_i) \sum_{j\neq i} w(\phi_j)\,\text{d}{\cal{S}}, \end{eqnarray} \tag{ 3 }$$

with non-dimensional capillary and interaction number, Ca and In, respectively. The first is a Cahn-Hilliard energy, with $W(\phi_i) = \frac{1}{4}(\phi_i^2 - 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(\phi_i) = \frac{\phi_i + 1}{2}$, a simple shift of $\phi_i$ to values in $[0,1]$ and a cell-cell interaction potential

$$\begin{equation} w(\phi_j) = 1 - (a+1)\left(\frac{\phi_j -1}{2}\right)^2 + a \left(\frac{\phi_j -1}{2}\right)^4 \end{equation} \tag{ 4 }$$

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 [57]. We here consider a = 2 which models repulsive and attractive interactions. For the definition of $\mathbf{v}_i$ in eq. (1), we follow the simplest possible approach, which can be considered as a generalization of active Brownian particles to deformable objects [50]. In this approach the specified propulsion speed v₀ is the same for each cell, but the specified direction of motion, determined by the angle $\theta_i$ is controlled by rotational noise ${\rm d}\theta_i(t) = \sqrt{2 D_r} {\rm d} W_i(t)$, with diffusivity D_r and a Wiener process W_i , which results in $\mathbf{v}_i = (\cos{\theta_i}, \sin{\theta_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. In [51] the effect of this definition for $\mathbf{v}_i$ on emerging macroscopic properties is compared with various other proposed definitions. While quantitative differences occur, all considered models lead to a qualitatively similar solid-liquid phase diagram, similar number of neighbor distributions and the emergence of nematic liquid crystal order. We therefore expect similar behavior also on a sphere and only consider the described definition of $\mathbf{v}_i$ to analyze the effect of curvature.

The resulting system is solved by surface finite elements [58,59] within the toolbox AMDiS [60,61]. 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 [51]. Following [62] we include a degenerate term in the Cahn-Hilliard energy. This modification does not affect the asymptotic behaviour as $\epsilon \to 0$ but helps to ensure $\phi_i \in [-1,1]$ and to increase accuracy, see [62,63]. 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 [64]. 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 $\epsilon = 0.006$, $D_r = 0.005$, $In = 6.0$ and vary $Ca \in [150,450]$ and $v_0 \in [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. Different cell numbers thus correspond to different curvatures. For $N = 32$ the radius of the sphere R = 1.

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 [65] 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 [66,67]. Figure 1 shows the corresponding cell configurations, with an appropriate color coding indicating the number of neighbors. These results correspond 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}\| \leq 0.06$ is used to identify the solid and the liquid state. This criterion 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 criterion the mean square displacement (msd) of the centres of mass of all cells as used in [50]. 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 \leq 0.02$ to distinguish between solid and liquid. As known from flat space both activity v₀ 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 indicates 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 plot the pairwise angles between the centres of mass of all cells over time, see fig. 3. They are exemplary shown for $v_0 = 1.5$ and $Ca = 150$ and one cell fixed.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

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 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 [48,51] 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 [68]. The different behavior might also result from different sizes of the systems. As pointed out in [37–39] using active Brownian particles velocity correlations emerge. This defines a mesoscopic length scale which briefly is determined by $\xi = \sqrt{\mu/(D_r \zeta)}$, where μ is the shear modulus or stiffness of the solid phase, and ζ the substrate friction. The derivation is generic, and it should hold as long as there is no active feedback between passive and active mechanics. Such correlations have also been found in the multiphase field models [48,51]. Assuming this to be true also on a sphere, rotating states would emerge for $\xi \leq R$, and disappear for larger R, which is consistent with our finding. Another property common for many active systems is also shown in fig. 3. We plot the mean absolute velocity as a function of v₀. 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₀ corresponds to the transition from rotation to less ordered movement in the liquid phase. Within the rotating state we observe a decrease in slope with increasing N and thus decreasing curvature.

We are next interested in the distribution of the number of neighbors as a function of v₀. We compute this quantity in each time step and average over all times. Figure 4 shows the resulting distributions. While the deviations 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₀ is also visible by looking at snapshots. Figure 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₀.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

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) / (\overline{AR} - 1.0)$ with $\overline{AR}$ the mean of the aspect ratio of the cells. In [69] 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} / \Gamma(k)$, with Legendre-Gamma function Γ. Figure 6 shows the obtained results for 92 cells with $Ca = 150$ and various values for v₀. According to [69] diverse epithelial systems, including Madin-Darby canine kidney (MDCK) cells, Human bronchial epithelial 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.

Zoom In Zoom Out Reset image size

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 to the presence of topological defects, see fig. 7. This coarse-graining allows to establish a connection to surface active gel models, similar to that considered for the flat case, where epithelial tissues 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. For elongated cells 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 tangential Q-tensor field $\bm{q}$. Defects in this field are identified by maxima of $||\nabla \bm{q} ||^2$. To determine the topological charge we project $\bm{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 [70], which has been successfully used for multiphase field models in [51,56], we use the sign of $\delta=\frac{\partial q_{11}}{\partial x}\frac{\partial q_{12}}{\partial y}-\frac{\partial q_{11}}{\partial y}\frac{\partial q_{12}}{\partial x}$ to distinguish between $+1/2$ and $-1/2$ defects, where q₁₁ and q₁₂ denote the two independent entries of a 2-dimensional Q-tensor. If cells are not elongated the coarse-graining leads to a smoothed random orientation field and the defect identification might fail. We track the trajectories of the identified defects using the software trackpy 0.5.0 [71] which implements the Crocker-Grier linking algorithm [72]. The resulting data can now be compared with results of surface active gel models on a sphere, see [25,73], and other modeling approaches and experimental data.

Zoom In Zoom Out Reset image size

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. Due to the possible failure of the defect identification algorithm if cells are isotropic, this is not established for each time step, but fulfilled on average, see fig. 8(a). It shows the number of defects averaged over all time steps and for various v₀. On average we have roughly six $+1/2$ defects and two $-1/2$ defects, the numbers slightly increase with v₀, 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 [74], 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 [25,41–43,73,74]. At least the continuous approaches [25,73] 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 epithelial tissue are not know. Figure 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₀. In agreement with results in flat space [51] no significant difference can be seen for the velocity, neither with respect to v₀, nor with respect to defect type. This is in contrast to properties of active nematics on a sphere [25,41–43,73].

Zoom In Zoom Out Reset image size

We have considered a minimal multiphase field model for epithelial 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 which are similar to the flat case. 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 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 epithelial tissue. Within a broader perspective solid-liquid phase transition, coordinated cellular movement and the creation and annihilation of defects are essential in morphogenesis [75]. But even if in principle possible to address with the considered numerical approach, such shape changes are beyond the scope of this work.

AV 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.

Data availability statement: The data that support the findings of this study are available upon reasonable request from the authors.