Effective mechanical potential of cell–cell interaction in tissues harboring cavity and in cell sheet toward morphogenesis

Measuring mechanical forces of cell–cell interactions is important for studying morphogenesis in multicellular organisms. We previously reported an image-based statistical method for inferring effective mechanical potentials of pairwise cell–cell interactions by fitting cell tracking data with a theoretical model. However, whether this method is applicable to tissues with non-cellular components such as cavities remains elusive. Here we evaluated the applicability of the method to cavity-harboring tissues. Using synthetic data generated by simulations, we found that the effect of expanding cavities was added to the pregiven potentials used in the simulations, resulting in the inferred effective potentials having an additional repulsive component derived from the expanding cavities. Interestingly, simulations by using the effective potentials reproduced the cavity-harboring structures. Then, we applied our method to the mouse blastocysts, and found that the inferred effective potentials can reproduce the cavity-harboring structures. Pairwise potentials with additional repulsive components were also detected in two-dimensional cell sheets, by which curved sheets including tubes and cups were simulated. We conclude that our inference method is applicable to tissues harboring cavities and cell sheets, and the resultant effective potentials are useful to simulate the morphologies.

1 Introduction

In multicellular living systems, various three-dimensional morphologies are emerged, which is thought to be primarily dependent on the mechanical properties of the constituent cells (Steinberg, 2007; Maître et al., 2012; Maître et al., 2015; Guillot and Lecuit, 2013; Fletcher et al., 2014; Bazellières et al., 2015; Petridou and Heisenberg, 2019). Measurement of cellular mechanical parameters is inevitable for understanding morphogenetic events. On the other hand, morphogenetic events are also influenced by non-cellular structures such as cavities, extracellular matrix, etc (Style et al., 2014; Wang et al., 2018). Therefore, it is challenging to identify all components contributing to morphogenetic events, and then, to measure their mechanical forces especially under three-dimensional situations.

Image-based methods for inferring mechanical forces in tissues have been developed, where mathematical models are assumed which are fitted to microscopic images (Miyoshi et al., 2006; Brodland et al., 2010; Ishihara and Sugimura, 2012; Koyama et al., 2012; Koyama et al., 2023; Merkel and Manning, 2017; Kondo et al., 2018). These methods include inference of cell–cell junction tensions in two-dimensional epithelia, cell surface tensions or bending rigidities, mechanical potential of pairwise cell–cell interactions, etc. Basically, these methods do not assume non-cellular structures in their mathematical models, and, it is unknown whether these methods are applicable to tissues with non-cellular components.

We previously reported a method for inferring mechanical potential of pairwise cell–cell interactions (Koyama et al., 2023). The pairwise potentials are composed of attractive and repulsive forces of cell–cell interactions. The attractive forces are at least partially derived from cell–cell adhesion forces mediated by cell adhesion molecules (Figures 1A,B), whereas the repulsive forces are from excluded volume effect of cells (Drasdo et al., 1995; Forgacs and Newman, 2005; Steinberg, 2007; Maître et al., 2012; Guillot and Lecuit, 2013; Yong et al., 2015; Merkel and Manning, 2017). The summation of these forces would be dependent on cell–cell distance (Figure 1C), by which the pairwise potential energy is calculated. We previously showed that, by applying our method to the mouse and C. elegans early embryos which did not still harbor clear cavities or thick extracellular matrix, we obtained mechanical potential of pairwise interactions (Koyama et al., 2023). The profiles of the pairwise potentials were quantitatively different among the tissues, by which the morphologic features of the corresponding tissues were reproduced.

Figure 1

Figure 1. Overview of strategy for inferring the effective potential of cell–cell interactions. (A) Microscopic forces contributing to cell–cell interactions are exemplified. (B) The microscopic forces are approximated as attractive/repulsive forces where cells are described as particles with the mean size. (C) Both the repulsive and attractive forces are provided as cell–cell distance-dependent functions, and the summation (black broken line) is the distance–force curve of the two particles. The relationship between the curve and the mean diameters of the cells is shown. (D) A cavity-harboring tissue is illustrated with liquid flux into the cavity.

In the present study, we evaluated whether our inference method is applicable to cavity-harboring tissues such as the mouse blastocysts (Figure 1D). First, we theoretically analyzed the effect of liquid cavities by using synthetic data which were generated by simulations under pregiven cell–cell forces and cavities. In the case where the cavities are expanding, the inferred effective pairwise forces of cell–cell interactions became significantly different from the pregiven forces, where an additional repulsive component was detected. But we showed that the repulsive component was derived from the expanding cavities, and that the inferred effective pairwise forces were the summation of both the pregiven forces and the effect of the expanding cavities. Therefore, simulations using inferred pairwise forces successfully reproduced the cavity-harboring structures. Application of our method to the mouse blastocysts also resulted in the same conclusion. Because the effect of cavities on the inferred effective forces is additive, we may evaluate differentials of cell–cell interaction forces between two tissues such as wild-type and mutant ones, if one can assume that the effect of the cavities is equivalent between the two tissues. Finally, we systematically analyzed what kind of profiles of the effective pairwise potentials are critical for forming cavity-harboring structures, which were presented as phase diagrams, and expanded the analyses to curved sheets.

2 Theory; force inference procedure with simulation model

2.1 Overview

In our inference method, a cell-particle model is assumed where attractive/repulsive forces of cell–cell interaction is considered (Figures 1B,C), and the mathematical model is fitted to nuclear tracking data (i.e., xyz-coordinates of each nucleus with the time series) as described previously (Koyama et al., 2023). Both the modeled particles and the nuclei were assumed as point particles. Briefly, the fitting procedure is based on a minimization of the differences between the xyz-coordinates of point particles in the mathematical model and in the nuclear tracking data, with an additional cost function for cell–cell distance-dependent weight. The fitting parameters are the force values for each pair of two interacting cells for each time frame, and then, by binned averaging the inferred forces along cell–cell distances (D), a distance–force curve [F_C–C(D)] is calculated. In the model to be fitted, we do not explicitly consider non-cellular structures exerting forces on cells. It can be argued that, if tissues harbor non-cellular structures, the model should assume the non-cellular structures. But, due to high degrees of freedom in such models, the fitting problems would become difficult to be solved in general. We intended to analyze whether valuable information can be obtained when our inference method is applied to tissues with cavities (Figure 1D). To this aim, we first generated simulation data considering both cell–cell interaction forces and cavities, to which we applied our inference method.

2.2 Simulation model

In our simulations based on particulate cells, the cell–cell interaction forces were set to be provided by the Lennard-Jones (LJ) potential (i.e., the pregiven forces) in a similar manner to our previous work (Koyama et al., 2023). Briefly, due to the low Reynolds number in the case of cellular-level phenomena, the inertial force is negligible. Viscous drag force or frictional forces are provided from surrounding liquid medium or tissues (Forgacs and Newman, 2005; Basan et al., 2011; Fletcher et al., 2014; Li et al., 2017). If the frictional forces are provided from non-moving surrounding medium or substrates, the equation of a particle motion may be written as V_C|p = F_C|p/γ for each particle, where V_C|p is the velocity of the pth particle, F_C|p is the summation of cell–cell interaction forces exerted on the pth particle, and γ is the coefficient of viscous drag and frictional forces provided from surrounding medium or substrates; we call it “Model-A”. This assumption has been used in both particle-based models and vertex models (a standard model for epithelial cells) including simulations for the blastocysts (Honda et al., 2008a; Krupinski et al., 2011; Fletcher et al., 2014; Li et al., 2017; Nissen et al., 2017; 2018). On the other hand, if viscous frictional forces are dominantly provided from neighboring particles in motion, relative velocities between particles become important, and the equation of a particle motion may be written as follows;

$F_{\left.\mathrm{C}\right|p}={\displaystyle \sum _{m=1}^N}\left[\left\{{\omega }_{\mathrm{g}}\left(D_{pm}\right)\gamma \right\}\left\{\left(V_{\left.\mathrm{C}\right|p}-V_{\left.\mathrm{C}\right|m}\right)·e_{pm}\right\}{e}_{pm}\right],$

where m is the ID of particles interacting with the pth particle, V_C|m is the velocity of the mth particle, and N is the total number of the interacting particles. D_pm is the distance between the pth and mth particles, ω_g (D_pm) is a distance-dependent weight for γ, and e_pm is the unit vector from the xyz-coordinates of mth to pth particles. ω_g (D_pm) decays along D_pm. We call it “Model-B”. This assumption has been used in particle-based models especially for three-dimensionally packed cell aggregates (Nikunen et al., 2003; Basan et al., 2011; Podewitz et al., 2015). A similar assumption based on relative velocities was also adopted in a three-dimensional vertex model (Okuda et al., 2015b). Note that, in our inference method, equivalent force values were obtained between the above two assumptions (Koyama et al., 2023).

