A soft active matter models explains spiral epithelial cell migration on in-vivo corneas

Reviewer #1 (Public review):

Summary:

The manuscript by Kostanjevec et al. investigates the mechanism behind spiral pattern formation in the cornea. The authors demonstrate that the spiral motion pattern on the mammalian corneal surface emerges from the interaction between the limbus position, cell division, extrusion, and collective cell migration. Using LacZ mosaic murine corneas, they reveal a tightening spiral flow pattern and show that their cell-based, in silico model accurately reproduces these patterns without global guidance cues. Additionally, they present a continuum model that extends the XYZ hypothesis to describe cell flux on the cornea, offering a quantitative explanation for tissue-scale processes on curved surfaces.

Strengths:

The manuscript is well-written, with a systematic approach that clearly explains experimental setups, model construction, assumptions, parameter selection, and predictions. The discussion also provides insightful perspectives on the broader implications of the results for both physics and biology.

Weaknesses:

The authors emphasize polar alignment as a key feature of the spiral pattern based on simulation results. However, they do not provide experimental evidence for this polar alignment.

Reviewer #2 (Public review):

In K. Kostanjevec et al., the authors study a possible mechanism for the formation of spiral patterns in the cornea. First the authors analyze an inferred velocity field, which is deduced from images of fixed corneas, and then determine the position-dependent spiral angle of this velocity fields. Next, the authors analysed two possible markers of cell polarity: the direction of the centrosome-nuclei and the axis of mitosis. Then the authors introduce a stochastic agent-based model of self-propelled particles with over-damped dynamics and with aligning interactions to the orientation of the nearest neighbors and to the particle's velocity. The authors claim to be able to reproduce the equal-time autocorrelation function and the velocity Fourier spectrum. Then the authors introduce the geometry of the cornea by constraining the dynamics on a spherical cap and show that their model can reproduce a typical trajectory in experiments. Finally, the authors produce a phase diagram of the states at a fixed time point as a function of the spherical cap radius and the strength of the coupling aligning constant. Finally, the authors propose an interpretation of the cell fluxes based on the equation of mass conservation.

Author response:

The following is the authors’ response to the original reviews.

Public Reviews:

Reviewer #1 (Public review):

Summary

The manuscript by Kostanjevec et al. investigates the mechanism behind spiral pattern formation in the cornea. The authors demonstrate that the spiral motion pattern on the mammalian corneal surface emerges from the interaction between the limbus position, cell division, extrusion, and collective cell migration. Using LacZ mosaic murine corneas, they reveal a tightening spiral flow pattern and show that their cell-based, in silico model accurately reproduces these patterns without global guidance cues. Additionally, they present a continuum model that extends the XYZ hypothesis to describe cell flux on the cornea, offering a quantitative explanation for tissue-scale processes on curved surfaces.

Strengths

The manuscript is well-written, with a systematic approach that clearly explains experimental setups, model construction, assumptions, parameter selection, and predictions. The discussion also provides insightful perspectives on the broader implications of the results for both physics and biology.

We thank the reviewer for their positive assessment of the manuscript. We are pleased that the reviewer found the work well written and systematic, and that the experimental design, model construction, assumptions, parameter selection, predictions and broader discussion were clearly presented. We have aimed to preserve these strengths in the revised manuscript while substantially expanding the discussion and analysis in response to the reviewer’s concerns.

Weaknesses

The central premise of the manuscript, that the spiral patterning of epithelial corneal cells occurs without guidance cues, is not fully supported. The authors overlook the potential role of axons in guiding epithelial cells, despite clear evidence of spiral axon patterns in their own Fig. 1b. Previous literature indicates that axon patterning precedes epithelial cell patterning, suggesting that epithelial migration might be influenced by pre-existing neural structures (e.g., Leiper et al. 2002, IOVS 2013). The authors need to address this point, possibly by exploring whether axonal patterns serve as a template for epithelial cell migration, or by providing experimental evidence to rule out axon-based guidance.

The reviewer raises an important point that we now address in the revised manuscript. We respectfully disagree with the assertion that the premise of our work is faulty. At no point do we claim that global guidance cues are absent or ignore the fact that the corneal nerves swirl; rather, our results show that such global or contact-mediated cues are not required to explain the observed swirling patterns of epithelial cell radial migration in the adult cornea.

