Numerical integration

We solved the governing equation (Eq 8) using the Euler-Maruyama method [37] to obtain the motion of all cells. Simulations were run using a timestep Δt = 10⁻³τ_M, much smaller than all relevant timescales in the system, namely τ_M and τ_D. The discretized version of Eq 8 that we integrated numerically reads

Simulations were initialized with the undivided egg (first cell) at the origin for unconfined simulations or to have Gaussian distributed initial positions with variance b/10 around the origin for confined simulations. Simulations either ended at the timestep before cells at the 4-cell stage would divide again in non-equilibrium simulations, when 4 cells reached a tetrahedron in simulations searching for the equilibrium relaxation time, or after 8000τ_M to determine equilibrium configurations in embryos with confining shells.

Cell-cell interaction potential

The cell-cell interaction potential must have a repulsive core (accounting for volume exclusion) and an attractive part that accounts for cell-cell adhesion. In addition, since cells have a finite size, the interaction potential must vanish for distances r larger than the cell size. We built a potential U_c(r) = S(r)U_p(r) with these characteristics by making it the product of a support function S(r) = 1/(1 + f(r_ij)), with f(r_ij)≡exp[(r_ij − (R_i + R_j))/a] (R_k being the radius of cell k and a being a length scale characterizing the transition from 1 to 0 of the support function), and a Lennard-Jones like potential U_p(r) = A/r^α − B/r^β. Imposing that the minimum of the potential has a magnitude −U₀ (U₀ being the scale of the potential) and is located at $$ r=r_{ij}^{*}$$ ($$ r_{ij}^{*}$$ being the equilibrium distance between cells i and j; this is equivalent to normalizing lengths with $$ r_{ij}^{*}$$), we obtain the values of A and B and the potential U_c(r) reads

Confining shell

The shape of an axisymmetric ellipsoidal shell, and the associated ellipsoidal level set, are given implicitly by

Since cells cannot penetrate the shell, the confining potential must diverge at the positions where the shell is located (c = 1). Moreover, the potential must vanish when the cell can no longer be in contact with the shell, which occurs when a cell is located at a distance larger than R from the shell. With this in mind, we define the confining potential of a repulsive shell U_s (Fig 1F) as

In the case of a sticky shell (Fig 1F), we use the same shell-cell interaction potential as the interaction potential between two cells, albeit with different adhesion strength. In this case, the equilibrium distance is changed to $$ r_i^{*}$$ instead of $$ r_i^{*}+r_j^{*}$$ since there is only one cell interacting with the shell. In these conditions, the interaction potential for a sticky shell reads

Topology inference

Because we simulated particles in a confined volume, it was necessary to move beyond a simple distance metric to determine if two cells were neighbors. We used a 3D Voronoi partitioning as an extra constraint in addition to distances. The package we used, Voro++ [38], determines the 3D Voronoi polytope around each cell by starting with a large 3D volume and then cutting it using the midplanes from the current cell to each of its neighbors. In this work, each cell starts with a dodecahedral volume that surrounds an inscribed sphere with the cell radius R. When a pair of cells are close enough, their dodecahedral volumes are cut by the weighted midplane between them (adjusted from the midpoint by their respective radii) (Fig 1G). Each new cut face of the Voronoi polytope is identified with a neighboring cell allowing all neighbors to be identified. Two cells are defined to be in contact when they are within a distance R_i + R_j (R_i and R_j being the radii of cells i and j) of each other and their Voronoi polytopes share a face. An adjacency graph is created by defining each cell as a node and adding edges between each pair of cells found to be in contact. We determine the topology of each arrangement by checking if the adjacency graph is isomorphic to a reference adjacency graph for each type of topology (square, diamond, tetrahedron, T-shape, line) [22, 39].

Cell divisions