2.3 Force inference method

In our force inference method, the above model is fitted to nuclear tracking data. The cost function to be minimized is as follows; $G=G_{xyz}+G_{F_0}$, where G_xyz is the difference between xyz-coordinates of each particle of in vivo data and particle simulations, and $G_{F_0}$ is the additional cost function for cell–cell distance-dependent weight as previously defined (Koyama et al., 2023). The definition is described in Supplementary Material. Our inference method was validated by applying this method to simulation data generated under pregiven forces (Koyama et al., 2023).

3 Results

3.1 Inference of effective pairwise potential in synthetic data with cavity

To generate simulation data with expanding cavities, we modeled them as follows. In real tissues, cavities are generated by secreting liquid from cells. The cell layers surrounding the cavities are sealed by the tight junctions, and thus the pressure of the cavities push the cells outward (Figure 1D) (Moriwaki et al., 2007; Wang et al., 2018; Dumortier et al., 2019). In simulations, we implemented expanding spherical cavities as temporally-evolved spatial constraints; the particulate cells cannot penetrate into the cavities (Supplementary Figures S1A, B and Supplementary Movie S1). Although this constraint was not functional to prevent the particles from going outside of the cavities (Supplementary Figure S1B), the particles were almost constrained to move on the surfaces of the cavities.

Figure 2 is the simulation results based on Model-A, and Figures 2Ai shows the snapshot of the simulations where the spherical cavity was expanding from the radius = 9.6–12.6 μm. Figures 2Aii shows inferred effective distance–force (DF) curves under three different conditions of the expanding cavities; the inferred DF curves were not consistent with the DF curves from the LJ potential (Figures 2Aii, compare with the black broken line), where repulsive components were strongly detected especially at distant regions, which we call distant energy barrier (DEB). The strength of the repulsive components seemed to be positively related to the expansion rates. Under the higher expansion rate (Figure 2A; 9.7–>14.2), attractive components were diminished in the DF curves, meaning that the attractive forces derived from the LJ potential were effectively cancelled. In contrast, under the conditions of non-expanding cavities, the inferred DF curves were consistent with the LJ potential (Supplementary Figure S1C), probably because the non-expanding cavities did barely push the particles outwardly. Similar repulsive components were detected in systems with both larger cavities and larger numbers of particles (Supplementary Figure S1D; particle number = 150, radius of cavity = ∼50 μm). Moreover, our results in Figure 2A were scalable; i.e., if the sizes of the cavities, and both the distance and the force in the LJ potential were set to be multiplied by the same value in the simulations (e.g., ×5), inferred DF curves were also multiplied by the same value (Supplementary Figure S1E; e.g., radius of cavity = ∼40 μm).

Figure 2

Figure 2. Influence of external constraints to effective forces of cell–cell interaction Simulations were performed on the basis of a distance–force (DF) curve obtained from the Lennard-Jones potential (LJ). (A) Expanding cavities are assumed in the particle model. i) An example of simulations with the change in the radius of the cavity. The particles are of the mean size of cells. ii) Inferred effective DF curves with the distant energy barriers (DEBs). The pregiven DF curve (LJ) is shown in a broken black line. The numbers of the particles (=64) are the same among the different cavities. To exclude the possible influence of initial particle configurations on the force inference, particle tracking data from t = 210 min but not 0 min was applied to the force inference. (B) Expanding cavities are assumed whereas cell–cell interaction forces are not introduced. i) An example of simulations. ii) Inferred effective DF curves. (C) Reconstruction of cell–cell interaction forces under expanding cavities. i) All the data points in A and B under the radius of the cavity = 9.7 –> 14.2 μm are plotted on distance–force graph as a heat map, and the binned averages were calculated (yellow) which correspond to the DF curves in A and B. Each heat map is derived from forces inferred from one simulation. ii) DF curves obtained by subtracting the DF curves in B from those in A: i.e., A − B as exemplified in C-i. The definition of the heat map is described in Supplementary Material. (D) The subtraction performed in C is schematically explained in the relationship among (A–C) See the main text. Related figure: Supplementary Figure S1 (simulation procedures and inferred DF curves under both non-expanding cavities and larger expanding cavities with larger number of particles). Related movie: Supplementary Movie S1 (simulation).

We then examined whether this modification is meaningless or is a result of a correct reflection of the expanding cavities. By performing simulations under conditions with the expanding cavities but without the LJ potential (i.e., without both the attractive and repulsive forces of cell–cell interaction) (Figures 2Bi), we solely extracted the influence of the expanding cavities. Figures 2Bii shows that the inferred effective DF curves were exclusively composed of repulsive components; i.e., attractive components did not appear. Then, we checked the raw data where all data points were plotted on the distance–force space; in Figures 2Ci, “From B,” a heat map of the data points is shown. Force values calculated by binned averaging are also shown in yellow, which correspond to the DF curve in Figures 2Bii (red line). Interestingly, all the data points were closely located around the DF curve (“From B”), implying that the DF curve is a good approximation of the expanding cavities. Probably, due to the effect of the expanding cavities which causes the increases in distances among the particles, pairwise forces of particle–particle interactions became apparently repulsive.

Finally, by using the above information, we examined whether the correct DF curve (i.e., LJ potential) can be reconstructed. We subtracted the DF curves in Figure 2B from those in Figure 2A (Figures 2Ci, “From A–From B”), and the resultant DF curves were well consistent with the LJ potential under any expansion rates of the cavities (Figures 2Cii). These results suggest that the effects of the expanding cavities can be additive to effective DF curves.

In addition, the maximum values of the attractive forces in Figures 2C–2ii were 50–70% of that in the LJ potential. According to our previous study (Koyama et al., 2023), equivalent decreases in the inferred attractive forces were seen in systems without cavities, and sampling intervals also affected the values of attractive forces. Therefore, we think that the decreases are derived from a measurement error. The shorter the sampling interval, the more accurate the inferred forces (Koyama et al., 2023). However, in real three-dimensional tissues, we are usually unable to set very short sampling intervals due to technical difficulties, so we typically set it to be ≥3 min in the case of the mouse embryos as shown later. In Figure 2, we adopted a sampling interval equivalent to the case of the mouse embryos.

We interpret the above results as described in Figure 2D. A cell receives forces from cell–cell interactions (“LJ” in Figure 2D) and from the pressure of the cavity, where the vector of the pressure is perpendicular to the surface of the cavity (Figure 2D, top panel). The repulsive forces detected as cell–cell interactions in Figure 2A are described in Figure 2D, middle panel (light blue vectors), and the vector of the net force (a red vector) for a cell is perpendicular to the surface of the cavity, if the cells are almost uniformly distributed. Therefore, the force provided by the pressure of the cavity can be equivalent to the net force from the DEBs (Figure 2D, bottom panel), leading to the successful reconstruction of the LJ potential.

We conclude that the expanding cavities can be approximated as effective pairwise potentials. In general, effective pairwise forces/potentials are the summation of different components as explained in Figure 1C. External factors are one of the components to be summed up in molecular and colloidal sciences: e.g., by considering hydration of particles by solvents, effective pairwise potentials are modified so as to contain energy barriers around distant regions (Pettitt and Rossky, 1986; Israelachvili, 2011), and such effective potentials have been used to simulate behaviors of the systems. We will show the simulation outcome from the above inferred effective DF curve later.

3.2 Effective forces in mouse blastocyst contained DEB

We investigated whether a repulsive component appeared in effective forces inferred from the mouse blastocyst (Figure 3A, illustration). An inner cavity is formed through liquid secretion from the cells, and the embryos expand while maintaining their spherical shape (Moriwaki et al., 2007; Wang et al., 2018; Dumortier et al., 2019) (Figure 3A, bright field and cross-section). Trophectoderm (TE) cells form a single-cell-layered structure at the outermost surface of the embryo, whereas the cells of the inner cell mass (ICM) were tightly assembled, forming a cluster (Figure 3A, illustration). We performed cell tracking of both TE and ICM (Figure 3A, tracking) and inferred effective forces.

Figure 3