Although our work shows that a swirling prepattern of nerves is not required to obtain radial patterns of epithelial cell migration, we now address why nerve swirling occurs and how it affects interpretation of our data. Previous literature indicates that nerve swirling is visible from approximately 3 weeks, before epithelial striping patterns become apparent in transgenic LacZ and GFP reporter mosaics at about 5 weeks. However, both experimental observations and our modelling indicate that spiral epithelial cell migration proceeds for many days before reporter stripe patterns become visible. Thus, epithelial migration may already be underway before epithelial striping is detectable. If axons are following epithelial migration, they would therefore be expected to become radially aligned before the reporter epithelial stripes are evident.

We also considered the alternative possibility that epithelial cells follow an axonal prepattern. However, several experimental observations argue against the hypothesis that epithelial migration is primarily guided by axonal projections. In situations of genetic mutation or corneal injury, epithelial cell migration can progress independently of, or ahead of, corneal axon extension. In chimeric Pax6+/− LacZ⁺ ↔ Pax6^+/+ LacZ^- mice, where the normally disrupted radial migration of Pax6+/− epithelial cells is restored, the underlying nerves may continue to exhibit abnormal projection patterns. These findings suggest that axonal organization is at least partly dependent on epithelial behaviour rather than the reverse.

Importantly, we do not exclude the possibility that additional cues, including axonal contact, neurotrophins or other environmental signals, may modulate or refine epithelial migration. We have therefore expanded the revised manuscript to discuss these possibilities and the relevant literature more fully.

While the model is well-constructed, it currently falls short of its stated goal of elucidating the mechanisms of spiral formation. Key questions remain unanswered:

Is the curvature of the cornea necessary for spiral formation, or would a simpler disk geometry suffice?

What role do boundary conditions play?

How well do the model's predictions quantitatively match experimental data?

The current comparisons in Fig. 4c-f lack quantitative agreement, and this discrepancy should be discussed with possible explanations.

We thank the reviewer for identifying these points, which we have now addressed in the revised manuscript.

First, we have examined the role of geometry more systematically. The spiral pattern also appears in a disk geometry, and the same qualitative migration pattern is observed across a broader range of simulated geometries, including different curvatures, cap angles, a prolate ellipsoid, an oblate ellipsoid and a disk. The precise shape of the spiral depends on geometric features such as curvature and cap-angle opening, but spiral formation is robust across convex cornea-like geometries. These results are now discussed in the new section “Robustness of the spiral migration pattern and requirement for limbal stem cells” and shown in the revised Fig. 10.

Second, we have clarified the role of boundary conditions, particularly the role of limbal epithelial stem cell proliferation. Without limbal stem cells, and with all other parameters unchanged, the cornea fails to produce the radial striping pattern. At the alignment strength where robust spiral formation normally appears, the simulated tissue flow is disordered and resembles our in vitro calibration simulations. At higher alignment strengths, spiral formation is still not recovered; instead, defects become anchored to the boundary. These findings show that ordered influx from the limbus is important for promoting the spiral state. They are now discussed in the same new section and shown in revised Fig. 9.

Third, we have revised the quantitative comparison between model and experiment. We agree that the original comparison was limited. We have therefore reanalysed both experimental and simulation data, focusing on the time-averaged hydrodynamic velocity field rather than short-range fluctuations amplified by divisions in the numerical model. We also replaced the Fourier-space velocity correlation functions with spatial velocity correlation functions, which are more directly interpretable. The revised analysis shows that the experiments have mesoscale spatial and temporal correlations, of the order of 5-6 cell sizes in space and about one hour in time, and that these are well captured by the simulations for both plastic and explant substrates. We have also added representative experimental and simulation snapshots in revised Fig. 4g-j.

For the full cornea, we acknowledge that direct quantitative comparison remains limited by the available experimental data. We can compare with inferred migration direction fields and resurfacing timescales, but we do not yet have direct live measurements of the full corneal velocity field.

The authors emphasize polar alignment as a key feature of the spiral pattern based on simulation results. However, they do not provide experimental evidence for this polar alignment. The manuscript includes discussions of polar and nematic symmetries that, without supporting data, feel somewhat distracting. If direct experimental evidence for polar alignment is not available, the authors could instead quantify nematic alignment as the spiral forms. This would also allow them to explore potential crosstalk between nematic cell orientation and the polar alignment of self-propulsion, especially considering recent studies showing alternative mechanisms for vortex formation in similar systems.

We thank the reviewer for pointing out that the discussion of polar and nematic alignment was confusing. We have substantially revised this part of the manuscript.

We agree that we do not have direct experimental evidence for polar alignment. However, several observations support the interpretation that the system is dominated by substrate-based polar motility with weak polar alignment. We have now quantified nematic alignment of cell orientations and find no evidence of significant elongation or local nematic order in corneal epithelial cells, as shown in the new Fig. 3d.