Cell divisions occur at well defined intervals of time τ_D (Fig 1H), with all cells dividing simultaneously. The new daughter cells were placed at a daughter radius away from the mother cell in a single timestep. In the case of random divisions, the daughter cell divides in a random direction from the mother cell (with a check to ensure that the daughter is not placed within a distance that would cause it to substantially overlap with an already existing cell). For the case of ordered divisions, the egg first divides in the x direction, then both daughter cells divide in the y direction.

Since at the developmental stages considered herein (4-cell stage) total volume is conserved, daughter cells have half the volume of the mother cell, as described in the main text. However, this simulation framework allows also simulations where the volume changes. It is also possible to simulate changes in cell size following changes in physical quantities that evolve during the simulations. This can be done by providing a rule to decide how cell size evolves over time, how it changes upon division or establishing how some physical parameters in the simulations affect cell size. Then, the cell size is adapted for each cell and at each time step by changing the interaction potential cut-off (cell size) accordingly.

Volume adjustment

The overlap between blastomeres was determined by defining a sphere with radius R_i around cell i and then calculating its overlap volume V_o with neighboring spheres. This overlapping volume then added to cell i, making it larger. In particular, the radius of cell i is modified from R_i before the correction to $$ R_i^{\prime }$$ after it, with $$ 4\pi {\left(R_i^{\prime }\right)}^3/3=4\pi R_i^3/3+V_o$$ (Fig 1E). This adjustment is performed at every timestep, and the Voronoi dodecahedron is also scaled to surround a sphere of radius $$ R_i^{\prime }$$ after volume correction, at each timestep too.

Angular mean squared displacement

The angular mean squared displacement is defined as

Theoretical description

In order to simulate the 3D dynamics of blastomeres, accounting for the interactions between them, their divisions as well as embryo confinement, we use a minimal representation and describe each blastomere (cell) as a particle. In this particle-based representation, cells interact with each other through an interaction potential U(r_ij) that effectively accounts for the mechanical interactions between cells (adhesion, etc. [40]), with $$ r_{ij}=|{\overrightarrow{r}}_i-{\overrightarrow{r}}_j|$$ being the distance between two given cells located at positions $$ {\overrightarrow{r}}_i$$ and $$ {\overrightarrow{r}}_j$$. Cell-cell adhesion is represented by an attractive range in the potential, whereas a repulsive region ensures that cells do not interpenetrate when they become too close to each other (Fig 1C). To account for cell size in this particle description, we include a sharp cut-off of the potential at a distance R (Methods), which corresponds to the radius of an isolated blastomere. The balance of attractive and repulsive forces between two blastomeres occurs when they are separated by a distance 2r*. The ratio between this equilibrium distance between blastomeres and the blastomere size 2R corresponds to r*/R = cos θ, with θ being the contact angle θ between cells (Fig 1B). Since the contact angle is an easily measurable quantity that informs about the relative strength of adhesion and cortical tension [15, 30, 41] (Fig 1D), we use θ as control parameter instead of r*. Moreover, although it is not possible to enforce exact volume conservation in a particle-based description, we perform leading order corrections upon cell contact (Fig 1E; Methods); we have checked that the volume corrections are small and we have tested that our results do not qualitatively depend on them.

At the spatial and temporal scales of embryo development, the system is overdamped and inertia can be safely neglected [42]. In this case, force balance (momentum conservation) for a given blastomere reads

In order to properly determine what cells are in contact and can therefore apply forces on each other, we use Voronoi tessellation (Methods; Fig 1G). Previous particle-based simulations used distance-based metrics to determine the neighbors of each cell. However, distance-based metrics can give erroneous results for both cell-cell contacts and dynamics in the presence of confining shells. This is because when cells are highly confined, the distance between next-nearest neighbors can be smaller than the interaction potential range, thereby erroneously considering the forces of cells that are not in direct contact. Voronoi tessellation overcomes this problem and enables proper determination of cell-cell contact topology at each timestep of the simulation (Methods).