Figure 3. Inference of the effective force of cell–cell interaction in mouse blastocysts. (A) Blastocyst stages of mouse embryo are illustrated, and their confocal microscopic images are shown: bright field, maximum intensity projection (MIP) and cross-section (equatorial plane) of fluorescence of H2B-EGFP. Snapshots of nuclear tracking are also shown; blue spheres indicate the detected nuclei. ICM cells are located around the bottom region of the image, as indicated by *. Scale bars = 15 μm. (B) Inferred effective DF curve in the mouse blastocysts. In left panel, the blastocyst (#1) is shown by a red line. DEB is indicated. For comparison, an inferred effective DF curve from the compaction stage embryo before forming cavities is shown which we previously reported (green line) (Koyama et al., 2023). In the right panel, inferred effective DF curves from four blastocysts (#1–#4) are shown, whose bright field (BF) images are shown in Supplementary Figure S2A. (C) Inferred effective DP curves in the mouse blastocysts. Four embryos were examined. An enlarged view around longer distance is also shown in the inset. The cell numbers in the four embryos and their changes during imaging were as follows; 56–>58, 68–>76, 56–>65, and 58–>59. (D) Inferred effective DP curves in the mouse embryos where the cavitation was inhibited. S3226 whose target is Na⁺/H⁺ exchanger inhibits the cavitation of the blastocyst (Kawagishi et al., 2004). In the left panel, the microscopic images are shown. In the middle panel, the DF curves from three embryos are shown; the light blue line is from the embryo in the left panel. In the right panel, the DP curves are shown. The cell numbers in the three embryos were as follows; 32, 37, and 32 (the numbers were not changed during imaging). Scale bars = 15 μm. Related figure: Supplementary Figure S2 (inferred DF and DP curves based on Model-B) and Supplementary Figure S3 (experimental procedures of S3226 treatment).

Figure 3B shows the inferred effective DF curve from the blastocyst. The DF curve contained the repulsive component around distant regions (Figure 3B, DEB). For comparison, an inferred DF curve from an embryo before forming a cavity is shown, which did not contain DEB (Figure 3B “compaction” meaning the compacted embryo at the morula stage). In addition, Figure 3C shows the distance–potential energy (DP) curves [U_C–C(D)] which can be obtained by mathematically calculating the integral of the DF curves [F_C–C(D)] along the distance (D) (Koyama et al., 2023); DP and DF curves are mutually convertible. In the DP curves, the DEB was expressed as a long slope with a negative gradient.

To further examine whether the DEB was derived from the cavity, we inhibited the cavitation by using an inhibitor for Na⁺/H⁺ exchanger (S3226) (Kawagishi et al., 2004). Figure 3D shows the inhibitor-treated embryos; the experimental procedures are shown in Supplementary Figure S3. Figure 3D shows the inferred effective DP curves, where DEB was diminished, indicating that the cavity yielded the DEBs.

We also evaluated a model dependency. We inferred effective forces under the assumption of Model-B. As shown in Supplementary Figures S2B,C, the inferred DF and DP curves were similar to those in Figures 3B,C. Therefore, model choice did not significantly affect the inferred forces.

3.3 Effective potentials reproduced cavity-harboring structures

Effective DF/DP curves are the summation of the multiple components including the excluded volume effect, adhesive forces, and external factors, and are useful to study behaviors of systems in molecular and colloidal sciences. We evaluated whether inferred effective DF/DP curves from the cavity-harboring systems can reproduce their morphologies. The effective DF curves inferred from the mouse blastocyst, and the artificial tissue harboring a cavity (Figures 2, 3) were used for simulations, and the systems were relaxed to find stable states. When initial configurations of particles were set to be cavity-harboring structures (Figure 4A), the structures were maintained (Figures 4B,C, left panels). To examine whether the DEBs contribute to the formation of the cavities, we artificially eliminated the DEBs from the DF curves (i.e., the DF curves were truncated so as not to include the DEBs), and performed simulations. We found that cavities were not maintained but aggregates—filled with the particles—were formed (Figures 4B,C, “DF curve w/o DEB”). For comparison, we used the Lennard-Jones potential, and showed that the cavity-harboring structure was not maintained but became an aggregate (Figure 4D). In addition, when we used a DF curve inferred in cysts of Madin-Darby canine kidney cells, a cavity was also maintained (Supplementary Figure S4). In the case of the blastocysts, the particles were sparsely distributed (Figure 4B). This was caused by the very long range of DEB in the blastocysts. In the blastocysts, the mean diameters (defined in Figure 1C) were ∼10 μm (Figures 3B,C), and the ranges of DEBs were from ∼15 μm to ∼60 μm (Figure 3C and its inset), meaning 1.5∼6.0-fold of the mean diameters. On the other hand, in the Lennard-Jones simulation data, the mean diameters were ∼5 μm (Figure 2C), and the ranges of DEBs were from ∼7 to 15 μm (Figure 2A), meaning 1.4∼3.0-fold of the mean diameter. Note that these simulations explained the maintenance of cavities but not the generation of cavities from filled aggregates.

Figure 4

Figure 4. Outcomes of simulations based on inferred distance–force curves with DEBs. (A) Initial configurations of particles are exemplified, that harbor a cavity. A cross-section is also shown in the right panel. (B–D) Simulations results under the effective DF curves derived from the mouse blastocyst (B), and the synthetic tissue harboring an expanding cavity (C) in Figure 2 are shown. The DF curves used are from #1 in Figure 3B for the blastocyst, and from Figure 2A (8.6–>10.1). E is the simulation result under the LJ potential for comparison. In addition to the left panels (“Full DF curves”), the DF curves in which the DEBs were artificially eliminated were used in the right panel (“DF curves w/o DEB”); in C and D, analyzes for the right panel were not performed (“ND”). For each simulation result, cross-sections are also shown, and in the case of the mouse blastocyst, two parallel cross-sections (shown in blue and green) are merged. The diameters of the blue spheres are equivalent to the mean diameter.

3.4 Modeling of distance–force curves with two parameters

To comprehensively understand the effect of the differences of the profiles in the DF curves on morphogenesis, we modeled the DF curves as a simple mathematical equation. To analyze the contribution of DEBs, the frameworks based on the LJ and Morse potentials were not suitable, because they do not contain repulsive components around longer distances. Hence, we selected the following equation (Figure 5A, left panel, blue line):

$F\left(D\right)=\epsilon {\left(D-D_0\right)}^{-N}\cos \left\{\frac{2\pi \left(D-D_0\right)}{\delta }\right\},\text{for\hspace{0.17em}}D>D_0(1)$

Here, F is the force, D is the distance. N is the exponent of D, affecting the decay of the DF curves along D (X-axis) (Figure 5A, left panel, green broken line). To model DEB, we introduced a cosine function whose wavelength is affected by δ; the combination of the above exponential function (Figure 5A, left panel, green broken line) with the cosine function (Figure 5A, left panel, purple broken line) resulted in a function bearing DEB (Figure 5A, left panel, red line). We set a cut-off distance at the end of DEBs [i.e., F (D) = 0, if D > 5δ/4 + D₀ in Figure 5A, right panel]. D₀ can translate the DF curves along D, (Figure 5A, left panel, red vs blue lines). By this translation, the mean diameter (D₁) of particles is changed from δ/4 to δ/4 + D₀, and the location of DEB becomes (3δ/4 + D₀ ∼ 5δ/4 + D₀) where the center of DEB is at the distance of D₂ = δ + D₀ (Figure 5A, left panel). Therefore, this translation simultaneously modifies the relative positions (rD₃) of DEBs to the mean diameters (D₁) (Figure 5A, left panel, rD₃ = D₂/D₁). N also affects the strengths of DEBs (Figure 5A, right panel, F_rep/F_att is the relative strength of DEB). ε determines the scale of the forces, that does not contribute to morphologies under stable states.

Figure 5

Figure 5. Modeling of distance–force curves and diagrams of simulation outcomes; cavity-harboring structures. (A) Modeling of DF curves using four parameters: N, D₀, δ, and ε. In the left panel, a DF curve without the cosine term is always >0 (repulsive) (F = εD ^−N, green broken line; N = 3, ε = 1.0). Upon introduction of the cosine (purple broken line), the curves acquire regions with <0 (attractive) and a DEB (red line; N = 3, δ = 1.5, ε = 0.2). Upon introduction of D₀, the curves are translated along the distance (blue line; D₀ = 0.3). As shown in the blue line, the mean diameter of particles corresponds to δ/4 + D₀ (D₁), the DEB ranges from 3δ/4 + D₀ to 5δ/4 + D₀ where the center of the DEB is at δ + D₀ (D₂). In the right panel, N affects the gradient of the curves and the height of DEBs, and thus, the relative strength of DEB is modified. A cut-off distance was set: F (D) = 0, if D > 5δ/4 + D₀. (B) Outcomes of simulations under various values of N and D₀* in the presence of DEB. From the equation defined in A, DF curves from D = 0 to D at the end of DEBs (“w/DEB”) were used as exemplified in the graph of DF curve (w/DEB). The initial configuration is shown in the left bottom. On the diagram, some conditions were designated by symbols (circle, square, and triangle), and their simulation results are visualized on the bottom row (blue spheres) with the profiles of the DF curves. Cross-sections are visualized for some simulation results. The diameters of the blue spheres are equivalent to the mean diameter. For “cup” (orange rings, mouth of the cup; white dots, particles located on the ring; crosses, the particles located neither on the ring nor at the center of the ring). For “tube” (orange arrows, open ends of the tube). For “cavity with extra particles” (arrows, extra particles where multiple- but not single-layered). The definition of each phase is explained in Supplementary Material. Gray crosses, data points of all simulations for depicting the lines separating each phase (see Supplementary Material). (C) Simulation outcomes in the absence of DEBs. DF curves from D = 0 to D before starting DEB (w/o DEB) were used as exemplified in the graph of DF curve (“w/o DEB”). Related figures: Supplementary Figure S5 (phase diagrams under a different number of particles), Supplementary Figure S6 (maintenance of tubes was examined under DF curves with DEBs), Supplementary Figure S7 (smoothness of the surface of the cavities), Supplementary Figure S8 (phase diagram for conditions of D₀* <1.2). Related movie: Supplementary Movie S2 (the cavity-harboring structure was stably maintained under fluctuations).

All parameters and variables related to lengths/distances were normalized by (δ/4), leading to generation of dimensionless parameters and variables such as D₀^* = D₀/(δ/4). Therefore, we can present simulation outcomes in two-dimensional space (D₀^* and N) in a manner similar to a phase diagram (Figure 5B) as shown in the next section. Although there can be other equations modeling DF curves, the advantage of this equation is that it is composed of just two parameters (D₀^* and N). All simulations were performed by using Model-A.

3.5 Appropriate position and strength of DEB were essential for cavity-harboring structures

First, we examined what kind of DF curves can maintain a cavity-harboring structures. In simulations of multicellular systems, the morphologies are determined by both DF curves and initial configurations, because energetic local minimum states (i.e., metastable state) as well as global minimum states are meaningful. A cavity-harboring structure was set as the initial configuration (Figure 5B, “Initial configuration”), and simulations were performed until the systems reached steady states including metastable states; i.e., the velocities of the particles became almost undetectable.

Figure 5B shows a phase diagram with representative images of the simulations, on which two regions of cavity-harboring structures were identified (“cavity” and “cavity with extra particles”). The “cavity” class absolutely depended on the DEBs, because DF curves whose DEBs were eliminated did not form the “cavity” as shown in Figure 5C (“w/o DEB”; truncated version of DF curves exemplified at the upper left). We confirmed that the cavity was not broken by slight positional fluctuations (Supplementary Movie S2 with the legend in Supplementary Material), indicating that the structures were stable. In Figure 5B, the “cavity” class was located around smaller N values (less than ∼2.0), otherwise the “aggregate” class was formed, suggesting that the strengths of the DEBs, defined in Figure 5A (right panel), may be effective. Furthermore, larger D₀ values were required for the “cavity” class, otherwise the “multiple clusters/disorganized” class was formed, suggesting that the relative positions of DEBs (rD₃ = D₂/D₁ in Figure 5A, left panel) should be close to the mean diameters. Under a different number of particles, a similar trend was observed (Supplementary Figure S5; particle number = 64 and 150). The parameter values leading to the “cavity” class were almost equivalent between Figure 5B and Supplementary Figure S5, and thus there was no clear relationship between the features of DEBs and the sizes of the cavities.

In addition, the “cups” and “tubes” classes were found between the “cavity” and the “multiple clusters/disorganized” classes in Figure 5B but not in Figure 5C, indicating that the DEBs also contribute to the formation of these structures. Because particle models are discretized systems, symmetry breaking for generating the cups and tubes can occur even under the directionally uniform DF curves. The maintenance of tubes was confirmed to be dependent on DEBs (Supplementary Figure S6). Detailed discussion about the “cup,” “tube,” “multiple clusters/disorganized,” and “cavity with extra particles” classes is described in Supplementary Material.

The simulation outcomes of the “cavity” class showed slightly deformed spheres (Figure 5B and Supplementary Figure S5A). We evaluated the smoothness of the surface of the cavity by calculating the distances of each particle from the center of the cavity (Supplementary Figure S7A). The values of the standard deviation (SD) of the distances were 3–8% of the radii of the cavities (Supplementary Figure S7B, C). Therefore, absolutely smooth cavities were not generated by the DF curves with DEBs, which is a limitation of DEBs for describing spherical cavities. The SD values depended on both D₀ and N; smaller values of D₀ and N resulted in smooth cavities (Supplementary Figure S7B, C)).