Our computational model begins from uncorrelated, substrate-based polar active cell migration. We tested whether polar alignment is necessary and found that the in vitro data are inconsistent with the complete absence of alignment: the flocking order parameter is too high, and the spatial and temporal correlations are larger than expected from persistent driving alone.

We also discuss that cell-cell stress patterns in the epithelium may remain nematic, but any such effect must be sufficiently weak not to dominate the observed axon motion. We have further revised the discussion to include recent related work showing spiral formation in substrate-based cell migration models with polar dynamics, as well as recent work indicating that nematic-like phenomenology can arise from minimal ingredients such as uncorrelated polar activity and cell deformability. This revision is intended to make the interpretation clearer and to avoid overemphasising unsupported claims.

Reviewer #2 (Public review):

In K. Kostanjevec et al, the authors study a possible mechanism for the formation of spiral patterns in the cornea. First the authors analyze an inferred velocity field, which is deduced from images of fixed corneas, and then determine the position-dependent spiral angle of this velocity fields. Next, the authors analysed two possible markers of cell polarity: the direction of the centrosome-nuclei and the axis of mitosis. Then the authors introduce a stochastic agent-based model of self-propelled particles with over-damped dynamics and with aligning interactions to the orientation of the nearest neighbors and to the particle's velocity. The authors claim to be able to reproduce the equal-time autocorrelation function and the velocity Fourier spectrum. Then the authors introduce the geometry of the cornea by constraining the dynamics on a spherical cap and show that their model can reproduce a typical trajectory in experiments. Finally, the authors produce a phase diagram of the states at a fixed time point as a function of the spherical cap radius and the strength of the coupling aligning constant. Finally, the authors propose an interpretation of the cell fluxes based on the equation of mass conservation.

We thank the reviewer for their careful assessment of the manuscript and for recognising the work as a solid theoretical study. We have revised the manuscript substantially in response to the reviewer’s major concerns, particularly regarding the terminology of topological defects and stagnation points, the comparison with experiments, and the role of corneal geometry.

Regarding the terminology of topological defects, we have clarified the distinction between a velocity-field stagnation point and the topological classification of the velocity direction field. Stagnation points can be assigned a topological index when one considers the direction of the vector field away from the core, while ignoring the magnitude. We agree that the physical origin of interactions in a velocity field differs from that in a director field, and we have revised the text to avoid confusion. The revised manuscript now includes Box 1, which summarises the relevant topological concepts and caveats.

Regarding the comparison with experiments, we have expanded and clarified the validation of the inferred velocity field. The LacZ reporter system was designed for lineage tracing and therefore reports coarse-grained cell motion and growth patterns. Direct live imaging of the full cornea remains technically difficult because of the macroscopic size, curved surface and long resurfacing time. However, live fluorescent reporter systems are consistent with our in vivo model and support the interpretation that inferred velocity fields recapitulate epithelial migration in vivo. We have also clarified the interpolation procedure using the XY model: approximately 30% of the corneal surface is covered by directly estimated velocity directions from stripe edges, rising to more than 50% near the central spiral. These measured directions provide sufficient boundary conditions for the annealing procedure to converge to a slowly varying field consistent with the observed stripe geometry.

We have also revised the comparison between simulation and experiment. For in vitro data, we now compare spatial velocity correlation functions of the time-averaged hydrodynamic velocity field, rather than relying on Fourier-space correlations. The revised comparison shows good agreement between experiments and simulations for both plastic and explant substrates. For the full cornea, we acknowledge that the available quantitative data are limited to inferred migration direction fields and resurfacing timescales.

Regarding the role of geometry, we have now simulated a wider range of substrate geometries, including spherical caps with different curvatures and cap angles, prolate and oblate ellipsoids, and a flat disk. The spiral migration pattern is robust across these convex geometries, although the detailed spiral shape depends on geometric features. We have also clarified that the boundary condition of inward limbal influx is crucial: without limbal stem cells, the radial striping pattern does not form, and defects may instead anchor to the boundary. These results are now shown in revised Figs. 9 and 10 and discussed in the new section “Robustness of the spiral migration pattern and requirement for limbal stem cells.”

Overall, these additions clarify that the proposed mechanism relies on the combination of polar motility, weak alignment, limbal influx, cell division and extrusion, and cornea-like confinement, while also acknowledging the current experimental limitations.

Recommendations for the authors:

