Learning developmental mode dynamics from singlecell trajectories
Abstract
Embryogenesis is a multiscale process during which developmental symmetry breaking transitions give rise to complex multicellular organisms. Recent advances in highresolution livecell microscopy provide unprecedented insights into the collective cell dynamics at various stages of embryonic development. This rapid experimental progress poses the theoretical challenge of translating highdimensional imaging data into predictive lowdimensional models that capture the essential ordering principles governing developmental cell migration in complex geometries. Here, we combine mode decomposition ideas that have proved successful in condensed matter physics and turbulence theory with recent advances in sparse dynamical systems inference to realize a computational framework for learning quantitative continuum models from singlecell imaging data. Considering panembryo cell migration during early gastrulation in zebrafish as a widely studied example, we show how cell trajectory data on a curved surface can be coarsegrained and compressed with suitable harmonic basis functions. The resulting lowdimensional representation of the collective cell dynamics enables a compact characterization of developmental symmetry breaking and the direct inference of an interpretable hydrodynamic model, which reveals similarities between panembryo cell migration and active Brownian particle dynamics on curved surfaces. Due to its generic conceptual foundation, we expect that modebased model learning can help advance the quantitative biophysical understanding of a wide range of developmental structure formation processes.
Editor's evaluation
This work proposes a method to obtain a reduced description of the collective dynamics of thousands of cells moving together during zebrafish gastrulation as a few fundamental modes, and to derive effective dynamics for these modes. This wellwritten work enables a simplified picture of the key features of cellular collective motion, that will be useful to physicists and biologists looking for a quantitative understanding of morphogenesis.
https://doi.org/10.7554/eLife.68679.sa0Introduction
Embryogenesis, the development of a multicellular organism from a single fertilized egg cell, requires coordinated collective motions of thousands of cells across a wide range of length and time scales (Gilbert and Barresi, 2016; SolnicaKrezel, 2005). Understanding how a highly reproducible and robust tissue organization arises from the dynamics and interactions of individual cells presents a major interdisciplinary challenge (Collinet and Lecuit, 2021). Recent advances in highresolution live imaging make it possible to track the internal biological states and physical movements of many individual cells on panembryonic scales throughout various stages of development (Stelzer, 2015; Power and Huisken, 2017; Hartmann et al., 2019; Shah et al., 2019). This unprecedented wealth of data poses two intertwined compression problems of equal practical and conceptual importance. The first concerns the efficient reduction of highdimensional tracking data without loss of relevant information; the second relates to inferring predictive lowdimensional models for the developmental dynamics. Mathematical solutions to the first problem are aided by taking into account the geometry and symmetries of the developing embryo, which suggest suitable basis functions for a coarsegrained and sparse mode representation of raw data (Levy, 2006). Efficient algorithmic approaches tackling the second problem appear within reach thanks to recent advances in the direct inference of dynamical systems equations from data (Brunton et al., 2016; Rackauckas et al., 2021). Building on these ideas, we construct and demonstrate here a computational framework that translates developmental singlecell trajectory data on curved surfaces into quantitative models for the dominant hydrodynamic modes.
Widely applied in physics (Kac, 1966; Goldenfeld and Woese, 2011; Kantsler and Goldstein, 2012; Bhaduri et al., 2020), engineering (Soong and Grigoriu, 1993; Heydari et al., 2021), and spectral computing (Driscoll et al., 2014; Burns et al., 2020; Fortunato et al., 2021), mode representations (Schmid, 2010; Tu et al., 2014) provide a powerful tool to decompose and study system dynamics at and across different energetic, spatial and temporal scales. In quantum systems, for example, mode representations in the form of carefully constructed eigenstates are used to characterize essential energetic system properties (Slater and Koster, 1954; Jaynes and Cummings, 1963). Similarly, turbulence theory has seen significant progress by studying the coupling between Fourier modes that represent dynamics at different length scales. This approach enabled a better understanding of energy cascades (Kolmogorov, 1941; Wang et al., 2021) and provided insights into the nature of turbulence in nonliving (Kraichnan and Montgomery, 1980; Pope, 2000) and in living (Dunkel et al., 2013; Bratanov et al., 2015; Ramaswamy and Jülicher, 2016; Alert et al., 2020) systems. Additionally, the multiscale nature of many biological processes make them particularly amenable to a representation in terms of spatial and temporal modes (Marchetti et al., 2013). Despite this fact, however, mode representations are not yet widely used to characterize and compress cell tracking data, or to infer dynamic models from such data.
To demonstrate the practical potential of mode representations for the description of multicellular developmental processes, we develop here a computational framework that takes cell tracking data as inputs, translates these data into a sparse mode representation by exploiting symmetries of the biological system, and utilizes recently developed ODE inference techniques (Rackauckas et al., 2021) to infer a predictive dynamical model. The model will be specified in terms of a learned Green’s function that propagates initial cell density and flux data forward in time. To validate the approach, we demonstrate that it correctly recovers the hydrodynamic equations for active Brownian particle (ABP) dynamics on curved surfaces. Subsequently, as a first example application to experimental singlecell tracking data, we consider the panembryonic cell migration during early gastrulation in zebrafish (Shah et al., 2019), an important vertebrate model system for studying various morphogenetic events (SolnicaKrezel, 2005; Krieg et al., 2008; Morita et al., 2017). During gastrulation, complex migratory cell movements organize several thousand initially undifferentiated cells into different germlayers that lay out the primary body plan (Rohde and Heisenberg, 2007). The underlying highdimensional singlecell data make this process a prototypical test problem for illustrating how spatiotemporal information can be efficiently compressed to analyze and model biological structure formation.
Results
Broadly, our goal is to translate experimentally measured singlecell trajectories on a curved surface into a quantitative model of collective cell migration dynamics. As a specific example, we consider recently published lightsheet microscopy data that captures the individual movements of thousands of cells during early zebrafish development from epiboly onset at 4 hours postfertilization (hpf) to about 18 hpf (Shah et al., 2019). This developmental period is characterized by a collective symmetry breaking event during which cells collectively migrate over the yolk cell surface (Rohde and Heisenberg, 2007). Namely, they rearrange from an initial localization around the animal pole (AP) (Figure 1A, left) into a more elongated configuration that already indicates the basic geometry of the fully developed zebrafish larva (Figure 1A, right). Working with a twodimensional (2D) sphere projection of the experimental data, we first describe a coarsegraining approach that faithfully captures cellmass transport on a curved surface. We then construct a sparse mode representation of the resulting hydrodynamic fields in terms of scalar and vector spherical harmonic basis functions, discuss mode signatures of morphogenetic symmetry breaking events, and connect them to the dynamics of topological defects in the cellular flux. We validate this mode representation framework and the subsequent model inference using synthetic data of ABPs on a sphere, for which coarsegrained fields and learned models can be directly compared against analytical predictions. Finally, we infer a linear model for the mode dynamics of the experimental zebrafish data, which enables us to study the characteristics of cell interactions through kernels that couple cell density and flux and compare their features with the hydrodynamic meanfield signatures of ABPs on a sphere.
Coarsegraining of cellular dynamics on a spherical surface
The experimentally observed cell motions are approximately twodimensional (2D): The radius of the yolk cell surface on which the dynamics takes place is much larger than the average height changes of the evolving cell mass (Shah et al., 2019). We therefore adopt a thin film approximation, in which the cellular motion is represented on an effective spherical midsurface (gray surface in Figure 1B); refined future models should aim to account for the full 3D dynamics. Focusing here on the inplane dynamics, we project all cell positions and velocities onto a spherical midsurface $S$ of radius $R}_{s}=300\phantom{\rule{thinmathspace}{0ex}}\mu \text{m$. On this spherical surface, each cell $\alpha =1,2,\mathrm{\dots},N$ has a position ${\mathbf{r}}_{\alpha}(t)$ and inplane velocity ${\mathbf{v}}_{\alpha}(t)=\mathrm{d}{\mathbf{r}}_{\alpha}/\mathrm{d}t$.
As a second processing step, a coarsegrained representation of the singlecell dynamics on a spherical surface is determined. To facilitate the applicability of our framework to a wide range of experimental inputs, we propose a coarsegraining approach that can flexibly integrate cell number variations stemming from cell divisions, but also those from experimental uncertainties in cell imaging and tracking. Consequently, we first consider an idealized scenario in which the total cell number is approximately constant. In this case, mass conservation informs the construction of selfconsistent coarsegraining kernels on a spherical surface. In a second step, we describe how this approach generalizes when there are variations in the total cell number.
Consistent coarsegraining of idealized microscopic data
Our specific aim is to translate microscopic cell positions ${\mathbf{r}}_{\alpha}(t)$ and velocities ${\mathbf{v}}_{\alpha}(t)$ into a continuous cell surface density $\rho (\mathbf{r},t)$ and an associated flux $\mathbf{J}(\mathbf{r},t)$ at any point $\mathbf{r}$ of the spherical midsurface. For an approximately constant total number of cells, the fields $\rho $ and $\mathbf{J}$ are related by the mass conservation equation
Here, ${\nabla}_{S}\cdot \mathbf{J}$ denotes the inplane divergence of the cell number flux. To convert cell position ${\mathbf{r}}_{\alpha}(t)$ and velocities ${\mathbf{v}}_{\alpha}(t)$ into a normalized cell surface density $\rho (\mathbf{r},t)$ and an associated normalized flux $\mathbf{J}(\mathbf{r},t)$, we consider a kernel coarsegraining of the form (Appendix 1)
where $N$ is the total number of cells and ${\overline{\mathbf{v}}}_{\alpha}={\mathbf{v}}_{\alpha}/{\mathbf{r}}_{\alpha}$ is the angular velocity of a given cell on a reference unit sphere (Appendix 1). The kernels $K(\mathbf{r},{\mathbf{r}}^{\prime})$ and $K(\mathbf{r},{\mathbf{r}}^{\prime})$ are given by a scalar and a matrixvalued function, respectively. The matrix kernel $K(\mathbf{r},{\mathbf{r}}^{\prime})$ takes into account contributions of a particle with velocity ${\mathbf{v}}_{\alpha}$ at ${\mathbf{r}}^{\prime}={\mathbf{r}}_{\alpha}$ to nearby points $\mathbf{r}$ on the sphere, which involves an additional projection to ensure that $\mathbf{J}(\mathbf{r},t)$ is everywhere tangent to the spherical surface (Appendix 1). Importantly, the mass conservation Equation 1 implies a nontrivial consistency relation between the kernels $K(\mathbf{r},{\mathbf{r}}^{\prime})$ and $K(\mathbf{r},{\mathbf{r}}^{\prime})$ in Equation 2a, Equation 2b. The kernels that obey this condition represent different coarsegraining length scales (Appendix 1—figure 2). Throughout, we fix an intermediate coarsegraining length scale to enable a sparse representation of the experimental data, while ensuring that spatial details of the dynamics remain sufficiently well resolved. The final surface density $\rho (\mathbf{r},t)$ and the associated normalized flux $\mathbf{J}(\mathbf{r},t)$, computed from Equation 2a and Equation 2b using a kernel with an effective greatcircle coarsegraining width of $\sim 70\mu \text{m}$, are shown in Figure 1C (see also Video 1).
Consequences of cell number variations in experimental data
Because cell divisions are essential to most developmental processes, total cell numbers will in many cases – including early zebrafish gastrulation (Kobitski et al., 2015) – vary over time. True cell numbers and cell number changes are often difficult to measure due to experimental uncertainties arising from singlecell imaging and tracking within dense cellular aggregates. We therefore merely assume here that single cells are tracked in a representative fashion so that local relative surface densities found from Equation 2a reflect the probability that cells are present at a given point $\mathbf{r}$. In the absence of further information on cell deaths and cell divisions, we additionally make the more restrictive assumption that cell appearances or disappearances are everywhere proportional to the local cell density. With these assumptions, we can define a cell number surface density $\stackrel{~}{\rho}(\mathbf{r},t)=N(t)\rho (\mathbf{r},t)$, where $N(t)$ is the cell number at time $t$ and $\rho (\mathbf{r},t)$ is the normalized surface density given in Equation 2a. Similarly, a cell number flux is given by $\stackrel{~}{\mathbf{J}}(\mathbf{r},t)=N(t)\mathbf{J}(\mathbf{r},t)$, where the flux $\mathbf{J}(\mathbf{r},t)$ is computed from the data as described by Equation 2b. Using these definitions in Equation 1, we find that the fields $\stackrel{~}{\rho}(\mathbf{r},t)$ and $\stackrel{~}{\mathbf{J}}(\mathbf{r},t)$ obey a continuity equation
where $k(t)=\dot{N}(t)/N(t)$ denotes a timedependent effective growth rate. Importantly, under the two above assumptions, Equation 3 encodes for any timedependent total cell number $N(t)\S gt;0$ the same information as Equation 1 for coarsegrained normalized surface density $\rho (\mathbf{r},t)$ and associated flux $\mathbf{J}(\mathbf{r},t)$ given by Equation 2a and Equation 2b, respectively. In the following analysis, we hence focus on these normalized fields.
Spatial mode representation on a spherical surface
To obtain a sparse mode representation of the hydrodynamic fields $\rho (\mathbf{r},t)$ and $\mathbf{J}(\mathbf{r},t)$ on the spherical surface, we expand them in terms of scalar and vector spherical harmonics (SHs) (Arfken et al., 2013; Sandberg, 1978) (Appendix 2.A). SHs are defined on points $\widehat{\mathbf{r}}=\mathbf{r}/{R}_{s}$ of the unit sphere, where ${R}_{s}=300$ μm is the midsurface radius. In this basis, the scalar density field is represented as
which conveniently separates the time and spacedependence of $\rho (\mathbf{r},t)$ into mode amplitudes ${\rho}_{lm}(t)$ and scalar harmonic functions ${Y}_{lm}(\widehat{\mathbf{r}})$, respectively. The maximal mode number ${l}_{\text{max}}$ is a proxy for the maximal spatial resolution at which $\rho (\mathbf{r},t)$ is faithfully represented. Similarly, the vectorvalued flux $\mathbf{J}(\mathbf{r},t)$ can be decomposed into timedependent mode amplitudes ${j}_{lm}^{(1)}(t)$ and ${j}_{lm}^{(2)}(t)$, while its spatial dependence is described by vector SHs ${\mathbf{\Psi}}_{lm}(\hat{\mathbf{r}})$ and ${\mathbf{\Phi}}_{lm}(\hat{\mathbf{r}})$ (Sandberg, 1978) (Appendix 2, Video 2).
Besides the inplane divergence ${\nabla}_{S}\cdot \mathbf{J}$ that leads to local density changes (see Equation 1), the cell number flux $\mathbf{J}(\mathbf{r},t)$ also contains an inplane curl component ${\nabla}_{S}\times \mathbf{J}$ that is associated with locally rotational cell flux. The two sets of vector SHs $\{{\mathbf{\Psi}}_{lm}\}$ and $\{{\mathbf{\Phi}}_{lm}\}$ conveniently decompose the flux into these contributions: Because ${\mathrm{\nabla}}_{\mathcal{S}}\cdot {\mathbf{\Phi}}_{lm}={\mathrm{\nabla}}_{\mathcal{S}}\times {\mathbf{\Psi}}_{lm}=0$, and $\hat{\mathbf{r}}\cdot \left({\mathrm{\nabla}}_{\mathcal{S}}\times {\mathbf{\Phi}}_{lm}\right)={\mathrm{\nabla}}_{\mathcal{S}}\cdot {\mathbf{\Psi}}_{lm}=l(l+1){Y}_{lm}/{R}_{s}$ (Sandberg, 1978), we see from Equation 5 that ${j}_{lm}^{(1)}(t)$ corresponds to modes that drive density changes and ${j}_{lm}^{(2)}(t)$ represents modes of local rotational cell motion that change relative cell positions but do not change local density. Indeed, using harmonic mode representations of the cell number density Equation 4 and the cell number flux Equation 5 directly in the continuity Equation 1, we find a set of ordinary differential equation in mode space
where $l=0,1,\mathrm{\dots},{l}_{\text{max}}$ and for each value of $l$, $m=l,l+1,\mathrm{\dots},l1,l$. Equation 6 offers an alternative way of determining the modes ${j}_{lm}^{(1)}(t)$ directly from the modes ${\rho}_{lm}(t)$ of the coarsegrained cell number density (see Equation 4 and Equation 2a), while ensuring that the resulting fields obey mass conservation exactly. In practice, the modes ${j}_{lm}^{(1)}(t)$ found from a vector harmonic representation of the coarsegrained cell number flux (Equation 2b) will often deviate from modes ${j}_{lm}^{(1)}(t)$ determined from Equation 6, even if cell numbers are expected to be conserved. This can be, for example, due to limited accuracy in determining velocities ${\mathbf{v}}_{\alpha}(t)$ from noisy singlecell trajectories ${\mathbf{r}}_{\alpha}(t)$, or due to spatially inhomogeneous appearances and disappearances of cells in tracking data. Consistent with our simplifying assumption that cell number changes in the data can be sufficiently well approximated by a globally homogeneous growth rate (compare Equation 1 with Equation 3), the subsequent analysis uses the modes ${j}_{lm}^{(1)}(t)$ as determined from the density modes ${\rho}_{lm}(t)$ via Equation 6, together with modes ${j}_{lm}^{(2)}(t)$ from the explicit velocity coarsegraining Equation 2b. The complete construction is detailed in Appendix 2 and the full coarsegrained dynamics is shown in Video 1.
The representation of $\rho (\mathbf{r},t)$ and $\mathbf{J}(\mathbf{r},t)$ in terms of spherical harmonic modes with $l\le {l}_{\text{max}}$ leads in total to $3{({l}_{\text{max}}+1)}^{2}$ mode amplitude trajectories, displaying only a few dominant contributions (Figure 1D) with almost no signal remaining for $l\ge 5$ (Figure 1—figure supplement 1, Video 2). This demonstrates that the underlying coarsegrained experimental data is sufficiently smooth and implies that a spectral representations is indeed meaningful. Thus, the coarsegraining approach outlined above provides a sparse spectral representation of highdimensional microscopic singlecell data. The associated harmonic basis functions and vectors have an intuitive physical meaning, convenient algebraic properties and, as we will see, encode information about the length scales and symmetries of the collective dynamics.
Temporal mode representation
We further compress the dynamical information by representing the time series of the modes in terms of Chebyshev polynomial basis functions ${T}_{n}(t)$ (Driscoll et al., 2014; Mason and Handscomb, 2002). To simplify notation, we define a dynamic mode vector $\mathbf{a}(t)={[{\rho}_{lm}(t),{j}_{lm}^{(1)}(t),{j}_{lm}^{(2)}(t)]}^{\top}$ that collects all the modes up to $l={l}_{\text{max}}$ determined in the previous section and consider an expansion
in terms of the spatiotemporal mode coefficients ${\widehat{\mathbf{a}}}_{n}$ with temporal mode number $n$ (Appendix 2). This compression allows us to accurately evaluate time derivatives of the mode amplitudes (Supekar et al., 2021), an important step when using Equation 6 to determine flux modes ${j}_{lm}^{(1)}(t)$ directly from density modes ${\rho}_{lm}$. Fixing ${l}_{\text{max}}=4$ and ${n}_{\text{max}}=30$ in the remainder, the initial singlecell data set of about 1.4 million recorded cell position entries, or 4.2 million degrees of freedom, has thus been reduced to 2250 mode coefficients, corresponding to a compression ratio $\gtrsim 1800$.
Characterization of the developmental mode dynamics
A harmonic mode decomposition naturally integrates the geometry of the underlying domain and simultaneously provides useful insights into spatial scales and symmetries of the dynamics. For each mode $(lm)$ in the sets of SHs $\{{Y}_{lm}\}$, $\{{\mathbf{\Psi}}_{lm}\}$ and $\{{\mathbf{\Phi}}_{lm}\}$, the integer index $l$ indicates the spatial scale of the harmonic, with $l=0$ being a constant and larger $l$ indicating progressively finer spatial scales. The second index $m\in \{l,l+1,\mathrm{\dots},l\}$ provides additional information about the orientation of the harmonic scalar function or vector field. The modes $l=1$ and $l=2$ are particularly useful for characterizing the symmetry of spatial patterns on a spherical surface (Mietke et al., 2019; Scholich et al., 2020): Modes with $l=1$ indicate patterns with a global polar symmetry, whereas modes with $l=2$ represent spatial patterns with a global nematic symmetry. We now exploit these features for a detailed characterization of the symmetry breaking that takes place during cellular rearrangements and to study the properties of the cellular flux in more detail. To this end, we discuss spatial averages
of different realspace observables $O(\mathbf{r},t)$ over the midsurface $S$.
Mode signatures of developmental symmetry breaking
To study how different developmental stages and their associated symmetry breaking events are reflected in the mode representation, we first consider the average cell surface density fluctuations
For each mode $l$, the power spectrum ${P}_{\rho ,l}(t)={\sum}_{m=l}^{l}{\rho}_{lm}^{2}(t)$ in Equation 9 provides a rotationally invariant quantity (Çetingül et al., 2012; Schwab et al., 2013) that can effectively serve as an order parameter to characterize the symmetry of cell density patterns on the spherical surface. The dynamics of the density fluctuations given in Equation 9 broken down into contributions ${P}_{\rho ,l}(t)$ from each mode $l\le {l}_{\text{max}}=4$ is shown in Figure 2B. Several features of this representation are particularly striking and can be directly related to specific developmental stages. First, patterns of cell surface density fluctuations evolve from a dominantly polar symmetry ($l=1$) into density patterns with a prominent nematic symmetry ($l=2$). These mode signatures intuitively reflect the essential symmetry breaking that takes place when cells collectively reorganize from an initially localized cell dome (Figure 1B, 52 min) into an elongated shape that wraps in an open ringlike pattern around the yolk cell (Figure 1B, 760 min). Second, during this transition at around 300 min (9 hpf) (black triangle in Figure 2B), the cell surface density is most homogeneous as fluctuations become minimal for all modes $l$. Interestingly, this time point approximately marks the completion of epiboly, when the different cell layers have fully engulfed the yolk. Finally, although in a less pronounced manner, the power spectrum of the mode $l=4$ also exhibits an increased amplitude towards later times, indicating the formation of structures at finer spatial scales as development progresses. We find that mode signatures of the symmetry breaking and progression through developmental stages are robust (Figure 2—figure supplement 1B, D), illustrating that modebased analysis can provide a systematic and meaningful characterization of developmental symmetry breaking events.
Mode signatures of emergent topological defects in cellular flux
The vectorial nature of the cell number flux $\mathbf{J}(\mathbf{r},t)$ on a spherical surface implies the presence of topological defects (colored circles in Figure 2A, see Materials and methods) (Kamien, 2002). Several recent experimental results pertaining to the selforganization of multicellular systems suggest an important role of such topological defects in organizing morphogenetic events (Doostmohammadi et al., 2016; Saw et al., 2017; Guillamat et al., 2020; Copenhagen et al., 2020; Meacock et al., 2020; MaroudasSacks et al., 2020). We therefore analyze how defects within the cell number flux $\mathbf{J}(\mathbf{r},t)$ are dynamically organized during early zebrafish gastrulation and if signatures of defect formation and annihilation are present in the mode representation Equation 5. We first consider the average squared divergence and curl of the cell number flux given by
which are shown in Figure 2C (top). The two contributions to the collective cellular dynamics – locally compressible, divergent flux quantified by the divergence ${\nabla}_{S}\cdot \mathbf{J}$ and locally incompressible, rotational cell motion characterized by the curl ${\nabla}_{S}\times \mathbf{J}$ – are independently determined by the modes ${j}_{lm}^{(1)}(t)$ and ${j}_{lm}^{(2)}(t)$. Therefore, each contribution can be evaluated conveniently and with high accuracy from a representation of $\mathbf{J}(\mathbf{r},t)$ in terms of vector SHs. From Figure 2C (top), we see that the most significant divergent flux (blue curve) occurs around 300 min at the transition from epiboly towards the convergence and extension stage. A quantification of the incompressible rotational flux relative to the total cell number flux is shown in Figure 2C (bottom), where we plotted the relative curl amplitude
This measure suggests a correlation between incompressible rotational cell motion and the occurrence of topological defects (circles in Figure 2A) in the cell flux $\mathbf{J}(\mathbf{r},t)$. The total number of topological defects present at any time point is depicted in Figure 2C (bottom, blue curve). Because the vectorvalued flux is defined on a sphere, we observe that the total topological charge always sums to +2 (Kamien, 2002), while additional defect pairs with opposite charge (red +1 and white 1 circles in Figure 2A) can be created, resulting in total defect numbers greater than two (see Figure 2C, bottom). Interestingly, the relative curl amplitude ${S}_{\text{curl}}$ defined in Equation 11 indicates that increased contributions from incompressible rotational flux are associated with the formation of topological defects in the cell number flux, a feature that is robustly identified by our framework (Figure 2—figure supplement 1A, C, Figure 2—figure supplement 3, Figure 2—figure supplement 4). The appearance of additional defects at the end of epiboly, when the developing embryo begins to extrude more significantly in the radial direction, suggests that topological defects in the 2D projected cellular flux fields could signal the start of the formation of more complex structures in three dimensions.
Learning a linear hydrodynamic model of the developmental mode dynamics
The results in Figure 2 confirm that a lowdimensional mode representation can capture essential characteristics of developmental symmetry breaking processes. The mode representation therefore provides a natural starting point for the inference of hydrodynamic models from coarsegrained celltracking data. For a given timedependent mode vector $\mathbf{a}(t)={[{\rho}_{lm}(t),{j}_{lm}^{(1)}(t),{j}_{lm}^{(2)}(t)]}^{\top}$ that contains all modes up to $l={l}_{\text{max}}$, the simplest hydrodynamic model corresponds to the linear dynamical equation
where the constant coefficient matrix $M$ encodes the couplings between different modes. Intuitively, Equation 12 aims to describe an experimentally observed density and flux dynamics in terms of a relaxation process, starting from inhomogeneous initial conditions represented by $\mathbf{a}(0)$. The mathematical learning problem is then to find a coefficient matrix $M$ such that the linear model Equation 12 holds for the mode vector time series $\mathbf{a}(t)$ that was determined from the coarsegraining procedure described in the previous sections.
Validation of the learning framework using active Brownian particle dynamics
Before applying the combined coarsegraining and inference framework to experimental data, we illustrate and validate the learning approach on synthetic data for which coarsegraining results and hydrodynamic meanfield equations are analytically tractable. To this end, we consider the stochastic dynamics of noninteracting active Brownian particles (ABPs) on the unit sphere of radius ${R}_{0}=1$ (Sknepnek and Henkes, 2015; Fily et al., 2016; CastroVillarreal and Sevilla, 2018). Similar to a migrating cell, an ABP at position $\mathbf{x}(t)$ moves across the unit sphere at constant speed v_{0} in the direction of its fluctuating orientation unit vector $\mathbf{u}(t)$. The strength of the orientational Gaussian white noise is characterized by a rotational diffusion constant ${D}_{r}$ (Figure 3A, Appendix 3).
Compared with conventional passive Brownian motion, selfpropulsion of an ABP along its orientation direction $\mathbf{u}$ introduces a persistence to the particle’s motion that is reduced as rotational noise ${D}_{r}$ is increased. Additionally, the topology of the spherical surface implies that in the lownoise regime, ${R}_{0}{D}_{r}/{v}_{0}\S lt;1$, particles are expected to return to the vicinity of their starting points after a duration $\mathrm{\Delta}t\approx 2\pi {R}_{0}/{v}_{0}$. The conjunction of persistent motion and topology then leads to oscillatory dynamics in the positional correlation $\u27e8\mathbf{x}(t)\cdot \mathbf{x}(0)\u27e9$ (blue dots in Figure 3B, Appendix 3). Comparing correlations from stochastic ABP simulations in different noise regimes with theoretical predictions (solid lines in Figure 3B) validates our numerical ABP simulation scheme.
To generate a test data set for our coarsegraining and inference framework, we simulated noninteracting ABPs in both the lownoise (${R}_{0}{D}_{r}/{v}_{0}\S lt;1$) and the highnoise (${R}_{0}{D}_{r}/{v}_{0}\S gt;1$) regime with initial positions drawn from the experimental data shown in Figure 1. Specifically, at each cell position present in the data, we generated 60 particles with random orientation, amounting to approximately $1.2\times {10}^{5}$ particles in total, and simulated their dynamics on a unit sphere. The resulting trajectory data were coarsegrained following the procedure outlined in the previous sections, yielding dynamic density fields $\rho (\mathbf{r},t)$ and fluxes $\mathbf{J}(\mathbf{r},t)$ (Video 3), together with their mode representations ${\rho}_{lm}(t),{j}_{lm}^{(1)}(t)$ and ${j}_{lm}^{(2)}(t)$.
In the second ‘learning’ step, we infer a sparse mode coupling matrix $M$ that approximates the dynamics Equation 12 for the dynamical mode vectors $\mathbf{a}(t)={[{\rho}_{lm},{j}_{lm}^{(1)},{j}_{lm}^{(2)}]}^{\top}$ obtained from the coarsegrained simulated ABP data. Our inference algorithm combines adjoint techniques (Rackauckas et al., 2021) and a multistep sequential thresholding approach inspired by the Sparse Identification of Nonlinear Dynamics (SINDy) algorithm introduced by Brunton et al., 2016. The full algorithm is detailed in Appendix 4 and illustrated in the summary flowchart Appendix 4—figure 1. Importantly, we perform the sparse regression using dynamical mode vectors $\mathbf{a}(t)$ rescaled by their median absolute deviation (MAD) to compensate for substantial scale variations between different modes. The final output matrix $M$ of this learning algorithm is shown in the right panel of Figure 3C and can be compared against the analytically coarsegrained dynamics of ABPs on curved surfaces (Fily et al., 2016; CastroVillarreal and Sevilla, 2018). Under suitable closure assumptions (Appendix 3), the meanfield dynamics of ABPs on a unit sphere is given in harmonic mode space by
from which we can read off the mode coupling matrix $M$ shown in the left panel of Figure 3C. A direct comparison between the theoretical and the inferred matrices shows that our framework recovers both the structure and the quantitative values of $M$ with good accuracy. Due to the finite number of ABPs used to determine the coarsegrained fields, we do not expect that the theoretically predicted coupling matrix is recovered perfectly from the data. Instead, some mode couplings suggested by Equation 13a may not be present or modified in the particular realization of the ABP dynamics that was coarsegrained. Indeed, direct simulation of the learned model projected in real space (Figure 3E) reveals a density and flux dynamics that agrees very well with the dynamics of the the coarsegrained input data (Figure 3D). Altogether, these results demonstrate that the proposed inference framework enables us to to faithfully recover expected meanfield dynamics from coarsegrained fields of noisy particlebased data.
Learning developmental mode dynamics from experimental data
The same inference framework can now be directly applied to the coarsegrained experimental zebrafish embryo data shown in Figure 1C and D, yielding a sparse coefficient matrix $M$ (Figure 4A and B) that encodes the dynamics of the developmental mode vector $\mathbf{a}(t)={[{\rho}_{lm}(t),{j}_{lm}^{(1)}(t),{j}_{lm}^{(2)}(t)]}^{\top}$ according to Equation 12. The inferred coupling between the time derivative of density modes ${\rho}_{lm}$ and flux modes ${j}_{lm}^{(1)}$ faithfully recovers mass conservation (Figure 4C; see Equation 6). Overall, the learned matrix $M$ has 395 nonzero elements, effectively providing further compression of the experimental data, which required 2,250 spatiotemporal mode coefficients collected in ${\widehat{\mathbf{a}}}_{n}$ (see Equation 7) for its representation. Using the mode vector $\mathbf{a}(t=0)$ of the first experimental time point as the initial condition, the inferred minimal model Equation 12 with $M$ shown in (Figure 4A and B) faithfully recovers both the mode and realspace dynamics seen in the coarsegrained fields of the experimental input data (Figure 4E–G, Video 4).
It is instructive to analyze the inferred matrix $M$ and the linear model it encodes in more detail. Comparing the MADrescaled matrix (see Appendix 4) learned for the experimental zebrafish data (Figure 4B) with the nondimensionalized matrix learned for the active Brownian particle dynamics (Figure 3C), we find similar patterns of prominent diagonal and blockdiagonal couplings. Consistent with the analysis of single cell trajectories (Shah et al., 2019), this suggests that a random, but persistent movement of cells akin to ABPs moving on a sphere partially contributes to the early gastrulation process in zebrafish. This is complemented in the minimal model of the experimental dynamics by significant offdiagonal contributions (Figure 4B), which are absent in the noninteracting ABP model. Such offdiagonal contributions represent effective linear approximations of cellcell interactions, environmental influences or other external stimuli reflected in the experimental timeseries data. Ultimately, such contributions to the mode coupling matrix $M$ help realize the symmetry breaking process observed in the underlying experimental data (Figure 2).
The inferred mode coupling matrix $M$ shown in Figure 4B together with Equation 12 provides a highly robust minimal model. Specifically, despite being linear, it is numerically stable over a period approximately four times as long as the input data from which the matrix $M$ was learned. Furthermore, simulations with modified initial conditions (see Figure 4—figure supplement 1) still exhibit a characteristic symmetry breaking and lead to the emergence of density and flux patterns similar to those seen in Figure 4F and G. For example, simulating Equation 12 using the initial condition of a different experimental data set (Figure 2—figure supplement 1) leads to final patterns with the same symmetry as in the original training data, further corroborating that the observed symmetry breaking is directly encoded in the interactions represented by the matrix $M$. A similar robustness is observed under moderate perturbations of the initial condition, such as a rotation of initial cell density patterns relative to the coordinate system in which $M$ was inferred, or a local depletion of the initial density, emulating a partial removal of cells as experimentally realized in Morita et al., 2017. Taken together, these numerical experiments demonstrate that the inferred mode coupling matrix $M$ meaningfully captures the dynamics and interactions of cells that facilitate the symmetry breaking observed during early zebrafish development.
Green’s function representation of learned models in real space
To characterize the inferred spatial interactions in more detail, we can analyze the realspace representation of the learned mode coupling matrix $M$. While the density dynamics represented by $M$ (the first row in Figure 4AB) simply reflects mass conservation Equation 1 in real space, the dynamics of the flux (the second and third row in Figure 4A and B) corresponds in real space to the integral equation (Appendix 4)
where $d{\mathrm{\Omega}}^{\prime}=\mathrm{sin}{\theta}^{\prime}d{\theta}^{\prime}d{\varphi}^{\prime}$ is the spherical surface area element. The vectorvalued kernel ${\mathbf{m}}^{\rho}(\mathbf{r},{\mathbf{r}}^{\prime})$ in Equation 14 connects the distribution of cell density $\rho $ across the surface to dynamic changes of the flux $\mathbf{J}$ at a given point $\mathbf{r}$. Similarly, the matrixvalued kernel ${M}^{J}(\mathbf{r},{\mathbf{r}}^{\prime})$ describes how the distribution of cell fluxes at ${\mathbf{r}}^{\prime}$ affects temporal changes of the flux at $\mathbf{r}$.
To analyze the spatial range of interactions between points $\mathbf{r}$ and ${\mathbf{r}}^{\prime}$, we use the fact that the matrixvalued kernel ${M}^{J}(\mathbf{r},{\mathbf{r}}^{\prime})$ has only one nonzero eigenvalue (Appendix 4—figure 2). Consequently, the trace $\text{tr}({M}^{J})$ serves as a proxy for the distancedependent interaction strength mediated by ${M}^{J}$. Averages of $\text{tr}({M}^{J})$ over pointpairs with the same angular distance $\omega =\text{acos}(\mathbf{r}\cdot {\mathbf{r}}^{\prime})$ are shown for the ABP dynamics and for the minimal model inferred from experimental data in Figure 4D. Note that to make the models amenable to comparison, we compute ${M}^{J}(\mathbf{r},{\mathbf{r}}^{\prime})$ from the known meanfield model of ABPs (Equation 13a) using the same finite number of modes as used to represent the ABP and the zebrafish data (${l}_{\text{max}}=4$). In theory, one expects for the ABP dynamics a highly localized, homogeneous kernel $\text{tr}({M}^{J})\sim \delta (\mathbf{r}{\mathbf{r}}^{\prime})$, so that an exact spectral representation would require an infinite number of modes (see Appendix 4). In practice, using a finite number of modes leads to a wider kernel range (Figure 4D ’ABP theory’) and introduces an apparent spatial inhomogeneity, as indicated by the nonzero standard deviation of $\text{tr}({M}^{J})$ at fixed distance $\omega $ (blue shades). Both the quantitative profile of $\text{tr}({M}^{J})$ and its variation are successfully recovered by applying the inference framework to stochastic simulations of ABPs (Figure 4D ’ABP simulation’) where ${M}^{J}(\mathbf{r},{\mathbf{r}}^{\prime})$ was computed from the learned mode coupling matrix $M$ shown in Figure 3C. For the inferred minimal model of the cell dynamics (Figure 4D ’Zebrafish experiment’), we find a similar shortranged fluxflux coupling mediated by ${M}^{J}$. However, the increased variability of $\text{tr}({M}^{J})$ at fixed distances $\omega $ indicates more substantial spatial inhomogeneities of the corresponding interactions. These inhomogeneities are absent in a noninteracting system of ABPs and represent an interpretable realspace signature of the symmetrybreaking mechanisms built into the underlying mode coupling matrix $M$.
A similar analysis can be performed for the kernel ${\mathbf{m}}^{\rho}(\mathbf{r},{\mathbf{r}}^{\prime})$ that couples the density at position ${\mathbf{r}}^{\prime}$ to dynamics of fluxes at position $\mathbf{r}$ (see Equation 14), where we average the magnitude ${\mathbf{m}}^{\rho}(\mathbf{r},{\mathbf{r}}^{\prime})$ over pairs ($\mathbf{r}$, ${\mathbf{r}}^{\prime}$) with the same angular distance $\omega $ (Figure 4D insets). Using a finite number of modes to compute this kernel in the different scenarios again introduces apparent spatial inhomogeneities in all cases. Additionally, all kernel profiles exhibit a distinct maximum at short range, indicating a coupling between density gradients and the flux dynamics that emerges microscopically from a persistent ABP and cell motion (see Appendix 3 and 4) – an observations that is consistent with the similar blockdiagonal structure of both inferred matrices $M$ (compare Figures 3C and 4B).
In conclusion, the realspace analysis and comparison of inferred interaction kernels further highlights potential ABPlike contributions to the collective cellular organization during early zebrafish development and reveals an effectively nonlocal coupling between density and flux dynamics. The latter could result, for example, from unresolved fastevolving morphogens (Hannezo and Heisenberg, 2019), through mechanical interactions with the surrounding material (Münster et al., 2019) or due to other relevant degrees of freedom that are not explicitly captured in this linear hydrodynamic model. More generally, a realspace representation of kernels provides an alternative interpretable way to study the interactions and symmetrybreaking mechanisms encoded by models directly learned in mode space.
Discussion
Leveraging a sparse mode representation of collective cellular dynamics on a curved surface, we have presented a learning framework that translates singlecell trajectories into quantitative hydrodynamic models. This work complements traditional approaches to find quantitative continuum models of complex multicellular processes (Etournay et al., 2015; Hannezo et al., 2015; Morita et al., 2017; Streichan et al., 2018; Münster et al., 2019) that match problemspecific constitutive relations of active materials in realspace with experimental observations. We have demonstrated here that length scales and symmetries associated with a mode representation can directly inform about the character of symmetry breaking transitions and topological features of collective cellular motion even before a model is specified. The successful applications to synthetic ABP simulation data and experimental zebrafish embryo data show that model learning in mode space provides a promising and computationally feasible approach to infer quantitative interpretable models in complex geometries.
The learned linear minimal model for cell migration during early zebrafish morphogenesis quantitatively recapitulates the spatiotemporal dynamics of a complex developmental process (Figure 4F and G), and highlights similarities between collective cell migration and analytically tractable ABP dynamics on a curved surface. An extension to nonlinear modecoupling models or an integration of additional, experimentally measured degrees of freedom, such as concentration fields of morphogens involved in mechanochemical feedbacks (Hannezo and Heisenberg, 2019), is in principle straightforward by including nonlinear terms in Equation 12. Furthermore, the above framework could be generalized to describe the dynamics within a spherical shell of finite height by complementing the surface vector SHs used in this work by their radial counterpart (Barrera et al., 1985).
To provide a concrete example, we focused here on applying the model learning framework to singlecell tracking data of early zebrafish morphogenesis. However, the essentially spherical organization of cells during gastrulation observed in zebrafish is shared by many species whose early development occurs through a similar discoidal cleavage (Gilbert and Barresi, 2016), and the framework introduced here is directly applicable once tracking data becomes available for these systems. More generally, as novel imaging technologies are being developed (Keller et al., 2010; Royer et al., 2016; Shah et al., 2019), we expect that even larger and more detailed imaging data will further facilitate the exploration of finer scales and lengthscale bridging processes (Lenne et al., 2021) through learning approaches that directly built on modebased data representations.
Materials and methods
Data preprocessing
Request a detailed protocolWe obtained two singlecell tracking data sets from the experiments described in Shah et al., 2019. These data consist of the Cartesian coordinates of each cell together with a tracking ID. Some of the data is accessible at https://idr.openmicroscopy.org with ID number idr0068. We first denoised each cell trajectory using MATLAB’s (Matlab, 2019) wavelet denoiser function wdenoise, and centered the cloud of cells by leastsquares fitting a spherical surface through it and shifting the origin at each time to coincide with the center of this sphere. We then computed the velocity of each cell by using Tikhonovregularized differentiation as described in Knowles and Renka, 2014 and implemented in the MATLAB thirdparty module rdiff (Wagner, 2020). After examination of the cells’ velocity distribution, we further removed outlier cells whose speed is in the 95th percentile or above and verified that this operation only removes aberrant cells. Finally, we rotated the data to align the animal pole of the embryo with the $z$axis, as determined by the direction of the center of mass of the initial cell distribution. The resulting single cell data are shown as point clouds in Figure 1B and in Video 1.
Topological defect tracking
Request a detailed protocolWe have developed a defect tracker that identifies topological defects in vector fields tangent to a spherical surface via integration along suitable Burger circuits. The corresponding software is available at (https://github.com/NicoRomeo/surfvecdefects; Romeo, 2022, copy archived at swh:1:rev:6dc742c376b0d085e19ece65f932ac6935342aba).
Appendix 1
Consistent coarsegraining on curved surfaces
We describe the derivation of selfconsistent coarsegraining kernels that are used in the main text to convert single cell information into a continuous density field and its associated fluxes on a spherical surface. We first motivate this problem for a flat surface and then proceed with a detailed derivation for the case of a spherical surface.
Kernel consistency in Euclidean space
It is instructive to first consider a set of particles $\alpha =1,2,3,\mathrm{\dots}$ at positions ${\mathbf{X}}_{\alpha}(t)$ moving with velocities ${\mathbf{V}}_{\alpha}(t)=\mathrm{d}{\mathbf{X}}_{\alpha}/\mathrm{d}t$, where capitalized vectors indicate position and velocity in Euclidean space, e.g. particles move on a flat surface or within some threedimensional volume. A coarsegrained density $\rho (\mathbf{X},t)$ and a mass flux $\mathbf{J}(\mathbf{X},t)$ can be defined from this microscopic information by
where ${K}_{e}(\mathbf{X},{\mathbf{X}}^{\prime})$ and ${K}_{e}(\mathbf{X},{\mathbf{X}}^{\prime})$ represent a scalarvalued and a matrixvalued kernel function, respectively. At the same time, in a system with a constant number of particles, mass conservation implies, in general,
relating the density $\rho (\mathbf{X},t)$ and the mass flux $\mathbf{J}(\mathbf{X},t)$ of particles. Using the coarsegraining prescriptions Equation 15a and Equation 15b directly in Equation 16 and assuming the resulting relation must hold for any set of particle trajectories, one finds a general kernel consistency relation
This condition is automatically satisfied for any translationally invariant and isotropic pair of kernels ${K}_{e}(\mathbf{X},{\mathbf{X}}^{\prime})={K}_{e}(\mathbf{X}{\mathbf{X}}^{\prime})$ and ${K}_{e}(\mathbf{X},{\mathbf{X}}^{\prime})={K}_{e}(\mathbf{X}{\mathbf{X}}^{\prime})I$, where $I$ is the unit matrix. Coarsegraining with such kernels is frequently employed in practice: Positions and velocities can be, for example, simply convolved with a Gaussian function of mean zero (Supekar et al., 2021).
Kernel consistency on a curved surface
For a surface parameterized by $\mathbf{r}({s}^{1},{s}^{2})\in {\mathbb{R}}^{3}$ with generalized coordinates ${s}^{1},{s}^{2}$, two tangential basis vectors are defined by ${\mathbf{e}}_{i}=\partial \mathbf{r}/\partial {s}^{i}$ ($i=1,2$). Partial derivatives are, in the following, denoted ${\partial}_{i}:=\partial /\partial {s}^{i}$. The metric tensor is given by ${g}_{ij}={\mathbf{e}}_{i}\cdot {\mathbf{e}}_{j}$. The mean curvature is defined by $H\mathbf{n}={\nabla}_{i}{\mathbf{e}}^{i}/2$, where $\mathbf{n}={\mathbf{e}}_{1}\times {\mathbf{e}}_{2}/{\mathbf{e}}_{1}\times {\mathbf{e}}_{2}$ denotes the unit surface normal and the Einstein summation convention is used. The covariant form of mass conservation Equation 1 (main text) on a curved surface reads
with ${J}^{i}={\mathbf{e}}^{i}\cdot \mathbf{J}$ and ${\nabla}_{i}$ denotes the covariant derivative. In general, we are interested in describing an effective dynamics for cell positions and velocities that are projected onto a common reference sphere of radius ${R}_{s}$. Such a description can be found by first formulating the coarsegraining approach for a unit sphere, on which particle positions and velocities are fully determined by angular coordinates and corresponding angular velocities, and finally rescaling the density and flux fields by suitable factors of ${R}_{s}$. The corresponding coarsegraining Equation 2b (main text) of inplane angular velocities ${\overline{\mathbf{v}}}_{\alpha}(t)={\mathbf{v}}_{\alpha}(t)/{\mathbf{r}}_{\alpha}(t)$ for particles $\alpha $ on a unit sphere reads covariantly
where ${\overline{v}}_{\alpha}^{i}={\mathbf{e}}^{i}\cdot {\overline{\mathbf{v}}}_{\alpha}$ and we drop the dependence on time to simplify the notation. The twopoint kernel tensor $K{(\mathbf{r},{\mathbf{r}}^{\prime})}_{i{j}^{\prime}}$ (a ‘bitensor’) is evaluated in the tangent space of $\mathbf{r}$ for its first index and in the tangent space of ${\mathbf{r}}^{\prime}$ at the second, primed index (Appendix 1—figure 1). Mass conservation on a curved surface, Equation 18, together with the coarsegraining prescriptions Equation 2a (main text) and Equation 19 then implies a covariant kernel consistency relation
Solving the kernel consistency relation on a sphere
We solve Equation 20 in the following on the unit sphere, such that $\mathbf{r}=\mathbf{n}$ corresponds to the surface normal. The final result can simply be rescaled to any spherical surface of radius ${R}_{s}$. Furthermore, we note that the parameter
provides a measure for the great circle distance $\omega (x)=\text{acos}(x)$ between two points on a sphere. Hence, we consider an ansatz for the kernels in Equation 20 of the form
with two unknown scalar functions $f(x)$ and $g(x)$. The relevant derivatives of the ansatz Equation 22a and Equation 22b can readily be evaluated to
Here, ⊗ denotes a dyadic product and we use ${\partial}_{i}x={\mathbf{r}}^{\prime}\cdot {\mathbf{e}}_{i}$ and ${\partial}_{{i}^{\prime}}x=\mathbf{r}\cdot {\mathbf{e}}_{{i}^{\prime}}$, which follows from Equation 21, as well as ${\nabla}_{i}{\mathbf{e}}^{i}=2\mathbf{r}$ in the second equation, which holds on a unit sphere and follows from the definition of the mean curvature. We then use the expansion of the identity matrix in ${\mathbb{R}}^{3}$ on the spherical basis $I={\mathbf{e}}_{i}\otimes {\mathbf{e}}^{i}+\mathbf{n}\otimes \mathbf{n}$, such that in our case with $\mathbf{r}=\mathbf{n}$ we have ${\mathbf{e}}_{i}\otimes {\mathbf{e}}^{i}=I\mathbf{r}\otimes \mathbf{r}$. Hence, Equation 23b becomes
Using Equation 24 and Equation 23a in the kernel consistency relation Equation 20 and dividing by $\mathbf{r}\cdot {\mathbf{e}}_{{j}^{\prime}}$ (at $\mathbf{r}={\mathbf{r}}^{\prime}$, for which $\mathbf{r}\cdot {\mathbf{e}}_{{j}^{\prime}}=0$, Equation 20 is obeyed for any $f(x)$, $g(x)$), we find that the scalar functions in the kernel ansatz Equation 22a and Equation 22b have to obey
Hence, the general covariant consistency relation Equation 20 implies for the kernel ansatz Equation 22a and Equation 22b that the weighting functions $g(x)$ and $f(x)$ must be related by
Kernel functions with compact support
In the last step, we determine a family of kernel functions $g(x)$ and $f(x)$ defined on the interval $x\in [1,1]$ that satisfy Equation 25, along with the requirements:
$f(x)$ and $g(x)$ must be ${C}^{1}$ regular on $[1,1]$
$f\ge 0$ on $[1,1]$
$f$ is normalized to one on the unit sphere.
Recalling $x=\mathrm{cos}[\omega (\mathbf{r},{\mathbf{r}}^{\prime})]$ with angular distance $\omega $ between $\mathbf{r}$ and ${\mathbf{r}}^{\prime}$, a family of functions fulfilling these conditions is given by
where $\mathbf{1}}_{\{\mathrm{cos}\omega >0\}$ is an indicator function that is one if $\mathrm{cos}\omega \S gt;0$ and vanishes otherwise (Appendix 1—figure 2). In this work, we have chosen the kernels Equation 22a and Equation 22b with $f={f}_{k}$ and $g={g}_{k}$ for $k=6$. For these kernels derived here, densities $\rho (\mathbf{r},t)$ and associated fluxes $\mathbf{J}(\mathbf{r},t)$ that are coarsegrained on a unit sphere can be converted into effective densities and fluxes on a spherical surface of radius ${R}_{s}$ through the rescaling $\rho \to \rho /{R}_{s}^{2}$ and $\mathbf{J}\to \mathbf{J}/{R}_{s}$. Equivalently, rescaled kernels $K(\mathbf{r},{\mathbf{r}}^{\prime})\to K(\mathbf{r},{\mathbf{r}}^{\prime})/{R}_{s}^{2}$ and $K{(\mathbf{r},{\mathbf{r}}^{\prime})}_{i{j}^{\prime}}\to K{(\mathbf{r},{\mathbf{r}}^{\prime})}_{i{j}^{\prime}}/{R}_{s}$ can be used directly, as was done in Equation 2a and Equation 2b of the main text to generate the data shown in Figure 1 (main text).
Appendix 2
Spatiotemporal mode decomposition
In this section, we provide explicit expressions for the scalar and spherical harmonic basis functions (‘spatial modes’), as well as for the Chebyshev basis functions (‘temporal modes’) used in this work. Additionally, we describe a systematic approach to determine the minimal number of modes needed to describe the coarsegrained data, while preserving a high level of accuracy in the representation.
Spatial basis: spherical harmonics
In this work, we use the real spherical harmonics defined in spherical coordinates $(\theta ,\varphi )$ by Arfken et al., 2013 as
where ${P}_{l}^{m}(x)$ is the associated Legendre polynomial of degree $l$ and order $m$, and
Vector spherical harmonics can be defined and expressed as vector fields in 3D or covariantly as (Sandberg, 1978; Mietke et al., 2019)
where ${\nabla}_{S}={\mathbf{e}}_{\theta}{\partial}_{\theta}+{\mathbf{e}}_{\varphi}{\mathrm{sin}}^{1}\theta {\partial}_{\varphi}$ denotes the gradient operator an a unit sphere, ${\u03f5}_{ij}$ is the covariant LeviCivita tensor, and ${g}_{ij}$ the metric tensor. Scalar harmonics ${Y}_{lm}$ and either vector harmonic $\mathbf{\Lambda}}_{lm}\in \{{\mathbf{\Psi}}_{lm},{\mathbf{\Phi}}_{lm}\$ are orthogonal:
where $\mathrm{d}\mathrm{\Omega}=\mathrm{sin}\theta \mathrm{d}\theta \mathrm{d}\varphi$. The increasing complexity of patterns and accuracy of reconstruction with larger $l$ is illustrated in Appendix 2—figure 1 and Video 2.
Temporal basis: Chebyshev polynomials
Chebyshev polynomials of the first kind ${T}_{n}$ are defined by Arfken et al., 2013.
Chebyshev polynomials form an orthogonal basis of continuous functions on the interval $[1,1]$, such that an expansion
uniformly converges as ${n}_{\text{max}}\to \mathrm{\infty}$ (Driscoll et al., 2014). This representation also allows computing derivatives spectrally from
Information loss through coarsegraining
Coarsegraining microscopic data into smooth fields is an irreversible operation, during which some of the original particle information is irretrievably lost. The choice of coarsegraining scale is thus dictated by a tradeoff between smoothness and information content  choosing larger coarsegraining scales leads to smoother fields but blurs finer scale structures which may be of interest. To inform our choice of coarsegraining scale, we quantify the loss of information incurred by the coarsegraining operation.
The measure we introduce to quantify information loss is based on the the wellknown relationship between the smoothness of functions in real space and Fourier space (Stein and Shakarchi, 2011): A smooth function in real space should have a peaked, quickly decaying spectrum in Fourier space while a collection of pointlike objects such as delta functions should have a uniform nondecaying spectrum. Specifically, we describe a uniformly sampled field as a $M\times N$ matrix with components being the field values ${X}_{i,j}=X({\theta}_{i},{\varphi}_{j})$. In our case, ${X}_{i,j}$ represents either the density field $\rho $ or any of the Cartesian components of the flux vector field $\mathbf{J}$ at a given time point. We find the complex discrete Fourier spectrum ${\widehat{X}}_{i,j}$ of this matrix using the twodimensional fast Fourier transform. We then calculate the power spectral density (PSD) of the Fourier spectrum as ${R}_{i,j}={{\widehat{X}}_{i,j}}^{2}$ and interpret the normalized PSD.
as a discrete probability distribution. The spectral entropy $S$ characterizing the information content of the field $X$ is then defined by
Smooth fields are sharply peaked in Fourier space and have a low spectral entropy, whereas fields that resolve discrete single particle information are rather flat in Fourier space and have a large spectral entropy. The difference in entropy between particle data and smoothed fields then measures the information eliminated by the coarsegraining procedure. If we additionally normalize by the entropy of the spectral entropy ${S}_{0}(X)$ of the raw particle data, we finally obtain a relative measure of the information that is lost in the coarsegraining process. In general, a measure as given in Equation 34 can be defined for any transform with the property that smoothness in real space leads to a fast decaying spectrum in transform space.
We compute the spectral entropy of density and flux component fields at a representative time point and for varying coarsegraining length scales (Appendix 2—figure 2). Specifically, we coarsegrain density and flux through the procedure described in the main text and in Appendix 1 for different values of the kernel parameter $k$ (see Equation 26a). Large values of $k$ correspond to small coarsegraining length scales, with the effective halfwidth at halfmaximum (HWHM) of the kernels given in Equation 22a and Equation 22b with weight functions Equation 26a and Equation 26b scaling as HWHM $=\mathrm{arccos}({2}^{1/k})$. Normalized spectral entropies $S(X)/{S}_{0}(X)$ with $X\in \{\rho ,\mathbf{J}\}$ are then computed using Equation 34. For the flux field, we define $S(\mathbf{J}):=S({J}_{x})+S({J}_{y})+S({J}_{z})$ ("Flux sum" in Appendix 2—figure 2) and interpret the sum of these three contributions ("Flux x", "Flux y", "Flux z" in Appendix 2—figure 2) as the total information contained in the flux field. We find that the spectral entropies of all fields show similar features. In particular, an increasing coarsegraining width first results in a sharp loss of information as individual particle positions are blurred, followed by less steep information loss as continuous fields progressively lose details of finer structures. In this work, we use an intermediate value of the coarsegraining parameter $k=6$ (yellow data in Appendix 2—figure 2).
Optimal compression in space and time
Spectral representations are exact in the limit of an infinite number of modes. In practice, we choose a maximal harmonic mode number ${l}_{\text{max}}$ and maximal Chebyshev mode number ${n}_{\text{max}}$. A too large value of ${l}_{\text{max}}$ and ${n}_{\text{max}}$ provides little compression benefit, while too small values suffer accuracy penalties. Hence, there is a compressionaccuracy tradeoff that we seek to optimize. To evaluate the tradeoff quantitatively, we define a heuristic compression metric $C$ by
where ${N}_{t}$ is the number of sampled time steps and ${N}_{s}$ is the number of spatial grid points used for coarsegraining. Larger values of $C$ correspond to a higher compression factor. To define accuracy metrics, we consider the norm
where the sum runs over ${N}_{t}$ regularly sampled time points t_{i}. We denote a particular mode representation $\{{\stackrel{~}{\rho}}_{lm}(t),{\stackrel{~}{j}}_{lm}^{(1)}(t),{\stackrel{~}{j}}_{lm}^{(2)}(t)\}$ of the data that was coarsegrained via Equation 2a and Equation 2b (main text) for $l=0,\mathrm{\dots},{l}_{\text{max}}^{\text{ref}}=20$ as the ‘uncompressed’ reference. A measure to characterize the accuracy of a modetruncated ‘compressed’ data representation is then given by a relative average mode reconstruction error
This measure compares the compressed mode representation $\{{\rho}_{lm}(t),{j}_{lm}^{(2)}(t)\}$, truncated at maximal Chebychev mode number ${n}_{\text{max}}$ (temporal representation Equation 32, Appendix 2) and maximal harmonic mode number ${l}_{\text{max}}$ (spatial representation, Equations 4; 5, main text) with the reference modes $\{{\stackrel{~}{\rho}}_{lm}(t),{\stackrel{~}{j}}_{lm}^{(2)}(t)\}$. To find a compromise between accuracy, characterized by ${E}_{\text{modes}}({n}_{\text{max}},{l}_{\text{max}})$, and compression $C$ defined in Equation 35, the aim is to find a pair $({n}_{\text{max}},{l}_{\text{max}})$ on the Pareto front (Jin and Sendhoff, 2008) of ${E}_{\text{modes}}$ vs. $1/C$ (red dots in Appendix 2—figure 3).
Note that the modes ${\stackrel{~}{j}}_{lm}^{(1)}(t)$ and ${j}_{lm}^{(1)}(t)$ are so far omitted from this analysis, because the latter are in practice found directly from density modes via Equation 6 (main text). However, taking temporal derivatives of ${\rho}_{lm}(t)$ using Equation 33 to determine ${j}_{lm}^{(1)}(t)$ introduces undesirable oscillations for too large Chebychev cutoffs ${n}_{\text{max}}$. This implies an additional tradeoff between the need for accuracy (higher ${n}_{\text{max}}$) and stability (lower ${n}_{\text{max}}$). In practice, we wish to find values of $({n}_{\text{max}},{l}_{\text{max}})$ such that relative amplitudes of pairs $({\stackrel{~}{j}}_{lm}^{(1)},{\stackrel{~}{j}}_{lm}^{(2)})$ and $({j}_{lm}^{(1)},{j}_{lm}^{(2)})$ are preserved by the compression. This can be achieved by comparing the relative curl amplitude
to the analog quantity ${\stackrel{~}{S}}_{\text{curl}}(t)$ computed from the reference modes $\{{\stackrel{~}{j}}_{lm}^{(1)}(t),{\stackrel{~}{j}}_{lm}^{(2)}(t)\}$ and analyzing the curl reconstruction error
as a function of ${n}_{\text{max}}$ and ${l}_{\text{max}}$ (Appendix 2—figure 4). From this, we find a region of low error around ${l}_{\text{max}}=4,{n}_{\text{max}}=30$, which also is on the Pareto front of the accuracy vs. compression tradeoff (orange circles in Appendix 2—figures 3 and 4) and represents the final values used throughout this work.
Appendix 3
Active Brownian particles on the sphere
In this section, we describe the stochastic dynamics of noninteracting, active Brownian particles (ABPs) (Romanczuk et al., 2012) on curved surfaces and derive analytically coarsegrained meanfield equations, as well as a kernel representation of ABP dynamics. These results are used in the main text to validate our coarsegraining and inference framework.
We consider active Brownian particles at position $\mathbf{x}\in {\mathbb{R}}^{3}$ that move with speed v_{0} on the surface of a unit sphere (radius ${R}_{0}=1$) in the direction of their unit orientation vector $\mathbf{u}\in {\mathbb{R}}^{3}$. Since $\mathbf{x}=1$ at all times, we can interpret v_{0} as the particle’s angular speed on the unit sphere. The orientation vector is at all times tangential to the surface, but is subject to random inplane fluctuations characterized by a rotational diffusion coefficient ${D}_{r}$. The corresponding dynamics of $\mathbf{x}(t)$ and $\mathbf{u}(t)$ is given by the stochastic differential equations (in units ${R}_{0}=1$)
where the stochastic differential Equation 38b is interpreted in the Stratonovich sense, as denoted by the symbol "o" (Braumann, 2007). It follows from Equation 38a that $\mathbf{x}(t)$ and $\mathbf{u}(t)$ are normalized at all times. In the absence of rotational diffusion (${D}_{r}=0$), the vectors $\mathbf{x}$ and $\mathbf{u}$ rotate over time by an angle ${v}_{0}t$ around the axis $\mathbf{u}\times \mathbf{x}$. Consequently, particle trajectories in the absence of noise trace out great circles in the plane defined by $(\mathbf{u}\times \mathbf{x})$.
Spatial correlation of APBs on a sphere
$C(t)=\u27e8\mathbf{x}(t)\cdot \mathbf{x}(0)\u27e9$ To illustrate how ABPs on a sphere differ from ABPs in Euclidean space, we study first the correlation function , where the angled brackets denote a Gaussian whitenoise average. To this end, we rewrite the ABP dynamics Equation 38a in their equivalent Itô form given by
In the Itô formulation any smooth function $f(\mathbf{x},\mathbf{u})$ obeys $\u27e8f(\mathbf{x},\mathbf{u})\mathrm{d}\xi \u27e9=0$, such that (Winkler et al., 2015)
and
which yields a damped harmonic oscillator equation for the correlation function
Normalization and orthogonality of $\mathbf{x}(t)$ and $\mathbf{u}(t)$ imply the initial conditions $C=1$ and $\mathrm{d}C/\mathrm{d}t=0$ at $t=0$. The behavior of solutions of Equation 40 is a function of the rotational Péclet number ${\mathrm{Pe}}_{r}:={v}_{0}/{D}_{r}$ that quantifies the ratio between active motion and orientational diffusion. For ${\mathrm{Pe}}_{r}\S lt;1$, (‘highnoise regime’), the position correlation function $C(t)=\u27e8\mathbf{x}(t)\cdot \mathbf{x}(0)\u27e9$ decays according to Equation 40 monotonically to zero. For ${\mathrm{Pe}}_{r}\S gt;1$, (‘low noise regime’) position correlations exhibit damped oscillations. To validate our simulation method (described in the following section), analytic predictions for $C(t)$ are in Figure 3B (main text) compared against the ensemble average $\u27e8\mathbf{x}(t)\cdot \mathbf{x}(0)\u27e9$ over $3\times {10}^{4}$ simulated ABPs.
Stochastic simulation of active Brownian particles on the sphere
To ensure a numerically exact normalization of the particle’s position and orientation vectors on the unit sphere, we simulated the dynamics
We numerically solve the Itô formulation of this system using the EulerMayurama scheme (Higham., 2001), and confirm that this system reproduces the correlation dynamics predicted by Equation 40 (Figure 3B, main text).
FokkerPlanck equation
To study the continuum dynamics of a large number of noninteracting ABPs on a sphere, we determine the dynamics of the probability density $p(\mathbf{x},\mathbf{u},t)$ of particle positions $\mathbf{x}$ and orientations $\mathbf{u}$ at time $t$. To do so, it is convenient to express particle positions in terms of a parameterisation $\mathbf{x}(t)=\mathbf{x}[{x}^{1}(t),{x}^{2}(t)]$ that defines tangential basis vectors by ${\mathbf{e}}_{i}=\partial \mathbf{x}/\partial {x}^{i}$ ($i=1,2$) and a metric tensor ${g}_{ij}={\mathbf{e}}_{i}\cdot {\mathbf{e}}_{j}$. By definition, we have $\mathrm{d}\mathbf{x}={\mathbf{e}}_{i}\mathrm{d}{x}^{i}$ and Equation 38a can be rewritten as
General tangential vectors on the surface can be written as $\mathbf{u}={u}^{i}{\mathbf{e}}_{i}$ and on a unit sphere the surface normal can be identified with particle positions $\mathbf{n}={\mathbf{e}}_{1}\times {\mathbf{e}}_{2}/{\mathbf{e}}_{1}\times {\mathbf{e}}_{2}=\mathbf{x}$. Hence, on the unit sphere the GaussWeingarten relation reads ${\partial}_{i}{\mathbf{e}}_{j}={C}_{ij}\mathbf{x}+{\mathrm{\Gamma}}_{ij}^{k}{\mathbf{e}}_{k}$, where ${\mathrm{\Gamma}}_{ij}^{k}$ denote Christoffel symbols and ${C}_{ij}$ is the curvature tensor. This implies together with Equation 42 the geometric relation
Comparing this identity with the stochastic dynamics $\mathrm{d}\mathbf{u}$ in Equation 38b and using that ${C}_{ij}{u}^{i}{u}^{j}={g}_{ij}{u}^{i}{u}^{j}={\mathbf{u}}^{2}=1$ for unit vectors $\mathbf{u}$ on the unit sphere, we find the covariant stochastic differential equation
In Equation 43, ${\u03f5}_{ij}=\mathbf{x}\cdot ({\mathbf{e}}_{i}\times {\mathbf{e}}_{j})$ denotes the LeviCivita tensor on the unit sphere.
In this covariant basis, we define the scalar probability density
where $\delta (x)$ denotes a Dirac function. Combining Equations 42; 43, standard methods (Fily et al., 2016; CastroVillarreal and Sevilla, 2018) allow us to obtain the FokkerPlanck equation for $p(\mathbf{x},\mathbf{u},t)$ as
Using the identity ${\u03f5}_{k}^{i}{\u03f5}_{l}^{j}={g}^{ij}{g}_{kl}{\delta}_{l}^{i}{\delta}_{k}^{j}$, the dynamics of the probability density is finally given by
which agrees with the result in CastroVillarreal and Sevilla, 2018.
Hydrodynamic expansion
To connect the FokkerPlanck dynamics given in Equation 46 to hydrodynamic fields, we define (probability) density and fluxes by $\rho (\mathbf{x},t)=\int {\mathrm{d}}^{2}\mathbf{u}p(\mathbf{x},\mathbf{u},t)$, and ${J}^{i}(\mathbf{x},t)={v}_{0}\int {\mathrm{d}}^{2}\mathbf{u}{u}^{i}p(\mathbf{x},\mathbf{u},t)$. Their dynamics on the unit sphere is given by CastroVillarreal and Sevilla, 2018.
where couplings to higher order fields are neglected, as they vanish at shorter timescales due to the presence of rotational noise. Expressing Equation 47a and Equation 47b in terms of scalar and vector spherical harmonics (see Appendix 2) for an arbitrary sphere radius R_{0} yields the mode dynamics Equation 13a, Equation 13b and Equation 13c the main text.
Appendix 4
Learning and interpreting the linear model
We describe details about the inference procedure used to learn the linear ordinary differential equation (ODE) model considered in the main text. We then discuss how the matrix $M$ found by this procedure can be further studied in terms of its realspace kernel representation and derive this kernel for the ABP dynamics introduced in Appendix 4.
Inference of the dynamical mode coupling matrix M
Given a dynamical mode vector $\mathbf{a}(t)={[{\rho}_{lm}(t),{j}_{lm}^{(1)}(t),{j}_{lm}^{(2)}(t)]}^{\top}$, the goal is to learn a linear minimal model
of the mode dynamics. Here, $M$ is an unknown $n\times n$ mode coupling matrix, where generally $n=3{({l}_{\mathrm{max}}+1)}^{2}2$. In systems with global mass conservation, as considered in this work, one can additionally use that the mode ${\rho}_{00}$ is constant and eliminate the corresponding couplings from $M$.
To describe the algorithm that was used to infer the mode coupling matrix $M$, we parameterize $M$ by a vector $\mathbf{p}$ that contains all nonzero entries and introduce a function $\mathcal{M}$ that represents the underlying matrix structure. Together, they generate the explicit form $M=\mathcal{M}(\mathbf{p})$ of the mode coupling matrix. Imposing structure on the matrix, such as rank constraints, or sparsity leads to a shorter vector $\mathbf{p}$ and modifies the definition of $\mathcal{M}$ accordingly. Denoting $\mathbf{A}(t;\mathcal{M},\mathbf{p},{\mathbf{a}}_{0})$ as the result of numerically integrating the system of ODEs Equation 48 up to time $t$ from initial condition ${\mathbf{a}}_{0}$ with $M=\mathcal{M}(\mathbf{p})$, we define the loss function
where the t_{i} are time points in an interval $[{t}_{I},{t}_{N}]$ at which the data and the ODE solution are sampled. Using the ODE solvers and optimization functions provided by the Julia modules DifferentialEquations.jl and DiffEqFlux.jl (Rackauckas et al., 2021), we can differentiate through the ODE solver to calculate derivatives of the loss function Equation 49 with respect to parameters $\mathbf{p}$ and subsequently apply gradientbased optimization algorithms. The loss function is minimized using the ADAM algorithm (Kingma and Ba, 2017), followed by the BroydenFletcherGoldfarbShannon (BFGS) algorithm (Nocedal and Wright, 2006). To increase the robustness of the optimization and promote sparsity, we use a sequentially thresholded algorithm (Supekar et al., 2021; Brunton et al., 2016; Reinbold, 2020). A complete overview of this procedure is shown in Appendix 4—figure 1 and the details of the specific design decisions made in the algorithm are discussed in the following:
1. To account for the variation in scale between the different modes in the data $\mathbf{a}(t)$, each mode is normalized by its median absolute deviation (MAD) across the full timespan in which the data are available. Specifically, we scale each mode by
where ${\overline{a}}_{i}={\text{median}}_{k}[{a}_{i}({t}_{k})]$ and the median is taken over all timepoints, giving rise to a scaled mode vector $\stackrel{~}{\mathbf{a}}(t)$. Losses analogous to Equation 49 that are computed using scaled data are denoted in the following by $\stackrel{~}{L}$.
2. To prevent overfitting, we divide the data into two regions, a learning region from ${t}_{I}$ to ${t}_{N}$ and a validation region from ${t}_{N}$ to ${t}_{F}$. Only data from the learning region is used in the optimization of the loss function Equation 49. However, the model is integrated into the validation region, and a corresponding validation loss using only the data in the validation region is calculated. During each optimization run, we choose the model with the lowest loss in the validation region, lowering the likelihood of overfitting to the specific data in the learning region.
3. To prevent the optimization from getting stuck in local minima, we incrementally increase the timespan of the data included in the optimization objective (blue box in Appendix 4—figure 1). We increase the time window backward from a fixed endpoint ${t}_{1}={t}_{F}$, choosing each iteration an earlier initial condition at time $t}_{i}\S lt;{t}_{i1$. The advantage of stepping backward rather than forward from a fixed initial condition is twofold: first, the validation region stays unchanged throughout the optimization, making comparisons of the validation loss easy. Second, because the initial condition changes with each run, the learned matrix tends to be more robust to fluctuations in the initial condition.
4. After the optimization step, sparsity is promoted by thresholding the elements in the matrix (Brunton et al., 2016), removing small magnitude elements that do not noticeably contribute to the mode dynamics (purple box in Appendix 4—figure 1). The optimization procedure is then repeated until the thresholding converges. The threshold is chosen to generate a sparse matrix that still reproduces the dynamics faithfully.
5. Once the sparsity pattern is obtained from the sequential thresholding and optimization procedure a final run of the optimization is performed on the unscaled mode data to find the final dynamical matrix $M$, which removes any potential slight bias the MAD scaling might have introduced in the parameter values $\mathbf{p}$.
Finally, the numerical stability of the model can be checked by examining the eigenvalues of the learned matrix. For the ABP test data, we learn a matrix $M$ for which the largest real part of its eigenvalues is at machine precision. For the experimental data, the largest real part in the eigenvalues is $7.4\times {10}^{4}$, which corresponds to a time scale of around 675 mins. While the corresponding dynamics will eventually become unstable, solutions remain bound over a period of approximately 45 hours, which is four times as long as the input data from which the mode coupling matrix was learned.
Learning and validation regions used in this work
For the ABP data, the first 15 frames are excluded, so that – consistent with coarsegraining assumptions (see Appendix 3, Equation 47a, Equation 47b) any remnants of higher orientational order introduced by the initial conditions have decayed. The subsequent 140 frames are used as the learning region, followed by a validation region of 20 frames. Each frame corresponds to a time interval of approximately 0.06 in units of ${R}_{0}/{v}_{0}=1$. We exclude the first and last 10 frames of the experimental zebrafish data and split the remaining data into a learning region of 360 frames, with the remaining 40 frames used for validation. Each frame corresponds to a time interval of 2 min.
Green’s function representation of the learned matrix
The learned matrix $M$ consists of 9 blocks each with $[{({l}_{\text{max}}+1)}^{2}1]\times [{({l}_{\text{max}}+1)}^{2}1]$ entries. Each block relates a mode family to time derivatives of another and we write
We denote the components of each block by ${\left({M}^{{m}_{1}{m}_{2}}\right)}_{lm,{l}^{\prime}{m}^{\prime}}\equiv {M}_{\alpha \beta}^{{m}_{1}{m}_{2}}$, where ${m}_{1},{m}_{2}\in \{\rho ,1,2\}$, and $\alpha $, $\beta $ are multiindices that represent the harmonic modes $(lm)$. Using the mode representation Equation 5 and the form of the linear minimal model Equation 48, we find
Using Equation 30a, Equation 51 can be cast into the dynamic kernel Equation 14 given in the main text, where we defined the vector kernel
and the matrix kernel
where ⊗ denotes a dyadic product. The matrix ${M}^{J}(\mathbf{r},{\mathbf{r}}^{\prime})$ has a 0 eigenvalue with right eigenvector ${\widehat{\mathbf{r}}}^{\prime}$ and left eigenvector $\widehat{\mathbf{r}}$, which implies $det\left({M}^{J}\right)=0$. Numerical analysis of the matrix invariants shows that a second eigenvalue is 0 (Appendix 4—figure 2), leaving only a single nonzero eigenvalue that can be conveniently found from $\text{tr}\left[{M}^{J}(\mathbf{r},{\mathbf{r}}^{\prime})\right]$ and is shown in Figure 4D (main text).
Realspace kernels of active Brownian particle dynamics
In the following we determine a realspace kernel representation in the form Equation 14 for the flux dynamics of ABPs given in Equation 47b. We can read off the kernel coefficients in Equation 52 and in Equation 53 from the coarsegrained ABP dynamics in mode space, given in Equation 13b and Equation 13c. For the kernel ${\mathbf{m}}^{\rho}(\mathbf{r},{\mathbf{r}}^{\prime})$, we have ${M}_{\alpha \beta}^{1\rho}=\frac{{v}_{0}^{2}}{2}{\delta}_{\alpha \beta}$ and ${M}_{\alpha \beta}^{2\rho}=0$ ($\alpha ,\beta =(lm)$), such that Equation 52 becomes
Here, we have used in the first step the definition of ${\mathbf{\Psi}}_{lm}(\hat{\mathbf{r}})$ given in Equation 29a and in the second step the completeness of the spherical harmonic basis functions ${Y}_{lm}(\widehat{\mathbf{r}})$, where $\delta (\mathbf{r}{\mathbf{r}}^{\prime})=\delta (\varphi {\varphi}^{\prime})\delta (\mathrm{cos}\theta \mathrm{cos}{\theta}^{\prime})$ denotes the delta function on a sphere. Note that a unit sphere was considered throughout this analysis, such that $\mathbf{r}=\widehat{\mathbf{r}}$. Similarly, Equation 13b and Equation 13c imply for the kernel coefficients in Equation 53 that ${M}_{\alpha \beta}^{11}={M}_{\alpha \beta}^{22}={D}_{r}{\delta}_{\alpha \beta}$ and ${M}_{\alpha \beta}^{12}={M}_{\alpha \beta}^{21}=0$. Consequently, we have
where ${P}_{\parallel}=I\mathbf{r}\otimes \mathbf{r}$ is the tangential projector on the unit sphere. The hydrodynamic flux Equation 47b of ABPs on a sphere can therefore be written in the equivalent integral kernel form
To make analytic kernel properties comparable to practical inference scenarios in which we work with a finite number of harmonic modes, we computed the sums in Equations 54; 55 up to a maximum mode number ${l}_{\text{max}}=4$. The resulting kernels – depicted in Figure 4D (main text) – approximate the Dirac delta function $\delta (\mathbf{r}{\mathbf{r}}^{\prime})$ and its derivative, leading to the finite range of $\text{tr}({M}^{J})$ with amplitude maximum at $\omega =0$, while ${\mathbf{m}}^{\rho}$ vanishes at and peaks away from $\omega =0$. Additionally, finite mode representations introduce an apparent kernel inhomogeneity across the spherical surface as evident from the nonzero standard deviation depicted in Figure 4D of the main text (blue shades).
Data availability
Raw data used in this study can be obtained at https://doi.org/10.1038/s41467019136250 and https://imbdev.gitlab.io/cellflownavigator/.
References

Universal scaling of active nematic turbulenceNature Physics 16:682–688.https://doi.org/10.1038/s4156702008544

BookMathematical Methods for Physicists: A Comprehensive GuideAmsterdam, Netherlands: Elsevier Science.https://doi.org/10.1016/C20090306297

Vector spherical harmonics and their application to magnetostaticsEuropean Journal of Physics 6:287–294.https://doi.org/10.1088/01430807/6/4/014

Anomalous refraction of optical spacetime wave packetsNature Photonics 14:416–421.https://doi.org/10.1038/s4156602006456

New class of turbulence in active fluidsPNAS 112:15048–15053.https://doi.org/10.1073/pnas.1509304112

Itô versus Stratonovich calculus in random population growthMathematical Biosciences 206:81–107.https://doi.org/10.1016/j.mbs.2004.09.002

Dedalus: A flexible framework for numerical simulations with spectral methodsPhysical Review Research 2:023068.https://doi.org/10.1103/PhysRevResearch.2.023068

Active motion on curved surfacesPhysical Review. E 97:052605.https://doi.org/10.1103/PhysRevE.97.052605

Conference9th IEEE International Symposium on Biomedical ImagingAn algebraic solution to rotation recovery in HARDI from correspondences of orientation distribution functions. pp. 38–41.https://doi.org/10.1109/ISBI.2012.6235478

Programmed and selforganized flow of information during morphogenesisNature Reviews. Molecular Cell Biology 22:245–265.https://doi.org/10.1038/s41580020003186

DefectMediated Morphologies in Growing Cell ColoniesPhysical Review Letters 117:048102.https://doi.org/10.1103/PhysRevLett.117.048102

Fluid dynamics of bacterial turbulencePhysical Review Letters 110:228102.https://doi.org/10.1103/PhysRevLett.110.228102

The ultraspherical spectral element methodJournal of Computational Physics 436:110087.https://doi.org/10.1016/j.jcp.2020.110087

Life is Physics: Evolution as a Collective Phenomenon Far From EquilibriumAnnual Review of Condensed Matter Physics 2:375–399.https://doi.org/10.1146/annurevconmatphys062910140509

Modal Analysis of Turbulent Flow near an Inclined Bank–Longitudinal Structure JunctionJournal of Hydraulic Engineering 147:04020100.https://doi.org/10.1061/(ASCE)HY.19437900.0001856

Comparison of quantum and semiclassical radiation theories with application to the beam maserProceedings of the IEEE 51:89–109.https://doi.org/10.1109/PROC.1963.1664

ConferenceIEEE Transactions on Systems, Man, and Cybernetics, Part CParetoBased Multiobjective Machine Learning: An Overview and Case Studies. pp. 397–415.https://doi.org/10.1109/TSMCC.2008.919172

Can One Hear the Shape of a Drum?The American Mathematical Monthly 73:1.https://doi.org/10.2307/2313748

The geometry of soft materials: a primerReviews of Modern Physics 74:953–971.https://doi.org/10.1103/RevModPhys.74.953

The local structure of turbulence in incompressible viscous fluid for very large Reynolds numbersC R (Dokl) Acad Sci URSS 30:301–305.https://doi.org/10.1098/rspa.1991.0075

Twodimensional turbulenceReports on Progress in Physics 43:547–619.https://doi.org/10.1088/00344885/43/5/001

Tensile forces govern germlayer organization in zebrafishNature Cell Biology 10:429–436.https://doi.org/10.1038/ncb1705

ConferenceIEEE International Conference on Shape Modeling and Applications 2006LaplaceBeltrami Eigenfunctions Towards an Algorithm That “Understands” Geometry. pp. 13–19.https://doi.org/10.1109/SMI.2006.21

Hydrodynamics of soft active matterReviews of Modern Physics 85:1143–1189.https://doi.org/10.1103/RevModPhys.85.1143

Bacteria solve the problem of crowding by moving slowlyNature Physics 17:205–210.https://doi.org/10.1038/s41567020010706

Minimal Model of Cellular Symmetry BreakingPhysical Review Letters 123:188101.https://doi.org/10.1103/PhysRevLett.123.188101

BookNumerical OptimizationBerlin, Germany: Springer Science & Business Media.https://doi.org/10.1007/9780387400655

BookThe scales of turbulent motionIn: Pope SB, editors. Turbulent Flows. Cambridge: Cambridge University Press. pp. 182–263.https://doi.org/10.1017/CBO9780511840531

Using noisy or incomplete data to discover models of spatiotemporal dynamicsPhysical Review. E 101:010203(R).https://doi.org/10.1103/PhysRevE.101.010203

Conference2018 IEEE High Performance Extreme Computing ConferenceInteractive Supercomputing on 40,000 Cores for Machine Learning and Data Analysis. pp. 1–6.

Zebrafish gastrulation: cell movements, signals, and mechanismsInternational Review of Cytology 261:159–192.https://doi.org/10.1016/S00747696(07)610043

Active Brownian particlesThe European Physical Journal Special Topics 202:1–162.https://doi.org/10.1140/epjst/e201201529y

Softwaresurfvecdefects, version swh:1:rev:6dc742c376b0d085e19ece65f932ac6935342abaSoftware Heritage.

Adaptive lightsheet microscopy for longterm, highresolution imaging in living organismsNature Biotechnology 34:1267–1278.https://doi.org/10.1038/nbt.3708

Tensor spherical harmonics on S2 and S3 as eigenvalue problemsJournal of Mathematical Physics 19:2441–2446.https://doi.org/10.1063/1.523649

Dynamic mode decomposition of numerical and experimental dataJournal of Fluid Mechanics 656:5–28.https://doi.org/10.1017/S0022112010001217

Quantification of nematic cell polarity in threedimensional tissuesPLOS Computational Biology 16:e1008412.https://doi.org/10.1371/journal.pcbi.1008412

Rotation invariant features for HARDIInformation Processing in Medical Imaging 23:705–717.https://doi.org/10.1007/9783642388682_59

Active swarms on a spherePhysical Review. E, Statistical, Nonlinear, and Soft Matter Physics 91:022306.https://doi.org/10.1103/PhysRevE.91.022306

Simplified LCAO Method for the Periodic Potential ProblemPhysical Review 94:1498–1524.https://doi.org/10.1103/PhysRev.94.1498

Conserved patterns of cell movements during vertebrate gastrulationCurrent Biology 15:R213–R228.https://doi.org/10.1016/j.cub.2005.03.016

BookRandom Vibration of Mechanical and Structural SystemsNASA STI/Recon Technical Report A.

BookFourier Analysis: An IntroductionNew Jersey, United States: Princeton University Press.

On dynamic mode decomposition: Theory and applicationsJournal of Computational Dynamics 1:391–421.https://doi.org/10.3934/jcd.2014.1.391

BookEmergence of puffs, weak and strong slugs from a stochastic predatorprey model for transitional turbulence with streamwise shear interactionsIn: Shih HY, editors. Bulletin of the American Physical Society. New York, United States: Simons Foundation. pp. 65–69.
Decision letter

Raymond E GoldsteinReviewing Editor; University of Cambridge, United Kingdom

Naama BarkaiSenior Editor; Weizmann Institute of Science, Israel

Sebastian FürthauerReviewer; Flatiron Institute, United States
In the interests of transparency, eLife publishes the most substantive revision requests and the accompanying author responses.
Decision letter after peer review:
Thank you for submitting your article "Learning developmental mode dynamics from singlecell trajectories" for consideration by eLife. Your article has been reviewed by 2 peer reviewers, and the evaluation has been overseen by a Reviewing Editor and Naama Barkai as the Senior Editor. The following individual involved in review of your submission has agreed to reveal their identity: Sebastian Fürthauer (Reviewer #2).
The reviewers have discussed their reviews with one another, and the Reviewing Editor has drafted this to help you prepare a revised submission.
Essential revisions:
1) How good is the mode decomposition?
– It is not a priori clear that the mode decomposition suggested is a good representation of the data. It would be helpful if the authors quantified and discussed the error between the continuous low dimensional representation of the data and the raw data.
– Along the same lines, are there a conditions under which mode decomposition is expected to fail? The authors should give a clear idea on when the ideas presented here are applicable and when not.
2) How trustworthy are the inferred summary statistics?
The summary statistic ( defect numbers / fluctuations ) need to be tested against some ground truth. If this is not possible from the data used here it would be illuminating if the authors could test against artificial data.
3) How predictive is the inferred forward model?
– In Figure 3 the authors compare the inferred dynamics and the real data. This is very appealing and the strong point of the paper. We are, however, confused whether the matrix M is inferred using the full time course of the experiment, or just parts of the time course. In other words is M just a way of rewriting the coefficients in Equation 7?
– As it stands we are unclear of what is being learned in M since there is no clear separation between a 'training' and a 'test' set of samples.
4) The paper could greatly benefit if the authors used their method to infer the known dynamics of some artificial data. Then they could validate the inferred model against time courses generated from different initial conditions and against unseen data.
5) We understand that the integral of the field ρ on the surface is 1? We would suggest to mention this already in Equation 1. We also assume that with this definition, J is not actually a cell number flux, but a normalised cell number flux?
6) Because the effective Equations 2a and 2b are written for normalised densities, we understand that cell division and death are effectively absorbed in a flux. So if cells were undergoing cell division in some region of space, and cell apoptosis in another region, this would be interpreted as a net flux between the two regions. Isn't that a potential issue of the method? Can the authors further comment on this?
7) We would argue that one central conclusion of the manuscript is in Figure S6, showing that local rules of interaction cannot explain the dynamics of the system. That is also one potential weakness of the approach: what is exactly learned by the effective dynamics? There is a complex set of interactions between modes suggested by Figure 3D but what is their meaning, given that Figure S6 also indicates, as the authors also point out, that there are missing unresolved players in the equations? The whole approach would be more convincing if the authors could show that the effective learned model have some predictive value. For instance, would changing the initial condition in the integration of the effective system, gives a result that makes sense in regard to the developmental dynamics; indicating that the effective dynamics is not tied to one particular realisation of the biological system?
Reviewer #1:
The authors present ideas for obtaining a low dimensional representations of the complex and rich data obtained in cell tracking experiments on zebra fish gastrulation. They argue that such a low dimensional representation can be used to infer a dynamical model from the data. These ideas are potentially very important. However, the methods presented here need to be more explicitly validated.
Reviewer #2:
Romeo et al. describe a method for the analysis of the collective flow of cells and its decomposition into modes of motion. The authors start by defining a local density and a local flux by coarsegraining the discrete data in a consistent way. The density and flux fields are then projected onto scalar and vectorial spherical harmonics, and the authors show that a relatively small number of modes is sufficient to account for the flow field; the evolution of some of these modes with lower symmetry reflects major changes in the cell density pattern in the embryo. The authors then look for a linear mode coupling, first order dynamics model for the modes of lowest order they have kept in their analysis. This procedure defines an effective model for the dynamics of the system. Looking in real space at the spatial couplings shows that these effective couplings do not decay spatially, but are longranged; indicating that local rules for fluxes cannot capture the dynamics of the system.
The authors presents an elegant and solid work which is done with a very high level of care and rigour. The excellent level of clarity of the manuscript, which presents its concepts in a wellstructured manner, makes it really nice to read. The method described by the authors, notably, gives a way to analyse cellular trajectories in a reduced dimensional space which is helpful to characterise the complex cellular flows.
https://doi.org/10.7554/eLife.68679.sa1Author response
Essential revisions:
1) How good is the mode decomposition?
– It is not a priori clear that the mode decomposition suggested is a good representation of the data. It would be helpful if the authors quantified and discussed the error between the continuous low dimensional representation of the data and the raw data.
We thank the reviewer for raising this important point, which has been addressed in detail in the extended Appendix 2 (see also Appendix 2 – Figure 2).
First, we note that it is generally not straightforward to define a meaningful error between the raw particle data, which are restricted to a discrete point set, and the coarsegrained hydrodynamic fields, which are defined on a continuous domain (here given by the surface of the sphere). A theoretically meaningful and practically feasible approach is to quantify the information loss associated with different coarsegraining scales. In the new section of Appendix 2, “Information loss through coarsegraining”, we present measurements of the spectral entropy, which quantifies the variation in information content as the coarsegraining scale is varied. As the coarsegraining scale (kernel width) is increased, one finds that finer structures in the data become rapidly blurred before reaching a slower decay of information (Appendix 2 – Figure 2). Guided by this analysis, we selected a coarsegraining scale (yellow) in the domain where the information loss begins to decay more slowly. This choice ensures that key features of the cell density and fluxes remain wellresolved while also ensuring robustness against variations of the coarsegrain scale.
– Along the same lines, are there a conditions under which mode decomposition is expected to fail? The authors should give a clear idea on when the ideas presented here are applicable and when not.
We thank the reviewer for this comment. As is generally the case for hydrodynamic continuum models, mode decomposition succeeds when gradients in continuous fields are sufficiently smooth. This approach fails when fields are dominated by finescale structures or sharp features, such as shock fronts that lead to Gibbs ringing when represented spectrally. In those cases, the mode decomposition has long tails and fails to converge rapidly so that a direct real space representation of data may be more advantageous. In intermediate cases, where convergence is still fast, but the mode representation lacks sparsity, there is no practical benefit in switching from a realspace “pixel basis” to a mode representation. From a computational perspective, a mode representation is advantageous when spectral codes become faster than particlebased or finiteelementbased simulation techniques.
For the simulated and experimental systems considered in our study, the density modes decay rapidly with wavenumber and the raw flux modes derived from coarsegrained numerical derivatives of noisy data rapidly reach a noiselimited plateau at small mode amplitudes. We have added a supplementary figure to show this more explicitly (Figure 1—figure supplement 1). Therefore, our approach yields a satisfactory reconstruction of the deterministic part of the overall signal in terms of longwavelength and lowfrequency modes that facilitate the inference of hydrodynamic models and are clearly distinguishable from noise. To clarify this further, the revised manuscript contains a newly added Appendix 2 – Figure 1 that complements Video S2 and illustrates that a small number of modes is sufficient to represent the coarsegrained data well. Finally, we note that Appendix 2 provides a detailed description of the thresholdselection procedure (see also Appendix 2 – Figures 3 and 4), where we select thresholding parameters on the Pareto front to optimize the tradeoff between reconstruction accuracy and sparsity of the mode representation.
2) How trustworthy are the inferred summary statistics?
The summary statistic ( defect numbers / fluctuations ) need to be tested against some ground truth. If this is not possible from the data used here it would be illuminating if the authors could test against artificial data.
We thank the reviewer for this suggestion. We have now added two additional analyses to validate the accuracy and robustness of the defect tracking. First, we successfully tested the defect tracking on simple example harmonic vector modes for which defects can easily be identified by eye (see new Figure 2—figure supplement 2). Second, to verify that experimental defect statistics robustly relate to the underlying microscopic data, we performed the coarsegraining and mode representation with two other kernel length scales and found very similar defect dynamics (new Figure 2—figure supplement 3). We also note that the defects tracked in the experimental data obey at all times the nontrivial condition that defect charges add up to 2, confirming that the extracted information is indeed accurate and robust. Because the defect tracking algorithm works reliably on single snapshots of cell fluxes, dynamic measurements are expected to be similarly robust.
3) How predictive is the inferred forward model?
– In Figure 3 the authors compare the inferred dynamics and the real data. This is very appealing and the strong point of the paper. We are, however, confused whether the matrix M is inferred using the full time course of the experiment, or just parts of the time course. In other words is M just a way of rewriting the coefficients in Equation 7?
– As it stands we are unclear of what is being learned in M since there is no clear separation between a 'training' and a 'test' set of samples.
We thank the reviewer for this question and have substantially expanded the description of the inference framework (main text pages 1012 and Appendix 4; see in particular also Appendix 4 – Figure 1 for a detailed flowchart of the learning procedure).
In addition, we now validate the learning framework on synthetic test data generated from the analytically tractable active Brownian particle (ABP) model (pages 1012 and new Figure 3 in the main text; see also new Appendix 3 for additional details and derivations). For this test case, we follow traditional approaches and divide the data into learning, validation, and prediction segments (new Figure 3D,E). We find that our inferred model fits the learning and validation segments well and has predictive power.
For the experimental data, we divide the data into learning and validation segments. In this case, the linear model Equation 12 (Equation 7) in the initially submitted manuscript with the learned mode coupling matrix M can indeed be understood as a compressed representation of the system dynamics. Given the complexity of the experimental data, we require the full range of the data to learn a model that displays good agreement with the training data.
The inferred model can, however, still be used to make predictions: we show in the new Figure 4—figure supplement 1 the output of the model for different initial conditions, such as another experimental zebrafish sample, a rotated initial condition, and an initial condition where 10 percent of the initial density at the animal pole has been numerically “scraped off''. All of these tests show dynamics that reproduce density patterns with the characteristic, elongated “fish shape”, and suggest that indeed the learned model can be used in a predictive context.
4) The paper could greatly benefit if the authors used their method to infer the known dynamics of some artificial data. Then they could validate the inferred model against time courses generated from different initial conditions and against unseen data.
We thank the reviewers for this excellent suggestion. We have now applied our framework to artificial ABP data sets with known ground truth that emulate the key challenges posed by the experimental data. As described in the revised manuscript on pages 1012 and in Appendix 3, the particle data were generated by simulating the stochastic dynamics of ABPs on a sphere for which an analytic coarsegrained linear continuum model can be derived (Equation 47 in Appendix 3 and Equation 13 in the main text of the revised manuscript). We then applied our inference pipeline consisting of coarsegraining, mode projection, and model inference steps to this particle data. We found that the model inference approach learns a close approximation to the theoretical prediction (new Figure 3), even though the output from the particle data is obscured by noise and spatial information is blurred by coarsegraining. Remarkably, even modes that overlap little with the tested initial condition and thus have a low signaltonoise ratio were found to follow the model with good agreement.
5) We understand that the integral of the field ρ on the surface is 1? We would suggest to mention this already in Equation 1. We also assume that with this definition, J is not actually a cell number flux, but a normalised cell number flux?
The reviewers are correct that J represents a normalized cell number flux and we thank them for pointing out that this was not explained clearly enough. We have now added this information directly to the text following Equation 1 in line 124 and amended the revised manuscript in line 134.
6) Because the effective Equations 2a and 2b are written for normalised densities, we understand that cell division and death are effectively absorbed in a flux. So if cells were undergoing cell division in some region of space, and cell apoptosis in another region, this would be interpreted as a net flux between the two regions. Isn't that a potential issue of the method? Can the authors further comment on this?
Cell division and death are indeed an important aspect of developmental processes.
As detailed in the main text subsection “Consequences of cell number variations in experimental data”, we show that under an assumption of a spatiallyinvariant and densityindependent – but possibly timedependent – growth rate, one can absorb effects of varying cell numbers into a normalization of cell density and cell flux. This assumption is motivated by the technical limitation of being unable to distinguish between “true” cell division/death events and imaging or processing artifacts leading to loss of cells, an issue that we discussed with the acknowledged experimental groups that kindly shared and explained their data. As the reviewers correctly point out, inhomogeneous growth would then contribute to effective fluxes. To compensate for those effective inhomogeneous fluxes, we enforce Equation 6 of the main text to relate density and flux modes. This correction only applies to the divergent part of the flux – more details on how we chose spectral projection parameters that optimize reconstruction of fluxderived metrics to heuristically preserve the implicit relationship between divergent and rotational flux components are presented in the section “Optimal compression in space and time” of Appendix 2 (Appendix III B in the previously submitted manuscript).
The overall agreement between the learned model and the experimentally observed dynamics suggests that this procedure provides a good approximation to the exact dynamics of the zebrafish system. We have made limitations arising from currently available data more explicit in the corresponding main text section “Consequences of cell number variations in experimental data”. It is important to emphasize, however, that this is not a fundamental limitation of the mode representation and inference framework. Indeed, as more accurate tracking data is likely to become available in the near future, one can expect that it will soon be possible to distinguish between tracking defects and mitosis or apoptosis (by, say, using suitable mitotic or apoptotic markers). This will make it possible to account for inhomogeneous cell growth in the modeling framework through an extra “growth” field which can be directly included in the local mass balance.
7) We would argue that one central conclusion of the manuscript is in Figure S6, showing that local rules of interaction can not explain the dynamics of the system. That is also one potential weakness of the approach: what is exactly learned by the effective dynamics?
We agree that the learned model in the previous manuscript version while describing the experimental data well, lacked interpretability. To address this point, and guided by the newly added analytical tractable case of ABP dynamics (pages 1012 and new Figure 3 in the main text), we have substantially expanded the inference framework (new Appendix 4 section “Inference of the dynamical mode coupling matrix M” and flow chart Appendix 4 – Figure 1).
The computational improvements, which exploit recent advances in automatic differentiation of ODE solutions, (i) increased the numerical robustness of the inference scheme, (ii) identified a sparser model with higher predictive capability (Figure 4—figure supplement 1), and (iii) enabled through the comparisons with the ABP reference model a direct interpretation of the mode coupling matrix learned for the cellular dynamics during early zebrafish development.
Specifically, the improved model reveals interesting similarities and differences between models of ABPs on the sphere and the dynamics of migrating cells during zebrafish gastrulation. The improved model shows that fluxflux couplings are qualitatively similar for ABPs and zebrafish, while “offdiagonal” densityflux couplings represent longerrange and more inhomogeneous interactions in the experimental data (new Figure 4D, see also new Figure 3C vs. new Figure 4A,B). This suggests that one can interpret the experimentally observed collective cell migration in the zebrafish embryos as ABP dynamics with additional interactions.
These new results suggest an intriguing direction for future research: It will be interesting to study through theory and simulations how ABP models with different types of short, intermediate and longerrange particles interaction give rise to “offdiagonal” modecouplings in effective linear models, and if (and how) such models could be improved by including nonlinear coupling terms in Equation 12, similar to triad interactions in nonlinear pattern formation and turbulence models.
There is a complex set of interactions between modes suggested by Figure 3D but what is their meaning, given that Figure S6 also indicates, as the authors also point out, that there are missing unresolved players in the equations?
We agree with the reviewer that Figure 3D, as presented in the previously submitted manuscript, was not ideal for conveying the structure of the learned model. We have therefore replaced this panel with a detailed comparison of the models learned for ABP and zebrafish data, respectively. We believe that, combined with new ABP results in Figure 3, this change leads to much better interpretability of the learned models. As stated in the response to the preceding comment/question, an interesting avenue for future theoretical analysis is to understand how different types of interactions give rise to the “offdiagonal” linear modecouplings as well as to nonlinear mode couplings in generalized models.
The whole approach would be more convincing if the authors could show that the effective learned model have some predictive value. For instance, would changing the initial condition in the integration of the effective system, gives a result that makes sense in regard to the developmental dynamics; indicating that the effective dynamics is not tied to one particular realisation of the biological system?
We are grateful for this question, which led to the major theoretical and algorithmic improvements included in the revised manuscript. As shown in the new Figure 4—figure supplement 1, the sparser model identified by the enhanced inference framework can be integrated for different initial conditions.
Interestingly, we found that flux and density patterns similar to those seen in the coarsegrained experimental data robustly emerge under various perturbations of the initial conditions. Test examples included in the new Figure 4—figure supplement 1 correspond to an initial condition from the alternative experimental sample, a rotated initial condition, and a removal of cells near the animal pole of the original sample. It would be very interesting to see how the predictions of this model would compare to corresponding celldisplacement experiments in which a subset of cells is ‘surgically’ removed and/or ectopically transplanted during the early stages of zebrafish gastrulation.
https://doi.org/10.7554/eLife.68679.sa2Article and author information
Author details
Funding
European Molecular Biology Organization (ALTF 5282019)
 Alexander Mietke
