Introduction

The transparent cornea, the most anterior structure of the vertebrate eye, consists of a hypocellular, collagenous stroma covered by a stratified epithelium, resting on an inner endothelial monolayer (Forrester et al., 2020). The corneal epithelium, which is contiguous with the conjunctival epithelium that covers the rest of the anterior of the eye and the inner edge of the eyelids, is a paradigm model of tissue maintenance by stem cells (Cotsarelis et al., 1989; Lavker et al., 1991; Lehrer et al., 1998). Cells are constantly lost from the corneal surface by abrasion and must be replaced. Multiple lines of evidence have shown that the stem cells that maintain the corneal epithelium reside in a narrow ring at the edge of the cornea, at the boundary between the corneal and conjunctival stem cells – the limbus (Cotsarelis et al., 1989; Richardson et al., 2017). Limbal epithelial stem cells themselves are organised as an outer ring of slow-cycling, label-retaining cells that most likely contribute to corneal regeneration after wounding, and an inner layer of rapidly dividing cells that contribute to normal corneal homeostasis (Altshuler et al., 2021; Farrelly et al., 2021). Stem cells in the basal layer of the limbal epithelium divide and produce proliferative, undifferentiated ‘transit-amplifying’ (TA) cells that migrate into the corneal epithelium to replace those lost during normal life. Studies of human ‘hurricane’ keratopathies and mosaic patterns of cell staining in genetically tagged rodents, frogs and fish have shown that the tracks of migrating cells show a remarkable, regular, radial pattern in the uninjured corneal epithelium, usually with a central vortex (Bron and Tripathi, 1973; Dua et al., 1993; Collinson et al., 2002, 2004; Nagasaki and Zhao, 2003; Iannaccone et al., 2012; Di Girolamo et al., 2015; Nejad et al., 2015). The pinwheel patterns of cell staining observed suggest that migration is highly regular and potentially globally coordinated. The vortices of cell motion patterns observed at the corneal epithelial centre resemble a logarithmic spiral and lend themselves to mathematical analysis (Iannaccone et al., 2012). Modelling suggests the spiralling may correlate with tension vectors due to intraocular pressure (Nejad et al., 2015), but there are no empirical studies suggesting how it might occur in vivo.

The drivers for directed cell migration across the ocular surface are poorly known. Various mechanisms have been suggested, including chemorepulsion from vascular components at the limbus, population pressure from the highly proliferative peripheral corneal epithelium, biased cell loss near the corneal centre, electrical currents caused by ion leakage through abrasion, and contact-mediated cues (Bron and Tripathi, 1973; Buck, 1985; Sharma and Coles, 1989; Lemp and Mathers, 1989; Lavker et al., 1991; Wolosin et al., 2000; McCaig et al., 2005; Foster et al., 2014). Patterns of cell migration are known to be disrupted, at least temporarily, by wounding, and in disease states (Collinson et al., 2004; Mort et al., 2009; Richardson et al., 2018). The planar polarity protein, Vangl2, is required for normal in vivo patterns of cell migration, suggesting aligned cell-cell interactions (Findlay et al., 2016). Corneal epithelial cell migration can be directed in a contact-mediated manner by substrate topography, and an integrin-mediated, cAMP-dependent interaction with their basement membrane may maintain some degree of radial orientation of migrating cells in vitro but is not sufficient to push all cells in a centripetal direction (Walczysko et al., 2016). A major centripetal driver may be a step-change increase in substrate compliance at the limbal-corneal boundary (from soft to hard substrate) causing limbal cells to make an initial durotactic movement into the cornea (Walczysko et al., 2016; Gouveia et al., 2019).

However, it has not been shown experimentally whether these genetic couplings or mechanical and extrinsic directional cues would be enough to produce stable patterns of long-distance migration in vivo. Understanding the coordinated, highly patterned, constrained collective migration patterns of the corneal epithelial cell sheet in vivo is a complex three-dimensional biological problem requiring an integrated, systems-based resolution that goes beyond the individual contributing elements.

Such an in silico model of in vivo ocular surface maintenance that produces radial striping is not conceptually difficult, and coordinate-based models based on cell proliferation dynamics and/or cell density have been proposed (Lobo et al., 2016; Moraki et al., 2019). They, however, have unsatisfactory aspects from a biological and predictive perspective. A good model needs to explain, not simply describe, both the formation and maintenance of stripes and the central spiral of cell movement, within a confluent cell sheet. It needs to be formulated in three dimensions to account for the curvature of the cornea, to be quantitative and scaleable, in both size and simulation times, and to apply to corneas of different shapes; importantly, it needs realistic local dynamics of physical cell-cell interactions, incorporate active cell migration and planar alignment, and to make no assumptions about cell behaviour or guidance that are not rooted in observations.

The physics of active or living matter seeks to understand the collective properties of individually motile agents, such as birds, fish, bacteria or eukaryotic cells (Vicsek and Zafeiris, 2012; Marchetti et al., 2013; Needleman and Dogic, 2017). One of its hallmarks is a state of collective migration, or ‘flocking’ that generally appears when agents locally coordinate their motion, without the need for a leader or guidance cue (Vicsek et al., 1995; Giardina, 2008). Such flocking collectives where the directions of motion are all aligned with each other are known as polar active materials. The apparent flocking behaviour of cells across the corneal surface, together with known behaviours of individual cells, such as self-propulsion, local coordination, contact-mediated cell proliferation and differentiation, leads us to hypothesise that a collective migration model rooted in the physics of soft active matter (Alert and Trepat, 2020) can explain the observed in vivo migration patterns.

Within the framework of active matter, it has become clear that the system’s topology plays a prominent role. For example, a sphere is a closed surface with no holes, with an associated topological invariant, the Euler characteristic χ = 2, loosely speaking, a number that describes an object’s shape regardless of how it is bent. A mathematical consequence is that the velocity profile of any flow confined to the surface of a sphere has to vanish at least at one point. Flocks, therefore, cannot maintain fully aligned velocities everywhere on the surface, and instead create points with locally mismatched directions. Near such ‘singular points’, if one loops over a contour encompassing it, the flow field winds an integer number of times (Do Carmo, 2016). Such points are known as topological defects with topological charge c (Mermin, 1979). Pictorially, +1 defects are circular objects corresponding to inward, outward or spiral flows around the defect core, while −1 defects have four-leaf clover symmetry, with flow inwards along two directions and outwards along the orthogonal ones (Alexander et al., 2012). The Poincaré-Hopf theorem requires that the total topological charge equals the Euler number, i.e. ∑_i c_i = χ (Do Carmo, 2016). Flocks of active particles on a sphere then typically generate patterns consisting of a circling band and two vortex-like +1 topological defects at the poles (Sknepnek and Henkes, 2015; Shankar et al., 2017; Hsu et al., 2022).

The cornea, with the limbal boundary, is a spherical cap, i.e. it has the same topology as a disk. In this topology, with Euler characteristic χ = 1, if the flows close to the boundary are uniform as for the cornea, we expect at least one topological defect with charge +1 (Stein, 1979). We propose that this is the vortex pattern corresponding to the central spiral of the cornea. Indeed, in an in vitro collective cell system with a disk topology, such a spiral defect was described (Guillamat et al., 2022). Pictorially, the migration pattern for the cornea is akin to water flowing in from a source and then circling down a sink, except that we can create and extrude cells everywhere within the sheet. Additionally, biological development can also mould the surfaces themselves and then generate deformations that couple to topological defects – as has recently been shown for hydra (Maroudas-Sacks et al., 2021), and argued to play an important role for organoids (Rozman et al., 2020).

In this paper, we build a quantitative model of spiral flocking on the cornea. We previously showed that the swirling migration patterns of confluent corneal epithelial cell monolayers in vitro can be understood using a model of active matter at high densities (Henkes et al., 2020), and that one can include mechanical coordination between cells as polar alignment (Saraswathibhatla et al., 2021). The model treats each cell as a soft particle, able to self-propel, align its planar polarity with, and exert elastic forces on its neighbours.