Reviewing Editor:

The authors could strongly improve the manuscript by following the recommendations given below.

Reviewer #1 (Recommendations for the authors):

There are, however, substantial shortcomings that authors need to address to make their claims supported by enough evidence:

Major concerns:

(1) Neglect of Potential Axon Guidance

As it stands the premise of the paper is unfortunately faulty. The authors overlook the potential role of axons in guiding epithelial cells, despite clear evidence of spiral axon patterns in their own Fig. 1b. Previous literature indicates that axon patterning precedes epithelial cell patterning, suggesting that epithelial migration might be influenced by preexisting neural structures (e.g., Leiper et al. 2002, IOVS 2013). The authors need to address this point, possibly by exploring whether axonal patterns serve as a template for epithelial cell migration, or by providing experimental evidence to rule out axon-based guidance.

See for example:

from: https://iovs.arvojournals.org/article.aspx?articleid=2126612:

It is therefore not clear why the authors clearly show the axon vortex in Figure 1, then talk about prepatterning - and never mention the axon vortex again.

The reviewer raises an important point that we now address in the revised manuscript. We respectfully disagree with the assertion that the premise of our work is faulty. At no point do we claim that global guidance cues are absent or ignore the fact that the corneal nerves swirl; however, our results show that such global or contact-mediated cues are not required to explain the observed swirling patterns of epithelial cell radial migration in the adult cornea. Although our work shows that a swirling ‘prepattern’ of nerves is not required to obtain radial patterns of epithelial cell migration, we address below why it occurs and how it affects our data.

An intuitive interpretation of corneal anatomy would be that sensory axons follow the path of least resistance between migrating epithelial cells. We believe this is the most likely explanation in the normal wild-type cornea; however, the reviewer correctly notes that previous literature indicates that axon patterning precedes epithelial cell patterning. Nerve swirling is observed from approximately 3 weeks, earlier than the epithelial striping patterns that become apparent in transgenic LacZ and GFP reporter mosaics at about 5 weeks (e.g. Collinson et al., 2002 [PMID: 12203735]; Iannaccone et al., 2012 [PMID: 22347498]; McKenna and Lwigale 2011 [PMID: 20811061]). We now address this in the revised manuscript and below.

The early appearance of radial axonal projections can be readily explained even if axons are following the epithelial cells. Experimental observations (Collinson et al., 2002 [PMID: 12203735]) and our modelling (Fig. 6a,b) both indicate that spiral epithelial cell migration proceeds for many days before stripe patterns become visible in reporter mosaics. Thus, epithelial migration is already underway before the striping pattern becomes detectable. If axons are following the epithelial migration, they would be expected to become radially aligned before the reporter epithelial stripes were evident. The apparent precedence of radial axonal projections over epithelial striping is therefore fully consistent with the biological scenario in which axons follow the migrating epithelial cells. Corneal epithelial basal cells have been shown to wrap around individual and grouped subbasal axons, acting as surrogate glia (Stepp et al., 2016, Investigative Ophthalmology & Visual Science 57, 1292), which represents a plausible mechanism to allow migrating epithelial cells to shepherd axons in their direction of movement.

We considered that epithelial cells may follow an axonal prepattern, but several experimental observations argue against the hypothesis that epithelial migration is guided by axonal projections. In situations of genetic mutation or corneal injury, epithelial cell migration can progress independently of, or ahead of, the extension of corneal axons (Leiper et al., 2009 [PMID: 19029029]; Song et al., 2004 [PMID: 14744881]). Furthermore, in chimeric Pax6^+/− LacZ⁺ ↔ Pax6^+/+ LacZ^- mice, where the normally disrupted radial migration of Pax6^+/− epithelial cells is restored, the underlying nerves may continue to exhibit abnormal projection patterns (Leiper et al., 2009 [PMID: 19029029]). These results suggest that axonal organization is at least partially dependent on epithelial behaviour rather than the reverse.

Importantly, none of this evidence excludes the possibility that additional cues (including axonal contact or other environmental signals) may modulate or refine epithelial migration. For example, Walczysko et al. 2016 [PMID: 27563231] showed that isolated epithelial cells cultured on de-epithelialised and de-nervated corneal stroma still migrate with a small but significant radial bias, indicating that epithelial cells can respond to physical features of their environment. Our model likewise includes a component of alignment with environmental structure. Since the central radial striping in many of our simulations is somewhat less ordered than in vivo, it is plausible that additional biological guidance cues or axonal neurotrophins help refine the pattern.

