# Spectral decomposition unlocks ascidian morphogenesis

Joel Dokmegang, Emmanuel Faure, Patrick Lemaire, Ed Munro, Madhav Mani

<figure><img src="image_1.png" /></figure>

## To cite this version:

Joel Dokmegang, Emmanuel Faure, Patrick Lemaire, Ed Munro, Madhav Mani. Spectral decomposition unlocks ascidian morphogenesis. 2023. hal-04310643

HAL Id: hal-04310643

https://cnrs.hal.science/hal-04310643v1

Preprint submitted on 27 Nov 2023

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
# Spectral decomposition unlocks ascidian morphogenesis

Joel Dokmegang<fnref n="1" />,<fnref n="2" />*, Emmanuel Faure<fnref n="4" />,<fnref n="6" />, Patrick Lemaire<fnref n="4" />,<fnref n="5" />, Ed Munro<fnref n="3" />,<fnref n="4" />, and Madhav Mani<fnref n="1" />,<fnref n="2" />*

## ABSTRACT

Describing morphogenesis generally consists in aggregating the multiple high resolution spatiotemporal processes involved into repeatable low dimensional morphological processes consistent across individuals of the same species or group. In order to achieve this goal, biologists often have to submit movies issued from live imaging of developing embryos either to a qualitative analysis or to basic statistical analysis. These approaches, however, present noticeable drawbacks, as they can be time consuming, hence unfit for scale, and often lack standardisation and a firm foundation. In this work, we leverage the power of a continuum mechanics approach and flexibility of spectral decompositions to propose a standardised framework for automatic detection and timing of morphological processes. First, we quantify whole-embryo scale shape changes in developing ascidian embryos by statistically estimating the strain-rate tensor field of its time-evolving surface without the requirement of cellular segmentation and tracking. We then apply to this data spectral decomposition in space using spherical harmonics and and in time using wavelets transforms. These transformations result in the identification of the principal dynamical modes of ascidian embryogenesis and the automatic unveiling of its blueprint in the form of spectograms that tell the story of development in ascidian embryos.

## Introduction

Morphogenesis, the emergence of shape in living systems, is a continuous process littered with spatiotemporal dynamics at various timescales and lengthscales and significance. Developmental biology consists in the identification, localisation, and timing of these processes. Once this work is carried out in a given species, embryogenesis can then be described as a series of stages delineated in space and time by the identified landmarks (1; 2; 3; 4; 5; 6; 7; 8; 9). In order to rigorously define development landmarks, biologists have mostly had to submit imaged embryos either to qualitative analyses, or to rudimentary statistical analysis. These methods however present major drawbacks. On the one hand, they can be time consuming, hence unfit for scale. On the other hand, since morphogenetic processes tend to be unique to a species, these simple methods often lack a general language and framework that permit comparative analyses. For instance, whereas the analysis of cell counts can inform about the proliferation dynamics in a tissue, it does not reveal anything about the shape of the system. For this purpose, other measurements such as length, width, height, aspect ratios or curvatures would be more suitable. Although efforts have been made to automate the staging of development in living systems (10), these methods still rely on preliminary examination using traditional methods.

A standardized method able to identify key milestones in development and lay out the blueprint of morphogenesis in a given system is henceforth needed. Recent breakthroughs in microscopy technology have propelled the resolution of live imaging data to the sub-cellular scale, allowing for the uncovering of precise cell and tissue shape dynamics (11; 12; 13). These advances have created an unprecedented opportunity for the leveraging of computational methods in the study of morphogenesis (9; 11; 14; 15; 16). The rigorous and physically motivated framework of continuum mechanics accommodates itself well to the flow-like dynamics of biological tissues (17; 18; 19; 20; 21). Within this framework, strain-rate fields, which measure the rate at which the shape of a system changes with time, are suited to characterise the dynamical behaviour of the system. Moreover, mounting evidence have informed of the requirement for embryo-wide approaches in the study of morphogenetic flows (21; 22). However, although the evaluation of such global fields across the spatial and temporal domains spanned by a system of interest may reveal valuable insights into its dynamical workings (23; 24; 25), their sole determination might not be sufficient for a holistic description of the behaviour of the system.

<footnote n="1">¹Northwestern University</footnote>

<footnote n="2">²NSF-Simons Center for Quantitative Biology</footnote>

<footnote n="3">³University of Chicago</footnote>

<footnote n="4">⁴University of Montpellier</footnote>

<footnote n="5">⁵CNRS</footnote>

<footnote n="6">⁶LIRMM</footnote>
This is especially true when it comes to morphogenesis (26; 27; 28). The processes involved in development are inherently multiscale, both in the spatial and temporal domains, and may interact or overlap (29; 30; 31). As it is the case with several species (29), ascidian gastrulation is a playground featuring important displays of cellular divisions and tissue mechanics, most notably endoderm invagination (5). The local behaviors captured by indicators such as the strain-rate field might therefore arise from a non-trivial superposition of these dynamical modes, essentially making these measurements complex to interpret without further analysis. Spectral decomposition, whereby a signal is broken down into its canonical components, is well suited to the study of systems that exhibit multimodal behaviors (26; 27). The benefits are at least two-folds: (i) individual constituents may represent distinct dynamical processes, thereby enabling the decoupling of physical processes entangled in the data; (ii) only a handful of components may significantly contribute to the original function, resulting in a compressed, lower dimensional representation that capture the main features of the studied process. The canonical components usually take the form of well known families of functions whose linear combination can reconstitute the original field.