These features are accounted for by a set of equations of motion that are evolved in time to simulate swirling spatiotemporally correlated patterns of cell motion observed in in vitro experiments. The model is quantitative, and its parameters can be directly derived from experimental observation. Here, we extend this model of corneal epithelial cell migration to a three-dimensional constrained geometry that recapitulates the mouse cornea, with stem cells at the limbus and proliferating TA cells within the cornea. The model produces radial striping flow patterns with a central spiral representing a +1 topological defect that replicates the spiral flow patterns of cell migration observed in vivo. We link flocking spiral formation to effective resurfacing of the cornea, and reformulate the XYZ hypothesis of corneal maintenance (Thoft and Friend, 1983) in physical terms. We are thus able to demonstrate the central role of confinement and curved geometry in corneal patterning. This model represents the first comprehensive, predictive tool for understanding the maintenance of the ocular surface in health and disease.

Results

Experimental flow fields

Figure 1a shows a representative X-linked LacZ-transgenic adult female mouse cornea, with a pin-wheel spiral pattern of β-galactosidase mosaicism due to the flow of cells from patches of blue and white stem cells at the periphery (limbus) to the corneal centre throughout adult life (Collinson et al., 2002). The patterns of cell migration are reflected by the projection of sensory axons between the cells in the basal layer, which largely follow the flow of epithelial cells (Figure 1b). Both axons and cells form a smooth flow pattern with a vortex at the corneal centre. This vortex is an example of a topological defect (i.e. a point where the orientational order of cell velocity vectors vanishes) with a +1 topological charge. We hypothesise that the flow pattern and the central topological defect can be explained as an emergent behaviour of a collective of self-propelled epithelial cells closely packed on a hemisphere. As a differential hypothesis, epithelial cells on a substrate have been described as an active nematic material (Saw et al., 2017; Balasubramaniam et al., 2021), where elongated cells or the stresses of the cells pulling on each other generate topological defects with ±1/2 charges (topological defects with ‘half-integer’ charges are allowed if the field has nematic, i.e. head-tail symmetry). Despite the compatible elongated shape of the axons, we do not observe such defects, which for a +1/2 defect resemble a fold in a material made of long, thin objects. Instead, in rare occurrences, additional pairs of +1 and −1 defects occur in axon patterning (Figure 1c). We infer that here the polar, crawling motility of the cells dominates over cell-cell stresses with nematic symmetry. Therefore, cell motion in a cornea is an example of a dense polar soft active matter system (Marchetti et al., 2013; Alert and Trepat, 2020).

Figure 1.

To test this hypothesis, we reconstructed the epithelial cell flow fields across the corneal surface, by dissecting, flattening and imaging corneas as described in Materials and Methods. We manually traced boundaries between blue and white stripes from images of n = 8 dissected wild-type adult XLacZ corneas in Inkscape (Figure 9). We wrote custom analysis software in MATLAB to digitally reconstruct the cornea surface in 3D and then inferred the local cell flow directions from the stripe edges. Finally, we minimised an XY model of polar alignment to interpolate between known directions. This enabled us to infer the direction of velocity vectors everywhere on the cornea. In Figure 2a, we show a two-dimensional projection of the flow pattern along the curved corneal surface inferred from one of the corneas (experiment 1). For all corneas, striping patterns are shown in Figure 2 – Figure supplement 1, and inferred flow patterns are shown in Figure 2 – Figure supplement 2.

Figure 2.

Figure 2–Figure supplement 1.

Figure 2–Figure supplement 2.

We quantified the flow direction by the flow angle α, determined as follows. Owing to the nearly spherical shape of the eye, we used spherical coordinates to identify points on the cornea. We placed the origin at the centre of the eye, and the location of the +1 topological defect represents the ‘north pole’, such that the z−axis passes through it (Figure 2c). The angle θ with θ = 0° being the spiral centre then measured the angular distance from the pole (like an inverse latitude), where θ ≈ 70° corresponds to the limbal region at the periphery of the cornea. We measured the azimuthal direction by the angle ϕ, with values between 0° and 360°, with an arbitrary origin as the cornea is radially symmetric.

The spherical coordinates allowed us to define local unit reference vectors ê_θ and ê_ϕ in the tangent plane of the corneal surface. Then we defined α as the angle between −ê_θ and the direction of the velocity vector v (Figure 2c). Here the minus sign ensures that |α| < 90° for an inward flow, as ê_θ points outwards by convention.

The flow pattern is clearly highly symmetric along the azimuthal direction. Therefore, we averaged α over the angle ϕ to extract the characteristic spiral angle flow profile α(θ), as a function of the polar angle θ (Figure 2b). Although both clockwise and counterclockwise patterns were approximately equally likely, we flipped the sign of angle α on clockwise-rotating corneas for visual clarity. This enabled us to identify a reproducible pattern for all corneas: Flow is straight inward near the corneal edge with α(θ = 70°) ≈ 0 and stays such until at an angle θ ≈ 15° − 25°, when α(θ) sharply increases, corresponding to a systematic tilt of v, i.e. the onset of a spiral. In all cases, the spiral angle further continuously tightened until reaching α(θ = 70°) ≈ 90° at the corneal centre, corresponding to flow in the azimuthal direction.

Finally, we identified topological defects in the flow field by computing angle differences over the loops defined by a Delaunay triangulation (Zapotocky et al., 1995; Henkes et al., 2018) (Figure 2d). We always found a single +1 topological defect surrounded by a spiral, located near the centre of the cornea (Figure 2a and Figure 2 – Figure Supplement 2).

Excluding prepatterning as the stripe formation mechanisms

Patterns of radial striping are potentially explicable if the corneal epithelial cells are aligned by a global patterning mechanism that, e.g. biases cell migration or aligns cell division axes to point towards the cornea centre.

To test whether the cell migration direction is biased, we used the centrosome. The position of the centrosome has been associated with cell polarity, motility, and division (Nigg and Raff, 2009). Intracellular movement of the centrosome has been shown both in vivo and in wound assays in vitro (Werner et al., 2017; Blitzer et al., 2011). Hence, the position of the centrosome relative to the nucleus can act as a proxy measurement for the cell polarity and, potentially, the direction of migration. To assay whether the corneal epithelial cells have an intrinsic directional polarity, we performed immunohistochemistry to visualise pericentrin in flat-mounted corneas and measured the angle of the centrosome-nuclear axis ϕ with respect to the corneal centre (6 eyes, 879 cells) (Figure 3a,b). Centrosomes were either behind the nucleus (0°−180°) (42±1%), in front of the nucleus (51 ± 2 %), or directly above or below the nucleus with no orientation (6 ± 1%). Statistically, there was a slight but significant tendency for centrosomes to locate in front of the nucleus (p = 0.0047), consistent with the observed mean flow towards the cornea centre. However, the overall scattered pattern of centrosomes with respect to the nucleus is more consistent with short-term changes in migration directions associated with the local cell environment than with any global patterning.

Figure 3.

To investigate whether radially-aligned cell division axes contributed to observed striping patterns, we identified mitotic figures in DAPI-stained flat-mounted corneas from adult mice as shown in Figure 3. We manually measured the angle ϕ of each mitotic division with respect to the centre of the cornea and binned it as ‘radial division’ (0° − 45°, 135° − 180°) or ‘tangential division’ (45° − 135°). Fig 3c shows that we did not observe any significant difference in the proportions of radial (47 ± 5 %) and tangential (53 ± 5%) divisions (9 eyes, 376 dividing cells). Therefore, we did not find any evidence for radial striping being caused by biased axes of cell division.

Soft active particle model

With the prepatterning mechanism being unlikely, we constructed an agent-based model that is able to explain the observed radial migration and spiral patterns as emergent phenomena rooted in the physics of soft active matter. To set the stage before analysing the effects of curvature, we first developed a planar model. As discussed in the next section, the planar model can be parametrised using corneal explants and monolayers of primary mouse corneal epithelial cells.