We now discuss these issues and the relevant experimental evidence more fully in the revised manuscript.

(2) Model Validation and Complexity

While the model is well-constructed, it currently falls short of its stated goal of elucidating the mechanisms of spiral formation. Key questions remain unanswered:

Is the curvature of the cornea necessary for spiral formation, or would a simpler disk geometry suffice?

To address the Reviewers’ concerns about the role of geometry, we first note that the spiral pattern also appears in a disk geometry. In fact, the geometric parameters of the cornea vary across mammals, including humans, yet the spiral pattern is preserved (Dua, et al., 1993 [PMID: 8325424]; Zander and Weddell, 1951 [PMID: 14814019]). Similar variations arise in disease states, e.g. in keratoconus the cornea becomes elongated.

On a spherical cap, and all shapes with the same topology, the boundary winding number fixes the interior index, so ongoing limbal influx maintains a total index of 1. To explore the role of geometry more systematically, we simulated a broader range of geometries, including different curvatures, cap angles, a prolate and an oblate ellipsoid, and finally a disk, and compared the resulting patterns with published data across mammals and with disease states.

For all of these shapes, we find the same qualitative migration pattern – a central spiral, although its precise shape depends on geometric features such as curvature and cap-angle opening. These results are discussed in a new section “Robustness of the spiral migration pattern and requirement for limbal stem cells” and shown in Fig. 10 of the revised manuscript.

What role do boundary conditions play?

In the revised manuscript, we address the role of boundary conditions, in particular the presence of limbal epithelial stem cell proliferation. Without limbal stem cells, and with all other parameters kept unchanged, the cornea fails to make the radial striping pattern. We observe two distinct changes: First, at J=0.1, the amount of alignment at which robust spiral formation appears normally, the simulated tissue flow is instead disordered, in fact very similar to our in vitro calibration simulations (revised Figure 4). This shows that the boundary cue of an ordered influx from the limbus promotes flocking when it otherwise would not (yet) appear. Second, when we increase the alignment to J=0.15 or J=0.2, we still do not observe spiral formation. Instead of the expected central vortex shape, we have anchoring of defects to the boundary, facilitated by the effective compressibility of the tissue because of the density feedback in the division/extrusion rates.

These findings are discussed in the new section “Robustness of the spiral migration pattern and requirement for limbal stem cells” and in Fig. 9 in the revised manuscript.

How well do the model's predictions quantitatively match experimental data?

The current comparisons in Fig. 4c-f lack quantitative agreement, and this discrepancy should be discussed with possible explanations.

Regarding the in vitro cell data and matching simulations: We are aware that we have limited data to work with, and our match is intended as a rough estimate of physical parameters.

For the revision, we have reanalysed both experimental and simulation data carefully. We realised that certain details of our numerical model amplify short-range fluctuations, namely the way cell divisions induce stress dipoles, and the fact that we did not include cell-cell friction forces. This is not merely a guess but emerged from related theoretical work by some of us (Keta and Henkes [PMID: 40556485]; Kammeraat et al, arXiv:2508.01046 (2025)).

We therefore compared the time-averaged velocity field excluding divisions to the experiment, the same quantity that we already introduced as hydrodynamic velocity for the corneal surface. We also carried out further simulations that included cell-cell friction forces for comparison. They led to very similar results, albeit with a transition to flocking at somewhat lower alignment strength J.

Instead of the Fourier-space velocity correlation functions that are hard to interpret, we switched to spatial velocity correlation functions. Our experiments have ‘swirly’ velocity fields with mesoscale spatial and temporal correlations, of the order of 5-6 cell sizes in space and one hour in time. As can be seen in the revised panels 4c,d for the spatial correlations of the hydrodynamic velocity, the match between experiment and simulations is in fact good for both plastic and explant substrates. There are also systematic changes in length and time scales between the two experiments that emerge without fine-tuning from simulation. We furthermore have included snapshots of both experimental in vitro conditions and matching simulations as new panels 4g-j, showing that the mesoscale correlations appear in the hydrodynamic velocity.

Ultimately the parameter values and length and time scales that we infer for our systems (see Table 1) are quantitatively consistent with three other estimates of in vitro epithelia (Henkes et al. 2020 [PMID: 32179745]; Saraswathibhatla et al, Extreme Mechanics Letters 48, 101438 (2021), Kammeraat et al, arXiv:2508.01046 (2025)).

For cell flows on the full cornea, regrettably we lack further quantitative data to compare to beyond the migration direction fields inferred in Figure 2, and the time scales of corneal resurfacing (Figure 6).