In this work, we take advantage of segmented cellular data to develop a generic computational framework able to identify and delineate the main features of morphogenesis. Our method takes as input precise segmentations of cellular geometries, although only the embryonic surface would suffice, and outputs spatiotemporal spectograms of ascidian development in the form of heatmaps that highlight key developmental processes and stages of ascidian gastrulation. By virtue of a novel meshing scheme derived from level-set methods, raw cell geometry data is first transformed into a single time-evolving embryonic surface on which the strain-rate tensor field can be computed. The accuracy of our inference of a strain-rate field relies on high-frequency temporal sampling, characterised by small deformations of the embryonic surface between subsequent time points. The morphomaps we present are a result of spectral analyses of the strain-rate fields, featuring spherical harmonics decomposition in the spatial domain and wavelet decomposition in the temporal domain. In summary, our method can identify and classify dynamical morphogenetic events. In particular, we are able to identify and distinguish the morphogenetic modes of gastrulation and neurulation phases in ascidian development, recover the characteristic two step sequence of endoderm invagination (6), and capture patterns of cellular divisions in ascidian development (2). Moreover, our method identifies a distinctive stage of ascidian gastrulation, ‘blastophore closure’, which follows endoderm invagination and precedes neurulation.

# Results

## Definition of Lagrangian markers on the surface of the embryo

In order to recover a continuum description of the dynamics in ascidian morphogenesis, we aim to examine the time evolution of strain-rate fields across the entire surface of developing embryos. This endeavour however presents at least two significant challenges. On the one hand, strain-rate computation requires the presence of fiducial markers on the surface of the embryo. Characteristically, this requirement is not always accounted for in the imaging of developing embryos. On the other hand, the outer layer of the embryo being constituted of single cell apical faces, even if such markers had been defined at an earlier time point, uncontrolled stochastic biological processes such as proliferation within the tissue might subsequently grossly uneven the distribution of these markers, thus rendering the computed mechanical indicators at best imprecise. Given the non-triviality of an experimental setup able to solve the described issues, a computational method is required. The goal of such a method would be to computationally discretise an embryo surface into a set of material particles whose trajectories can be tracked in small lapses of development. The positions of these markers over time can then be used to derive mechanical indicators of development dynamics (17).

To achieve this fit, we first take advantage of the level set scheme described in (32) to define static markers on the surface of the embryo at every timepoint of development. The gist of our method resides in the definition of a homeomorphic map between the surface of the embryo ($S_1(t)$) and a topologically equivalent mesh ($S_2(t)$) whose number of vertices (here, $N = 10242$), faces and edges remain fixed (fig. 2a). Conformal parametrisations of embryonic shape have been used in other systems (33). Using this map, the mesh can be deformed to match the surface of the embryo at each time point of development (fig. 2b). As in (32), the map is obtained by finding the positions of $S_2(t)$ vertices that minimize the distance between both surfaces (fig. 2a, right). At the initial time point, $S_2(t)$ is chosen to be a sphere enclosing the embryo. For various reasons, including computational efficiency, variants of this method can be defined such that, for instance, at subsequent time steps, linear combinations of the sphere and its deformations matching the embryo at preceding steps are used. Further details are given in the References section.

Next, via a numerical study, we assert that markers defined as such behave as Lagrangian particles in small increments of developmental dynamics. To support this point, we identify the positions of cellular junctions at the surface of the embryo and evaluate how well our virtual markers mimic their movements in small increments of development. We measure the relative displacement between a cellular junction (fig. 1c red points) and its corresponding vertex (fig. 1c green points) at consecutive time points. Despite gross approximations inherent to the nature of the dataset (geometric meshes) and the process of identifying cellular junctions (averaging the barycenters of closest triangles between three or more cells in contact), the displacement between cellular junctions and their markers remains on average relatively small (under 23% for 96%, fig. 1d).
This result can be seen in more favorable light on the method when considering that cellular junctions, the meeting point of three or more cells, are expected to exhibit more chaotic behaviour than material particles belonging to just one cell. Moreover, these numbers are skewed by large scale morphogenesis processes such as synchronised cellular divisions, as evidenced by spikes in fig. 1d, and fast-paced endoderm gastrulation, as highlighted by higher displacements at the vegetal pole of the embryo during this phase (red dots in fig. 1c, right).

## Strain rate Field describes ascidian morphogenesis