3.6 Profiles of DF curves for describing 2D sheet

To further investigate the capability of DF curves, we next applied a two-dimensional (2D) sheet as the initial configuration which corresponds to a monolayered cell sheet (Figure 6A, “Initial configuration”). Figure 6A shows a phase diagram with representative simulation images. The 2D sheets were maintained (the “sheet” class), which was largely dependent on the presence of DEBs (Figure 6A, w/DEBs vs 6B, w/o DEBs; truncated version of DF curves exemplified at the upper left in Figure 6B). We confirmed that the 2D sheets were not broken by slight positional fluctuations (Supplementary Movie S2 with the legend in Supplementary Material), indicating that the structures were stable. The region of the “sheet” class was located around larger values of D₀, which was a similar trend to the case of the cavity in Figure 5. The relationship between the sheet and cavity classes will be analyzed in detail later.

Figure 6

Figure 6. Diagrams of simulation outcomes; 2D sheets, spherical/non-spherical aggregates. (A) Outcomes of simulations starting from a two-dimensional (2D) sheet with slight positional fluctuations. The results are shown in a similar manner to Figure 5B, except that the diameters of the blue spheres are equivalent to 0.4× the mean diameter. D₁, mean diameter; W₁, width of attractive force regions; P₁, potential energy derived from attractive forces. Note that D₀* can have a minus sign; because the mean diameter (D₁) of a particle is δ/4 + D₀, D₀ should be > −δ/4 where D₀* becomes > −1. (B) Simulation outcomes without DEBs. Both Figures 6B, 5C were generated under the conditions of w/o DEBs, but the ranges of D₀* were different each other, resulting in that the “cavity” class was not appeared in Figure 5C. We did not divide the “aggregate” class in Figure 5C into the two “aggregate” classes as defined in Figure 6B. Related figure: Supplementary Figure S9 (possible morphological transitions). Related movie: Supplementary Movie S2 (2D sheets were stably maintained under fluctuations), and Supplementary Movie S3 (possible morphological transitions).

In addition, by expanding the parameter range in the phase diagram (i.e., D₀* = −0.2∼1.2), the cavity-harboring structures were newly generated even in the absence of the DEBs (Figure 6B); because the parameter values were different from those for the cavity class in the previous figure (Figure 5B), the cavity classes in Figures 5B, 6B were different each other. Then, by applying the similar parameter values to the phase diagram in Figure 5B, the cavity-harboring structures were also emerged (Supplementary Figure S8), indicating that the formation of this cavity class does not depend on initial configurations. Although we do not understand the precise mechanism of the de novo formation of the cavity, we guess that, if particles have a property to strongly assemble each other (e.g., due to strong cell–cell adhesion), the formation of a small cavity at the center can be energetically favorable; note that the sizes of the particles in Figure 6B were set to be 0.4× mean diameter of the particles, and thus, the size of the cavity should be quite small. By temporally changing N, D₀, and existence of DEB, we presented possible morphological transitions among the classes in Figures 5, 6 (Supplementary Figure S9 and Supplementary Movie S3).

3.7 Profiles of DF curves for describing spherical/non-spherical aggregates

In Figure 6, we found that both the “spherical” and “non-spherical aggregate” classes were generated w/and w/o DEBs. Smaller values of N and/or D₀ led to the spherical ones (Figure 6B). Under smaller values of D₀, because the width (W₁; see the bottom row in Figure 6) of the attractive force regions in DF curves became relatively long compared to the mean diameter (D₁, in Figure 6), the resultant increased attractive forces may cause the spherical aggregates. Under smaller values of N, the potential energies (P₁, in Figure 6) derived from the attractive forces are increased, possibly leading to the spherical aggregates.

3.8 Effective forces in 2D sheets included DEB

As shown in Figures 6, 2D sheets were maintained under the DF curves including DEB. Similarly, the cups and tubes in Figure 5 were also maintained under the DF curves including DEBs; these structures can be interpreted as a curved sheet. As we previously showed, DEB was detected in tissues with pressurized cavities (Figures 2, 3). But it is unlikely that 2D sheets, cups and tubes have such pressures. Is the relationship between these structures and DEB an artifact or are there any reasonable explanations?

To investigate this issue, we used synthetic 2D sheets generated by simulations (Figures 7Ai). The 2D simulations were performed under the Lennard-Jones (LJ) potential as the pregiven force, and the DF curves were inferred. As shown in Figures 7Ai, the inferred DF curve was similar to the LJ potential, but the curve contained a clear DEB. This was contrast to the inferred DF curve from the simulation of a 3D aggregate, where DEB was rarely detected (Figures 7Aii).

Figure 7