(3) Importance of polar alignment.

Polar alignment is put forward as one main feature for the observed patterns based on the simulation results. However, no experimental confirmation for such polar alignment is presented. The authors instead present various scattered discussions about polar versus nematic symmetry, which at times reads rather unnecessary and distracting from their main message.

If they are not providing direct experimental evidence on the polar alignment, at least they could quantify nematic alignment of cell orientation as the spiral forms. One possibility is that nematic alignment of cell orientations has a crosstalk with the polar alignment associated with the self-propulsion of the cells. This is important because, as acknowledged by authors, several recent works in the context of in vitro epithelial under disk confinement have revealed alternative mechanisms for spiral vortex formation.

We thank the Reviewer for pointing out the confusion, and we have therefore completely revised our discussion of alignment. While the evidence remains indirect, the following observations lead us to conclude that our system is dominated by substrate-based polar motility and weak polar alignment:

We have investigated nematic alignment of cell orientations. As can been seen in new Figure 3d, corneal epithelial cells do not show evidence of significant elongation in any direction, i.e. we find no local nematic order.

Our computational model starts from uncorrelated, substrate-based polar active cell migration. We have carefully investigated if polar alignment is a necessary ingredient and found that the in vitro data are inconsistent with the absence of alignment: the flocking order parameter is too high, and the spatial and temporal correlations are larger than those expected from persistent driving only.

Still, cell-cell stress patterns in the epithelium may remain nematic. While one cannot exclude anything, the effect must be sufficiently weak to not affect axon motion. Their naturally long, thin shapes would strongly react to nematic stresses, but we do not see ±1/2 defects in their growth patterns, only polar ±1 defects.

We were recently made aware the work of Lång et al. (2024) [PMID: 38630812], in which the authors also observe spiral formation in cell migration on a substrate, with +1 topological defects. They explain their observations using a polar model, and our results are consistent with their observations and model.

In the active matter community, the ‘active nematic cell sheet’ paradigm is currently undergoing a revision. Notably, nematic-like phenomenology, in particular ±1/2 defects, can also arise from the minimal ingredients of uncorrelated polar activity and cell deformability (Chiang et al. 2024 [PMID: 39302997]).

Minor comments:

- It would help the reader to see some representative images of the cells, explants, and the simulation (maybe with vectors overlaid), as it could be hard to imagine what exactly is happening and the other relevant properties to compare.

Please see new panels Fig. 4g-j, and new supplementary videos S3 and S4.

- Add a colorbar that represents direction to Fig. 1a.

We assume the Reviewer meant Fig. 2a. We have added a circular orientation colour chart to the image.

- When do additional +1,-1 defects appear? The authors just say it is unlikely, but it is not clear when they are observed. Disease state?

The reviewer is correct about pathology. We now cite (in conclusion) Collinson et al. 2004 [PMID: 15037575], which shows disruption and discontinuity in Pax6-mutant corneas with chronic corneal degeneration. We also have a publication in prep that shows extra +1 and -1 discontinuities occur during wound healing. We don’t want to include these data in this manuscript, but cite Sagga et al. (in prep).

Reviewer #2 (Recommendations for the authors):

Overall, the manuscript presents a solid theoretical work. However, I have a few major concerns on this manuscript. (1) The authors use the concept of topological defect to refer to stagnation points of the velocity field (2) The comparison to experiments remains qualitative and it is unclear whether the proposed mechanism is at work, (3) The role of geometry remains unclear.

Major points:

(1) On page 3 the authors claim that topological defects exist in velocity fields, however this statement is incorrect and can confuse readers. Unlike the director field of an ordered phase, their velocity field is a vector in R^2 with a norm that is not fixed, and therefore the velocity field has no topological defects. In fluid dynamics, these special points in the velocity field are called stagnation points, fluid sinks or sources. This distinction is important because the physical nature of the interaction forces between two "topological defects" in a velocity field is fundamentally different to the interaction forces between two topological defects in a director field. For these reasons, it can confuse readers to mix the two concepts (stagnation points vs topological defects). Note that if the authors address this concern, many parts of the main manuscript should be rewritten.

Stagnation points are topological defects since the topological classification ignores the magnitude and looks only at the direction of the vector field away from the core. There are some caveats related to the role of boundary conditions, since in the case of a fluid, if boundary conditions are not fixed, one can eliminate the defect by setting the flow field to zero everywhere. Mathematically, a pair of point sources or sinks in an incompressible potential flow interacts through the same Green’s function that gives the elastic interaction between two-point disclinations in a 2D nematic director field, so their long-range pair potentials are formally identical (~ln(r) in 2D). The physical mechanism behind these forces is, as correctly pointed out by the reviewer, quite different. In addition, our coarse-grained velocity field is effectively compressible, which has consequences for the hydrodynamic equations we (can) write, see below.