The corneal epithelium is stratified, but the active migration is performed solely by the basal layer of cells, with the more apical layers essentially riding passively on the basal flow. For simplicity, therefore we modelled the corneal epithelium as a single sheet of cells, themselves modelled as soft disks. Each disk has a soft repulsive core, which sets the cell radius σ and accounts for elastic repulsion between overlapping neighbours. In addition, cell-cell adhesion with their neighbours is modelled via a short-range attractive potential. We introduced cell migration on the basement membrane covering the stroma through an active self-propulsion force that randomly changes its direction over time, reorienting the direction of a cell’s migration. We also allowed cells to divide or be extruded with rates d_z and a, respectively. In this case, ‘extrusion’ represents the differentiation of basal epithelial cells and their movement into the apical layers of the epithelium in vivo. We modelled cell division as density-dependent, where the rate of cell division d_z = d₀(1 − z/z_max) depends on the local cell density through the number of nearest neighbours, z, and allows the tissue to reach homeostasis when the rate of production of cells by division equals the rate of extrusion, and d₀ and z_max are parameters that set the actual numerical value of d_z (Matoz-Fernandez et al., 2017). On the other hand, we assumed the extrusion rate a to be fixed and independent of the local cell density.

Furthermore, cells can coordinate their migration direction through a combination of chemical or mechanical cues. For corneal cells, the detailed mechanochemistry is unknown. However, for the extensively studied MDCK epithelial cell lines, many types of coordination have been proposed (Alert and Trepat, 2020), such as plithotaxis, i.e. reorientation of the migration direction with the mechanical stress that it experiences (Trepat and Fredberg, 2011). Here, in the absence of information about cellular stresses and the observation of a directional global flow, we assumed two simultaneously acting mechanisms by which a cell can reorient its direction of migration, i.e. reorientation in response to the local force feedback (Henkes et al., 2011; Petrolli et al., 2019; Baconnier et al., 2024) and alignment with its nearest neighbours (Saraswathibhatla et al., 2021). If sufficiently strong, both can cause cells to ‘flock’, i.e. to collectively migrate in the same direction akin to a flock of birds (Giardina, 2008; Vicsek and Zafeiris, 2012).

We describe the active crawling force F_act of the cell labelled by index i with a magnitude |F_act| = ζ ν₀ and a unit-length direction vector , where ϑ_i is the angle between and the x−axis of the cell sheet plane, spanned by unit-length vectors ê_x and ê_y. The equations of motion for the position vector r_i and the angle ϑ_i are

where ν₀ is the migration speed, ζ is the friction coefficient with the substrate, and the pair forces F_ij derive from our short-range attractive-repulsive potential with stiffness k, and J as the alignment strength (see Materials and Methods). The dynamics of the active migration direction is determined by alignment with the cell’s velocity angle Ψ_i (first term in square brackets) and pair alignment with the cell’s nearest neighbours (second term in square brackets), with added rotational white noise η_i with time correlation ⟨η_i(t)η_j(t^′) ⟩ = 2D_rδ_ij δ(t − t^′). D_r is the rotation diffusion constant and it sets the cell motion persistence time τ = 1/D_r, i.e. the typical time a cell takes to reorient its migration direction. Please see Materials and Methods and (Sknepnek and Henkes, 2015; Henkes et al., 2020) for full model details; we used our simulation package (SAMoS, 2024) to carry out the computations and our matching analysis package (SAMoSA, 2024) to analyse the data. The definition of parameters and their values are given in Table 1.

Table 1.

We briefly comment on the choice of a model with polar symmetry. Epithelial cells can generate mechanical stresses, which enable them to exert forces on the substrate and their neighbours. In many cases, those forces are dipolar and thus result from stresses with nematic (i.e. head-tail) symmetry. Among other things, the nematic symmetry leads to the appearance of ±1/2 topological defects in the cell-shape texture and active turbulence (Shankar et al., 2017; Blanch-Mercader et al., 2018; Duclos et al., 2018; Mueller et al., 2019), phenomena that can be successfully explained using models for active nematics (Doostmohammadi et al., 2018). Dipolar active force cannot cause a single cell to migrate, and the flow is an emergent collective phenomenon. Epithelial cells can, however, also generate polar (i.e. directional) self-propulsion forces that can cause a single cell to crawl. Forces with polar symmetry can explain, e.g. observed velocity correlations in confluent epithelial sheets (Henkes et al., 2020). This highlights that cells can generate and transmit complex mechanical stresses (Armengol-Collado et al., 2023) that can lead to rich and only partly understood collective behaviours.

Although we do not exclude the possibility that active nematic stresses are present in corneal cells, as mentioned already, the absence of half-integer topological defects suggests that predominant mechanical forces have polar symmetry. This makes a model with self-propelled agents subject to polar alignment appropriate to describe cell migration of corneal cells. We stress that the term ‘polar’ refers to the symmetry of active mechanical forces and not to the internal cell polarity.

Extracting model parameters from experiment

To parametrise the model, we experimentally measured cell sizes, velocities, flow patterns, and cell division orientation for mouse corneal epithelial cells.

We determined cell migration speeds by time-lapse analysis in cultured mouse corneal explants (labelled ‘explants’), and also in monolayers of primary mouse corneal epithelial cells cultured on a substrate as described in (Hazlett et al., 1996) (labelled ‘plastic’) (see Materials and Methods), with the results shown in Figure 4. The in vivo epithelium moves over a basement membrane on a collagenous stroma, a substrate with different mechanical properties compared to tissue culture plastic. This is consistent with the emerging understanding that there is a trade-off of the active stresses arising from the cytoskeleton between contractility at the apical surface and traction with the substrate (Ladoux and Mège, 2017). Therefore, in this study, we adjusted the main size parameter, cell radius and migration speed ν₀ of modelled cells to match the explant data derived empirically (summarised in Table 1).

Figure 4.

We estimated the mean planar radius of basal corneal epithelial cells in vivo by microscopy in flat-mount and tissue sections to be σ ≈ 5 − 6 µm. Cells in the tissue culture on ‘plastic’ were larger, with an estimated radius of σ ≈ 11 µm. Therefore, we used σ = 5 µm for simulations of explants and, subsequently, the full cornea simulations, and σ = 11 µm for simulating monolayers grown on ‘plastic’.

Next, we investigated the cell migration patterns using particle-image velocimetry (PIV) (Adrian, 1991) on phase-contrast images of the confluent cell layers. In comparison to our previous PIV analysis of HCE-S cell monolayers that found the root mean squared velocity to be 12± 2 µm/h (Henkes et al., 2020), ν_rms was higher (15 − 25 µm/h) for mouse monolayers in tissue culture (‘plastic’), but lower (2 − 6 µm/h) for mouse epithelial cells in whole corneal explants moving on the endogenous stroma (‘explant’), see Figure 4a. In combination with previous measurements of corneal epithelial wound healing migration of up to 60 µm/h in vitro and in vivo (Leiper et al., 2006; Garcia et al., 2015; Walczysko et al., 2016), these data set the range of cell velocities biologically permissible for the model. Examples of this collective cell motion together with the velocity field inferred from PIV are shown in supplementary videos S1 and S2 for explant and plastic conditions, respectively.

The mean cell cycle time in basal corneal mouse epithelial cells has been measured previously to be ≈ 3.5 − 5.0 days (Sagga et al., 2018), with individual basal epithelial cells capable of division once every 2.5 days (Lehrer et al., 1998), compared to 1 - 1.5 days in HCE-S cells (Henkes et al., 2020). To create a quantitative model of a homeostatic system where the rate of cell loss matches the rate of cell division, we first set the cell extrusion rate a to the inverse of the observed cell cycle time, resulting in a = 0.02 h⁻¹ (i.e., 50 hours between divisions). This rate was then applied in both ‘plastic’ and ‘explant’ matching simulations. We then chose d_z to be compatible with the observed confluent tissue densities (Matoz-Fernandez et al., 2017). The cell-cell interaction potential strength k and the cell-substrate friction coefficient ζ are difficult to estimate individually, but fortunately, they enter the model only through their ratio. We, therefore, used k/ζ = 55 h⁻¹ extracted in vitro for HCE-S cells (Henkes et al., 2020) as the starting value of the fitting process. To match the experimental field of view, we simulated a system of size 867 × 662 µm². We estimated the cell migration speed in isolation ν₀ and the alignment strength J through measuring ν_rms and the flocking order parameter (Vicsek et al., 1995), which varies between 0 (random migration directions) and 1 (full flock migration). In Figure 4a-b, we show ν_rms and Φ_v for experiments and simulations. As mentioned above, cells on plastic are considerably faster than explants, and we used ν₀ = 30 µm/h (ν₀ = 3 µm/h) for simulations of experiments on ‘plastic’ (‘explants’). Qualitatively, the observed difference in speed is due to explant cells being tall and narrow with few focal adhesions, while cells on plastic are spread out with strong cytoskeletal coupling to the substrate. We found that the observed values of the flocking parameter, Φ_v = 0.3 − 0.7, were incompatible with J = 0 h⁻¹ used to model the HCE-S data in (Henkes et al., 2020). Instead, we found a good simultaneous match to ν_rms and Φ_v in the range J = 0.1 − 0.2 h⁻¹, as befits the considerable variance between experimental repeats, and consistent with results for MDCK sheets (Saraswathibhatla et al., 2021). This range corresponds to the region of the flocking phase transition in our model, a first hint that this transition provides the mechanism responsible for the stripe formation process in vivo.