Figure 7. Radial distribution functions (RDF) in 2D and 3D situations. (A) Comparison of inferred DF curves in 2D i) and 3D ii) situations. The DF curves are shown with snapshots of simulations and with the pregiven potentials (LJ). DEB is indicated. The condition of the simulations is (sampling interval = 3 min, SD value of force fluctuation = 1,000%, and persistency of force fluctuation = 1 min). (B) The relationship between RDF and effective pairwise potentials is described. An RDF for the central particle (dark blue) is illustrated. Densities of particles are shown as the function of distances from the central particle. Higher density at a distance means an energetically favorable state, and therefore, the potential is expected to be lower at the distance. k_T is the effective temperature as an indicator of non-equilibrium fluctuations in multicellular systems, in analogy to k_BT where k_B is the Boltzmann constant and T is the temperature. (C) RDFs in 2D and 3D situations were calculated in simulations under the LJ potential. The snapshots of the simulations and the RDFs are shown. In these simulations, fluctuations of the forces given by the LJ potential were introduced with the persistency. Under longer persistency, the particles tend to easily change their positions. The condition is (SD value of force fluctuation = 300%, and persistency of force fluctuation = 1 min, sampling interval = 500 min). The numbers of particles are 100 (2D) and 500 (3D). Related figure: Supplementary Figure S10 (distribution of particles under the closest packing situations under 2D and 3D conditions).

We hypothesized that the radial distribution functions (RDF) between 2D and 3D situations were different, because RDFs are related to the profile of distance–potential curves in atomic/molecular/colloidal sciences (Israelachvili, 2011). RDFs mean the density of particles at a certain distance from a particle of interest (Figure 7B). The density may be reflected as the effective potential energy at the distance; high density corresponds to low potential energy (Figure 7B). In other words, by using RDFs, effective pairwise potentials are calculated if the systems are under equilibrium conditions.

We calculated RDFs in 2D and 3D simulations performed under the LJ potential (Figure 7C). In these simulations, fluctuations of the forces of cell–cell interactions were assumed with the persistency of the fluctuations (the details are described in Supplementary Material). Under the 2D situations, clear peaks were found where the particle density was high. Importantly, between the peaks, the density was very low (Figures 7Ci, black arrow). By contrast, under the 3D simulations, relatively unclear peaks were detected, and the regions with low particle densities almost disappeared (Figures 7Cii). These results may be reasonably interpreted by considering the closest packing situations under 2D and 3D conditions, where the peaks of high densities were frequently detected along the distances under the 3D condition compared with 2D (Supplementary Figure S10). Additionally, even under the 2D simulations, the longer persistency caused disappearance of the regions with low particle densities, because the long persistency promotes the movements of the particles. These results suggest that 2D sheets spontaneously yield regions with low particle density in RDF, probably due to the geometric reasons, and that this feature can be reflected to the effective pairwise potentials as DEB.

3.9 Expanding cavity was found before forming multiple clusters

We have been presented the two “cavity” classes in the phase diagram (Figures 5B, 6). The formation of “cavity” class in Figure 6 did not depend on DEB, and the “cavity” class in Figure 5B may be derived from the geometric reason of cell sheets (i.e., curved sheets) because the location of the “cavity” class in Figure 5B overlapped that of the “sheet” class in Figure 6. In the case of the simulation of the blastocyst (Figure 3B), the cavity was largely expanded, resulting in the sparse distribution of the particles. We failed to find such sparsely distributed pattern in the phase diagrams in Figures 5, 6 and Supplementary Figure S5. Then, we carefully observe the dynamics of the simulations in Figure 5 and Supplementary Figure S5 whose initial configurations were “cavity,” and we found that, before forming the “multiple-clusters/disorganized” class, expanding cavities really existed. Figure 8 shows the temporal dynamics of the simulations in Supplementary Figure S5 with the particle number = 64 (D₀* = 1.33 and N = 0.2). The initial cavity was rapidly expanded, while maintaining its spherical shape. After it stopped expansion, the multiple clusters were formed. Therefore, the process of the expansion can be simulated by considering the DEB. Note that the boundary of the two origins of DEB may correspond to D₀* ∼2.0, around which there is a risk to mix the two origins.

Figure 8

Figure 8. Simulation of expanding cavity. An outcome of a simulation under a condition of N = 0.2 and D₀* = 1.333 is shown with its temporal dynamics. Before forming multiple-clusters, the cavity was expanded. The particle number and simulation condition were the same as that in Supplementary Figure S5A.

4 Discussion

In this study, by inferring effective forces of cell–cell interactions in tissues harboring expanding cavities, we showed that the effect of the cavities was addible to the effective forces of cell–cell interactions as a repulsive component (DEB). The resultant distance–force curves reproduced the tissue structures harboring cavities in simulations. A similar repulsive component was also detected in two-dimension (2D) cell sheets, which may be derived from geometric reasons of 2D systems. Finally, by assuming distance–force curves with the repulsive components, tissue structures of 2D sheets, cups, and tubes were reproduced in simulations; it is conceivable that, because 2D sheets, cups, and tubes may not harbor pressurized cavities, the origin of their DEBs is the geometric reasons of 2D or quasi-2D systems. These results suggest that external components and geometric constraints can be correctly added to the effective forces of cell–cell interactions. This framework considering the repulsive components may be expanded for simulating cell sheets with complicatedly curved surfaces such as epithelial folding as well as cavity-harboring structures, cups, and tubes (Wang et al., 1990a; Montesano et al., 1991; Honda et al., 2008b; Eiraku et al., 2011; Solnica-krezel and Sepich, 2012; Belmonte et al., 2016; Koyama et al., 2016; Nissen et al., 2018; Rivron et al., 2018; Serra et al., 2019).

4.1 Influence of external factors on effective potentials

Our results suggest that expanding cavities and geometric features of 2D sheet can be incorporated into the effective pairwise potentials. Another external factor in the mouse blastocyst is the Zona Pellucida (ZP) surrounding the embryo. In the late blastocyst, the ZP is outwardly pushed by the pressure of the cavity through the TE cells when the blastocysts are inflating, resulting in that the ZP is stretched. Conversely, the stretched ZP inwardly pushes the TE cells. Therefore, we think that the ZP should effectively reduce the magnitude of the vector from the pressure in Figure 2D, “Pressure from cavity (red vector),” meaning that the effect of the ZP is also incorporated into the inferred effective potentials. In atomic/molecular/colloidal sciences, it is known that some external factors are additive to effective pairwise potentials as mentioned in the Result section; e.g., hydration of ions, external pressures which can modify RDF (Israelachvili, 2011). In multicellular systems, it is unknown what kind of external factors are additive to the effective pairwise potentials. We speculate that, if external factors uniformly exert forces on cells, the effect may be correctly incorporated into the effective potentials.

In our present study, we focused on tissues with a single cavity of relatively larger sizes. On the other hand, at the beginning of the cavitation in the mouse embryos, multiple cavities with sizes much smaller than the cell size are generated (Dumortier et al., 2019). Smaller cysts are also known in the MDCK cells and other tissues (Wang et al., 1990b; Jewett and Prekeris, 2018). In such cases, particle model-based descriptions do not have sufficient spatial resolutions to simulate their morphologies, whereas other high-resolution models including the phase field models are suitable (Nonomura, 2012). On the other hand, because such small cavities are generated on cell–cell junctions by partially breaking cell–cell adhesion, it might be possible that, in our inference method, the generation of small cavities could be reflected as repulsive components of cell–cell interactions.

In the blastocysts and MDCK cysts, hydrostatic pressure of the cavities largely contributes to the cavitations. On the other hand, there can be other type of cavitation/lumenogenesis where hydrostatic pressure may be lesser (Jewett and Prekeris, 2018; Asleh et al., 2023). We do not know whether DF curves in such tissues become similar to those in the blastocysts, the MDCK cysts, or “cavity class” in Figure 6.

According to our result that the effect of external factors is additive to the effective potentials, it is speculated that, if one can assume that external factors are equivalent between two tissues (e.g., wild-type and mutant tissues), differentials of the inferred effective DF curves between the two tissues correspond to the differences of cell–cell interaction properties. Therefore, even in the case that the precise effect of the external factors is not experimentally measurable, one can evaluate the differences of cell–cell interaction properties which may involve cell–cell adhesion. Similarly, we can also predict what will happen if cell-cell interaction forces are modified under the condition of non-cellular components. For instance, under the assumption of the fixed condition of non-cellular components (e.g., fixed expansion rate of cavities), we can simulate and predict cell behaviors in the case that cell-cell interaction forces such as adhesive forces are modified.

4.2 Relationship of the effective potential to cell shape, polarities, temporal changes, and cell–cell junction tension

