An integrated modelling framework from cells to organism based on a cohort of digital embryos

We conducted a quantitative comparison of developing sea urchin embryos based on the analysis of five digital specimens obtained by automatic processing of in toto 3D+ time image data. These measurements served the reconstruction of a prototypical cell lineage tree able to predict the spatiotemporal cellular organisation of a normal sea urchin blastula. The reconstruction was achieved by designing and tuning a multi-level probabilistic model that reproduced embryo-level dynamics from a small number of statistical parameters characterising cell proliferation, cell surface area and cell volume evolution along the cell lineage. Our resulting artificial prototype was embedded in 3D space by biomechanical agent-based modelling and simulation, which allowed a systematic exploration and optimisation of free parameters to fit the experimental data and test biological hypotheses. The spherical monolayered blastula and the spatial arrangement of its different cell types appeared tightly constrained by cell stiffness, cell-adhesion parameters and blastocoel turgor pressure.

. c) Affine transform parameters (α e , β e ) used in the previous top curves. d) Total cellular volumes W e (t). e) Temporally and spatially rescaled total cellular volumes, (γ e ) 3 W e (α e t + β e ) (top), and field cellular volumes (γ e ) 3 W k e (α e t + β e ) (same order as b), using the same parameters (α e , β e ) as before. f) Coefficients γ e used in the previous top curves. g) Total cellular surface areas Z e (t). h) Temporally and spatially rescaled total cellular surface areas, (γ e ) 2 Z e (α e t + β e ) (top), and field cellular surface areas (γ e ) 2 Z k e (α e t + β e ) (same order as b), using the same parameters (α e , β e ) and γ e as before. i) Temporally rescaled total numbers of contacts C e (α e t + β e ) and field numbers of contacts C k e (α e t + β e ).
Supplementary Figure 3: Mean and standard deviation bars representing the normal and log-normal approximations of the cell feature distributions in various cell groups L g . For each feature (one per frame), the cell groups g = (e, n, k) cover all four cell types k (except cell cycle length), one or several generations n ∈ [6,10] per type, and most of the five embryos e per generation-type combination (using the same colours blue, green, orange, cyan, pink for each embryo as Supplementary Fig. 2 with additional bars (black, bold) for the prototype).
Supplementary Figure 4: Same as Supplementary Fig. 3 c to f without applying the log function to the volume and surface area. This figure is intended to highlight the value 0.5 in the daughter/mother volume ratio (here in c).
Supplementary Figure 5: Top Panel (above the red line) -The 12 cell feature distributions (two per variable type) with the highest p-value, i.e. the best cases of fit by normal or log-normal curves with respect to the chi-squared test. Left to right, top to bottom: the two X g (cell cycle length) and two M g (division time) distributions closest to normal curves; the two V g (average volume), two S g (average surface area), two A g (daughter/mother volume ratio) and two B g (daughter/mother surface area ratio) distributions closest to lognormal curves. Titles indicate the embryo e, cell type k, generation n and p-value. Bottom Panel (below the red line) -The 12 worst-fit cases (p-value < 0.01) of cell feature distributions, where the assumption of normal or log-normal curves could not be retained with respect to the chi-squared goodness-of-fit test. They fall into four n, k categories: 8, Mac; 8, Mes; 9, Mac; and 9, Mes. Most of them concern X (cell cycle length) and M (division time). -See Supplementary Note 2.2.1 Supplementary Figure 6: Histogram of all p-values obtained by chi-squared goodness-of-fit test on 124 distributions of cell features. The six variables considered here were X g (cell cycle length), M g (division time), logV g (log average volume), log A g (log average surface area), log S g (log daughter/mother volume ratio) and log B g (log daughter/mother surface area ratio) across the five embryos e, four cell types k, and all generations n within the period during which the cell feature was observable.  Figure 11: Simulation of artificial cell lineages and comparison with real data. The three rows correspond to embryos 1, 3, and 5. The three columns represent the total number of cells N(t), the total cellular volume W (t) and surface area Z(t). Inside each chart, the top curves show the embryo-level quantities (N e , W e , or Z e ), while the four groups of lower curves show their type-specific components (N k e , W k e , or Z k e ), from top to bottom: Mes, Mac, LMic, SMic. In black: mean of 300 simulated cell lineages; in gray: their standard deviation; in colour: measured values in the real specimen (all curves obtained by temporal and spatial rescaling). Shaded areas: incomplete cycles, for which an accurate estimation was not possible.
Supplementary Figure 12: Representation of cells by cylindrical particles. Each cell i or j is defined by a lateral radius R i , a half-height R ⊥ i and a normal vector U i (R j , R ⊥ j , U j respectively).
Supplementary Figure 13: Comparison of the attraction-repulsion force used here with alternative classical forces. The equilibrium distance is r eq i j = c eq (R i + R j ) with R i = 0.1 and R j = 0.15. The solid line represents the attraction/repulsion potential F i j with w adh = w rep = 500. The dashed line represents a force derived from the Morse potential: F M (r) = 2Dke −k(r−r eq i j ) − 2Dke −2k(r−r eq i j ) with D = 0.265 and k = 10. The dot-dashed line is a force derived from the Lennard-Jones potential: with ε = 0.117 and σ = 2 −1/6 r eq i j . Parameters were selected to make the equilibrium distance and maximum value of the three curves approximately match on the plot.
Supplementary Figure 14: Diagram of the planar rigidity force between neighbouring epithelial cells. Here, cells i, j and k are represented in 2D, with their respective normal axes U i (yellow arrow), U j (red arrow) and U k (blue arrow). The planar rigidity forces exerted between neighbour cells i and j, and between i and k are aligned with the bisectors of the normal vectors n i j (orange), and n ik (green) respectively. Forces are represented by thicker arrows. The intensity of the force is less between i and k because the relative position vectors u ik and n ik are nearly orthogonal. As a result of these forces, neighbour cells are always pulled back to the lateral domain of each cell, thus maintaining the planarity of the tissue.   Cell group L g with g = (e, n, k) with g − 1 = (e, n − 1, k) and g = (e, n < n, k) Table 4: Summary of the relationships between the different random variables in the multilevel probabilistic model. The random variables considered are the division time M, the cell cycle length X, the average volume V , the average surface area S, the daughter/mother volume ratio A, the daughter/mother surface area ratio B. See details in Supplementary Note 2.3.
Cells j ∈ L g−1 and i ∈ L g i is one of j's two daughters b i is a realization of B g s j is a realization of S g−1 φ g is the deterministic surface area microdynamics u is the percentage of elapsed cycle length of cell i Supplementary  Table 6: Exploration range and increments, with chosen values, for the parameters of the spatial model. w adh,e is the heterotypic adhesion coefficient, w adh,o is the homotypic adhesion coefficient, w rep is the repulsion coefficient, c max is used to define the maximum cell deformation in the surface plane, these four coefficients are involved in the computation of the attraction-repulsion force. k rig is the rigidity coefficient used to compute the planar rigidity force. α gab is a coefficient used to define the neighbourhood of each cell through a 3D generalization of the Gabriel criterion. λ 0 is used to compute the damping coefficient, required to compute the equation of motion. The full description of the spatial model can be found in Supplementary Note 3.1 Supplementary Video 1: Visualization of reconstructed embryo 3 visualized from 4.56 hpf (32-cell stage) to 8.65 hpf (hatching blastula) with Mov-IT software: A) Raw nuclei, coloured dots are detected nuclear centres (mesomeres in blue, macromeres in red, large micromeres in pink, small micromeres in purple). Volume rendering of raw images from two-photon laser scanning microscopy with nuclear staining by H2B-mCherry. Note that at the onset of imaging, membrane signal is leaking in the nucleus channel. Lateral view animal to the top. B) Raw membranes, the coloured dots are the same as in A). Three orthoslices (xy imaging plane, xz and xy planes orthogonal to the imaging plane) showing raw images from two-photon laser scanning microscopy with membrane staining by farnesylated eGFP. The animal orientation can be inferred from the organization of the cell populations. C) Surface rendering of segmented cell membranes.
(mesomeres in blue, macromeres in red, large micromeres in pink, small micromeres in purple). Lateral view animal to the top. D) Same as C). The embryo is cut along a radial plane and shown from the blastocoel.
Supplementary Video 2: Modelling of sea urchin prototypical early development with the MecaGen platform and simulation with four different sets of mechanical parameters: A) Simulation of a monolayered, spherical embryo with realistic cell population borders. This specimen falls into the best-fit domain. It was obtained with the following parameters: w adh,o = 600, w adh,e = 110, w rep = 100, k rig = 10000, α gab = 1, λ 0 = 3000, and c max = 2. The embryo is cut along a radial plane and shown from the blastocoel. Black lines materialize cell-cell contacts. B) Simulation of a monolayered, spherical embryo with nonrealistic cell population borders. Although the global shape is correct, this simulated embryo shows cell population borders that differ from the live specimens with nonrealistic cell mixing (pink area in Fig. 4c). It was obtained with the following parameters: w adh,o = 100, w adh,e = 200, w rep = 100, k rig = 10000, α gab = 1, λ 0 = 3000, and c max = 2. The embryo is cut along a sagittal plane and shown from the blastocoel. Black lines materialize cell-cell contacts. C) Simulation of a monolayered, aspherical embryo with nonrealistic cell population borders. The simulated embryo still displays an epithelial monolayer but not the appropriate sphericity (in addition to cell population borders that differ from the live specimens with nonrealistic cell mixing). Since the embryo collapsed, it was not possible to define a radial plane. It is displayed in toto from the outside with the whole scaffold of cell-cell contacts (black lines) and was obtained with the following parameters: w adh,o = 400, w adh,e = 220, w rep = 100, k rig = 10000, α gab = 1, λ 0 = 3000, and c max = 2. D) Simulation of a multilayered, aspherical embryo with nonrealistic cell population borders. This simulated specimen cumulates all three defects and crumpled upon itself. It was not possible to define a radial plane in this case. The embryo is cut along a longitudinal plane and shown from the blastocoel. Black lines materialize cell-cell contacts. It was obtained with the following parameters: w adh,o = 700, w adh,e = 700, w rep = 100, k rig = 10000, α gab = 1, λ 0 = 3000, and c max = 2.