As we show in videos S1 and S2, cells in the sheet moved coherently over many cell sizes. We finished the analysis by matching cell velocity-velocity correlations in space and time. The time correlation of cell velocities is defined as ⟨v(t) ⋅ v(0) ⟩, where ⟨… ⟩ is the average over all cells and all time points along cell trajectories that are time t apart, and it measures how long it takes cells to change their migration direction, as shown in Figure 4e-f for ‘explants’ (top) and ‘plastic’ (bottom). In both cases, while variable, the experimental correlation decays over several hours, with longer correlation times for the tissue culture plastic. As before, we chose a simulated value of the cell persistence time τ = 2 h, which together with J = 0.1 − 0.2 h⁻¹ brackets the values extracted from experiments. The simulations above the flocking threshold showed a long-time correlation plateau as do most experiments.

Finally, we obtained the Fourier velocity correlation function ⟨ |v(q) |²⟩ as average of the square modulus of the spatial Fourier transform of the velocity field, . As shown in Figure 4c-d, it measures which wavelengths contribute to cell motion, and is characterised by a decay beyond a spatial correlation length , where μ is the shear modulus, which is set by the particle elastic stiffness through k/ζ (Henkes et al., 2020).

For both ‘plastic’ and ‘explant’ data, we obtained good agreement for the low modes q ≲ 2π/ξ where the continuum approach is valid, giving a final value of k/ζ = 30 h⁻¹, similar to other epithelial sheets (Henkes et al., 2020; Saraswathibhatla et al., 2021). Supplementary videos S3 and S4 show our matched simulations at J = 0.15 h⁻¹ for explant and plastic conditions, respectively.

Corneal surface model

To combine active cell dynamics with the realistic three-dimensional corneal shape, we modelled the cornea as a spherical cap of radius R = 100 − 1500 µm and a 70° opening, consistent with experimental measurements. The corneal epithelial cells, i.e. the TA cells derived from stem cell divisions in the limbus, were modelled as soft particles able to migrate on the curved surface, divide, and be removed. To ensure that the cells remain confined to the surface throughout the simulation, we projected the velocity vectors onto the local tangent planes at each time step. The equations of motion, therefore, contain the same term as in the flat case, eqn. (1), with additional projection operators (Sknepnek and Henkes, 2015) (see Materials and Methods). The model parameters, listed as ‘explants’ in Table 1, were determined as discussed in the previous section. At the bottom boundary, we modelled the limbal stem cell niche as a single layer of continuously dividing cells, a simplification that takes into account both dividing stem cells and a population of TA cells that may divide before leaving the limbus. In the mouse limbal epithelium, slow-cycling stem cells, comprising approximately 25% of the cell population, were previously determined to divide once every 10 - 20 days with the other 75% of rapidly dividing cells having a cell cycle time of 3 - 3.5 days (Sagga et al., 2018). To match the above and the empirically determined corneal resurfacing time of 14 days at full-size, we chose a density-independent division rate of d_L = 0.4 h⁻¹ for the limbal ring.

Finally, we included a double ring of tightly spaced stiff cells as a barrier to constrain motion to the interior of the cornea area. This reflects experimental observations that the progeny of limbal epithelial stem cells never migrate outwards into the conjunctiva. We recapitulated the lineage tracing provided by XLacZ staining by separating the cells into ‘blue’ and ‘white’ types. Each ‘blue’ limbal stem cell gives rise to ‘blue’ TA cells, which can also only divide further into more ‘blue’ cells, as in vivo. We used the same dynamics for the ‘white’ cells. Blue and white cells were identical in their interactions and they could mix freely. We set the base for the stripes to 12 alternating regions of fixed stem cells along the limbal stem cell line (Figure 5). For the initial set-up of simulations, we randomly seeded cells with a 50:50 ‘blue’:’white’ mixture on the corneal surface. For the largest cornea radii R_c = 1000 µm and R_c = 1500 µm, due to computational limitations, we started with an empty corneal surface. As above, we implemented the model through custom configuration files and initialisation scripts in our package (SAMoS, 2024). We simulated between 1400 cells at R_c = 100 µm, and 241 000 cells at R_c = 1500 µm. We ran the simulations using an Euler-Maruyama algorithm for 500 000 simulation time steps at dt = 0.001 h, corresponding to 20.8 days of development. Please see supplementary video S5 for a R = 100 simulation corresponding to Figure 5.

Figure 5.

Emergence of the spiral pattern

The corneal spiral is not present from birth (Figure 6a): Young XLacZ mice exhibit a random patch-work arrangement of blue and white clones from mosaic LacZ expression of cells. Cell loss stimulates the immigration of cells into the cornea from the limbal epithelium, such that by adulthood the patterns resolve as the spiral radial arrangement of clones. The transition from a primarily patchwork to a primarily striped pattern takes 3 - 4 weeks, with intermediate patterns with stripes emerging from the periphery, retaining a patchwork centre. We were able to successfully reproduce both this transition and the spiral striping pattern using numerical simulations of the corneal surface model. In Figure 6b-e, we show the time evolution of a simulated cornea of radius R = 700 µm. As can be seen in the second row (Figure 6b and accompanying supplementary video S6), the initially randomly distributed ‘blue’ and ‘white’ cell pattern was gradually replaced by a striped pattern made of cells derived from the limbus. This was not simply a case of cells being replaced since the change happened over several generations of cells, with the overwhelming majority of striping pattern cells being derived from other corneal epithelial TA cells. By day 20, a cornea with a stable central spiral emerged. This pattern of replacement, including the time scale over which it occurred, was very similar to the corneal evolution observed in neonatal mice shown in Figure 6a. We can understand the emergence of the spiral flow pattern by examining the temporally averaged (or hydrodynamic) velocity field (Figure 6c and supplementary video S7), where we tracked the displacement Δr_i(t) = r_i(t + Δ/2) − r_i(t − Δ/2) of each surviving cell i over Δt = 2.5 h, giving . This is distinct from the instantaneous cell migration speed v_i which was more random, consistent with the nearly random cell orientations seen in experiments (cf. Figure 3b). Figure 10 in the Materials and Methods shows how the smooth hydrodynamic flow pattern emerges as we increase the averaging time Δt. The hydrodynamic velocity field gradually forms a spiral pattern with flow parallel to the blue/white stripe edges of the spiral. This confirmed our experimental approach of inferring the flow field direction from the stripe edges.

Figure 6.