Since shells of many species have spherical or ellipsoidal shapes, we consider only these cases in what follows. We approximate the shell surrounding the embryo by an axisymmetric ellipsoid with major and minor axes a and b, respectively, with volume $$ V_s=\frac{4\pi }{3}a{b}^2$$ and aspect ratio a/b. Since blastomeres cannot penetrate the shell, we use confining potential forms that diverge at the shell boundary (Methods; Fig 1F). Moreover, the confining potential vanishes for distances larger than R from the shell, as these distances are not within the reach of cells.

Beyond physical interactions among cells and with the confining shell, blastomeres in early embryos divide at regular intervals, with a time τ_D between division events (Fig 1H). We simulate division events accounting for the change in volume of the cells upon division (Methods). Since the volume of the daughter cells is half of cell volume before division, the cell radius R changes after each division cycle to R_n = R₀/2^n/3, where R₀ is the radius of the initial egg (and $$ V_c=4\pi {R_0}^3/3$$ is the initial egg volume) and n is the number of divisions that have occurred (Fig 1H). Finally, in order to study the role of division rules, we control the spatial direction along which cell division occurs, which corresponds to the direction perpendicular to the mitotic plane. While division rules are known to exist [33, 36, 43], the specific rules and the parameters that control them are still under debate, especially for different species. As a consequence, to study the role of division rules, we consider two limiting cases: (1) Ordered divisions, in which we impose representative division rules at early developmental stages (division axis is perpendicular to the division axis in the two previous division cycles; for first division, perpendicular to previous division), and (2) Random divisions, in which there are no division rules and we randomize the direction of cell divisions for each cell and division cycle (Fig 1I).

Normalizing all lengths by the initial egg radius R₀, all forces with U₀/R₀ and time with the mechanical relaxation time τ_M, which is given by $$ {\tau }_M\equiv \mu R_0^2/U_0$$ and represents the characteristic timescale over which mechanical disturbances relax to equilibrium, we obtain the relevant dimensionless parameters in the problem (Table 1).

Table 1

Definition of the relevant dimensionless parameters in the problem.

Dimensionless Parameters
Parameter	Description
θ	Contact angle: Specifies the relative strength of adhesion between cells to the cortical tension of each cell.
τD/τM	Ratio of characteristic timescales: Time between divisions relative to the mechanical relaxation time.
Vs/Vc	Shell to cells volume ratio: Ratio of confining shell volume and the total volume of cells.
$$	Cell-shell to cell-cell adhesion strength: Ratio of the attractive energy scale for cell-shell interactions $$ and cell-cell interactions U0.
a/b	Shell aspect ratio: Ratio of ellipsoid major axis length a to minor axis length b.

In what follows, we simulate the stochastic movements of the multiple interacting blastomeres using Langevin dynamics (Eq 8; Methods) in different conditions.

Unconfined cellular packings

To understand the packing configurations at the 4-cell stage that arise from the system dynamics, we first simulate the cellular dynamics upon divisions in the absence of embryo confinement (Fig 2A). We define the 4-cell stage packing configurations (Fig 2B) as the cellular arrangement just before cells at the 4-cell stage undergo the next division cycle. At equilibrium, the minimal energy configuration of 4 blastomeres in contact with each other is a tetrahedron, as already established both theoretically and experimentally for clusters of four particles with attractive interactions [21–23, 26]. However, if blastomeres divide much faster than the time required for cells to undergo mechanical relaxation (τ_D/τ_M < 1), cells do not have time to reach mechanical equilibrium in between divisions and the cellular packings do not coincide with the equilibrium packing configuration, as expected. When divisions are fast compared to mechanical relaxation (τ_D/τ_M = 0.1; Fig 2C) and cells divide following ordered divisions, either squared or diamond configurations are observed, with diamond configurations being more prevalent as the cell contact angle increases. If cells divide in random directions, other packing configurations appear and depend on the contact angle, but both squared and tetrahedral packings are missing. In contrast to the case where division occur fast, if blastomeres take much longer than the mechanical relaxation time to divide, equilibrium packings are expected because cells should have enough time to reach mechanical equilibrium between divisions. However, our simulations show that for τ_D ≫ τ_M (specifically, τ_D/τ_M = 10; Fig 2C), the expected tetrahedral configurations are not observed for ordered divisions (only diamond configurations are observed) and barely observed for random divisions.