Once a mesh representing the surface has been constructed for the embryo surface, we proceed with the computation of the strain-rate fields across the surface of the embryo and throughout development timeline. Thanks to the Lagrangian nature of mesh vertices, a velocity field can be defined on the mesh. Although particles at every given time point live on the 2D surface of the embryo, their trajectories in time involve greater degrees of freedom in the 3D space. A correct parametrisation of the velocity field at every position henceforth requires three coordinates $v(\mathbf{x}) = (v_x(\mathbf{x}), v_y(\mathbf{x}), v_z(\mathbf{x}))^T$. The strain rate field is derived as the symmetric part of the gradient of the velocity field. Intuitively, the strain rate evaluated on a given mesh vertex measures how the velocity vector varies in the neighbourhood this point (34; 35).

$$D(t, \mathbf{x}) = \mathrm{sym}(\nabla v(t, \mathbf{x})) \qquad (1)$$

The mathematical construction of the strain rate eq. 1 implies that its algebraic representation takes the form of a second order tensor which can be written down as a $3 \times 3$ matrix $(D)_{ij}$. The diagonal elements of this matrix capture the linear strain rate in the $x, y, z$ axes, depicting the change in length per unit length. The non-diagonal elements stand for shear strain rates in the $xy, xz$ and $yz$ directions. Because $D$ is symmetric, there exists an alternative representation which holds stronger local geometric meaning. This representation is obtained by computing the eigenvectors and eigenvalues of the strain rate tensor. Eigenvectors in this sense stand for orthogonal spatial directions that are not rotated, but only stretched, by the application of the strain rate tensor. They define the principal axes of a coordinate system in which the strain rate tensor would be solely composed of maximal linear strain rates (fig. 2a). From this decomposition, we derive a scalar field that is computed at every mesh particle as the square root of the sum of the squared eigenvalues of the strain-rate (fig. 2b). Intuitively, this field describes the magnitude of the rate of change underwent by a particle at the surface of the embryo in the three orthogonal spatial directions of most significant rate of change.

In order to minimize undesirable artifacts that may arise from numerical inefficiencies, we apply a Gaussian filter to the strain rate tensor field before deriving the scalar field. At each particle location, we apply convolution with a Gaussian mask spanning the 2-ring neighbourhood of the particle. A similar smoothing process is also used in the time domain. Interestingly, this strain-rate derived scalar field remarkably mirrors well-known features of ascidian development. A clear correlation between the spatiotemporal distribution of morphogenesis processes described in the literature and heatmaps of this field on the evolving embryo surface emerges. On the one hand, wider spatial bands of yellow and red depicting higher morphological activity portray the spatiotemporal locations of endoderm invagination in the embryonic vegetal pole (fig. 2b, center-left), synchronised rounds of division in the animal pole, and zippering in the notochord (fig. 2b, center-right, right). On the other hand, known spatiotemporal locations of low morphological significance in the embryo exhibit stronger concentration of mechanical activity on cell boundaries, with the corollary that cellular identities are mostly preserved (fig. 2b, center-left). A notable by-product of this scalar field is the evidencing of the duality of the embryo as both a sum of parts constituted of cells and an emerging entity in itself: the strain rate field clearly discriminates between spatiotemporal locations where isolated single cell behaviours are preponderant and those where coordinated cell behaviours dominate.

This brief overview already demonstrates the riches in a quantitative, spatially global and not event-driven approach to study morphogenesis. It also sets the stage for further analysis of morphogenesis dynamics in the ascidian embryo.

## Spectral decomposition in space: Spherical harmonics reveal the main modes of ascidian morphogenesis

In order to capture relevant features of the strain rate field in the spatial domain, we conduct a spectral analysis of the scalar strain rate field. The family of spherical harmonic functions stands out as a de-facto standard for the study of signals defined on a unit sphere, and by extension on surfaces homeomorphic to the sphere. Spherical harmonics form an infinite orthonormal basis of functions defined on the surface of the sphere and represent a generalisation of the Fourier series for functions of two variables (36). Unsurprisingly, these functions play an important role in many branches of science including quantum mechanics and geophysics (36; 37). Spherical harmonics have recently been used in studies of morphogenesis in zebrafish and mouse (26; 27).

Spherical harmonic basis functions are indexed by two parameters $(l, m)$, such that $l \ge 0$, $|m| \le l$ representing respectively the degree and order of the harmonic. A signal defined on the sphere can be written as a linear combination of these functions. Decomposing a function into spherical harmonics hence amounts to finding the coefficients $f_{lm}$ of this weighted sum. In the
<figure><img src="image_3.png" /></figure>

<figure><img src="image_4.png" /></figure>

### C. Method validation: Tracking cellular junctions