In the revised version, we clarify these points by adding Box 1, which summarizes the idea of topological defects.

(2.1) How did the authors check that the velocity field that is inferred from the stripe edges matched the coarse-grained cell velocity field in a life sample? How did the authors validated the interpolation of the inferred velocity field using an XY model? What is the scale of the inferred velocity field? Can the authors clarify also this point?

The murine cornea LacZ reporter system was specifically designed to allow for lineage tracing, i.e. following coarse-grained cell motion and growth patterns. The combination of macroscopic size (3.6 mm diameter), two-week resurfacing time and the curved surface however stymied our early attempts to directly measure the cell velocity field on the cornea using confocal time-lapse microscopy. However live sample fluorescent reporter systems are fully consistent with our in vivo model and shown conclusively that inferred velocity fields are recapitulated in vivo (Park et al., 2019 [PMID: 31843909]).

For the XY model inference: We first note that approximately 30% of the corneal surface is covered by directly estimated velocity directions from the stripe edges (red arrows, Fig. 11f). This rises to more than 50% near the central spiral (red arrows, Fig. 11h) due to the way the stripes narrow near the centre due to cell extrusion. Therefore, the inferred areas are only slightly more than half of the cornea, and in regions where we expect the flow field to be largely uniform with no defects. The red arrows provide sufficient boundary conditions that a simple annealing simulation (or equivalently an energy minimisation) of the XY model rapidly converges to a slowly varying field consistent with those boundary conditions.

Furthermore, in simulation we observe that stripe edges and the macroscopic velocity field correlate strongly with each other once the spiral has fully formed (Fig. 6b-c).

The scale of the inferred velocity field is the distance between arrows in our digital version of the cornea, approximately 40 concentric rings over a 70° cone angle for a R = 1800 μm micron cornea, resulting in a spacing of 55 μm between velocity arrows. This is about at the scale of the in vitro velocity correlations. We are not able to obtain velocity magnitudes using this procedure.

Furthermore, the spiral angle profile reported in Fig. 7b-c appear to be different to that found in experiments Fig. 2b. Can the authors clarify if their theoretical framework reproduces the spiral angle profiles?

First, we note that empirically, simulations with the largest two radii (R = 1000 μm, R = 1500 μm), approaching the full experimental size, match the observed angle profile best. They both consist of a radially inward profile α(0) = 0° at the edges, only increasing, corresponding to a tightening spiral, below about θ = 20°. Note that very near the corneal centre, few cells contribute to data, and additionally the central defect position fluctuates somewhat. Therefore, the angle profile not reaching α(0) = 90° in the simulations is due to fluctuations and lack of statistics. As these simulations were run with our best fit experimentally matched parameters, this convergence is meaningful and cannot be scaled out.

Second, our partial model (eq. 4, reproduced below) links the spiral profile angle α(θ) with the macroscopic velocity magnitude v(θ) and the net cell loss rate A(θ), and it depends explicitly on the corneal radius R.

That is one equation for three radial fields. If we make the reasonable assumption that simulations at fixed 𝐽 that differ only in R have the same constitutive law A(ρ), and would follow the same continuum velocity equation that ultimately sets v(ρ), the radius R still appears explicitly in the equation. This is consistent with the different radial profiles for different R that we observe in Fig. 7c. Furthermore, in Fig. 7b we observe that above the flocking threshold, different J lead to very similar profiles. This indicates that the system enters a fully polar phase.

The same is true in experiment: We expect that cell mechanics and planar cell polarisation coordination are local effects that will set J, A(ρ), and v(ρ). Thus, we do predict that the radius R will explicitly affect the spiral profile consistent with the equation above, but we would need more direct measurements of J, A(ρ), and v(ρ) to go any further.

Properly answering this question would first require simulations of a non-dimensionalised model with different radii, boundary influx, alignment strengths and constitutive laws for the division / extrusion dynamics. Then one would want to construct a matching equation for the polarisation and / or velocity field, going beyond the Malthusian flock approximations of constant magnitude v(θ) and net cell loss rate A(θ) = 0. This is well beyond the scope of this publication, and there is certainly no experimental data to compare to.