Supplementary Note 1 Digital embryo reconstruction
As anticipated in [1], we provide the first complete digital reconstruction of the sea urchin embryo's early development, comprising the position and shape of every cell over time. The methods and tools used here are part of the BioEmergences platform developed in our laboratory [2]. Cell positions were detected at the local maxima of cell nucleus images filtered by a difference of Gaussians (DoG), while cell shapes were reconstructed from cell membrane images using the "subjective surface" algorithm (SubSurf) [3], following earlier results [4]. To reconstruct cell trajectories, cell movements were tracked by minimizing a heuristic functional with simulated annealing (SimAnn) [5]. The choice of parameters for DoG filtering was supervised, and all cell positions and trajectories were manually checked and validated through our dedicated Mov-IT software interface [6,2].

Supplementary Note 1.1 Nucleus center detection
In order to detect cell positions, we assumed that cells could be identified by the spatial coordinates of nuclear centres. Raw imaging data was preprocessed to simultaneously smooth images and preserve only the spatial information relevant to the size of cell nuclei. Filtering was carried out by a convolution of the original image with two Gaussian functions G 1 and G 2 of variances σ 1 < σ 2 , and calculating the difference

Supplementary Note 1.2 Cell tracking
The objective of cell tracking is to detect the invariant identity of each cell by relating its positions at two consecutive time steps, from its birth at the last mitosis to its disappearance at the next mitosis. Starting with the initial cell population, temporal links create together a lineage tree whose branches split at cell divisions. Considering that sea urchin embryos were immobilized and cells did not move fast at the observed early stages of development, our tracking strategy essentially relied on a nearest-neighbour criterion to link detected nuclear centres. First, tree branches were constructed by connecting nearest cells on consecutive images when that proximity criterion was reciprocal. Then, mitoses were detected by linking cells left over without a predecessor to their closest neighbour at the previous time step. Finally, the resulting lineage tree was refined via minimization of a functional E( f ). This energy term assumes that the embryo behaves like an elastic mass, therefore cells that are neighbours at time t remain close to each other at time t + 1: where i, j represent the cells, f is the link between cells at time t and cells at time t + 1, and r i = (x i , y i , z i ) is the cell position. The minimization of the functional was performed by a simulated annealing algorithm [5], searching for the lowest-energy configuration by randomly changing a link, removing/adding a link or removing/adding a cell. Finally, we refined the automated tracking by removing cells with a very short cell cycle (e.g. a few time steps), as they were most likely to be false positives.

Supplementary Note 1.3 Digital specimen validation and correction
All cell positions and trajectories were checked by at least two experts through visual inspection comparing the digital reconstruction with the raw data. All detectable errors (false positive and false negative nuclei, false links, missing links) were corrected to build error-free lineages-hereafter called "gold standard". All five lineage trees were fully validated for all five datasets in this manner. Moreover, each cell was unambiguously identified by a unique ID in our database. Both operations were a prerequisite to perform the next steps comprising cell segmentation, which uses cell positions for initialization and cell IDs for labelling, and data analysis.

Supplementary Note 1.4 Cell segmentation
Three-dimensional segmentation of in vivo imaging presents a difficult challenge, especially since fast and deep two-photon laser scanning microscopy produces a low signal-to-noise ratio and incomplete membrane contours, thus dedicated methods have to be used for this task. To this purpose, cell shapes were obtained by applying a generalized version of the SubSurf method [4] for its ability to reconstruct missing boundaries, making it particularly suitable for this type of segmentation. To remove the noise and smooth the images while faithfully preserving edge information, the data was first filtered by geodesic mean curvature flow (GMCF) [7]. Segmentation of cell membranes was then performed on the filtered data using the manually validated cell positions as starting points. Each membrane was segmented independently on a region of interest located around the initialization points and ranging from ±25µm to ±16µm in every direction depending on the number of cells in the processed volume. Finally, the whole embryo shape was recomposed from the union of the segmented cells and each region labelled according to the ID of the corresponding cell centre. Information was exchanged between adjacent membranes to avoid overlapping. In particular, superimposed regions were split depending on the distance from the nearest cell centre and the ID of the obtained subregions set accordingly. Segmentation results were exported in two formats: voxel-based image data to support the calculation of cell volumes and the detection of neighbour cells; and triangulated surfaces to support the calculation of cell surface areas. The GMCF and generalized SubSurf methods are briefly summarized in the next paragraphs. Numerical implementation of both filtering and segmentation algorithms was performed via finite difference schemes [8,3]. Further details about the methods employed can be found in [7] and [4].

Supplementary Note 1.4.1 Image filtering
We model the 3D input image by a real positive function I 0 (x, y, z), in some spatial rectangular domain Ω ⊂ R 3 , and its low-level local features by an edge indicator function g (see below). According to GMCF, a smooth-edge enhanced version of the original image can be obtained by extending the image function in time, creating a family of "level sets" I(x, y, z,t), and solving the following partial differential equation: where G σ is a Gaussian kernel with a small variance σ . This is accompanied by the initial condition I(x, y, z, 0) = I 0 (x, y, z) and a null Neumann boundary condition (i.e. zero gradients) on the domain edges ∂ Ω. Selecting a final time t = T , the solution I(x, y, z, T ) represents the filtered image I f . In the model, the mean curvature motion of the level sets of function I is determined by the edge indicator function g(s) = 1/(1 + Ks 2 ), K ≥ 0, applied to the image intensity gradient pre-smoothed by G σ .