<figure><img src="image_5.png" /><figcaption>**Figure 1.** Level-sets inspired Lagragian markers. a) Left Schematics of the level set method. Right Fundamentals of the numerical scheme that shapes $S_2$ into $S_1$. b) Illustration of the method in action. Left Raw data consisting of geometric meshes of single cells spatially organised into the embryo. Center Embryo surface mesh resulting from the application of the level set scheme. Right Rendering of the embryonic surface. c) Tracking of cellular junctions. Left Identification of cellular junctions (red dots). Center Corresponding markers (green dots, defined as vertices on the computed embryonic surface closest to the junctions. Right Relative displacement between junctions and their markers at consecutive timepoints. d) Plot over time of the relative displacement between cellular junctions and their markers.</figcaption></figure>
<figure><img src="image_6.png" /><figcaption>**Figure 2.** Strain-rate field describes morphogenesis. The strain-rate tensor field measures the rate at which morphological changes occur in the embryo as a function of time. The strain-rate tensor field is locally represented as a 3x3 symmetric matrix and is completely determined by its eigenvector fields. a) Heatmap of the eigenvector fields of the strain rate tensor. Each row represents a vector field distinguished by a distinct root color (yellow, pink, white). The gradient from the root color to red represents increasing levels of morphological activity. Top Spatiotemporal dynamics of the first eigenvector field. Middle Spatiotemporal dynamics of the second eigenvector field. Bottom Spatiotemporal dynamics of the third eigenvector field. b) Heatmap of the scalar strain rate field. The gradient from yellow to red depicts regions of increasing morphological activity, while black stands for areas of low morphological activity. The heatmaps show high morphological activity in the invaginating endoderm and zippering notochord, but also across the embryonic animal during rounds of synchronized division.</figcaption></figure>

case of our spatiotemporal scalar strain rate field, the coefficients $f_{lm}$ are also a function of time and can be obtained as shown in equation 2.

$$f_{lm}(t) = \oint f(\theta, \phi, t) Y_{lm}^*(\theta, \phi) dA \quad (2)$$

Here, $Y_{lm}^*$ stands for the complex conjugate of the spherical harmonic $Y_{lm}$. Moreover, for a given degree $l$, each of the $(2l+1)$ spherical harmonics $(Y_{lm})_{|m|<l}$ spatially partitions the unit sphere into as many spatial domains indicating when a signal is positive, negative or null. Figure 3a illustrate the projections of the scalar strain rate field to spherical harmonics $(Y_{lm})_{l \le 2, |m|<l}$ at $t = 4.74$ hpf, and their mapping onto the surface of the embryo. These plots reveal for instance that while there is no embryo-wide dominant morphogenesis process at this time ($l = 0, m = 0$), smaller regions, notably the vegetal pole are experiencing significant morphological activity ($l = 1, m = 0$).
The contributions of each spherical harmonic to the global signal can be assessed more rigorously, and interpreted in the light of biology. To this effect, we observe the temporal dynamics of the coefficients $f_{lm}(t)$ associated with each spherical harmonic. In analogy to *Principal Components Analysis*, we measure the average variance ratio over time of each harmonic with respect to the original signal (Fig. 3b). With a variance ratio of 64.4%, the spherical harmonic $Y_{00}$, capturing embryo-wide morphological activity, contributes the most to ascidian morphogenesis. Spherical harmonic $Y_{10}$ is the next contributor, coming second with a variance ratio of 4.1%. This observation is warranted, as $Y_{10}$ maps to the animal and vegetal poles of the embryo, which are the epicenters of synchronised cellular divisions and endoderm invagination respectively (3; 5). Interestingly, variances in the directions of $Y_{00}$ and $Y_{10}$ evolve in an antiphased pattern, most notably in earlier parts of the plot, with $Y_{00}$ contributing maximally (and $Y_{10}$ minimally), during periods of high proliferation, before relinquishing some variance shares to $Y_{10}$, which then peaks. This suggests that while sporadic deformations induced by cellular divisions often dominate the landscape of morphological activity, an observation consistent with studies in other species (38), other localised, slower processes are at play in the embryo. The described pattern tends to fade out a little in the later parts of the plot, suggesting a shift in development dynamics.

Furthermore, by observing the time dynamics of the coefficients themselves (Fig. 3c), one can easily identify which parts of the embryo are concerned by the morphological changes depicted. For instance, the positive peaks in $f_{10}(t)$ (Fig. 3c) indicate that the morphological processes at hand take place in the northern hemisphere of the sphere. Remarkably, these coincide with rapid growth in cell population and thus synchronous cellular divisions, which are known to be restricted to the animal pole of the ascidian embryo in its gastrulation phase (3). In addition, most of the dynamics captured by $f_{10}(t)$ are in the negative spectrum, pointing to the lower hemisphere of the embryo, the foyer of several morphological mechanisms at play in ascidian early development.

The sporadic short-time scale proliferation events in the animal pole coexist with numerous other features of morphogenesis, most notably, the larger scale continuous deformation process in endorderm invagination at the embryonic vegetal pole. Beside the peaks on the plots of the time series, it is not a trivial task to identify what other rich insights may be hidden in this data. A simple observation of the oscillatory patterns of these main modes hence paints an incomplete picture of ascidian morphogenesis. Extracting the footprint of all morphogenesis processes in these time series requires further analysis.

## Spectral decomposition in time: Wavelets analysis of spherical harmonic signals unveils the blueprint of morphogenesis