Fig 2

Packing configurations and dynamics of 4-cell stage unconfined embryos.

A, Example of time evolution of cellular packings, showing how cell volume decreases by half upon division. B, Definition of possible topological arrangements at the 4-cell stage. C, Frequency of packing arrangements for slow and fast divisions (τ_D/τ_M = 0.1, 10, respectively), different contact angles θ and both ordered and random division rules (n = 10³ simulation runs for each parameter set). D, Histogram of the time to reach the equilibrium tetrahedral configuration (n = 10⁵ simulation runs for each condition). E, Frequency of packing configurations as the system relaxes to equilibrium (t_E being the time to reach equilibrium) for both ordered (left) and random (right) divisions (n = 10² simulation runs for each condition). F, Top view (with cells) and cross section (without cells) showing the equipotential surface (orange tones) caused by three fixed cells on a fourth cell. G, Angular mean squared displacement (MSD) of a cell moving in the potential generated by three cells fixed in a triangle (n = 8 simulation runs), showing its diffusive nature (fit, gray line).

To understand why the expected tetrahedral configurations are not observed, we characterized the time necessary to reach the tetrahedral equilibrium configurations at the 4-cell stage by preventing the next division round. For both ordered and random divisions, we find that cells require times three orders of magnitude longer than τ_M to reach equilibrium (Fig 2D). By monitoring the time evolution of packing configurations at the 4-cell stage as the system relaxes to equilibrium, we found that for ordered divisions cells are always in a diamond configuration before reaching the tetrahedral configuration, whereas for random divisions cells mostly evolve towards the diamond configuration and stay in that configuration until reaching the equilibrium packing (Fig 2E). This results indicate that it takes a long time for the cluster to leave the diamond configuration, suggesting that the transition between diamond and tetrahedral configurations may involve the rotational diffusion of a blastomere. Indeed, the equipotential surface that a blastomere perceives when the three other blastomeres form a triangle indicates that in order to transit from diamond to tetrahedral configurations, one blastomere needs to traverse a flat region of the potential (Fig 2F). The angular mean squared displacement of the movements of such blastomere scales linearly with time (Fig 2G), showing that the transition between diamond and tetrahedral configurations occurs via rotational diffusion and explaining the long times necessary to reach the tetrahedral state.

The consequence of this floppy mode in the dynamics of the blastomeres is that it imposes an extraordinarily long time for the system to reach mechanical equilibrium, effectively leading to a degeneracy in the packing configurations at the 4-cell stage for normal division times, with degenerate packings being strongly dependent of division rules and adhesion strength between blastomeres. Such large degeneracy in the packing configurations and their strong dependency of multiple parameters is not an adequate strategy to robustly specify cellular packings. Since tetrahedral packings at the 4-cell stage are observed in embryos of several species, our results suggest that another mechanism must contribute to establishing the observed tetrahedral packings, as otherwise cell divisions would need to be extraordinarily slow (τ_D∼ 24h) to allow blastomeres to reach equilibrium between divisions.

Cellular packings under spherical confinement