We next quantified our observations by tracking the +1 and −1 topological defects of v^H, shown in white and black, respectively in Figure 6c. As discussed above, topological defects correspond to regions where the flow field is discontinuous, with +1 defects being circular objects with spiral flows around the defect core, while −1 defects have four-fold inward/outward symmetry. We computed topological defect charges c using the same discrete sum around Delaunay triangles of a tesselation constructed from the particle positions as for the experiment (Figure 2d), and further merged the charges of defects recursively within 15 µm, to avoid regions of small-scale disordered flow dominating our measurements (Henkes et al., 2018; Prathyusha et al., 2018). Please see Materials and Methods for full details; both the computation of v^H and the topological defects were carried out in our analysis package (SAMoSA, 2024). Since the cornea is topologically a disk, i.e. a surface with a boundary, in general, the total topological charge is influenced by that boundary. However, due to flows being uniformly inward close to the limbal region, the sum of the topological charges is expected to be +1, i.e. N₊₁ − N₋₁ = 1 (Stein, 1979). We indeed observed that as the simulated cornea matures, +1 and −1 defects were initially very numerous and then reduced until only one central, +1 topological spiral defect remained, which corresponded to the centre of the blue-white spiral. In Figure 6d, we show the polar angle of the ±1 defects as a function of time. One can clearly see the mutual annihilation of +1 and −1 defects such that by day 8, only the central +1 defect remained. Here the top graph corresponds to the defects of the flow field shown in Figure 6b, and the bottom graph corresponds to the director field, i.e. the field of directions of migration of individual cells. Figure 6e shows the spatial trajectory of the central defect superimposed on the flow pattern, again for velocity and director defects. Both show that the central defect is almost stationary and nearly identical for velocity and director defects. These measures all show that the flow field pattern is a direct result of the polar alignment of the cells, i.e. the spiral pattern indeed emerges from the interaction between flocking and hemispherical geometry.

Key parameters of simulated spiral corneas

We next sought to understand the properties of the spiral flow field and how it compared to the experimental observations. In Figure 7, we show the dependence of the flow patterns on cornea radius R and polar alignment strength J. Figure 7 – Figure supplement 1 and 2, respectively, show simulation snapshots of all simulated values of J and R.

Figure 7.

Figure 7–Figure supplement 1.

Figure 7–Figure supplement 2.

From these simulations, we made a key observation: polar alignment is necessary for spiral formation, and the cornea is not effectively resurfaced in its absence. For example, in the bottom grey oval in Figure 7a, we show the flow field (left) and a snapshot of cell configuration (right) for J = 0.025 h⁻¹ and R = 500 µm after 20 days. The flow field remained disordered with numerous topological defects and no emerging spiral. Supplementary videos S8 and S9 show the time evolution (or lack theoreof) of this process. Conversely, for J = 0.2 h⁻¹ (top grey oval) we obtained a tighter spiral with a well-developed spiral flow field. Supplementary videos S10 and S11 show the time evolution, including the more rapid appearance of the single central defect.

We quantified these observations by computing azimuthally averaged profiles using spherical coordinates with the central defect along the z axis, using the same method as for the experiment (Figure 2c). In Figures 7b and 7d we, respectively, show the spiral angle α and the mean flow velocity ν^H as a function of polar angle θ. The former is the same measure as in Figure 2b obtained in experiments. We found a transition between a disordered state with small ν^H in the interior (and thus no effective resurfacing) and no developed spiral for J ≤ 0.05 h⁻¹ and increasingly effective resurfacing for J ≥ 0.075 h⁻¹, in the same region as the flocking transition in the flat case. Both spiral shape and flow velocity saturated for J ≈ 0.2 h⁻¹, in the limit where cells fully flock together. Thus, we predict that the emergence of the topological spiral flow due to polar alignment is essential for a healthy cornea.

We observed that a spiral emerges for all simulated cornea radii (R = 100 − 1500 µm), showing the robustness of the flocking mechanism. However, the spiral shape is strongly radius-dependent. In Figure 7c, we show the spiral angle α with the inward direction as a function of polar angle θ for different values of R. We see that there is a systematic change from a pattern of a wide spiral (large α at most θ) at small values of R to a flux that is nearly straight inwards on the outer areas of the cornea (small α at large θ) followed by an increasingly sharper spiral pattern near the centre. In all cases, we had α(0) → 90°, i.e. a spiral that becomes a vortex at the centre. The flow field ν^H (Figure 7e) shows that all corneas except the smallest reached a steady flow velocity of approximately 2 µm/h on the corneal surface. Near the centre, however, ν^H drops to 0 approximately linearly, with an increasingly sharper gradient as the corneal size increases.

The XYZ model revisited

A mainstay of understanding of ocular surface maintenance has been the XYZ hypothesis which states that the amount of cells flowing into the cornea from the limbal stem cell reservoir is exactly balanced by the difference between the amounts of TA cells being born and leaving the cornea (Thoft and Friend, 1983). In other words, the XYZ hypothesis posits that the cornea is in a steady state with a constant but nonzero flux of cells from the limbus inward to the tip of the cornea. We can test and expand on this hypothesis by writing a continuum conservation law for cell numbers, and apply it to the spiral flow profiles, linking cell flux v^H, the spiral angle of the cornea α, and the difference between birth and extrusion rates.

The XYZ hypothesis states that the number of cells in a patch of the cornea is, on average, constant. If we denote the number of cells in a patch as N, N_in cells flowing into the patch, N_out cells flowing out of the patch, N_d cells dividing, and N_a cells that are extruded, we obtain the following balance relation:

This equation could be reformulated in terms of the density of cells, ρ = N/S, where S is the patch area. The flux J = ρv^H is the mean number of cells with velocity v^H flowing through the boundary of S per unit time. Similarly, we can define d(ρ) and a(ρ), the division and death rates per cell. Eqn. (2) can then be reformulated as a flux conservation or continuity equation (see Materials and Methods) that states that the divergence of the flux at each point needs to be balanced by the difference in cells being extruded and those that are dividing at the same point, i.e.

It is possible to link this equation with the XYZ hypothesis as follows. We first assume that the flow has a radially symmetric spiral pattern, i.e. v^H = v(θ) cos α(θ)e_θ +sin α(θ)e_ϕ, where v(θ) = |v^H| is its magnitude and α(θ) is the angle between v_H and −e_θ. The difference between division and extrusion can be combined into the density-dependent net loss rate A(ρ) = [a(ρ) − d(ρ)]. We further assume that the deviations of ρ from its mean are small, though we retain their influence on A(θ) ≡ A(ρ(θ)). Written in spherical coordinates with the defect along z again as in Figure 2d, eqn. (3) becomes (see Materials and Methods),

This equation links the spiral angle profile α(θ), shown in Figures 7b and 7c with its velocity profile v(θ) (Figures 7d and 7e) through the net loss rate A(θ) which we show for the simulations in Figures 8b-d. We immediately see that A(θ) has a distinct emerging angular profile, despite cell division and extrusion in the model incorporating no spatial information. In particular, we note that for small values of the alignment strength J, where the system cannot form a spiral, there is a strongly positive A at large θ, near the limbus. This corresponds to a ‘traffic jam’ of cells stuck near their place of birth that are extruded before being able to migrate up the surface. In contrast, for spiral forming configurations, A < 0 near the limbus, signalling an enhanced birth rate and lower extrusion rate is balanced by a higher extrusion rate near the centre. Additionally, these A profiles are also explicitly dependent on the corneal radius.

Figure 8.

Figure 9.

Figure 10.

We write an equation for the inward flux F (θ) = v(θ) cos(α(θ)), i.e. the radially inward component of the velocity. Integrating over the polar angle θ from the centre outwards then gives a relationship between the flux at angle θ and the integrated loss rate inwards of this point,

In Figures 8a and 8c, we show the predicted F (θ)/R (dashed lines) compared to the empirically computed flux in the simulations, showing a very good fit. We can clearly see that the flux is also a nontrivial function of both θ and cornea radius R, explaining the change in spiral shape that we observed as a function of angle θ and radius R. This dependence on R is both explicit in eqn. 5 and due to the nontrivial form of A(θ), a quantity that we cannot currently measure, but that may at some point become accessible experimentally.

Since we have established that our model is a polar active material, it should in principle be possible to derive continuum hydrodynamic equations for it, in a similar approach taken for in vitro systems (Guillamat et al., 2022). The delicate density-dependent balance that links division, death and velocity magnitude in our system however complicates this: the appropriate symmetry class are the ‘Malthusian’ flocks (Toner, 2012), which include division and death. However, the strong approximations inherent in this approach, most notably that A = 0 and that v is constant, independent of cell density, make them insufficiently detailed to explain the corneal flux. A future fully quantitative hydrodynamic description of the cornea will require a more sophisticated approach.