Wavelets have been put forward as effective multi-resolution tool able to strike the right balance between resolution in time and resolution in frequency (39; 40). Although they have been taken advantage of in the broader context of biology, most notably in the analysis of brain and heart signals (41), they have so far been under-used in studies of biological development. The reason might be found in the reality that morphogenesis data is often not understood in term of time series. Our spherical harmonics decomposition of morphogenesis, inspired by similar endeavours in other fields (36; 37), offers an unprecedented opportunity to leverage the existing rich signal processing toolbox in development biology. In particular, enlisting the help of wavelet transforms in unlocking the complex entanglements of the multiple morphological process at play during ascidian early development. We proceed to apply the Ricker wavelet transform to our spherical harmonics time series, normalised by mean and standard deviation in different time windows of interest. The result is a set of spectrograms which decompose the signals into canonical components organised in timelines that reveal the blueprint ascidian development.

First, we apply the wavelet transform on the timeseries $f_{00}(t)$ to the entire time range covered by our dataset, comprising both gastrulation and neurulation (fig. 4a). Mirroring this timeseries, the high frequency events depicted by yellow blob at the top of the heatmap represent periods of synchronized division across the embryo. The dark band in the middle separating two large regions depicts a short transition phase delimiting two phases of ascidian early development. The timing of these stages as reflected in the scalogram matches the timeline of gastrulation and neurulation. Within both phases, the concentric gradients from red to yellow culminating in dense yellow spots in the center of both regions portray increasing morphodynamics. Moreover, the brighter blob on the right translates the fact that the second phase of gastrulation characterised by lateral shortening happens at a faster pace than the first apical constriction driven phase.

To better understand the specifics of ascidian gastrulation, we restricted the wavelet transform to the gastrulation period ($t \in [4, 6.3]$ hpf). The resulting spectrogram (fig. 4b) shows that ascidian gastrulation unfold itself in two major phases, delineated on the scalogram by the dark region at the center of the heatmap. The timeline of these events, strengthened by an analysis of topological holes in the embryo (supplementary fig. 4a) support the hypothesis that these phases correspond to endoderm invagination followed by the near-closing of the resulting gut, a process initiated by the collective motion of lateral mesoderm cells that we coin *blastophore closure*. Both the timeseries (mostly in the negative spectrum) and the scalogram of $f_{10}(t)$ (fig. 4c) adds another layer of validity to this conclusion: the large yellow blob occupying the majority the plot surface highlights that fact that regions of the embryo covered by spherical harmonic $Y_{10}$, hence endoderm and mesoderm cells, are subject to intense and prolonged morphological processes.
<figure><img src="image_7.png" /><figcaption>b. Dominant modes of morphogenesis</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>c. Timeseries of f₀₀ and f₁₀</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>**Figure 3.** Spherical harmonics decomposition of morphogenesis. a) Spherical harmonics decomposition of the scalar strain rate field at *t* = 4.74 hpf. b) Time evolution of the variance ratios of the main modes of ascidian early morphogenesis (*Y*₀₀ and *Y*₁₀). The cell population dynamic is also included in the plot for clarity. c) Time evolution of the coefficients *f*₀₀ and *f*₁₀ associated with spherical harmonics (*Y*₀₀ and *Y*₁₀). The cell population dynamic is also included in the plot for clarity.</figcaption></figure>

The first of these two phases, namely endoderm invagination, has been thoroughly investigated in literature. Most notably, it was identified that endoderm invagination was driven by two distinct mechanisms of endoderm single cells (5): on the one hand, cells constricted apically by reducing the surface area of the apices, thus flattening the convex vegetal pole of the embryo setting the stage for invagination. This was followed by a collective drastic decrease of their lateral faces, increasing the pace of endoderm invagination. The wavelet transform restricted to the period of endoderm invagination applied to *f*₂₀(*t*), whose corresponding spherical harmonic function *Y*₂₀(*t*) maps more precisely to the endoderm, beautifully captures this two-steps process (fig. 4d). The timing revealed by this scalogram is in accordance with an analysis endoderm cell shape ratios
(supplementary fig. 4b).

<figure><img src="image_10.png" /><figcaption>**Figure 4.** Wavelet analysis highlights multi-timescale modes of morphogenesis. **a)** Scalogram resulting from the ricker wavelet transform applied to $f_{00}(t)$ over the whole period covered by the dataset $t \in [4, 8.6]$ hpf. **b)** Scalogram resulting from the ricker wavelet transform applied to $f_{00}(t)$ restricted to the gastrulation period $t \in [4, 6.3]$ hpf. The high frequency events highlighted here represent time points of synchronized division across the embryo. The dark band in the middle separating two large red regions indicates that there are two phases of invagination characterized by large deformations and a relatively calm transition phase in between. **c)** Scalogram resulting from the ricker wavelet transform applied to $f_{10}(t)$ restricted to the gastrulation period $t \in [4, 6.3]$ hpf. Similar to **b)**, the high frequency events indicate synchronized division in the embryo. **d)** Scalogram resulting from the ricker wavelet transform applied to $f_{20}(t)$ restricted to endoderm invagination $t \in [4, 5.1]$ hpf.</figcaption></figure>
## Spectral decomposition of morphogenesis in mutant embryo