Deutsche Forschungsgemeinschaft (431144836)
 Alexander Mietke
James S. McDonnell Foundation (Complex Systems Scholar Award)
 Jörn Dunkel
Alfred P. Sloan Foundation (G202116758)
 Jörn Dunkel
MathWorks
 Nicolas Romeo
 Alasdair Hastewell
Robert E Collins Distinguished Scholarship Fund
 Jörn Dunkel
The funders had no role in study design, data collection and interpretation, or the decision to submit the work for publication.
Acknowledgements
We thank Nico Scherf and Ghopi Shah for providing singlecell tracking data and for giving helpful advice on zebrafish development and we thank Paul Matsudaira for discussions. We thank the MIT SuperCloud (Reuther et al., 2018) for providing us access to HPC resources. This work was supported by a MathWorks Science Fellowship (NR and ADH), a Longterm Fellowship from the European Molecular Biology Organization (ALTF 528–2019, AM), a Postdoctoral Research Fellowship from the Deutsche Forschungsgemeinschaft (Project 431144836, AM), a Complex Systems Scholar Award from the James S McDonnell Foundation (JD), the Robert E Collins Distinguished Scholarship Fund (JD) and the Alfred P Sloan Foundation (G2021–16758, JD).
Senior Editor
 Naama Barkai, Weizmann Institute of Science, Israel