Many embryos displaying tetrahedral packings at the 4-cell stage feature a spherical confining shell [8, 30]. To understand the potential role of embryo confinement on cellular packings at the 4-cell stage, we simulate the dynamics of blastomeres in the presence of (repulsive) spherical confinement (Fig 1F). We focus on long division times (τ_D/τ_M = 10) as this was the limit in which tetrahedral packings were expected, but shown above to be missing due to the long times associated with rotational diffusion of the blastomeres. If the confining shell has a very large volume compared to the total volume of the cells (V_s/V_c ≫ 1; Fig 3A), the situation is similar to the unconstrained embryo. As the volume of the confining shell is decreased, the 4-cell stage packings start to change because cells start interacting with the repulsive shell. Finally, when the volume of the confining shell is comparable to the volume occupied by the cells (V_s/V_c ≃ 2), only tetrahedral configurations are observed (Fig 3A). In this case, we find that blastomeres robustly reach the tetrahedral packing at the 4-cell stage five orders of magnitude faster than in the absence of confinement (relaxation time < τ_M; Fig 3B) and irrespective of their contact angle or division rules (Fig 3C). These results indicate that the presence of spherical confinement removes the degeneracy of cellular packings and quickly imposes a tetrahedral blastomere configuration, overriding division rules.

Fig 3

Spherically confined cellular packings.

A, Frequency of different cell arrangements for embryos with spherical repulsive confinement as the ratio of the shell volume and total cells volume varies (random divisions; τ_D/τ_M = 10; n = 10⁴ simulation runs for each parameter set). B, Histogram of the time to reach the equilibrium tetrahedral configuration (n = 10³ simulation runs). C, Frequency of packing arrangements for embryos confined in a spherical repulsive shell (a/b = 1; V_s/V_c = 1.52), and both slow and fast divisions (τ_D/τ_M = 0.1, 10, respectively), different contact angles θ and both ordered and random division rules (n = 500 simulation runs for each parameter set).

Cellular packings under ellipsoidal confinement

While embryos of several species have spherical confining shells, other shell geometries are observed across species. Different nematode species display elongated axisymmetric shells of varying aspect ratios [8, 9] that can be approximated by an axisymmetric ellipsoidal geometry. Previous works have shown that the shape of the confining shell is important for cellular arrangements in nematodes [30], for which the 4-cell stage packing arrangements are critical for the survival of the embryo, as improper cell contacts lead to fatal developmental defects [44]. To understand the role of varying ellipsoidal shell geometries on the cellular packings across nematode species, we systematically studied the effect of confining shell volume and aspect ratio on the blastomere packing configurations at the 4-cell stage.

Similarly to spherical shell geometries, for all simulated aspect ratios of ellipsoidal shells (a/b = 1, …, 3) we find that when the confinement is negligible (V_s/V_c ≫ 1) and divisions are randomly oriented, packing configurations are strongly degenerate even if the time between divisions is much longer than the mechanical relaxation timescale (τ_D/τ_M ≫ 1; Fig 4). In this case, different packing configurations have similar relative frequencies, albeit with the diamond configuration always being predominant. Essentially, if the volume of the shell is sufficiently large compared to the total volume of the blastomeres (V_s/V_c ≫ 1), the observed packings at the 4-cell stage are similar to unconfined embryos (Figs 2C and 4), as expected. As the volume of the shell is decreased and the blastomeres start to feel the physical confinement, the relative frequencies of 4-cell stage packings start to change, removing some degeneracy in packing configurations, in a manner that depends on the shell aspect ratio. When the volume of the confining shell is comparable to the volume occupied by the cells (V_s ≃ 2V_c), the degeneracy in 4-cell stage packings largely disappears and different, but unique, packings exist for different aspect ratios (Fig 4).

Fig 4

Cell arrangements for repulsive ellipsoidal confinement.

Frequency of different cell arrangements for embryos with ellipsoidal repulsive confinement of varying aspect ratio (a/b) as the ratio of the shell volume and total cells volume varies (random divisions; τ_D/τ_M = 10; n = 10⁴ simulation runs for each parameter set). Snapshots of non-equilibrium packing configurations for large shell volumes (V_s/V_c = 6) are shown on the right. Equilibrium (t = 8000τ_M) packing configurations for each aspect ratio and V_s/V_c = 1.52 are shown on the left (n = 200 simulation runs for each parameter set).