To assess how our framework adapts to different phenotypes, we set out to conduct a spectral decomposition of morphogenesis in a mutant embryo. In this particular mutant, MEK kinase was inhibited, which resulted in a massive re-specification of vegetal cell fates, and a disruption of endoderm invagination. We applied to the mutant raw data (fig 5a *top*) each of the steps in our workflow. First an evolving mesh matching the shape of the embryo at every time point was obtained through the level set scheme. Then, a strain rate tensor field was computed over the surface of the embryo throughout development time (fig 5a *bottom*). A spatiotemporal spectral analysis was subsequently conducted using spherical harmonics on the mutant surface and wavelet analysis of the timeseries of the coefficients of the main hamonic modes. In order to meaningfully contrast the dynamics of the mutant development against those of the wild-type embryo, the analysis was restricted to development beyond the 64-cell stage.

Similar to the wildtype (WT) embryo, the main harmonic modes in the mutant development were $Y_{00}$ and $Y_{10}$ with respective variance ratios 73.68% and 1.65%, making the timeseries $f_{00}(t)$ and $f_{10}(t)$ the main focus of our examination. The temporal dynamics of these coefficients already serve massive discrepancies between the two strains of ascidians (fig 5b). On the one hand, the drop in the share of $Y_{10}$ is telling of the lower order of morphological activity at the hemisphere. On the other hand, the peaks and lows of $f_{10}(t)$ which coincide with growth in population size, are not restricted either to the negative or positive domain of the curve. This implies that, contrary to the WT, synchronous cellular divisions are not restricted to one hemisphere of the embryo. The fact the combined effects of native morphogenesis and synchronous cellular divisions at the mutant vegetal pole account for only 1.65% of morphological dynamics, a significant drop when compared to WT embryos which are deprived of synchronous cellular divisions at the vegetal pole, is an evidence of marginal morphological activity in the endoderm.

The wavelet transform applied to timeseries $f_{00}(t)$ and $f_{10}(t)$ yield spectograms which bring further insights into the mutant morphogenesis (fig 5c). On the one hand, the heatmap of $f_{00}(t)$ hints to two distinct phases of the mutant development during the observed time span (fig 5c *top*). The timing of the second phase seems to coincide with the advent of the most prolific cellular division in the embryo. On the other hand, plots of the scalogram of $f_{10}(t)$ appears to concur with the hypothesis of very low morphological activity in the vegetal pole. The constant red on this heatmap reflects unchanging levels of morphological activity at the vegetal pole of the embryo and confirms the absence of drastic cellular movements that usually drive invagination. This is in alignment with the perturbation induced by MEK kinase which prevents invagination for happening in the mutant endoderm.

## Discussion

Ray Keller's roadmap of morphogenesis studies establishes a clear path for understanding the biomechanical processes involved in development (4). In his proposed workflow, the first step is to determine when and where cells move. Identifying regions of significant morphological activity in space and time has usually followed a script consisting of observing via a microscope the developing system, formulating an hypothesis of what is happening in the system, and subsequently affirming or refuting the hypothesis using quantitative analysis. This method, which has remarkably propelled the field of developmental biology to its current heights, nevertheless has some limitations. Distinct morphogenesis events can overlap both in space and time, rendering eye observation vulnerable to misinterpretation. Second, these methods are not automated, hence do not scale.

In this work, we attempted to develop an alternative approach to probing development in living systems. Our approach takes advantage of the recent boom in the availability of single cell shape tracking data to propose a generic method able to identify interesting defining morphological processes through space and time in developing embryos. The method takes as input data consisting of evolving cellular geometries and outputs a series of spatial heatmaps showcasing in the time-frequency domains the most important traits of morphogenesis in the studied embryo. Our framework presents over the traditional eye test method multiple advantages. First, the workflow is fully automated, providing an unprecedented hands-off approach in preliminary studies of morphogenesis. Another outstanding advantage over traditional is that our workflow is able to compress the story of the development, such that, in a single image, one can grasp the essence of morphogenesis in a system of interest. In particular, our method has been able to neatly discriminate between the gastrulation and neurulation phases of ascidian early development, identify a second phase of gastrulation: *blastophore closure* which follows invagination, reconstitute the two-step process of endoderm invagination during the gastrulation phase, while clearly distinguishing between short scale division events and low frequency tissue-wide deformations.

In order to achieve this fit, raw cell shape data underwent a series of transformations including a level-sets driven homeomorphic map of the unit sphere to the developing embryo's surface, the computation of the strain rate field of embryo deformations through time using successive iterations of this map, a spherical harmonics decomposition of this strain-rate field, and wavelet decomposition of the most significant spherical harmonics time series. Each of these transformations comes with its own challenges, but also delivers new perspectives for the study of living systems. Our level set scheme excels at defining a homeomorphic map between the unit sphere and surface of the embryo. It goes without saying that in order for the deformed sphere to best match the shape of the embryo, a high sampling of points on the unit sphere is required. A compromise
<figure><img src="image_11.png" /><figcaption>b. Timeseries of $f_{00}$ and $f_{10}$</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>c. Wavelet transform on $f_{00}(t)$ and $f_{10}(t)$ for Mutant embryo</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>**Figure 5.** Spectral decomposition of morphogenesis in mutant embryo. **a)** *Top* Ascidian mutant morphogenesis. *Bottom* Spatiotemporal scalar strain rate field mapped to the mutant surface. **b)** Time evolution of the coefficients $f_{00}$ and $f_{10}$ associated with spherical harmonics ($Y_{00}$ and $Y_{10}$). The cell population dynamic is also included in the plot for clarity. **c)** Wavelet transform applied on $f_{00}(t)$ (top) and $f_{10}(t)$ (bottom).</figcaption></figure>