Reviewing Editor
 Raymond E Goldstein, University of Cambridge, United Kingdom
Reviewer
 Sebastian Fürthauer, Flatiron Institute, United States
Publication history
 Received: March 23, 2021
 Accepted: December 24, 2021
 Accepted Manuscript published: December 29, 2021 (version 1)
 Version of Record published: February 24, 2022 (version 2)
Copyright
© 2021, Romeo et al.
This article is distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use and redistribution provided that the original author and source are credited.
Metrics

 1,380
 Page views

 316
 Downloads

 0
 Citations
Article citation count generated by polling the highest count across the following sources: Crossref, PubMed Central, Scopus.
Download links
Downloads (link to download the article as PDF)
Open citations (links to open the citations from this article in various online reference manager services)
Cite this article (links to download the citations from this article in formats compatible with various reference manager tools)
Further reading

 Physics of Living Systems
Schooling in fish is linked to a number of factors such as increased foraging success, predator avoidance, and social interactions. In addition, a prevailing hypothesis is that swimming in groups provides energetic benefits through hydrodynamic interactions. Thrust wakes are frequently occurring flow structures in fish schools as they are shed behind swimming fish. Despite increased flow speeds in these wakes, recent modeling work has suggested that swimming directly inline behind an individual may lead to increased efficiency. However, only limited data are available on live fish interacting with thrust wakes. Here we designed a controlled experiment in which brook trout, Salvelinus fontinalis, interact with thrust wakes generated by a robotic mechanism that produces a fishlike wake. We show that trout swim in thrust wakes, reduce their tailbeat frequencies, and synchronize with the robotic flapping mechanism. Our flow and pressure field analysis revealed that the trout are interacting with oncoming vortices and that they exhibit reduced pressure drag at the head compared to swimming in isolation. Together, these experiments suggest that trout swim energetically more efficiently in thrust wakes and support the hypothesis that swimming in the wake of one another is an advantageous strategy to save energy in a school.

 Physics of Living Systems
When a fish beats its tail, it produces vortices in the water that other fish could take advantage of to save energy while swimming.