Discussion

Conclusion

We introduced an in silico agent-based model that describes the complex collective motion pattern observed in the corneal epithelium as an emergent phenomenon. The model only requires simple assumptions about local cell behaviours and neither relies on any organ-scale pre-patterning or regulation nor requires the imposition of any external guidance cues or any shear and tensile strains that may otherwise cause the sliding of epithelial cells in logarithmic vortices (Mohammad Nejad et al., 2015). In vivo, cell migration patterns can, therefore, be quantitatively explained by a simple agent-based model commonly studied in the physics of soft active matter systems. This is the first corneal model to successfully reproduce the central vortex motion of epithelial cells around the corneal centre in simulations. The key outcomes from this study are that to recapitulate experimentally observed patterns of cell migration in vivo, cells need to: 1) be able to actively migrate, 2) respond to mechanical pressure due to the population of other cells, 3) align their migration direction to that of their neighbours (i.e. ‘go with the flow’), 4) have density-dependent control of proliferation, and 5) maintain a high rate of limbal cell proliferation. Then the emerging pattern is robust to noise, i.e. stochastic components in cell motion and proliferation. Furthermore, we were able to determine model parameters from experiments, which allowed for quantitative comparison between the in silico model and the biological system. Changing model parameters results in markedly different cell migration patterns that may apply to other biological scenarios, underscoring the power of the agent-based active matter models in describing collective cell behaviours.

Notably, in the cornea, where motile cells are bounded by immobile stem cells, the model predicts that at no time do cells need guidance cues to direct them towards the corneal centre. Given the step into the cornea from the limbus, which we impose on the simulations and reflects the likely durotactic pressure on cells moving from the soft limbal stroma to the stiffer corneal stroma, the radial centripetal migration that has long fascinated biologists can then be explained as being generated by a biophysical model.

The high parameter of limbal cell division that is required to produce robust radial striping is superficially at variance with the paradigm that suggests limbal epithelial stem cells to be slow-cycling. However, within the inner limbus, cycling stem cells give rise to rapidly dividing TA cells. Although up to 25% of limbal cells may be slow cycling, we previously calculated that the mean limbal epithelial cell cycle time is lower than that observed in the corneal epithelium (Sagga et al., 2018). There is, therefore, no discrepancy between the stem cell parameters imposed on our model and the real biological situation.

Analysis of the agent-based simulations data shows that at the start of the simulations, the sheet is peppered with ±1 topological defects, but those collide with each other and cancel out until a single +1 defect persists. Across several simulations, the model also found multiple +1 and -1 defects during the first 5 days, but by about simulated day 10 invariably only one defect remained. During simulations, each topological defect moved stochastically, which enabled all but one +1 defect to collide and eventually cancel out. This ‘Highlander’ model, where only one defect can survive, is a direct result of simulations but also has strong relevance to our understanding of biology. The start of the simulations represents the young mouse corneal epithelium which consists of a random unoriented patchwork of ‘blue’ and ‘white’ cells, full of topological defects, with the development of the single central vortical swirl only manifesting after days or weeks of radial cell migration (Collinson et al., 2002). The simulation suggested that it may be possible that two or more +1 discontinuities survive long-term (providing the overall topology of the cornea remains +1) if one or more -1 and +1 discontinuities happen not to collide. An in vivo example of this is seen in Figure 1.

Outlook: Models of the cornea and active matter

The physics of active matter has emerged as a powerful tool to describe many collective processes in biological systems, revealing that many complex biological processes rely on common physical mechanisms (Alert and Trepat, 2020). The cornea studied here is one such example. What, however, sets it apart from most other active matter studies of collective cell migration is the presence of curvature. Curvature alters the common notions of straight and parallel lines, which affects how active agents interact with each other if confined to move on a curved surface. It can also lead to the appearance of topological defects (e.g. the +1 defect at the cornea top), which can have important implications for the development of the organism, e.g. as was recently argued to be the case for hydra (Maroudas-Sacks et al., 2021). Compared to the efforts put into understanding active behaviour in flat spaces, the current understanding of the interplay of curvature and activity is, however, very limited. We, therefore, hope that this work will motivate further research in this direction.

Previously, coordinate-based in silico models of ocular striping have attempted to model the cornea, but did not adequately recapitulate the biological state, due to not accounting for the correct geometry and curvature of the cornea (70° hemisphere) and often ignoring factors such as cell migration. In one study, the cornea was modelled on a flat surface with population pressure from neighbouring cells being solely responsible for mobility, with cells moving from higher to lower pressure (Lobo et al., 2016). In contrast, another study used a simple stochastic approach with a hemisphere geometry, in which patterning is driven by proliferation, ignoring migration (Moraki et al., 2019). Both studies explored the effects of manipulating the amount of cell loss in different regions from the surface of the cornea through higher terminal differentiation of the TA cells (e.g. mimicking UVR damage), employed lineage tracking and the ability of all the cell types to undergo asymmetric or symmetric divisions. However, neither model incorporated all the components of cell-cell interaction, persistence, correlation and noise, and neither was able to generate a spiral migration pattern, a key feature observed in vivo.

Outlook: The cornea and cell migration

An important feature of our model is its stability over scale and corneal shape, reflecting the fact that patterns of radial migration have been observed in multiple species, including humans, with very different corneal shapes and sizes. It also allows us to directly test hypotheses for corneal homeostasis and make predictions about the corneal response to genetic, pharmacological and physical manipulation.

Eye opening in mice, at around 10 days postnatally, predates the emergence of radial striping patterns in vivo. The increase in corneal abrasion and cell loss concomitant with eye opening is presumably the trigger for limbal stem cell activation, though this has never been experimentally demonstrated, and the radial orientation of cell movement and neuronal projection is observed relatively rapidly following eye opening (Iannaccone et al., 2012). We suggest that the increase in cell loss associated with eye opening is a wound mimic which drives the rapid initiation of radial stripe patterns in young mice, and this is supported by the in silico modelling. By changing the parameters of simulated cell behaviour, it will be possible to mimic the effect of mutations in genes controlling cell velocity (e.g. Pax6) and polar alignment (e.g. Vangl2), both of which have been shown to disrupt patterns of corneal epithelial migration in vivo.

Understanding corneal epithelial cell homeostasis is essential to improve wound healing after corrective eye surgery, and to develop treatments for corneal diseases. The cornea is a paradigm system of epithelial maintenance by a population of stem cells. It is a system of choice for interdisciplinary studies of the biomechanics of cell migration because it is simple (few cell types), transparent (thus accessible to in vivo imaging), nonessential for life (thus amenable for mutation), and accessible to modelling. However, the principles and mechanisms of the cell migratory response are likely to be similar to those that underlie cell migration patterns in other tissues that are morphologically or functionally less amenable to investigation. The findings here, that a dense active matter model can recapitulate and predict patterns of epithelial cell behaviour in vivo, likely have general relevance to multiple cell migration events throughout embryogenesis and wound healing. The morphogenetic movements and epithelial-mesenchymal interactions that produce swirls of hair follicles on the head or that underlie patterns of pigmentation commonly associated with ‘Blaschko’s lines’, may also be predicted by an active matter model (Plikus and Chuong, 2004; Kücken and Newell, 2005; Tenea, 2017).

Materials and methods

Mice

We obtained all animal tissue samples from the Medical Research Facility of the University of Aberdeen in accordance with the Animal (Scientific Procedures) Act 1986. H253 mice which carry an X-linked LacZ transgene, subject to X-inactivation (Tan et al., 1993), were maintained on a CBA/Ca genetic background, breeding transgenic males with WT females to derive hemizygous, LacZ-mosaic females. We stained eyes from adult female mice at least 15 weeks old for the flow pattern and cell migration work.

Whole corneal imaging

We collected, fixed and stained LacZ-mosaic eyes with X-Gal (Collinson et al., 2002). We dissected the corneas from the eyes and added eight radial cuts from the edge to near-center such that they could be flat-mounted on a microscope slide with glycerol and a coverslip applied. Following flat-mounting the corneas were imaged using a slide scanner.