Supplementary Note 1.4.2 Image segmentation
In the 3D case, the SubSurf method consists of minimizing the volume of a 3D manifold S embedded in a 4D Riemannian space with a metric constructed on the image features [4]. We consider the filtered image I f and its low-level local features given by the same edge detector function g f (x, y, z) = g(I f (x, y, z)) as for the GMCF. Let also the manifold be defined as S = (x, y, z, Φ), where Φ(x, y, z) is a positive distance function defined in the image domain Ω. The generalized version of the SubSurf method leading to the minimization of the volume of S can be expressed as a function of Φ as follows: where µ and ν are weight coefficients and τ is the value of the scale parameter t which corresponds to the steady-state condition of the first equation. The notation |∇Φ| ε = ε + |∇Φ| 2 represents a regularization of |∇Φ| [9] and H reads: where subscripts are shortcuts for derivatives: The first term on the right-hand side of equation 3 (top row) represents a mean curvature flow weighted by the edge indicator g f . The second term attracts the level sets of Φ in the direction of the image edges, given by the velocity field −∇g f . Different behaviors can be identified locally in the image regions according to one of these flows, and missing boundaries (the so-called "subjective contours"), which are a continuation of existing edge fragments, are completed by geodesics. The use of different weights between the regularization and the advective term (ν > µ) facilitates the control of the dynamic process, while the parameter ε is used to shift the model from the mean curvature flow of the level sets (ε = 0) to the mean curvature flow of the graph (ε = 1).