In our simulations (e.g., Figures 5, 6), neither cell shape nor polarities were explicitly considered. On the other hand, the TE cells in the blastocysts exhibit a highly stretched squamous epithelium and an apico-basal polarity (Latorre et al., 2018). In addition, polarized tissues can have the planar cell polarity (Koyama et al., 2019). Although a cell in particle models look like a spherical cell shape, cell shapes and polarities can be considered to some extent by assuming anisotropy of the pairwise potentials; i.e., direction-dependent profiles of the potentials (Nissen et al., 2018). Such anisotropic potentials can elevate the capability of describing complicated tissue morphologies (Nissen et al., 2018). From our inferred forces, we may extract anisotropic potentials in principle if the anisotropy exists; a potential perpendicular to apico-basal polarity and a potential along apico-basal polarity. In addition, the TE cells are gradually stretched during the expansion of the blastocysts, and thus, the distance–force curves may be temporally changed. In the case of epithelial cells, cell–cell junction tensions are often considered and measured (Ishihara and Sugimura, 2012). The relationship between the profiles of the effective potentials and the cell–cell junction tensions is not understood, and we can theoretically analyze it by applying our inference method to simulation outcomes generated by the epithelial cell models such as the vertex model where the cell–cell junction tensions are assumed.

4.3 Comparison of other inference methods with other models

To understand morphogenetic events, simulations based on mechanical parameters are essential. Models such as vertex, Cellular Potts, and phase field models are often utilized, especially for epithelial cells (Honda et al., 2008a; Nonomura, 2012; Camley et al., 2014; Fletcher et al., 2014; Belmonte et al., 2016). These models successfully recapitulated various morphogenetic events, including the formation of the mouse blastocyst (Honda et al., 2008a). In these models, sub-cellular structures such as cell–cell junctions and their line tensions containing adhesive forces are explicitly modeled. In particular, the vertex models recapitulated well the behaviors of epithelial cells exhibiting polygonal cell shapes and three-dimensional cells exhibiting polyhedrons (Okuda et al., 2015a). On the other hand, non-polygonal epithelial cells are also found where cell–cell junctions are curved or interdigitated (Ishimoto and Morishita, 2014; Miyazaki et al., 2023), and thus, corresponding simulation frameworks become required. Because particle-based models are coarse-grained models, they cannot explicitly express cell–cell junctions. Due to their lower costs for programing, particle models were applied to three-dimensional simulations including intricately curves cell sheets (Nissen et al., 2018). Moreover, particle models were also applicable to mesenchymal cells that should exhibit intricated shapes (Paramore et al., 2024).

For the simulation of multi-cell behavior, the acquisition of parameter values is important. However, it is still challenging to measure cellular and non-cellular parameters in three-dimensional systems with high spatiotemporal resolution. Image-based force inference methods are important techniques (Tambe et al., 2011; Ishihara and Sugimura, 2012). Our idea can reduce the degrees of freedom derived from non-cellular factors, and may be adapted to other image-based inference methods when non-cellular factors exist. In analogy to the discussion in 4.1, differentials of forces between two tissues may be evaluated in the above image-based force inference methods. Our goal is to develop a framework for explaining diverse morphogenesis in a minimal model. Approximations of some external factors as effective potentials of cell–cell interaction serve to reduce degrees of freedom.

5 Materials and methods

5.1 Mouse embryos

Embryos at the embryonic days 2.5 (E2.5) were obtained after mating homozygous R26-H2B-EGFP knock-in male mice which constitutively express EGFP (enhanced green fluorescent protein)-fused H2B (histone2B proteins) (Kurotaki et al., 2007) and ICR female mice (Japan SLC). The embryos were cultured in EmbryoMax KSOM + AA with D-Glucose (Millipore, United States) covered with mineral oil on a glass bottom dish (35-mm; 27-mm ϕ, Matsunami, Japan) at 37°C. Then, the E4.0 embryos were subjected to microscopic imaging. Four distinct embryos were analyzed. Note that the mouse blastocysts sometimes lose the integrity of the sealing conferred by the tight junctions in TE cells, resulting in reduction of the pressure of the cavities with the subsequent volume decrease of the cavities (Niimura, 2003; Chan et al., 2019). This event is called collapse. Our four embryos, however, did not experience the “collapse” events during our microscopic imaging.