Immunohistochemistry

We used mouse corneas for immunohistochemistry to visualise cells undergoing mitosis and determine the location of the centrosome. Mouse eyes were fixed in 4% paraformaldehyde for 2 hours, washed in phosphate-buffered saline (PBS), then we dissected the corneas off and post-fixed in methanol for 1 hour at -20°C. After 3, 20 minute washes in PBS, corneas were immersed in blocking buffer (PBS, 4% normal donkey serum, 0.3% bovine serum albumen, 0.1% Triton X-100) for 1 hour at room temperature. Primary antibody [rabbit anti-phosphohistone-H3 (06-570, Millipore) or rabbit anti-pericentrin (ab448, abcam)] was dissolved 1:200 in blocking buffer and added to corneas overnight. After 3, 20 minute washes in PBS, secondary antibody [Alexa-488-conjugated donkey anti-rabbit (Molecular probes)] was dissolved 1:200 in blocking buffer and added to corneas for 2 hours at room temperature in the dark. After 3, 20 minute washes in PBS, corneas were flat-mounted as above in a drop of Vectashield with DAPI (H-1200, Vector Labs) for fluorescence microscopy.

Estimating velocity fields from stripes

The mouse cornea is not a full hemisphere so we measured the spherical cap angle and eye diameter experimentally in nine adult corneas giving an average of 71.1° ± 0.3° SEM cap angle, and 3.60 ± 0.02 mm SEM in diameter, rounded to 70° and 3.6 mm, respectively for modelling.

In vivo, we inferred three-dimensional corneal epithelial flow patterns by the proxy of tracing the boundaries between LacZ-positive and negative radial stripes in XLacZ mosaic corneas. As shown in the in silico model, despite some small-scale mixing at the boundaries, the steady-state flow field is parallel to the stripe boundaries. We proceeded as follows to extract approximate flow patterns from XLacZ corneal epithelia: first, we captured scaled images of each eye (Figure 9a). After photography, we dissected the cornea from each eye, made 8 radial partial cuts then flat-mounted the cornea for high-resolution images (Figure 9b). Using Inkscape image processing software (Inkscape Project, 2020), we drew a circle corresponding to the limbus on the image in a first colour, together with the edges of the radial cuts in a second colour, generating a background image file. In a second layer, we drew the edges of the stripes. Because different clonal stripes were often in different shades of blue, it was usually also possible to resolve clonal borders within blue stripes, with as much detail and covering as possible (Figure 9c). Different colours were used to distinguish separate flow lines, and we saved this layer as a separate image file.

These digitised corneas were used to extract the flow fields using a custom MATLAB code: We first read the background image file and extracted the limbal circle as well as the edges of the cuts. We then placed equally spaced points in 40 concentric rings in only the interior of the irregular polygon of the actual cornea defined by this method (Figure 9d). On each ring of radius r, we defined n individual points with spacing dr = 2πr/n with coordinates r_k = (r cos α_k, r sin α_k), where α_k = 2πk/n; here the k corresponding to coordinates inside the cut areas were removed. The remaining points were where we then aimed to define the flow field.

We then read the image file containing flow lines and extracted the pixels corresponding to drawn lines, mapping each colour to a different label. We fitted directional unit vectors to pieces of lines using linear regression on groups of pixels from the same line. Here we chose the direction of the arrow so that arrows at the limbus pointed inward, and so that subsequent arrows on the same line were in the same direction. We applied an error correction step so that velocity vectors of adjacent ring points were required to point in the same direction, flipping individual directors if necessary. Finally, we averaged orientations over the different arrows within distance of a point to define the velocity unit vector on ring points. These inferred velocity vectors are shown as red arrows in Figure 9.

Since the cornea corresponds to a spherical cap with a 70° angle, we needed to close the gaps left by the radial cuts. To achieve this, we applied a virtual moulding technique to the image points starting from the inside and moving down ring by ring (Figure 9e): To move the tissue on either side of the cuts together, we introduced springs between all neighbouring pairs of points, and also between points between the ring and the previous (interior) one if they were within distance dr. We defined an energy function for our network of springs,

Using a Monte Carlo simulation with effective temperature kT = 0.01, we then sequentially equilibrated each ring for 1000 Monte Carlo sweeps, starting from the centre. The final angles were then equally spaced, closing the cuts. The final arrangements of points on the now curved cornea are shown in Figure 9f. We then inferred the flow direction in the region between stripe edges: Since neighbouring flow fields (or ‘spins’) are expected to be similar to each other, we introduce an aligning energy term akin to modelling magnets, i.e. the XY model, , where the sum is over neighbouring points only (black lines in Figure 9g). In spherical coordinates, (ê_r, ê_θ, ê_ϕ) as shown in Figure 2c, we could define the orientation of the flow direction through the deviation angle α from inward flow, i.e. in the local tangent plane to the sphere. One could then show using parallel transport (to first order), that the dot product corresponds to . Here α_j was shifted by the difference in longitude ϕ_j − ϕ_i between the two points, which is significant for direct neighbours close to the pole. We then again used a Monte Carlo approach with simulated annealing where we gradually reduced the effective temperature in steps from kT = 1 to kT = 0.001 to minimise the energy of our spin configuration while keeping the red arrows where we already knew the flow direction was fixed. As the directions are parameterised by declination α_i, all arrows remain tangent to the sphere, ensuring that our flow field is on the corneal surface, not into it or out from it. The resulting estimated flow field is shown in Figure 9h, and in main Figure 2a.

Defect finding algorithm and spherical coordinates

In order to identify the locations of topological defects in the flow field, we first constructed a Delaunay triangulation. For the experiment (Figure 2d), we first projected all points into the xy plane, while for the simulation, we derived the Delaunay triangulation from the convex hull of the particle positions on the sphere. The unit-length vectors along the direction of velocity can be written using the angle θ with an axis in the local tangent plane, i.e.. We located topological defects by summing over the angle differences around the loop of the projection plane if one traverses each triangle counterclockwise. The topological charge associated with each triangle is , (Kamien, 2002; Henkes et al., 2018). This is an approximation of the Burger’s loop integral ·dI, and allowed us to locate defects as the centres of loops with C ≠ 0. Generally, C = 0 (e.g. in loop A) for all triangles except the one at the centre of the spiral, where C = 1, corresponds to the central +1 topological defect (loop B). To obtain the spiral angle profile α(θ), we first located the topological defect(s) in the field. As shown in Figure 2a (experiment) and Figures 6c and 6d, we reliably found a single +1 topological defect near the centre for mature corneas. We then shifted our spherical coordinate system so that the defect was at the pole, and defined the corneal spiral angle (θ) through averaging over bins of latitude θ.

Time-lapse imaging

We dissected corneas, halved them, sterilised them in 1:10 penicillin/streptomycin solution for 10 minutes then washed them in PBS before attaching the halved sections, stromal side down, to uncoated 60 mm tissue culture dishes and flooding them with corneal culture medium following (Hazlett et al., 1996). Explant cultures were maintained at 37°C, 5% CO₂ until imaging. Cultures were imaged in HEPES-buffered media using phase-contrast optics on an inverted Leica (DM IRB) chamber microscope for 24 - 48 hours, taking images at 10-minute intervals for PIV analysis.

Division and centrosome angles

We used ImageJ to measure mitotic angles of division (phospho-histone-h3, DAPI) and centrosome position (pericentrin, DAPI) in relation to the centre of the eye (Figure 3).

Computational simulation details

Individual corneal epithelial TA cells obey the following equations of motion,

The first equation relates the time derivative of the position, r_i, of the i^th cell to the active self-propulsion velocity v₀ along a migration direction described by the unit-length vector . F_ij is the elastic force with the neighbouring cells labelled by indices j, and forces are summed over all neighbours within the interaction range. The constant ζ is the friction coefficient between cells and the substrate (i.e. the stroma). The second equation describes the rotational dynamics of the vector . In addition to the term ξ_i that accounts for random changes in the direction of , i.e. the rotational diffusion, there is a term that describes alignment with directions of neighbouring particles , and a term that models alignment of the particle’s direction with the direction of its velocity . The rotational diffusion is assumed to be the simple white noise with zero mean and variance 2D_r. Finally, the projection operator P_T (N, a) = a−(N⋅a)N enforces that the motion remains in the plane orthogonal to the local normal , the projection operator ⋅ a projects the rotation axis of to be along the local normal vector. In isolation, on the sphere, this model has a transition to a flocking ring state at sufficiently large J. A full analysis of this model is given in (Sknepnek and Henkes, 2015).