(3.1) It appears that the formation of the spiral velocity pattern results from a combination of a radial flow of cells due to cell division at the outer boundary and apoptosis at the geometrical center and azimuthal flow of cells due to flocking in confined geometries. Is this correct? Can the authors explain how does the 3d geometry of the spherical cap modify each of these flow fields?

The Reviewer is broadly correct; however, the true picture is more subtle. Almost all of the divisions and extrusions occur in TA cells, neither at the limbus or at the corneal centre (please see the A(θ) profiles in Figure 8). Flow is then a spiral flock that is partially radial and azimuthal, and with a variable velocity magnitude. It ultimately all has to follow the flux equation 4, in steady state.

For the modification of the flow field due to 3d geometry: Broadly, they simply change the amount of corneal surface available at different angles from the limbus when switching to isomorphic surfaces like, e.g. the disk. Thus, we still observe spirals, but of somewhat modified shapes. Please see the reply to Reviewer 1 above, and new Figure 10 for simulations of different geometries.

A full theory of radially symmetric corneal shapes would again need additional velocity and constitutive equations and is beyond the scope of this publication.

(3.2) In their model, it appears that the geometry is introduced by constraining the dynamics of agents, and it has not direct influence on the alignment of agents. Can the authors explain how the geometry of the spherical cap influences the emergent spiral states? Can the authors show whether their results are robust to changes in the substrate geometry? For example, by changing the substrate geometry from a spherical cap to another convex shape. Can the authors identify differences between the spiral patterns on a spherical cap vs that on a flat disk?

Please see the response to Reviewer 1, above, and new Figure 10. Briefly, the results are robust to changes in the corneal geometry as long as they are other convex shapes. There are differences in details of the spiral shape.

Minor comments:

- On page 3, the authors claim that the Euler characteristic of a spherical cap or a disk is 1. Unfortunately, this statement is incorrect. The Euler characteristic of a spherical cap or a disk is determined by the winding of the vector field around the open boundary. In their case the Euler characteristic is 1 because the velocity field is oriented towards the top of the spherical cap, which give a winding of +1. The authors explain this correctly on page 13.

The Euler characteristic is a topological invariant of the surface and does not depend on the tangent vector field. A disc or spherical cap has Euler characteristic 1. What depends on the boundary winding is the index formula for a vector field on a surface with boundary: the winding of the field along the boundary determines the corresponding boundary contribution, and hence the sum of interior indices. Thus, while the reviewer is correct about the role of boundary winding in computing the field’s index, this does not alter the Euler characteristic of the surface itself.

We note that the orientation of the velocity field does indeed matter in our simulations: In new Figure 10, we show that in the absence of the boundary condition of inward flux at the limbus, we do not observe a spiral robustly, and we do see anchoring of defects to the boundary, with winding numbers that are now different.

- On page 5, the authors claim the clockwise and counterclockwise -oriented spirals are equally likely, however no quantification is provided to support this claim. Can the authors clarify this point?

The approximately equal likelihood is as described in Collinson et al. (2002) [PMID: 12203735].

- On page 8 the authors state "Dipolar active force cannot cause a single cell to migrate, and the flow is an emergent collective phenomenon", here I was confused, because a bacteria swimming in a Newtonian fluid can self-propel by exerting a dipolar active force on the fluid. See for instance work by E. Lauga or I. Aronson. Can the authors clarify this point?

The Reviewer is correct about how bacteria can swim using dipolar active forces. However, and unfortunately, that language was straight ported to very different conditions, that of cells migrating on a frictional substrate with no induced flow, with the equation of motion . Unless there is a net force arising from the stress profile, the cell cannot move, and with microscopic models where the active stress is a single value per cell, that statement is always true. Of course, more detailed cell models with active stresses exist, but their motion is still due to the net force arising from them. We have clarified the statement in the revised manuscript.

- On page 16, the expression for the velocity seems to miss a parenthesis.

Fixed. Thanks!

- On page 16, the authors claim that the flux profiles in the simulations and in the XYZ model are in good agreement. However, there are clear differences for theta> 50 deg. Can the authors discuss the possible explanation for these differences? Note that one may expect a between agreement between two theoretical approaches.

These disagreements are due to imperfections in the observed spiral patterns in simulations. They are not perfectly radially symmetric, and the defect is not always in the dead centre of the cornea. Thus, the theoretical predictions are not quite accurate. Numerically, what happens for θ > 50° is that the radial bins sometimes include part of the limbal zone with strong proliferation, and the corneal edge itself. Both are not described by the flux prediction. We decided not to crop out this region and rather explain where it comes from in the text.