is however necessary between this sampling and, on the one hand, the overall spatial resolution of the raw data, on the other hand, the induced computational complexity. In its current form, the scheme produces approximations of Lagragian particles only under the assumption of small deformations in the embryo. Hence, the sampling rate during microscopy imaging is of critical importance: the shorter intervals between two successive frames of the movie, the more Lagragian-like the particles are expected to behave.

Given the provision of tracked surface particles meshed at every frame in a triangular network, the evaluation of the strain
rate field is straightforward, and enables, among others, a unified description of complex cell-level and tissue-level dynamics (19), such as drastic deformations and synchronised divisions. The accuracy of this field is affected, as previously, by both the spatial sampling of material points on the unit sphere and the timely sampling of morphogenesis frames. A visualisation of the eigenvalue field derived from this tensor field on the surface of the embryo can already highlight significant processes in morphogenesis. The decomposition of this field into spherical harmonics allows a better appreciation of the spatial patterns of morphological activity in the embryo, each harmonic mapping a region of space. Our spherical harmonics decomposition of morphogenesis results in a set of timeseries of coefficients associated with each harmonic, representing, to the best of our acknowledge, the first comprehensive timeseries-based description of morphogenesis. This transformation unleashes the full power of signal processing tools into studies of morphogenesis. The basis of spherical harmonics being infinite, a challenge here is to discriminate between harmonics that significantly contribute to the composed signal and those that do not. Furthermore, the filtration of principal harmonics modes enables the representation of morphogenesis in a significantly compressed form, in comparison to the initial datasets. This lower dimensional representation of morphogenesis can be helpful, among others, in modelling the physical dynamics of the system (27).

Despite describing canonical interactions in the space of spherical harmonic functions, our spherical timeseries still represent composed signals in time. We use the *Ricker* Wavelet as a mathematical microscope to zoom-in and zoom-out through these signals in order to identify their fundamentals components. This operation results in two-dimensional time-frequency heatmaps that showcase, for each time-series, the footprint of its canonical high and low frequency components, which can be mapped to biological processes. The sum of these tell the story of morphogenesis in the region corresponding to the spherical harmonic. Here, it is important to calibrate the wavelet in such a way that the undesirable effects of the uncertainty principle are avoided. It might also be useful to wisely target windows of time of interest. The wavelet heatmaps can be fed to analytic workflows such as deep neural networks for further studies. Example scenarios could include variational studies of morphogenesis processes in different wild-type or mutant embryo. Furthermore, the workflow presented in this paper can be applied to the examination of single cells morphological behaviour in development.

## References

1. Satoh, N. Cellular morphology and architecture during early morphogenesis of the ascidian egg: an sem study. *The Biol. Bull.* **155**, 608–614 (1978).

2. Nishida, H. Cell division pattern during gastrulation of the ascidian, halocynthia roretzi: (cell division pattern/gastrulation/neurulation/ascidian embryo). *Dev. growth & differentiation* **28**, 191–201 (1986).

3. Jeffery, W. R. A gastrulation center in the ascidian egg. (1992).

4. Keller, R., Davidson, L. A. & Shook, D. R. How we are shaped: the biomechanics of gastrulation. *Differ. ORIGINAL ARTICLE* **71**, 171–205 (2003).

5. Lemaire, P. Unfolding a chordate developmental program, one cell at a time: invariant cell lineages, short-range inductions and evolutionary plasticity in ascidians. *Dev. biology* **332**, 48–60 (2009).

6. Sherrard, K., Robin, F., Lemaire, P. & Munro, E. Sequential activation of apical and basolateral contractility drives ascidian endoderm invagination. *Curr. Biol.* **20**, 1499–1510 (2010).

7. Hashimoto, H., Robin, F. B., Sherrard, K. M. & Munro, E. M. Sequential contraction and exchange of apical junctions drives zippering and neural tube closure in a simple chordate. *Dev. cell* **32**, 241–255 (2015).

8. Hashimoto, H. & Munro, E. Dynamic interplay of cell fate, polarity and force generation in ascidian embryos. *Curr. opinion genetics & development* **51**, 67–77 (2018).

9. Guignard, L. et al. Contact area-dependent cell communication and the morphological invariance of ascidian embryogenesis. *Science* **369** (2020).

10. Jones, R. A., Renshaw, M. J., Barry, D. J. & Smith, J. C. Automated staging of zebrafish embryos using machine learning. *Wellcome Open Res.* **7**, 275 (2022).

11. Tassy, O., Daian, F., Hudson, C., Bertrand, V. & Lemaire, P. A quantitative approach to the study of cell shapes and interactions during early chordate embryogenesis. *Curr. Biol.* **16**, 345–358 (2006).