In the experiments using S3226 (Supplementary Figure S3), the E2.5 embryos were transferred to the medium with the drug. After 15 h culture, some embryos failed to form cavities (Supplementary Figure S3, #1 and #2), while others formed smaller cavities (Supplementary Figure S3, #3) than that under no-drug conditions (Supplementary Figure S3, No drugs). Live imaging was performed, and then a rescue process was carried out, resulting in that the embryos formed, at least, smaller cavities compared with that under no-drug conditions. The concentration of S3226 was 50 μM. We also tried Ouabain (Na⁺/K⁺ ATPase inhibitor) (Chan et al., 2019), but the cavitation was not drastically inhibited, and thus the embryos were not suitable for the cavity (−) control.

Animal care and experiments were conducted in accordance with National Institutes of Natural Sciences (NINS), the Guidelines of Animal Experimentation. The animal experiments were approved by The Institutional Animal Care and Use Committee of NINS (ref. #17A030, #18A026, #19A021, #20A014, and #21A039).

5.2 Microscopic live imaging

Fluorescence images were acquired on a Nikon A1 laser scanning confocal microscope (Nikon, Japan) equipped with a ×20 objective (PlanApo; Dry; NA = 0.75, Nikon, Japan). A stage-top CO₂ incubation system was used (INUG2-TIZ, Tokai Hit, Japan) on the inverted microscope. H2B-EGFP was excited using a 488-nm laser.

For imaging of mouse embryos, 50–100 frames were acquired at 3-min intervals in 30–40 z-slices separated by 2.5 µm.

The confocal images obtained from the mouse embryos were subjected to the procedure for nuclear detection and tracking. Using the Imaris software (Oxford instruments/Bitplane, UK), the nuclei were automatically detected, followed by manual corrections. The nuclear tracking was also performed automatically followed by manual corrections. Through the manual corrections, the accuracies of the nuclear detection and tracking became essentially 100%, to the best of our judgment.

5.3 Inference method and simulation procedures

The inference method is as we previously reported (Koyama et al., 2023), where handling of cell division and cell death were also defined. Simulation procedures are also as previously reported (Koyama et al., 2023), except that cavities are considered. Details of the modeling of the cavities are described in Supplementary Material.

Data availability statement

The datasets presented in this study can be found in online repositories. This data can be found here: https://doi.org/10.6084/m9.figshare.24125757.

Ethics statement

The animal study was approved by Institutional Animal Care and Use Committee of NINS. The study was conducted in accordance with the local legislation and institutional requirements.

Author contributions

HK: Conceptualization, Data curation, Formal Analysis, Funding acquisition, Investigation, Methodology, Project administration, Resources, Software, Validation, Visualization, Writing–original draft, Writing–review and editing. HO: Conceptualization, Investigation, Validation, Writing–review and editing. TO: Funding acquisition, Resources, Writing–review and editing. AI: Conceptualization, Investigation, Software, Writing–review and editing. KN: Formal Analysis, Investigation, Writing–review and editing. KK: Conceptualization, Software, Writing–review and editing. TF: Conceptualization, Supervision, Writing–review and editing.

Funding

The author(s) declare that financial support was received for the research, authorship, and/or publication of this article. This work was supported by following grants: Japan Ministry of Education, Culture, Sports, Science and Technology Grant-in-Aid for Scientific Research on Innovative Areas “Cross-talk between moving cells and microenvironment as a basis of emerging order” for HK, the National Institutes of Natural Sciences (NINS) program for cross-disciplinary science study for HK, a Japan Society for the Promotion of Science (JSPS) Grant-in-Aid for Young Scientists (B) for HK (17K15131), for Scientific Research (C) for TO (18K06234), and for Young Scientists (B) for TO (16K18544), and for Transformative Research Areas (A) for HK (24H01990, 24H01414), a MEXT/JSPS Grant-in-Aid for Scientific Research on Innovative Areas for TO (17H05627), the Inamori Foundation for TO, and the Takeda Science Foundation for TO.

Acknowledgments

We thank Drs. Hitoshi Niwa and Yayoi Toyooka for providing the piggy-bac plasmids and technical suggestions.

Conflict of interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher’s note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Supplementary material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fcell.2024.1414601/full#supplementary-material

References

Asleh, M. A., Zaher, M., Asleh, J., Jadon, J., Shaulov, L., Yelin, R., et al. (2023). A morphogenetic wave in the chick embryo lateral mesoderm generates mesenchymal-epithelial transition through a 3D-rosette intermediate. Dev. Cell. 58, 951–966. doi:10.1016/j.devcel.2023.03.017

PubMed Abstract | CrossRef Full Text | Google Scholar

Basan, M., Prost, J., Joanny, J.-F., and Elgeti, J. (2011). Dissipative particle dynamics simulations for biological tissues: rheology and competition. Phys. Biol. 8, 026014. doi:10.1088/1478-3975/8/2/026014

PubMed Abstract | CrossRef Full Text | Google Scholar

Bazellières, E., Conte, V., Elosegui-artola, A., Serra-picamal, X., Bintanel-morcillo, M., Roca-cusachs, P., et al. (2015). Control of cell–cell forces and collective cell dynamics by the intercellular adhesome. Nat. Cell. Biol. 17, 409–420. doi:10.1038/ncb3135

PubMed Abstract | CrossRef Full Text | Google Scholar

Belmonte, J. M., Clendenon, S. G., Oliveira, G. M., Swat, M. H., Greene, E. V., Jeyaraman, S., et al. (2016). Virtual-tissue computer simulations define the roles of cell adhesion and proliferation in the onset of kidney cystic disease. Mol. Biol. Cell. 27, 3673–3685. doi:10.1091/mbc.e16-01-0059

PubMed Abstract | CrossRef Full Text | Google Scholar

Brodland, G. W., Conte, V., Cranston, P. G., Veldhuis, J., Narasimhan, S., Hutson, M. S., et al. (2010). Video force microscopy reveals the mechanics of ventral furrow invagination in Drosophila. Proc. Natl. Acad. Sci. U. S. A. 107, 22111–22116. doi:10.1073/pnas.1006591107

PubMed Abstract | CrossRef Full Text | Google Scholar

Camley, B. A., Zhang, Y., Zhao, Y., Li, B., Ben-Jacob, E., Levine, H., et al. (2014). Polarity mechanisms such as contact inhibition of locomotion regulate persistent rotational motion of mammalian cells on micropatterns. Proc. Natl. Acad. Sci. U. S. A. 111, 14770–14775. doi:10.1073/pnas.1414498111

PubMed Abstract | CrossRef Full Text | Google Scholar

Chan, C. J., Costanzo, M., Ruiz-Herrero, T., Mönke, G., Petrie, R. J., Bergert, M., et al. (2019). Hydraulic control of mammalian embryo size and cell fate. Nature 571, 112–116. doi:10.1038/s41586-019-1309-x

PubMed Abstract | CrossRef Full Text | Google Scholar

Drasdo, D., Kree, R., and McCaskill, J. S. (1995). Monte Carlo approach to tissue-cell populations. Phys. Rev. E. 52, 6635–6657. doi:10.1103/physreve.52.6635

PubMed Abstract | CrossRef Full Text | Google Scholar

Dumortier, J. G., Verge-Serandour, M. L., Tortorelli, A. F., Mielke, A., de Plater, L., Turlier, H., et al. (2019). Hydraulic fracturing and active coarsening position the lumen of the mouse blastocyst. Science. 365, 465–468. doi:10.1126/science.aaw7709

PubMed Abstract | CrossRef Full Text | Google Scholar

Eiraku, M., Takata, N., Ishibashi, H., Kawada, M., Sakakura, E., Okuda, S., et al. (2011). Self-organizing optic-cup morphogenesis in three-dimensional culture. Nature 472, 51–56. doi:10.1038/nature09941

PubMed Abstract | CrossRef Full Text | Google Scholar

Fletcher, A. G., Osterfield, M., Baker, R. E., and Shvartsman, S. Y. (2014). Vertex models of epithelial morphogenesis. Biophys. J. 106, 2291–2304. doi:10.1016/j.bpj.2013.11.4498

PubMed Abstract | CrossRef Full Text | Google Scholar

Forgacs, G., and Newman, S. A. (2005). Biological physics of the developing embryo. Cambridge University Press.

Google Scholar

Guillot, C., and Lecuit, T. (2013). Mechanics of epithelial tissue homeostasis and morphogenesis. Science. 340, 1185–1189. doi:10.1126/science.1235249

PubMed Abstract | CrossRef Full Text | Google Scholar

Honda, H., Motosugi, N., Nagai, T., Tanemura, M., and Hiiragi, T. (2008a). Computer simulation of emerging asymmetry in the mouse blastocyst. Development. 135, 1407–1414. doi:10.1242/dev.014555

PubMed Abstract | CrossRef Full Text | Google Scholar

Honda, H., Nagai, T., and Tanemura, M. (2008b). Two different mechanisms of planar cell intercalation leading to tissue elongation. Dev. Dyn. 237, 1826–1836. doi:10.1002/dvdy.21609

PubMed Abstract | CrossRef Full Text | Google Scholar

Ishihara, S., and Sugimura, K. (2012). Bayesian inference of force dynamics during morphogenesis. J. Theor. Biol. 313, 201–211. doi:10.1016/j.jtbi.2012.08.017

PubMed Abstract | CrossRef Full Text | Google Scholar

Ishimoto, Y., and Morishita, Y. (2014). Bubbly vertex dynamics: a dynamical and geometrical model for epithelial tissues with curved cell shapes. Phys. Rev. E. 90, 1–23. doi:10.1103/PhysRevE.90.052711

PubMed Abstract | CrossRef Full Text | Google Scholar

Israelachvili, J. N. (2011). Intermolecular and surface forces. 3rd ed. Amsterdam: Elsevier Inc.

Google Scholar

Jewett, C. E., and Prekeris, R. (2018). Insane in the apical membrane: trafficking events mediating apicobasal epithelial polarity during tube morphogenesis. Traffic 19, 666–678. doi:10.1111/tra.12579

PubMed Abstract | CrossRef Full Text | Google Scholar

Kawagishi, R., Tahara, M., Sawada, K., Morishige, K., Sakata, M., Tasaka, K., et al. (2004). Na+/H+ exchanger-3 is involved in mouse blastocyst formation. J. Exp. Zool. Part A Comp. Exp. Biol. 301, 767–775. doi:10.1002/jez.a.90

PubMed Abstract | CrossRef Full Text | Google Scholar

Kondo, Y., Aoki, K., and Ishii, S. (2018). Inverse tissue mechanics of cell monolayer expansion. PLoS Comput. Biol. 14, 1–15. doi:10.1371/journal.pcbi.1006029

PubMed Abstract | CrossRef Full Text | Google Scholar

Koyama, H., Okumura, H., Ito, A. M., Nakamura, K., Otani, T., Kato, K., et al. (2023). Effective mechanical potential of cell–cell interaction explains three-dimensional morphologies during early embryogenesis. PLoS Comput. Biol. 19, e1011306. doi:10.1371/journal.pcbi.1011306

PubMed Abstract | CrossRef Full Text | Google Scholar

Koyama, H., Shi, D., and Fujimori, T. (2019). Biophysics in oviduct: planar cell polarity, cilia, epithelial fold and tube morphogenesis, egg dynamics. Biophys. Physicobiol. 16, 89–107. doi:10.2142/biophysico.16.0_89

PubMed Abstract | CrossRef Full Text | Google Scholar

Koyama, H., Shi, D., Suzuki, M., Ueno, N., Uemura, T., and Fujimori, T. (2016). Mechanical regulation of three-dimensional epithelial fold pattern formation in the mouse oviduct. Biophys. J. 111, 650–665. doi:10.1016/j.bpj.2016.06.032

PubMed Abstract | CrossRef Full Text | Google Scholar

Koyama, H., Umeda, T., Nakamura, K., Higuchi, T., and Kimura, A. (2012). A high-resolution shape fitting and simulation demonstrated equatorial cell surface softening during cytokinesis and its promotive role in cytokinesis. PLoS One 7, e31607. doi:10.1371/journal.pone.0031607

PubMed Abstract | CrossRef Full Text | Google Scholar

Krupinski, P., Chickarmane, V., and Peterson, C. (2011). Simulating the mammalian blastocyst - molecular and mechanical interactions pattern the embryo. PLoS. Comput. Biol. 7, e1001128. doi:10.1371/journal.pcbi.1001128

PubMed Abstract | CrossRef Full Text | Google Scholar

Kurotaki, Y., Hatta, K., Nakao, K., Nabeshima, Y.-I., and Fujimori, T. (2007). Blastocyst axis is specified independently of early cell lineage but aligns with the ZP shape. Science 316, 719–723. doi:10.1126/science.1138591

PubMed Abstract | CrossRef Full Text | Google Scholar

Latorre, E., Kale, S., Casares, L., Gómez-González, M., Uroz, M., Valon, L., et al. (2018). Active superelasticity in three-dimensional epithelia of controlled shape. Nature 563, 203–208. doi:10.1038/s41586-018-0671-4

PubMed Abstract | CrossRef Full Text | Google Scholar

Li, X., Balagam, R., He, T.-F., Lee, P. P., Igoshin, O. A., and Levine, H. (2017). On the mechanism of long-range orientational order of fibroblasts. Proc. Natl. Acad. Sci. 114, 8974–8979. doi:10.1073/pnas.1707210114

PubMed Abstract | CrossRef Full Text | Google Scholar

Maître, J.-L., Berthoumieux, H., Krens, S. F. G., Salbreux, G., Jülicher, F., Paluch, E., et al. (2012). Adhesion functions in cell sorting by mechanically coupling the cortices of adhering cells. Science. 338, 253–256. doi:10.1126/science.1225399

PubMed Abstract | CrossRef Full Text | Google Scholar

Maître, J.-L., Niwayama, R., Turlier, H., Nédélec, F., and Hiiragi, T. (2015). Pulsatile cell-autonomous contractility drives compaction in the mouse embryo. Nat. Cell. Biol. 17, 849–855. doi:10.1038/ncb3185

PubMed Abstract | CrossRef Full Text | Google Scholar

Merkel, M., and Manning, M. L. (2017). Using cell deformation and motion to predict forces and collective behavior in morphogenesis. Semin. Cell. Dev. Biol. 67, 161–169. doi:10.1016/j.semcdb.2016.07.029

PubMed Abstract | CrossRef Full Text | Google Scholar

Miyazaki, S., Otani, T., Sugihara, K., Fujimori, T., Furuse, M., and Miura, T. (2023). Mechanism of interdigitation formation at apical boundary of MDCK cell. iScience 26, 106594. doi:10.1016/j.isci.2023.106594

PubMed Abstract | CrossRef Full Text | Google Scholar

Miyoshi, H., Satoh, S. K., Yamada, E., and Hamaguchi, Y. (2006). Temporal change in local forces and total force all over the surface of the sea urchin egg during cytokinesis. Cell. Motil. Cytoskelet. 63, 208–221. doi:10.1002/cm.20118

PubMed Abstract | CrossRef Full Text | Google Scholar

Montesano, R., Schaller, G., and Orci, L. (1991). Induction of epithelial tubular morphogenesis in vitro by fibroblast-derived soluble factors. Cell. 66, 697–711. doi:10.1016/0092-8674(91)90115-f

PubMed Abstract | CrossRef Full Text | Google Scholar

Moriwaki, K., Tsukita, S., and Furuse, M. (2007). Tight junctions containing claudin 4 and 6 are essential for blastocyst formation in preimplantation mouse embryos. Dev. Biol. 312, 509–522. doi:10.1016/j.ydbio.2007.09.049

PubMed Abstract | CrossRef Full Text | Google Scholar

Niimura, S. (2003). Time-lapse videomicrographic analyses of contractions in mouse blastocysts. J. Reprod. Dev. 49, 413–423. doi:10.1262/jrd.49.413

PubMed Abstract | CrossRef Full Text | Google Scholar

Nikunen, P., Karttunen, M., and Vattulainen, I. (2003). How would you integrate the equations of motion in dissipative particle dynamics simulations? Comput. Phys. Commun. 153, 407–423. doi:10.1016/S0010-4655(03)00202-9

CrossRef Full Text | Google Scholar

Nissen, S. B., Perera, M., Gonzalez, J. M., Morgani, S. M., Jensen, M. H., Sneppen, K., et al. (2017). Four simple rules that are sufficient to generate the mammalian blastocyst. PLoS Biol. 15, e2000737. doi:10.1371/journal.pbio.2000737

PubMed Abstract | CrossRef Full Text | Google Scholar

Nissen, S. B., Rønhild, S., Trusina, A., and Sneppen, K. (2018). Theoretical tool bridging cell polarities with development of robust morphologies. Elife 7, e38407. doi:10.7554/elife.38407

PubMed Abstract | CrossRef Full Text | Google Scholar

Nonomura, M. (2012). Study on multicellular systems using a phase field model. PLoS One 7, 1–9. doi:10.1371/journal.pone.0033501

PubMed Abstract | CrossRef Full Text | Google Scholar

Okuda, S., Inoue, Y., and Adachi, T. (2015a). Three-dimensional vertex model for simulating multicellular morphogenesis. Biophys. Physicobiology 12, 13–20. doi:10.2142/biophysico.12.0_13

PubMed Abstract | CrossRef Full Text | Google Scholar

Okuda, S., Inoue, Y., Eiraku, M., Adachi, T., and Sasai, Y. (2015b). Vertex dynamics simulations of viscosity-dependent deformation during tissue morphogenesis. Biomech. Model. Mechanobiol. 14, 413–425. doi:10.1007/s10237-014-0613-5

PubMed Abstract | CrossRef Full Text | Google Scholar

Paramore, S. V., Trenado-Yuste, C., Sharan, R., Nelson, C. M., and Devenport, D. (2024). Vangl-dependent mesenchymal thinning shapes the distal lung during murine sacculation. Dev. Cell. 59, 1302–1315.e6. doi:10.1016/j.devcel.2024.03.010

PubMed Abstract | CrossRef Full Text | Google Scholar

Petridou, N. I., and Heisenberg, C.-P. (2019). Tissue rheology in embryonic organization. EMBO J. 38, e102497. doi:10.15252/embj.2019102497

PubMed Abstract | CrossRef Full Text | Google Scholar

Pettitt, B. M., and Rossky, P. J. (1986). Alkali halides in water: ion–solvent correlations and ion–ion potentials of mean force at infinite dilution. J. Chem. Phys. 84, 5836–5844. doi:10.1063/1.449894

CrossRef Full Text | Google Scholar

Podewitz, N., Delarue, M., and Elgeti, J. (2015). Tissue homeostasis: a tensile state. Epl. 109, 1–6. doi:10.1209/0295-5075/109/58005

CrossRef Full Text | Google Scholar

Rivron, N. C., Frias-Aldeguer, J., Vrij, E. J., Boisset, J. C., Korving, J., Vivié, J., et al. (2018). Blastocyst-like structures generated solely from stem cells. Nature. 557, 106–111. doi:10.1038/s41586-018-0051-0

PubMed Abstract | CrossRef Full Text | Google Scholar

Serra, D., Mayr, U., Boni, A., Lukonin, I., Rempfler, M., Meylan, L. C., et al. (2019). Self-organization and symmetry breaking in intestinal organoid development. Nature. 569, 66–72. doi:10.1038/s41586-019-1146-y

PubMed Abstract | CrossRef Full Text | Google Scholar

Solnica-krezel, L., and Sepich, D. S. (2012). Gastrulation: making and shaping germ layers. Annu. Rev. Cell. Dev. Biol. 28, 687–717. doi:10.1146/annurev-cellbio-092910-154043

PubMed Abstract | CrossRef Full Text | Google Scholar

Steinberg, M. S. (2007). Differential adhesion in morphogenesis: a modern view. Curr. Opin. Genet. Dev. 17, 1281–1286. doi:10.1016/j.gde.2007.05.002

PubMed Abstract | CrossRef Full Text | Google Scholar

Style, R. W., Boltyanskiy, R., German, G. K., Hyland, C., Macminn, C. W., Mertz, A. F., et al. (2014). Traction force microscopy in physics and biology. Soft Matter 10, 4047–4055. doi:10.1039/c4sm00264d

PubMed Abstract | CrossRef Full Text | Google Scholar

Tambe, D. T., Hardin, C. C., Angelini, T. E., Rajendran, K., Park, C. Y., Serra-Picamal, X., et al. (2011). Collective cell guidance by cooperative intercellular forces. Nat. Mat. 10, 469–475. doi:10.1038/nmat3025

PubMed Abstract | CrossRef Full Text | Google Scholar

Wang, A. Z., Ojakian, G. K., and Nelson, W. J. (1990a). Steps in the morphogenesis of a polarized epithelium I. Uncoupling the roles of cell-cell and cell-substratum contact in establishing plasma membrane polarity in multicellular epithelial (MDCK) cysts. J. Cell. Sci. 95, 137–151. doi:10.1242/jcs.95.1.137

PubMed Abstract | CrossRef Full Text | Google Scholar

Wang, A. Z., Ojakian, G. K., and Nelson, W. J. (1990b). Steps in the morphogenesis of a polarized epithelium II. Disassembly and assembly of plasma membrane domains during reversal of epithelial cell polarity in multicellular epithelial (MDCK) cysts. J. Cell. Sci. 95, 153–165. doi:10.1242/jcs.95.1.153

PubMed Abstract | CrossRef Full Text | Google Scholar

Wang, X., Zhang, Z., Tao, H., Liu, J., Hopyan, S., and Sun, Y. (2018). Characterizing inner pressure and stiffness of trophoblast and inner cell mass of blastocysts. Biophys. J. 115, 2443–2450. doi:10.1016/j.bpj.2018.11.008

PubMed Abstract | CrossRef Full Text | Google Scholar

Yong, J., Moeller, J., Hart, K. C., Ramallo, D., Vogel, V., Dunn, A. R., et al. (2015). Spatial distribution of cell–cell and cell–ECM adhesions regulates force balance while maintaining E-cadherin molecular tension in cell pairs. Mol. Biol. Cell. 26, 2456–2465. doi:10.1091/mbc.E14-12-1618

PubMed Abstract | CrossRef Full Text | Google Scholar