Supplementary Note 2 Multi-level probabilistic model of the cell lineage Supplementary Note 2.1 Multi-level measures and rescaling Supplementary Note 2.1.1 Individual cell features
The digital reconstruction of entire sea urchin specimens allowed us to extract a variety of features at the level of individual cells (Supplementary Note 1) [4]. For any cell i, we defined the following quantities (Supplementary Fig. 1): • its cell cycle length: x i , corresponding to the time between two consecutive divisions • a mitosis time: m i , denoting the moment at which the cell i divides into two daughter cells • its current volume: v i (t), with the average volume over its cell cycle length (using discrete time steps): • its current surface area: s i (t), with the average surface area over its cell cycle length: • a generation number: n i , denoting the number of past divisions on cell i's lineage branch, which includes itself and all its ancestors since fertilization • its current degree: d i (t), equal to the number of neighbours of cell i at time t, assuming that the spatial distribution of cells can be represented by an undirected graph of cellular contacts, in which cells correspond to nodes and cell-cell interactions to edges [11,12].
Each cell i also has a mother identified by j (a shorthand notation for the function j(i) mapping any cell index to its mother's index), which is used in the definition of two additional features: • a daughter/mother average volume ratio: a i = v i /v j • a daughter/mother average surface ratio: b i = s i /s j and three relationships: • the mitosis time of a cell is equal to the sum of its cell cycle length and the mitosis time of its mother: • the generation is incremented by one at each division: n i = n j + 1 • starting from the 32-cell stage, each cell type is conserved across divisions throughout the period considered here (formation of the blastula): k i = k j .

Supplementary Note 2.1.2 Cell populations and global state variables
The cell lineage, denoted by L , is the set of all cells, past and present, that the embryo contained during a given period of time. From there, various subsets were defined: • the current embryo, the set of all cells alive at time t: • cell populations, the sets of cells of a given type k across all generations: At the level of the entire embryo, and at each time step t, the global measures were: • the number of cells, cardinality of the current embryo: • the total cellular volume, sum of all individual cell volumes: • the total cellular surface area, sum of individual cell surface areas: • the total number of contacts, sum of the local degrees: These quantities were also calculated inside each cell population k, replacing L (t) by L k (t) in the above definitions and using the notations N k (t), W k (t), Z k (t), and C k (t).

Supplementary Note 2.1.3 Rescaled embryo dynamics
The number of cells, total cellular volume, and total cellular surface area were measured empirically in each embryo of the cohort E , indexed by e ∈ [1,5] and denoted N e (t), N k e (t), and so on. The evolution over time of these variables was plotted concurrently ( Supplementary Fig. 2a,d,g). We observed that the patterns of evolution across the various specimens were similar but did not overlap: some were globally faster, some slower; some were larger, other smaller. To filter out this variability and proceed to further interindividual comparison, we applied temporal and spatial rescaling functions.
Temporal rescaling For each observed cell number N e (t), we considered its inverse function t e (N) equal to the time at which embryo e contained a given number N of cells (symmetric of Supplementary Fig. 2a).
Since N e (t) is a monotonically increasing function (no cell death during the period considered), so is t e (N). However, because of the relative synchrony among cell cycles, N e can remain constant for several time steps (plateaus in Supplementary Fig. 2a). Therefore, to obtain a single-valued function, we defined t e (N) = min{t | N e (t) = N} over the size interval of e, [N e,min , N e,max ] (if N was not part of the observed discrete values, t e (N) was interpolated). Then, to minimize the discrepancy between the observed time of an embryo, t e (N), and the mean time of the five specimens, t (N), an affine transform was performed on the time axis using a cost function: where E N is the sublist of embryos e such that N ∈ [N e,min , N e,max ] (including by interpolation). For each embryo, the parameters adopted for the transform satisfied (α e , β e ) = arg min F e (α, β ). Finally, taking the inverse again, we obtained five rescaled total numbers N e (α e t + β e ) and four groups of five rescaled field numbers N k e (α e t + β e ), which are plotted in Supplementary Fig. 2b. The five parameter pairs (α e , β e ) are shown in Supplementary Fig. 2c.
Spatial rescaling Based on the same temporal rescaling, a linear transform was also performed along the spatial dimensions to minimize the discrepancy between the temporally rescaled cellular volume of an embryo, W e (α e t + β e ) and its mean value over all specimens: where E t is the sublist of embryos e such that (α e t + β e ) ∈ [t e,min ,t e,max ] (including by interpolation). For each embryo, γ e = arg min G e (γ) was then defined as the optimal coefficient. The original volumes W e (t) are plotted in Supplementary Fig. 2d. The five rescaled total volumes (γ e ) 3 W e (α e t + β e ) and four groups of five rescaled field volumes (γ e ) 3 W k e (α e t + β e ) each are plotted in Supplementary Fig. 2e, while coefficients γ e are shown in Supplementary Fig. 2f. Finally, we used the same coefficients squared for the rescaled cellular surface areas, (γ e ) 2 Z e (α e t + β e ) and (γ e ) 2 Z k e (α e t + β e ) (Supplementary Fig. 2h; the original curves Z e (t) are shown in Supplementary Fig. 2g).
neighbourhoods The total numbers of contacts in each embryo and each cell population did not require an extra spatial transform to highlight interindividual similarities. After carrying over the temporal rescaling parameters found above, the curves C e (α e t +β e ) and C k e (α e t +β e ) already presented a high degree of overlap among embryos e ∈ [1, 5] ( Supplementary Fig. 2i). This is consistent with the fact that neighbourhood relationships are topological features, thus independent from metric deformations.

Supplementary Note 2.1.4 Intermediate cell groups
Symmetries in the embryo, such as the rotational symmetry around the animal-vegetal (AV) axis, prevent the identification and matching of individual cells from one specimen to another. Unique identification of cells based on their morphological characteristics cannot be done without ambiguity either. To overcome these issues, a generic coarse-grained level of observation was needed. We chose here to define new subsets of L (adding to the list of Supplementary Note 2.1.2): • generational cell groups, or simply "cell groups", the subsets of generation-n cells of a given type k: Note that L n,k is not the same as L k (t): the latter is a snapshot of cell population k at time t and therefore may contain a mix of generations if some cells have divided more frequently than others. Conversely, cells in the former group may have to be taken at different times. These sets offer two useful viewpoints on the embryo, one focusing on its global state, the other on cell statistics. In any case, the greater n and t, the more distant these two sets are likely to be. Since each cell gives rise to two daughter cells at each cycle and there is no cell death, the number of cells of a given type continued doubling, i.e. |L n,k e | = 2 n−6 |L 6,k e | for n ≥ 6 during the period of interest. The six cell features considered here are the cell cycle length x i , mitosis time m i , average cell volume v i , average cell surface area s i , and the daughter/mother ratios of average volume a i and surface area b i . The dispersion of these features observed in each group L n,k e of the empirical data is represented by distributions, modelled by sequences of random variables. For example, (x 1 , x 2 , ..., x m ) denotes the sequence of the values that the first feature, cell cycle length, takes in the m cells of a given group L n,k e . Each x i is the realization of a random variable X i pertaining to cell i. Cells being indistinguishable within the same group, their order in the sequence is irrelevant. This property is referred to as "exchangeability" of the sequence of random variables and is formalized as follows: for any permutation π of the indices [1, m], the joint probability distribution of these variables is the same: P(X π(1) , X π(2) , ..., X π(m) ) = P(X 1 , X 2 , ..., X m ). 9 This property leads to de Finetti's theorem, which can be summarized by saying that an infinite sequence of exchangeable random variables is a "mixture" of independent and identically distributed (i.i.d.) sequences [13,14]. As for a finite exchangeable sequence, it can also be approximated by a mixture of i.i.d. sequences [15]. One consequence of this theorem is that the following empirical measure P m,X g constitutes a summary statistic for the probability density P of the X sequence: where I is an interval of values or union of intervals (i.e. an element of the "σ -algebra" on R) and δ x is the Dirac measure. If m could increase to ∞, P m,X g would converge asymptotically to the underlying theoretical probability distribution P of the sequence of exchangeable random variables.
The motivation of the probabilistic model presented here and in the next section was to provide an idealized representation of the dynamics of multicellular development relating the distributions observed in each cell group with the global dynamics of the cell lineage. To simplify notations, let us define an embryospecific cell group index g = (e, n, k) so that L n,k e can be equivalently written L g . With this, our goal was the following: • link the distribution in each group L g of cell cycle lengths: X g ∼ {x i | i ∈ L g } and the distribution of mitosis times: M g ∼ {m i | i ∈ L g } to the evolution of the total number of cells in the embryo N e (t) • link the distribution in each group L g of average volumes V g and the distribution of daughter/mother volume ratios A g to the evolution of the total cellular volume W e (t) • link the distribution in each group L g of average surface areas S g and the distribution of daughter/mother surface ratios B g to the evolution of the total cellular surface area Z e (t).

Supplementary Note 2.2 Observation and approximation of multi-level statistics
From the empirical distributions of individual cell features, we derived parameterized representations approximating their probability distributions in each cell group. Combining cell features along the cell lineage required considering simple models of inheritance. As explained below, we calculated correlations based on an assumption of linearity in intercellular dependencies and, obtaining low values, concluded that the individual features of a cell were largely independent from its mother's or sister's features. This independence was then taken as a founding hypothesis of our model.

Supplementary Note 2.2.1 Estimation of cell feature distributions in cell groups
To interpret the empirical feature distributions, capture their evolution and predict the developmental dynamics, it was best to describe them in terms of parametric models. First, however, because of the finite and relatively short duration of the period of observation, some groups L g were incompletely represented in the dataset. Therefore, a distribution of cell features in a group was considered to be significant only if the number of observed cells constituted at least 95% of the full cardinality |L g | calculated above. Moreover, certain features, such as the volume v i and surface area s i , required the observation of their evolution during whole cell cycles. Incomplete cell cycles were taken into account only if the mean period of observation of individual cell dynamics in a group was at least 20 min. Based on these criteria, we retained 252 distributions of individual cell features across all cell groups L g (corresponding to a rough average of 4 × 3 significant distribution-generation pairs in each one of the 5 × 4 embryo-cell types). Graphical assessment of these histograms led us to categorize them into two major types: normal distributions (i.e. Gaussian curves) for cell cycle lengths x i and mitosis times m i ; and log-normal distributions for average cell volumes v i , average cell surface areas s i , daughter/mother volume ratios a i and surface area ratios b i (where a random variable X is said to be log-normally distributed if the random variable Y = log X is normally distributed).
We performed a chi-squared goodness-of-fit test on each distribution to validate our assessment. This test evaluates the proximity of the empirical frequency with the theoretical one, i.e. whether the random variables are i.i.d. under a normal distribution or a log-normal distribution. To this aim, we first calculated the empirical mean and standard deviation of each random variable in each cell group L g using a classical maximum likelihood estimation [16]. For example, in the case of cell cycle lengths X g these two quantities are where N g = |L g |. The same formulas were applied to variables M g , logV g , log S g , log A g , or log B g , for each one of the 252 available distributions, yielding parameters ( µ M g , σ M g ), ..., ( µ B g , σ B g ). All 252 (µ, σ ) pairs are plotted in Supplementary Fig. 3. Then, based on the empirical mean and standard deviation, our categorization hypothesis was challenged in each distribution by calculating a p-value corresponding to the probability to find a good fit with a normal or log-normal curve. Taking again X g as an example, the distribution of cell cycle lengths {x 1 , x 2 , ..., x m } over the m cells of group L g was categorized into a fixed number l of discrete bins, (x 1 ,x 2 , ...,x l ), wherê . In this histogram, the observed distribution corresponds to the number of values falling in each bin and was given by the empirical measure P m,X g (I h ) (equation 10), where I h = [x h ,x h + ∆x) and I l is closed, while the expected distribution was given by P m, X g (I h ), where X g is the theoretical normal distribution calculated from ( µ X g , σ X g ). The chi-squared statistic computes the discrepancy between both distributions as follows: This statistic was assumed to follow a chi-squared distribution with κ = l − 3 degrees of freedom [17] (the number of frequencies, l, reduced by the number of parameters of the fitted distribution, µ and σ , and the constraint ∑ h P m,X g (I h ) = 1). This gave the p-value, which is the probability of observing a test statistic at least as extreme, i.e. 1 minus the cumulative distribution function of χ 2 g for κ degrees of freedom. Finally, our Gaussian-fit model was deemed adequate for a given feature in a group L g if its p-value was greater than 0.05 ( Supplementary Fig. 5).
Note that the chi-squared test of goodness-of-fit produces inaccurate results if the size of the sample is too small-i.e., by convention, less than 5 elements in a bin [17]. Reducing the number of frequencies l increases the number of elements in each bin, but also decreases the number of degrees of freedom kappa, which must remain greater than 3 for the test to be applicable here. Based on these constraints, 128 distributions out of 252 had to be excluded from the evaluation.
For most of the remaining 124 distributions, the validity of our assumptions about their type (normal or log-normal) could be retained. Only 20 of them had a p-value under 0.05, among which 12 were under 0.01 ( Supplementary Fig. 5). The complete ranges of p-values obtained were the following: [ Supplementary Fig. 6.
Ideally, the worst-fit cases could be used to subdivide cell groups L g into smaller classes (e.g., oral and aboral cells) reflecting the multimodal aspect of certain distributions, or to adopt different parametric models (e.g., exponential instead of Gaussian curves). Given the small size of the cohort, however, and the non-reproducibility of these outlier distributions in various specimens, it was difficult to rely on such deviations to define new cell types or new models. Moreover, biological relevance of the inferred cell clusters would have to be validated by correlation with the fate map and/or the gene expression in the form of a gene expression "atlas" [18], which is out of the scope of this study. This is why we preferred keeping the simplicity of normal and log-normal approximations as they captured the essential characteristics of statistical distributions (their mean and variance) even in marginal cases. For now, we left out any potential additional information for the sake of understandability and interpretability of the whole dataset. The description of distribution shapes can be refined in future work by using a larger cohort of specimens.

Supplementary Note 2.2.2 Cell volume and surface area dynamics
The volume and surface area of a cell are not constant during its cell cycle length because of various contractions and expansions. To highlight the fluctuations of the cellular geometry around its average, we define two other cell quantities: • the normalized volume: where u ∈ [0, 1] represents the cell's "normalized age", defined as the percentage of elapsed cell cycle length: (m j is the mother's mitosis time). From there, the microdynamics of cell geometry can be described by taking the empirical mean of these temporally realigned quantities within each group L g (Supplementary Figs. 7,8): making the assumption that the functions ω g and φ g are essentially deterministic, i.e. the variations in volume and surface rescaled in amplitude and time are the same for all cells i ∈ L g (hence the notations ω g , φ g instead of µ Ω g , µ Φ g ). Compared to the other six empirical means calculated previously (equation 11), the particularity of these values is their dependency on time, which is rescaled to align the signals. At the level of individual cells, by definition of a i and b i , the following relationships hold:

Supplementary Note 2.2.3 Independence along the lineage
Having defined parametric representations ( µ, σ ) for each cell feature distribution within each cell group L g , we wanted to investigate the relationships among these variables. Since cells are genealogically related, it was natural to hypothesize some degree of correlation between cells belonging to the same descent. With the goal to find a reduced number of parameters describing the phenomenology of the process, we assumed linear dependencies between the cell distributions of daughters and mothers, and among sisters. This is written Y = λ + νX, where the scalar parameters λ and ν can be estimated by linear regression based on a set of N realizations of the random variables X and Y : {(x i , y i )} i=1,...,N . The empirical optimal values λ and ν are the ones that minimize the residual sum of squares over the samples: where ε i denotes the residual error terms such that y i = λ + νx i + ε i and is assumed to be normally distributed around zero. Then, to characterize the accuracy of this linear estimation, we used the coefficient of determination, R 2 , which measures the percentage of variation of one variable explained linearly by the other [19,16]. If we define the total sum of squares by where y is the empirical mean of Y , then the percentage of variation that remains unexplained linearly is given by the ratio SS res /SS tot , and the coefficient of determination is equal to its complement: For a "simple" linear regression (i.e. one with a single explanatory variable), it can be shown that the coefficient of determination is equal to the square of the estimated Pearson's coefficient of correlation: The greater R 2 , the more X and Y tend to be linearly dependent, with a binary decision threshold generally set at 0.6. Note that R 2 does not measure nonlinear dependencies, thus can remain low even if X and Y are strongly related through quadratic or logarithmic functions, for example. Here, this coefficient was applied to six pairs of features calculated among the 252 distributions deemed sufficiently represented in the dataset (Supplementary Note 2.2.1, Supplementary Fig. 3). In the following, for a group index g = (e, n, k) we use the shorthand notation g − 1 = (e, n − 1, k) to refer to the previousgeneration cell groups of the same type in the same embryo.
Daugther's cell cycle length vs. mother's mitosis time As above, we assumed a simple model of linear dependency between the cell cycle length x i of a cell i ∈ L g and the mitosis time m j of its mother j ∈ L g−1 : x i = λ + νm j . The coefficients of determination R 2 = ρ 2 M g−1 ,X g were computed between the distributions X g ∼ {x i | i ∈ L g } and M g−1 ∼ {m j | j ∈ L g−1 } at several generations n. For this calculation, we used the 23 groups L g of Supplementary Fig. 3a in which all cell cycles were completely known, i.e. had observed start and end mitosis times, and their corresponding 23 groups L g−1 . In the end, the values of R 2 that we obtained fell in the range [3e-3,0.51], which indicated only weak or inexistent linear dependencies.
Sisters' cell cycle lengths For each pair of sister cells i 1 , i 2 ∈ L g , ordered such that x i 1 > x i 2 , we also assumed x i 1 = λ + νx i 2 , where the longer cell cycle length was included in a set X g,1 and the shorter one in another set X g,2 . The coefficient of determination R 2 = ρ 2 X g,1 ,X g,2 was then computed in the same 23 groups and the values obtained belonged to the interval [0.095,0.77], revealing some linear dependency among sisters' cell cycle lengths within certain groups. More precisely, R 2 coefficients scored above 0.6 in the following seven groups g = (e, n, k): (4,7,Mac), (5,7,Mac), (3,8,Mac), (4,8,Mac), (1,7,Mes), (3,7,Mes) and (3,9,Mes). In the other 16 groups, there was no clear linear dependency.
Volume and surface area Concerning the volume and surface area, several models of dependency can be explored. Usually, it is assumed that the volume of a mother is conserved in its two daughters through mitosis, and when volumes are measured immediately prior to division and after division, this relationship can even be used to assess segmentation accuracy [6]. In addition, as generations 6 to 10 are part of cleavage stages, the global cellular volume is expected to remain constant while individual cell volumes are decreasing by half at each generation, as confirmed on average by Supplementary Figs. 3, 4C. Thus we may expect average volumes to verify: v i 1 + v i 2 = v j . Dividing by v j , this is equivalent to a linear relationship a i 1 = 1 − a i 2 between the average daugther/mother volume ratios of the sister cells. To verify if there was such a dependency, we computed the coefficient of determination R 2 = ρ 2 A g,1 ,A g,2 as above between pairs of sister cells i 1 , i 2 ∈ L g ordered such that a i 1 > a i 2 . Over the 38 cell groups L g of Supplementary  Fig. 3c in which all volumes could be confidently measured, the coefficients of determination belonged to [1.4e-3,0.88], with 34 values below 0.61 and the other four values in [0.75,0.88] corresponding to cell groups (3,7,SMic), (3,7,Mac), (4,8,Mac) and (1,7,Mes). However, the λ and ν parameters for these groups were found in [−0.45,0.12] and [0.87,3.02] respectively, i.e. far from 1 and −1. Altogether, these results contradicted the hypothesis of average cell volume conservation across mitosis through consecutive generations, as we observed a variation in the range of ±20% around 0.5 for the average daughter/mother volume ratio.
Daughter's volume (or surface area) ratio vs. mother's volume (or surface area) We investigated whether daughter cells compensated for the mother's deviation from the average either in terms of volume or in terms of surface area. Such compensation would mean that the average daughter/mother volume ratio a i of a cell i ∈ L g and the average volume v j of its mother j ∈ L g−1 were correlated. However, the coefficients of determination R 2 computed in the same 38 groups as above ranged in [1.9e-4,0.43], indicating no linear relationship. The same tests were performed on surface area features, yielding similar negative results. All relationships are summarized in Table 3 and the histogram of all the coefficients of determination evaluated is shown on Supplementary Fig. 9.
However, although variability in volume and surface area at the individual cell level was not counterbalanced across consecutive generations, global uniformization could be observed at the macroscopic level, as shown in Fig. 3k (main article) for the evolution of the total cellular volume.
Overall, given the weak values of R 2 obtained for the above six features of individual features in the different cell groups (R 2 was greater than 0.6 in only 20 cases out of 198 investigated pairs of distributions), we adopted the viewpoint that the various random variables were independent between and within cell groups. Although some of these correlation values may be ascribed to underlying deterministic factors responsible for cell-to-cell variability [20], our goal in the present work was to find parameters that could summarize relationships between random variables, hence the poor statistical significance of the linear regression test led us to a global assumption of independence. The following section examines the relations between probability laws in the different cell groups and how they are combined in the cell lineage.

Supplementary Note 2.3 Multi-level probabilistic model
In this stage, we used the parameterized description of the individual feature distributions (cell cycle length, mitosis time, volume, surface area) obtained above in each cell group, along with the assumption of independence among the different groups, to design a probabilistic model. As explained below, our model accurately predicted the final distributions of individual cell features compared to inter-individual differences within the cohort. It reproduced the embryo-level dynamics in each specimen (number of cells, total cellular volume, total cellular surface area), and also accounted for the propagation of intra-individual variability along the cell lineage.

Supplementary Note 2.3.1 Proliferation dynamics
In the first step toward designing a probabilistic model of the dynamics and evolution of the cell lineage, based on our data analysis, we stated the following relationships between the mitosis times and cell cycle lengths: where N (µ, σ ) represents a normal distribution of mean µ and variance σ 2 [16]. These equations implement the assumption of independence between X g and M g−1 . They are also realized independently in each branch of the cell lineage by the usual relationship m i = m j + x i (where cell j ∈ L g−1 is the mother of cell i ∈ L g ). From there, we could relate the parameters governing these distributions as follows: where g = (e, n , k) denotes the "initial" group, i.e. n < n is the first generation at which cells of type k can be observed in embryo e, and we use the shorthand notation g + 1 = (e, n + 1, k). The last equation means that the variability of mitosis times, represented by the variance of their Gaussian distribution, is based on the sum of variabilities of cycle lengths in each preceding cell group along the lineage tree. Similarly, the mean mitosis time is the sum of mean cell cycle lengths over the preceding groups.
Proof. Equation 24 can be derived from equation 23 by developing the recurrence relation between the random variables: Recall that the observation window of the developing embryo in our data is limited. It begins at the earliest at the 32-cell stage, i.e. at the 6th cell cycle (n = 6) where cells can be found in each of the four subpopulations, and ends at the latest at the 408-cell stage, i.e. cycle n = 10 for Mes and Mac cells, cycle n = 9 for LMic cells, and cycle n = 7 for SMic cells. Since the recurrence could not be traced all the way back to the origin, we provided a boundary condition for the initial distribution of mitosis times in the form of a Gaussian distribution, as in equation 23: M g ∼ N (µ M g , σ M g ). Then, the probability distribution of mitosis time M g was derived from the property that the density of the sum of two independent random variables is equal to the convolution of their individual densities: where f Y denotes the probability density of Y and * is the convolution product. In the case of Gaussian distributions, the convolution of two distributions is another Gaussian whose mean and variance are the sums of their respective components [21]: Since each distribution of cell cycle lengths is a Gaussian distribution, f X g = f (x | µ X g , σ X g ) and f M g = f (m | µ M g , σ M g ), we obtained equation 24 by combining equations 26 and 28: This result guarantees that the knowledge of the distributions of cell cycle lengths in each cell group was sufficient to characterize the distributions of mitosis times.

Supplementary Note 2.3.2 Cell volume and surface area dynamics
In a second step, we modelled the stochastic component of the volume variations as follows (Supplementary Note 2.2.2): where Y ∼ log N (µ, σ ) denotes a log-normal distribution, i.e. log(Y ) is normally distributed with mean µ and variance σ 2 [22,16]. As above, these equations implement the assumption of independence between A g and V g−1 and are realized independently by v i = v j a i . From there, we also obtained: by a proof similar to the previous section, with the difference that we were dealing here with log-normal distributions, hence the recurrence relationship was: log(A h ). 32 Finally, to obtain v i from v i , we used equation 15 with a deterministic function ω g (the "microdynamics" model of the variation of volume occurring between two consecutive mitoses) estimated from the empirical curves of Supplementary Fig. 7. Note that there is no simple probabilistic model for V g directly, only for V g . All the above equations are the same for surface areas, replacing V g , A g , v i and ω g with S g , B g , s i and φ g , respectively.

Supplementary Note 2.3.3 Model summary
The model is summarized in Supplementary Table 4, which describes the relationships between random variables at the level of cell groups, and Supplementary Table 5, which recapitulates how the probabilistic variables are instantiated in each branch of the lineage tree. Note that, just like there was no easy analytical model for the non-renormalized volume and surface area random variables, V g and S g , it is also difficult or impossible to give an analytical expression for the probability distributions of the embryo-level quantities (defined in Supplementary Note 2.1.2 and measured in Supplementary Note 2.1.3): number of cells N k (t), total cellular volume W k (t), and total cellular surface area Z k (t). Although these quantities are themselves random variables, there is no simple theoretical model for the macroscopic dynamics of cell populations. They could be assessed, however, by numerical simulation via the generation of virtual cell lineages, as explained in Supplementary Note 2.3.5 below.

Supplementary Note 2.3.4 Model evaluation
To assess the relevance of our analytical model, we evaluated the difference between, on the one hand, the empirical distributions (M g , V g , S g ) measured directly on the embryos and, on the other hand, their counterparts (M * g , V * g , S * g ) predicted by the model from the first group g and the sums of (X h , A h , B h ) over the cycles. This was done in each "final" group g = g = (e, n , k), i.e. where n is the last generation at which the random variable can be observed on a sufficient number of cells of type k in embryo e, and should also verify n > n . For 5 embryos, 3 variables and 4 types, there is a total of 60 potential final groups, but in practice only 48 satisfied this condition and were retained.
To this aim, we relied on the "Kullback-Leibler (KL) divergence", which measures the discrepancy between two probability densities p and q [23,24]: With Gaussian distributions p = N (x | µ p , σ 2 p ) and q = N (x | µ q , σ 2 q ), the KL divergence becomes: and can be symmetrized into a "distance" as follows: Therefore, we calculated ∆ KL (M g , M * g ) by plugging in the above formula the following parameter values: where g = g , and same thing for ∆ KL (V g ,V * g ) and ∆ KL (S g , S * g ). An absolute distance, however, cannot be interpreted if it is not compared to other typical interindividual distances within the cluster of groups. This is why we replaced ∆ KL with a normalized version ∆ KL calculated relatively to each embryo of the cohort, for example in the case of M g in specimen 1: where g(1) = (1, n, k) and g(e) = (e, n, k). These normalized distances between the model and the data were computed for division times, volumes and surface areas, in each cell population k and relatively to each embryo e that verified n > n , i.e. displayed at least one measurable cycle. Results are summarized in the histogram of Supplementary Fig. 10. Among the 48 eligible groups, 75% showed a close proximity ( ∆ KL < 0.5 in 36 cases) between measured and predicted parameter sets, 21% were still in good agreement (0.5 ≤ ∆ KL < 1.5 in 10 cases), and 4% fell far away ( ∆ KL > 3 for V and S in the Mes cells of embryo 5). Therefore, except for these two outliers, our probabilistic model accurately predicted the dynamics of sea urchin development.

Supplementary Note 2.3.5 Simulation of artificial cell lineages
Based on the above model, realistic virtual cell lineages were generated and their statistics compared to the real data. Starting from the 32-cell stage (16 Mes,8 Mac,4 LMic, 4 SMic cells), a cell lineage was computed by letting each cell divide into two and implementing the relationships between microscopic features (Supplementary Table 5 For each embryo, 300 realizations of the cell lineage were produced by varying the random generator's seed. As Supplementary Fig. 1 illustrates, all cell lineages shared the same binary tree topology (if taking into account only the bifurcation nodes), whereas branches could be of variable lengths (along the temporal axis) depending on mitosis times. Another difference between embryos resided in their respective windows of observation, i.e. the specific intervals [n , n ] during which there were enough observable cells of type k.
Comparison between simulated specimens and empirical values is shown in Supplementary Fig. 11 on three macrosopic quantities: total number of cells N(t), total cellular volume W (t), and total cellular surface area Z(t). Given the finite window of observation, the last cell cycles were often incomplete. As a consequence, the volume and surface area microdynamics, ω g and φ g (representing the deterministic "carrier waves" in equations 15,16) could not be correctly estimated (shaded areas).

Supplementary Note 2.4 Prototype
In sum, our probabilistic model at the coarse-grained level of cell groups L g provides an invariant framework relating individual cell features to embryo-level dynamics. The whole cohort can be described by the same system of equations (based on normal distributions), while cell statistics allows us to extract specific parameter values for each variable of interest in each group.
In this final step, we wanted to obtain a unified representation, or prototype, of the sea urchin blastula's normal development. To this aim, individual measurements were aggregated over the embryo index e using the geometrical concept of centroid of a set of points, which generalizes the notion of average in highdimensional spaces such as probability distributions. The mathematical field of "information geometry" is especially relevant to treat these spaces, in particular equip them with a distance such as the "Bregman divergence" or the KL divergence [23,25]. Here, we kept with our use of the latter to compute centroid distributions for each generation-type pair (n, k). For example, in the case of variable X: where g(0) = (0, n, k) and 0 denotes the prototype's index. Extending the assumption of normal (or lognormal) distributions to the prototype, we write X g(0) ∼ N (µ X g(0) , σ X g(0) ), and therefore the prototype's mean and variance parameters can be found by solving the following minimization problem: Numerical results for the usual six variables, based on their empirical parameters in each cell group L g , are plotted in Supplementary Fig. 3. In addition, the prototypical volume and surface area microdynamics, which are deterministic functions, are simply obtained by averaging the group-specific variables over the cohort (Supplementary Figs. 7 and 8, in bold): . 42 Finally, as it was done for each embryo, 300 realizations of the cell lineage were also produced for the prototype by varying the random generator's seed, then compared with the empirical values ( Fig. 3j-l).

Supplementary Note 3 Spatial modelling
Based on the above digital reconstruction and probabilistic model, the last stage of our study consisted of a spatially explicit computational model and simulation of the sea urchin embryo, based on the MecaGen platform [26], [17].

Equation of motion
Newton's second law states that the acceleration a of a body is proportional to the sum of the forces F applied to it. The proportionality factor is the inverse of the mass of the body, a quantity assumed to remain constant over time: a = F/m. In particle-based physics, for each particle characterized by a mass m, acceleration a, location X, velocity v, and radius R, forces are a function of the last three quantities, and Newton's law reads At the scale of the biochemical world, however, this classical formulation does not reflect the emergent and rather counterintuitive phenomena that cells exhibit individually and collectively. Cells are so small and their interactions so "sticky" that the physics of their motion is radically different from the one to which we are accustomed at our scale. The main difference resides in the disappearance of inertial forces [27] and its consequence is that in a "low Reynolds number" environment, applied forces produce a velocity, not an acceleration. Therefore, displacement is proportional to instantaneous forces and inversely proportional to a damping coefficient λ : Generalizing to a multicellular system, the motion of each cell indexed by i is governed by its interactions with all other cells j belonging to its neighbourhood denoted by N i , according to Supplementary Note 3.

Cell properties
Each cell is represented by a single cylindrical particle, i.e. its relaxed shape is considered isotropic in the surface plane. The cylinder axis U is orthogonal to the epithelial surface. The radius R of the cylinder delimits the distance occupied by each cell in the surface plane, and its half-height denoted by R ⊥ characterizes the depth of the tissue. Throughout this study, we assume R ⊥ homogeneous and constant across the embryo: for each cell i and at any time t, R ⊥ i (t) = R ⊥ 0 . Therefore, the relation between a cell's current radius and volume V i (t) reads: Supplementary Note 3.

Spatial neighbourhood
A cell sustains forces exerted by the neighbouring cells in contact with it. At any point in time, the embryo is described by a set of cylindrical particles i, each of them characterized by a location X i = (x i , y i , z i ), a radius R i , a half-height R ⊥ i and an axis U i . We need to infer from these intrinsic parameters the neighbourhood links between cells. To this goal, we followed a two-step strategy: first, a preselection of potential neighbours that obey certain metric criteria, then a refinement of this list according to topological features. After that, we specify how the cell axis is updated from the new list.
Metric selection Two cells of radius R i , R j and locations X i , X j are considered neighbours according to a metric criterion if and only if their relative distance is less than the sum of their radii R i + R j multiplied by a constant factor c max . Thus the metric neighbourhood of cell i is defined as follows: Denoting the distance between the two cells by r i j and the cutoff value by r max i j , it means that j is a metric neighbour of i if and only if r i j ≤ r max i j . The radius r max i j sets the maximum distance of the cell deformation in the surface plane, i.e. the spheroidal domain of space inside which cells can potentially interact. We assume that this maximum distance is reached when cells deform by a factor c max = 2.0.
Topological selection A purely metric neighbourhood is not viable in our multicellular context as it often leads to cell volumes collapsing during simulation when the adhesion between interacting cells is high. This is why we filter the set N m i via a topological criteria to obtain a reduced set N t i . To do this, a 3D Delaunay triangulation could be employed but it is computationally too slow because it involves heavy data structures composed of multiple tetrahedra. Our own solution is inspired by a similar neighbourhood model, the Gabriel graph. In 2D, whereas the Delaunay triangulation imposes that no node be found inside the circumcircle of any triangle, the Gabriel graph method requires that no node be found inside the circle whose diameter is a valid neighbourhood edge.
Here, we propose to generalize the Gabriel criterion to the 3D case directly and add a parameter to control the size of the circle: two cells (i, j) are considered neighbours if no other cell k of the embryo is located inside the sphere whose origin is the centre of segment i j, which is ( X i + X j )/2, and whose diameter is α gab X i − X j , where : To our knowledge, the term Gabriel graph has not been formally defined in 3D. However, when it happens to be mentioned [28], triangles are used instead of segments to define circumspheres. This makes sense, as in n dimensions, the Delaunay triangulation uses n-simplices for its basic units (triangles in 2D, tetrahedra in 3D) while the Gabriel graph uses (n − 1)-simplices (edges in 2D, triangles in 3D). Here, however, we will continue using edges in 3D while still referring to our method as a Gabriel graph. Note that in the 2D case, the Gabriel graph is a subgraph of the Delaunay triangulation. Every pair of neighbours in a Gabriel graph is also a pair of neighbours according to the Delaunay criteria.
Cell axis Once the topological list of neighbours is determined, the surface orientation of the monolayered sea urchin can be expressed at each cell location by averaging the outward normal vectors of the n triangles formed by n topological neighbours (Fig. 4a). These triangles are obtained by sorting the list of neighbours in counterclockwise order around the current axis U i . Let's denote by a k with k = 0...n − 1 the neighbour indices in the sorted list. The surrounding triangles are X i X a k X a k+1 with k = 0...n − 2 and X i X a n−1 X a 0 . The updated cell axis is calculated by summing the surrounding outward vector axes: The cell axis is then normalized:

Supplementary Note 3.1.4 Forces
The interaction force F i j exerted by cell j over cell i leads the swarm of cells toward an equilibrium state. It is the sum of two components: • an attraction-repulsion interaction force F i j , which maintains the integrity of the cell's volume and controls both the stiffness and the adhesion of the interaction in the direction of the surface plane, and • a planarity conservation interaction force F ⊥ i j , which maintains the planarity or "monolayeredness" of the epithelium; this force can be modulated via a planar rigidity coefficient.
Moreover, the damping coefficient λ i is proportional to the surface of the cylindrical cell: . Therefore, the equation of motion used in this model becomes: Attraction-repulsion force This force is derived from an attractive/repulsive interaction potential between two neighbouring cells i, j, and it is colinear with the neighbourhood vector u i j = ( X i − X j )/ X i − X j (here, from j to i). The rationale of this potential is to maximize the contact surface area between any two neighbours in the swarm. For every pair of neighbouring cells (i, j) we focus on the attraction-repulsion interaction potential E i j from which F i j is derived, and model it from three juxtaposed domains: • a repulsion domain (decreasing E) at distances shorter than a certain equilibrium distance defined by r eq i j = c eq i R i + c eq j R j , where in the most general case c eq i is a coefficient that can vary from cell to cell • an attraction domain (increasing E) at distances greater than r eq i j • a neutral domain (constant E) beyond the maximum limit of the interaction field r max i j = c max (R i + R j ) Coefficient c eq i is established explicitly in this study by considering that the equilibrium state of the swarm is reached when the cells form the densest packing in the 2D plane, i.e. the hexagonal lattice if disks are identical. The ratio of the surface occupied by the disks to the surface containing the disks is equal to the ratio of the surface of a disk to the surface of an hexagon (its Voronoi domain), i.e. π/(2 √ 3) 0.9068997.... Therefore, in the perfect uniform arrangement here, the equilibrium distance between any cells (i, j) of equal radii R i = R j = R is constant and reads r eq i j = 2c eq R, with c eq = We now generalize the notion of equilibrium distance to pairs of unequally sized neighbour cells. In this case, the common radius of equilibrium depends on individual contributions from each cell, which is simply proportional to their radius, keeping the same constant coefficient c eq as above: Given these boundary conditions, we consider here the simplest form of schematic mechanical model that can realize both the short-range repulsion and long-range attraction parts, namely: spring-like forces derived from an elastic potential. It means that forces F are a linear function, and potentials E a quadratic function of (r − r eq i j ). Mechanical interactions between neighbouring cells have been extensively studied in biophysics, yet due the multiprotein and polymeric nature of the cytoskeleton, together with the membrane's flexibility and adhesiveness, it remains an extremely elusive topic and difficult to summarize by a definitive set of equations. Roughly, one consensus hypothesis is that a major driving mechanism can be characterized by intercellular surface tension. The two main components of this tension are cellular adhesion, which tends to lower it, and cellular cortical tension, which acts antagonistically. The balance between these two components determines the ultimate behavior of the interaction, i.e. attractive or repulsive.
Another fundamental mechanism must also be taken into account: volume preservation in cells. As we made the assumption that cells are cylindrical and cell interactions in the plane can be decoupled into two individual, additive contributions, we assume here that the volume conservation principle plays a role in the repulsive part of the potential only. This allow us to introduce a correction in the form of a discontinuity in the force derivative at the distance of equilibrium r eq i j . We will also globally refer as "adhesion" to the cumulated effect of true surface adhesion and its antagonist, cortical tension.
Accordingly, to modulate the intensity of adhesion, the expression of F is split into two parts, one below and one above the equilibrium distance r eq i j , corresponding to two different stiffness coefficients w. A low (resp. high) adhesion coefficient w will induce a weak (resp. strong) attraction: Cells i, j start to adhere only when their relative distance becomes less than r max i j . As the intensity of the attraction can not be suddenly maximal at that point, a half-bell shape between r eq i j and r max i j is better suited to model a more realistic cellular interaction. We obtain this shape by multiplying the linear force F ,lin i j by a term A(r i j , R i , R j ) scaling as the area of the surface between i and j. We decided to use a monomial of degree 2 to make it scale like a surface. If the distance between two double-size neighbouring cells is doubled, the area of the surface of contact should quadruple: In the whole interval [0, r max i j ], the depth of the well is controlled by a repulsion coefficient w rep (until r eq i j ) and an adhesion coefficient w adh (beyond r eq i j ; see Fig. 4b). Thus the final expression of the nonlinear version of the relaxation force that takes into account the contact surface area reads: Finally, to compare our custom attraction/repulsion force F i j with classical potentials, we show their superimposed plots in Supplementary Fig. 13.
Planar rigidity force Attraction-repulsion forces are not sufficient to ensure that the spatial configuration of the embryo will remain locally "planar", i.e. monolayered, during simulation, spherical. Indeed, such forces are colinear to the neighbourhood vectors u i j , causing the epithelium to eventually agglomerate into a multilayered tissue. To model a monolayered epithelial sheet, we introduce a specific planar rigidity force ( Supplementary Fig. 14), which aims at maintaining each neighbour cell in the lateral region, orthogonally to the normal vector. Between two neighbour cells i and j, this force is aligned with the bisector of their normal vectors: n i j = ( U i + U j )/ U i + U j , and is proportional to their distance, the dot product between u i j and n i j , and a planar rigidity coefficient, k rig : F ⊥ i j = −k rig r i j ( u i j . n i j ) n i j 57 Distance between distributions Given two empirical probability distributions p and q, considered on the bins X, such that ∑ x∈X p(x) = 1, we define the following measure of similarity [29]: which is always positive, equal to 1 when the two distributions do not overlap and 0 when they are equal.
Histogram resampling To overcome the problem of having different time samplings in the various embryos, we used a common time interval across all specimens and for the spatial simulation, ranging from 240 min post fertilization (mpf) to 702 mpf, with a sample every 3 min. We interpolated the value of the histograms used for the comparison at each time point of this time interval.
Metric We used the distribution over i of the numbers of neighbours b k→l i that a cell i shares at a border to compare simulations and digitally reconstruct empirical embryos using the distance defined above, averaged over the different time points, populations, and embryos. The comparison protocol is described in algorithm 1, where h k→l e (t) and h k→l s (t) denote the distributions of b k→l i (t) in an embryo e and a simulation s respectively. It leads to a metric D s whose value is comprised between 0 (best fit) and 1 (worst fit).

.2 Objective functions
In addition to the metric introduced above, we measured intrinsic features of the simulated embryos via objective functions. These functions are defined for a given time step and subsequently averaged over the entire time interval.
Sphericity objective function S s An important criteria to evaluate the model is the shape of the simulated embryos. For that purpose, we designed S s to be the ratio of the standard deviation over the average of the set of distances between the cell positions X i and the embryo centre X 0 = X i i . This function has a low value for a small standard deviation and a large average, and is null in the case of a spherical embryo: Planarity objective function P s For each cell, this function characterizes the unicity of the epithelium layer by measuring the ratio of cells belonging to the surrounding apicobasal neighbourhood. The neighbourhood N i is split into two complementary sets, one covering the lateral domain N i , the other the apical domain N ⊥ i : where η is a threshold value controlling the partition between the two sets, here η = cos(π/4). The planarity objective function is then written: For a monolayered epithelium, P s = 0, otherwise its value increases as cells agglomerate into a 3D tissue.

Supplementary Note 3.2.3 Validation in parameter space
The purpose of this exploration study is to identify the parameter sets that produce a realistic spatial folding of the sea urchin embryo during its development. The validation criteria are the sphericity of the global embryo shape, measured by S s , the maintenance of the monolayered epithelium, measured by P s , and the similarity of the interpopulation border shapes, measured by D s . Each criterion is explored separately and its optimal parametric region identified (main article, Fig. 4e-g). The common parameter space here is fourdimensional: it comprises the planar rigidity coefficient k rig , the Gabriel criteria coefficient α gab and the two adhesion coefficients specifically controlling the attractive part of the relaxation force F i j : w adh,o between homotypic cells and w adh,e between heterotypic cells. The other free parameters of the model (namely w rep , λ 0 , and c max ) are set to constant values (Supplementary Table 6). In Fig. 4 and Supplementary Videos 4-7, k rig and α gab are also set to constant values, while Fig. 4c-g shows only a restricted two-dimensional region of the explored domain in the (w adh,o , w adh,e ) plane. Finally, the time increment ∆t between two simulated time steps is 6 seconds. Note concerning the force amplitude coefficients w adh,e , w adh,o , w rep and k rig : these coefficients are all divided by the damping coefficient λ 0 in the master equation of motion. Thus, simulations will be strictly equivalent if all these parameters are multiplied by the same factor.
For the spatial simulation, an initial spatial state of the cells was needed in order to begin the simulation at the 32-cell stage: we chose the measured initial state of one of the five embryos of the cohort.

Supplementary Software 1 Supplementary Software 1
The compressed folder "prototype.zip" contains all the MATLAB files necessary to perform and visualize the complete statistical analysis of the cell lineage and to generate the prototype. It has been tested with Matlab R2011b and Matlab R2015a. A tutorial is provided in the file README.pdf contained in the folder.