For our simulations to match corneal ‘explants’ and tissue culture ‘plastic’ colonies, we used the same model on a flat plane, so that , simply the out-of-plane direction. If we define the polar angle ϕ_i through and the velocity angle through v_i = v_i (cos θ_ie_x + sin θ_ie_y), where e_x and e are the unit-length vectors, respectively, along the x and y axes of the simulation box, the model eqns. (7) and (8) reduce to a simpler form

We used a simple, short-range interaction potential that includes both repulsive and attractive parts, together with regulated cell division and removal. This model was introduced and characterised in (Matoz-Fernandez et al., 2017). Briefly, the resulting force between particles i and j with radii σ_i and σ_j and distance vector r_ij = r_j − r_i was given by

where σ_ij = σ_i +σ_j and the dimensionless number ϵ sets the attraction depth and width. Due to the soft potential, this model is unstable for ϵ ≥ 0.2, and we chose ϵ = 0.15 throughout, as well as 0.3 polydispersity of the cells. Particles (cells) were modelled to divide and disappear from the sheet according to a constant squeeze-out rate a but a division rate d = d₀(1 − z/z_max) that depends on the local density through z, the number of neighbours of the cell (defined as those in range for F_ij ≠ 0). As shown in (Matoz-Fernandez et al., 2017), there was a ‘self-melting’ dense liquid steady state in this model below a threshold a/d₀ ≈ 0.1, i.e. where extrusion is sufficiently rare for cells to divide to reach confluence at z ≈ 6 and so in turn inhibit division. The parameters we used for the simulations are summarised in Table 1. All simulations reported here were carried out in our simulation package (SAMoS, 2024) using custom configuration and input files. All analysis of simulation data, including coarse-grained fields, topological defects and angular averaging were carried out using our analysis package (SAMoSA, 2024).

Obtaining coarse-grained fields from simulation

Figure 10 shows the emergence of the hydrodynamic velocity field v^H upon time averaging.

The XYZ model revisited: details

Consider a small patch of area S of corneal tissue containing N cells at time t. One can write the change of number of cells between times t and t + Δt as

where the terms are given by the number of cells that flow into the patch (N_in), flow out of the patch (N_out), are born (N_d) and are squeezed out (N_a).

This is nothing but a form of a continuity equation which expresses that all cells need to be accounted for, and it can be written in a differential form. If ρ(r) is the local number density of cells, then N(t) = ∫_S ρ(r, t)dr. Other terms can be rewritten as functions of ρ: The total number of births during time span Δt is given by N_d = Δt ∫_S d(ρ)ρ(r, t)dr, where d(ρ) is the density-dependent division rate. Similarly, N_a = Δt ∫_S a(ρ)ρ(r, t)dr, where a(ρ) is the density-dependent extrusion rate. The flux of cells through the border of the patch can be expressed through a number flux current J = ρv, i.e. the number of particles with velocity v passing through a unit line in unit time, with units of inverse time and length,

where the integral is over the boundary of S and n is an outward pointing unit-length normal. Putting all of the pieces together, one obtains

Using the divergence theorem on the first term on the right-hand side, it can be rewritten as , followed by the limit Δt → 0. This leads to the continuity equation

Cornea spiral flux equations

We now explore the implications of the continuity equation for an existing spiral flux on the corneal surface. Consider a cornea that is on the surface of a sphere with radius R, with an opening angle at the limbus of θ_L = 70°. Using spherical coordinates as in Figure 2c, with θ as the angle from the cornea centre and ϕ the azimuthal angle, we assume that we have an established symmetric spiral flux ending at the cornea centre, making an angle α(θ) with the radially inward direction. We assume that the density is symmetric, i.e. that we have ρ(θ), and that the average cell velocity in the local polar coordinate system is given by . We also recall that in spherical coordinates, the divergence of vector is given by

Here, we have U = ρv, and we are always on the corneal surface with r = R_c and we eventually write

In steady-state, i.e. when the system becomes time-independent so that , eqn. (13) reduces to ∇ ⋅ (ρv) − A(ρ)ρ = 0, where the effective net loss rate A(ρ) = − [d(ρ) − a(ρ)]. We then find a time-independent, ordinary differential equation in the variable θ linking velocity, density and spiral angle:

We can make one more simplifying assumption, namely that ρ is approximately constant in the system. This is reasonable, first observationally as shown in Figure 11, and second, because the relaxation time of the potential is much much faster than the typical migration time. We will however retain two effects on the net loss rate A: First, A(ρ) is very sensitive to density, so we can not fully neglect its effect. Second, in our system TA cells have a typical lifetime, suppressing their death rate near the emergence from the limbus. We will therefore retain an angle dependence A = A(θ). Then the final continuity (XYZ) equation is given by

Figure 11.

Figures 6 and 7 show the numerically observed spiral angle α(θ), velocity magnitude v(θ) and the effective death rate A(θ). It becomes immediately clear that none of these functions is simple, and that there is a strong dependence on the radius of the simulated cornea, despite R only appearing in the scaling of the net death rate A in eqn. (15). A major contribution to the complexity is the form of A(θ): As can be seen, it is negative near the limbus where cells are recent descendants of stem cells, and then positive and approximately constant towards the inner regions of the cornea. Nevertheless, eqn. (15) is satisfied: the centrally inward (along ê^θ) part of the velocity is given by F (θ) = v(θ) cos(α(θ)), representing the inward flux of material. Undoing a number of the spherical coordinate derivatives, we can rewrite equation 15 as a function of the inward flux:

We supplement this equation with the following boundary conditions. First of all, the flux needs to vanish at the central defect, i.e. at θ = 0, F (0) = 0. This means that either v(0) = 0, i.e. the cells completely stop moving anywhere on average at the centre, or cos(α(0)) = 0, which means that the velocity is at an angle α(0) = ±90°, or completely orthogonal to the inward direction. We can then in fact find an integral solution for F (θ): First make the change of variables u(θ) = F (θ) sin θ, which when replaced in eqn. (16) leads us to

using the boundary condition u(0) = F (0) sin(0) = 0. Then we arrive at our final prediction for the flux on the cornea for a given death profile A(θ):

Acknowledgements

We would like to thank Luke Coburn for creating the initial version of the MATLAB stripe fitting algorithm. We thank Luiza Konkina for assistance with tracing corneal striping patterns. RS and SH would like to thank the Isaac Newton Institute for Mathematical Sciences, Cambridge, for support and hospitality during the programme New statistical physics in living matter, where part of the work on this paper was undertaken.

Additional information

Author contributions

Kaja Kostanjevec

Contributions: Experimental work. Experimental data analysis. Performed numerical simulations. Wrote and edited the manuscript.

Rastko Sknepnek

Contributions: Model and code development. Edited and commented on the manuscript.

Jon Martin Collinson

Contributions: Conceptualisation. Experiment setup and supervision. Experimental data analysis supervision. Wrote the manuscript.

Silke Henkes

Contributions: Conceptualisation. Model development. Performed numerical simulations. Experimental and simulation data analysis. Wrote the manuscript.

Funding

Martin Collinson - UK BBSRC (grant number BB/J015237/1)

Silke Henkes - UK BBSRC (grant number BB/N009150/1-2), and the support of the NWO sector plan science and Leiden University.

Kaja Kostanjevec - EASTBIO PhD studentship

Rastko Sknepnek - UK EPSRC (Award EP/W023946/1)

Additional files

Description of supplementary videos

Supplementary video 1

Supplementary video 2

Supplementary video 3

Supplementary video 4

Supplementary video 5

Supplementary video 6

Supplementary video 7

Supplementary video 8

Supplementary video 9

Supplementary video 10

Supplementary video 11