For some shell geometries we observe that even under confinement, two different packing geometries are possible, albeit with different relative frequencies (Fig 4). For aspect ratio of 1 (spherical limit), only tetrahedral packings are obtained, as described above and observed for nematode species with spherical shells [8, 30]. As the aspect ratio increases the relative frequency of the tetrahedral packing diminishes and the frequency of diamond packings increases, with only diamond configurations observed between a/b ∼ 1.5 and a/b ∼ 2.4. Increasing the aspect ratio above a/b ∼ 2.4 leads to the coexistence of diamond and linear configurations (Fig 4; a/b ∼ 2.5). For aspect ratios of a/b ∼ 3 and above, the only observed configuration is linear.

To check if the packing configurations obtained in confined embryos (V_s/V_c = 1.52; τ_D/τ_R = 10) correspond to the actual equilibrium packings for each specific shell geometry, as was the case for the spherical shell, we simulate the dynamics of blastomeres, preventing further cells divisions at the 4-cell stage and letting the system reach equilibrium. We find that for each value of the aspect ratio, the packing configurations observed for τ_D/τ_R = 10 (with random divisions) were the actual equilibrium configurations of the blastomeres at the same confining volume (V_s/V_c = 1.52) and aspect ratio. This indicates that the confining shell eliminates the degeneracy in packings and selects the equilibrium packing configurations for a given shell geometry. These results indicate that the geometry of the confining shell alone can determine the 4-cell stage blastomere arrangements regardless of the specific division rules, providing a robust mechanism to remove packing degeneracy and select the proper cellular packing.

Cellular packings in sticky shells

So far, we have only considered shells that confined the blastomeres by generating a repulsion force upon contact (repulsive shell, Fig 1F). However, in some species, the blastomeres can adhere to the confining shell, likely affecting blastomere packing configurations. In the case of sea urchin embryos (echinoderms), there is evidence of strong adhesion to the hyaline layer surrounding the blastomeres [10, 45–47] and the 4-cell stage blastomere packing configuration is a square [1], a configuration never observed in the cases described above. In the case of sea urchins, the geometry of the hyaline layer (confining shell) that surrounds the blastomeres is not exactly spherical and changes slightly over time. However, for the sake of simplicity, here we consider a spherical sticky confining shell (Fig 1F). Since echinoderms have stereotypical division rules (dividing perpendicular to the two previous divisions at early stages) we study the effect of shell-blastomere adhesion strength and confining volume on the 4-cell stage packing configurations for ordered divisions (Fig 1I).

Fig 5

Cell packing arrangements in a sticky confining shell.

A-C, Frequency of packing arrangements (A, diamond; B, square; C, tetrahedron) for varying values of the shell volume to cells volume ratio, V_s/V_c, and cell-shell to cell-cell adhesion strength, $$ U_s^0/U_0$$ (n = 5000 simulation runs for each parameter set).

For large shell volume compared to the total blastomere volume (V_s/V_c ≫ 1), the only observed configuration with finite blastomere-shell adhesion is the diamond configuration, with all blastomeres adhered to the shell (Fig 5A). When the confining volume becomes comparable to the blastomeres volume (V_s/V_c ≃ 1 − 1.5), the diamond configuration is suppressed and the tetrahedral and square packing configurations coexist at different frequencies depending on the relative strength of cell-cell adhesion and cell-shell adhesion (Fig 5B and 5C). In this case, if adhesion to the shell is very low, only tetrahedral configurations are observed, as expected in the limit of negligible shell adhesion (repulsive shell). When the adhesion to the shell dominates over cell-cell adhesion, we find that square and tetrahedral packing configurations are observed approximately 40% and 60% of the time, respectively.

The existence of strong adhesion to the confining shell introduces square packing configurations that were not observed in purely repulsive confining shells. However, the blastomere adhesion to the shell cannot be the only reason why 4-cell stage embryos of echinoderm species are square, since there is still a strong probability of tetrahedral packings even in the presence of adhesion.