12. Stelzer, E. H. Light-sheet fluorescence microscopy for quantitative biology. *Nat. methods* **12**, 23–26 (2015).
13. Power, R. M. & Huisken, J. A guide to light-sheet fluorescence microscopy for multiscale imaging. Nat. methods 14, 360–373 (2017).

14. Michelin, G. et al. Cell pairings for ascidian embryo registration. In 2015 IEEE 12th International Symposium on Biomedical Imaging (ISBI), 298–301 (IEEE, 2015).

15. Stegmaier, J. et al. Real-time three-dimensional cell segmentation in large-scale microscopy data of developing embryos. Dev. cell 36, 225–240 (2016).

16. Leggio, B. et al. Morphonet: an interactive online morphological browser to explore complex multi-scale data. Nat. communications 10, 1–8 (2019).

17. Humphrey, J. D. Continuum biomechanics of soft biological tissues. Proc. Royal Soc. London. Ser. A: Math. Phys. Eng. Sci. 459, 3–46 (2003).

18. Ambrosi, D. et al. Perspectives on biological growth and remodeling. J. Mech. Phys. Solids 59, 863–883 (2011).

19. Blanchard, G. B. et al. Tissue tectonics: morphogenetic strain rates, cell shape change and intercalation. Nat. methods 6, 458–464 (2009).

20. Humphrey, J. D. Cardiovascular solid mechanics: cells, tissues, and organs (Springer Science & Business Media, 2013).

21. Streichan, S. J., Lefebvre, M. F., Noll, N., Wieschaus, E. F. & Shraiman, B. I. Global morphogenetic flow is accurately predicted by the spatial distribution of myosin motors. Elife 7, e27454 (2018).

22. Mitchell, N. P. et al. Morphodynamic atlas for drosophila development. bioRxiv 2022–05 (2022).

23. Bar-Kochba, E., Toyjanova, J., Andrews, E., Kim, K.-S. & Franck, C. A fast iterative digital volume correlation algorithm for large deformations. Exp. Mech. 55, 261–274 (2015).

24. Stout, D. A. et al. Mean deformation metrics for quantifying 3d cell–matrix interactions without requiring information about matrix material properties. Proc. Natl. Acad. Sci. 113, 2898–2903 (2016).

25. Patel, M., Leggett, S. E., Landauer, A. K., Wong, I. Y. & Franck, C. Rapid, topology-based particle tracking for high-resolution measurements of large complex 3d motion fields. Sci. reports 8, 1–14 (2018).

26. Dalmasso, G. et al. 4d reconstruction of developmental trajectories using spherical harmonics. bioRxiv (2021).

27. Romeo, N., Hastewell, A., Mietke, A. & Dunkel, J. Learning developmental mode dynamics from single-cell trajectories. Elife 10, e68679 (2021).

28. Mitchell, N. P. & Cislo, D. J. Tubular: Tracking deformations of dynamic tissues and interfaces in 3d. bioRxiv (2022).

29. Godard, B. G. & Heisenberg, C.-P. Cell division and tissue mechanics. Curr. opinion cell biology 60, 114–120 (2019).

30. Dokmegang, J., Yap, M. H., Han, L., Cavaliere, M. & Doursat, R. Computational modelling unveils how epiblast remodelling and positioning rely on trophectoderm morphogenesis during mouse implantation. Plos one 16, e0254763 (2021).

31. Dokmegang, J. Modeling epiblast shape in implanting mammalian embryos. In Epiblast Stem Cells, 281–296 (Springer, 2022).

32. Zhao, H.-K., Osher, S., Merriman, B. & Kang, M. Implicit and nonparametric shape reconstruction from unorganized data using a variational level set method. Comput. Vis. Image Underst. 80, 295–314 (2000).

33. Alba, V., Carthew, J. E., Carthew, R. W. & Mani, M. Global constraints within the developmental program of the drosophila wing. Elife 10, e66750 (2021).

34. Mancinelli, C., Livesu, M. & Puppo, E. Gradient field estimation on triangle meshes. In STAG, 87–96 (2018).

35. De Goes, F., Butts, A. & Desbrun, M. Discrete differential operators on polygonal meshes. ACM Transactions on Graph. (TOG) 39, 110–1 (2020).
36. Knaack, R. & Stenflo, J. O. Spherical harmonic decomposition of solar magnetic fields. Astron. & Astrophys. 438, 349-363 (2005).

37. Dahlen, F. & Tromp, J. Theoretical global seismology. In Theoretical Global Seismology (Princeton university press, 2021).

38. Cislo, D. J. et al. Active cell divisions generate fourfold orientationally ordered phase in living tissue. Nat. Phys. 1-10 (2023).

39. Torrence, C. & Compo, G. P. A practical guide to wavelet analysis. Bull. Am. Meteorol. society 79, 61-78 (1998).

40. Lau, K.-M. & Weng, H. Climate signal detection using wavelet transform: How to make a time series sing. Bull. Am. meteorological society 76, 2391-2402 (1995).

41. Brunton, S. L. & Kutz, J. N. Data-driven science and engineering: Machine learning, dynamical systems, and control (Cambridge University Press, 2022).