Computer simulation of emerging asymmetry in the mouse blastocyst

The mechanisms underlying early mammalian development remain controversial(Hiiragi et al., 2006; Rossant and Tam, 2004). Blastocyst morphogenesis is unique to mammals, resulting in the formation of the extraembryonic structures. The mammalian blastocyst displays an obvious asymmetry, with the inner cell mass (ICM) at one end and the blastocyst cavity surrounded by the trophectoderm (TE) at the other. The position of the ICM in the blastocyst, i.e. the embryonic pole, provides a basis for further embryonic patterning in that the embryonic-abembryonic (Em-Ab) axis should in principle give rise to the dorsoventral axis of the mouse embryo(Beddington and Robertson,1999; Rossant and Tam,2004). Central questions include how this apparently initial asymmetry is established in the early mouse embryo, and whether the underlying mechanism in this process is comparable to that operating in non-mammalian`model' organisms, in which localized (`prepatterned') determinants play an essential role. Whereas recent reports(Gardner, 1997; Gardner, 2001; Piotrowska et al., 2001; Piotrowska and Zernicka-Goetz,2001; Piotrowska-Nitsche et al., 2005; Torres-Padilla et al., 2007) claim the presence of `prepatterning' in the mouse egg and preimplantation embryos reminiscent of non-mammals, our studies(Hiiragi and Solter, 2004; Motosugi et al., 2005) and those by others (Alarcon and Marikawa,2003; Alarcon and Marikawa,2005; Chroscicka et al.,2004; Kurotaki et al.,2007; Louvet-Vallee et al.,2005) provide evidence supporting a `regulative' and `mechanical constraint' (Alarcon and Marikawa,2003; Kurotaki et al.,2007; Motosugi et al.,2005) model of blastocyst morphogenesis. In this model, the first embryonic polarity, the Em-Ab axis, is not established until the blastocyst cavity is localized at one end (the abembryonic pole), and its eventual position is determined by the spatial constraints imposed by the zona pellucida (ZP) in conjunction with the epithelial seal in the outer cell layer. This model does not depend on any localized determinants, leading us to apply mathematical and physical modeling not only to examine the feasibility of the mechanical constraint model, but also to elucidate the essential cellular properties necessary to achieve blastocyst morphogenesis.

In order to reduce the parameter number without loss of generality, we introduce a new length unit R₀ and rewrite Eq. (8) using new dimensionless quantities r_i′,▿ _i′(=∂/∂r_i′), S_f′, S_m′, V_α′, V_std′ and z_α′ as follows:

Thus, Eq. (8) takes the form:

in which the new dimensionless quantities are defined as follows:

We take R₀(=V_std^⅓) =1, so that V_std′=1. Eq. (11) lacks explicit parameters corresponding to η and σ_S. Thus, without loss of generality, we can describe cell behaviors using the dimensionless parametersσ _O′, σ_I′,κ _V′, κ_VI′,κ _I′, w′_bound,λ′, a′, b′ and c′. Below, cell motions are measured in terms of the new length unit R₀=V_std^⅓ and the new time unit η/σ_S, which are the characteristic length scale and time scale of the cell aggregate, respectively. Hereafter, we omit primes(′) on the rescaled quantities in Eq.(11).

which is equivalent to Eq. (13) and Eq. (16) in the Results.

In addition to the equations of motion, our model involves an elementary process of reconnecting neighboring vertices(Honda et al., 2004). When the length of an edge connecting two neighboring vertices becomes short (less than a critical length δ), the relationship of the neighboring vertices changes and the neighbors reconnect themselves(Fig. 1B).

Cell boundary surface areas (S_k, O_k,I_k), cell volumes (V_α), and the volume of the blastocyst cavity (V_I) are calculated by the vertex coordinates, providing the potential (U) and its partial derivatives with respect to x_i, y_i, and z_i (▿_i U). The Runge-Kutta method(Ohno and Isoda, 1977) with step size h is applied to numerically solve the simultaneous equations of motion (Eq. 12). Thus, we obtain the movement of all vertices. After each Runge-Kutta time step, we apply the elementary process of reconnection (Honda et al.,2004) (Fig. 1B) to edges shorter than the critical length δ.

An initial structure of a cell aggregate consisting of 31 or 40 cells(n=31, 40) is produced according to the method described(Honda et al., 2004). To make an internal space (blastocyst cavity) in the cell aggregate, we choose a vertex in the cell aggregate and replace the vertex by a small tetrahedron(Fig. 1C). To simulate development to the expanded blastocyst, the volume of the tetrahedron(I_std) is made to increase to half the initial volume of the cell aggregate (I_std=15.5 or 20) during t=0-500 (Fig. 1D),while keeping the same conditions for the cell aggregate. In this way, the total volume of the cells slightly decreases (e.g. -11.8%), while the entire embryo volume significantly increases (e.g. +45.7%) accompanied by the cavity expansion (e.g. 0 to 19.09). This increase in the entire embryo volume during the blastocyst cavity formation is confirmed by time-lapse recording of preimplantation mouse embryos (data not shown) (see Motosugi et al., 2005). This period of linear increase in the cavity volume is followed by a plateau period(Fig. 1D; t=500-3500),which is technically necessary for computer simulation in order to obtain the stabilized structure with the minimum potential energy after the forced change in condition.

At present, accurate values for parameters related to the cell properties are not available. Thus, the present study investigates the general behavior of polyhedral cells in typical cell aggregates using Eq. (12) with dimensionless parameter values. The parameters σ_O,σ _I, κ_V, κ_VI,κ _I and w_bound are values relative to the coefficient of interface energy between the cells (cell-cell boundary energy),σ _S. The smaller the coefficient of the cell volume elasticity(κ_V), the easier the change in the cell volume and the more unstable the system becomes; when κ_V is larger, the shape of the cell aggregate hardly changes. The κ_V(κ_V=12) in this study allows both flexibility and rigidity. The coefficients of surface energy are chosen to be smaller(σ_o=σ_I=0.7) than that of interface energy between the cells, σ_S, so that the internal space smoothly expands. The coefficient of volume elasticity of the internal space,κ _VI, is considered to be smaller than the coefficient of cell volume elasticity (κ_V). However, when κ_VI is too small, the internal space does not expand although the standard volume of the internal space (I_std) is forced to increase in the simulations, leading us to choose κ_VI=6. In addition, we introduce the coefficient of surface elasticity of the cells facing the internal space, κ_I. The cell surfaces facing the internal space are expected to minimize the elastic energy (deviation from the relaxed surfaces). Without this elasticity, some vertices produce complicated tessellation patterns and the process freezes. We also confirm that the surface shows elastic properties when the relaxed area I_std is between 0.2 and 1. The process freezes when the elastic constant of the surface of polygons facing the blastocyst cavityκ _I is set at 0.2, suggesting that it is too large and leading us to choose I_std=0.433 and κ_I=0.1. The last term on the right side of Eq. (12) expresses hindrance of outward motion of a vertex from the ZP. Its parameters w_bound andλ allow a wide range of values. When w_bound is too large (e.g. w_bound=50), the space between the outer surface of the blastocyst and the ZP becomes too large. λ is a parameter indicating the slope of the logistic function, which imitates a step function. In this simulation we use w_bound=1 andλ=50. The critical length for reconnection, δ, is the lower cut-off length in the model and should be small for a detailed description,although a value too small in complicated tessellation patterns freezes the process. Here we use δ=0.05, where the relaxed cell volume(V_std) is 1. The step size h for the Runge-Kutta method is 0.005 to enable smooth changes in the process.

Computer programs are constructed in FORTRAN 90 running on a workstation(Octane-R12000, Silicon Graphics, USA) at Hyogo University (Kakogawa, Hyogo,Japan) or on a digital super computer (SGI Altix 3700, SGI, Japan) at the Institute of Statistical Mathematics (Tokyo, Japan). Mathematica (ver. 3.0,Wolfram Research) is used for analysis and drawing.

Collection of the mouse embryo and time-lapse recordings were carried out essentially as described (Hiiragi and Solter, 2004; Motosugi et al.,2005), except that embryos were photographed only once at the 2-cell stage, and were then time-lapse recorded from the morula to the blastocyst stage to minimize light exposure (see Fig. 3C and Movie 1 in the supplementary material). The ZP was partially digested by incubation with 0.5%pronase (Sigma, P8811) in H-KSOM-AA (KSOM with amino acids and 21 mM Hepes) at 37°C in 6% CO₂ for a few minutes with repeated microscopic observations, followed by several washes with H-KSOM-AA.

Our computer simulation is performed in order to obtain, under given conditions, the most stable pattern of the cell aggregate that is equivalent to the blastocyst in vivo, based on the following three assumptions. (Ass. 1)Dynamic cell behavior can be simulated by a 3D vertex dynamics model, as published previously (Honda et al.,2004). Briefly (see Materials and methods for details), an aggregate of multiple cells is considered as a 3D space tessellation composed of polyhedra without gaps or overlaps (Fig. 1A). The interfaces (cell surface) and volumes are expressed by 3D vertex coordinates, and the position of the vertex is rearranged according to an equation of motion so that the total potential energy of the cell aggregate, composed mainly of surface energy and elastic energy, should be minimal. (Ass. 2) All cells have equivalent mechanical properties. (Ass. 3) A cavity is created inside the cell aggregate by replacing one vertex with a small polyhedron with equivalent surface properties.

In this study, an aggregate of 40 cells is considered to simulate the morphology of the mid- to late-stage expanded blastocyst(Chisholm et al., 1985; Dietrich and Hiiragi, 2007; Gueth-Hallonet et al., 1993; Smith and McLaren, 1977). The number of cells is constant in this simulation, and the calculation is carried out to obtain the most stable (i.e. eventual) structure of the cell aggregate that has the minimum potential energy. An additional simulation is carried out under exactly the same conditions with an aggregate of 31 cells, equivalent to the nascent blastocyst stage (Chisholm et al., 1985; Dietrich and Hiiragi, 2007; Smith and McLaren, 1977), which in essence produces the same results (see Fig. S1 in the supplementary material) as described below. Therefore, the following study is presented only for the aggregate of 40 cells. A cavity is created in the aggregate by replacing one vertex with a small tetrahedron(Fig. 1C; the length of the six edges is 0.1). The volume of the tetrahedron (V_Istd),which subsequently forms the polyhedron, is forced to increase during t=0-500 (Fig. 1D) to half the total volume of the initial cell aggregate(V_Istd=20), reaching the expanded blastocyst stage. The simulation includes an additional phase, after the increase in the volume(V_Istd) stops (t=500), until the actual cavity reaches the maximal volume (which is not necessarily equivalent to that at the t=500 time point and is usually delayed owing to the elastic nature of the structure; see Fig. 4A)and the structure reaches a stable state, depending on the conditions. The blastocyst cavity is an intercellular space in the mammalian embryo, created by the directional fluid secretion of the outer layer cells(Aziz and Alexandre, 1991; Calarco and Brown, 1969; Wiley and Eglitis, 1981). This secretion most likely occurs in essentially all outer blastomeres with apicobasal epithelial polarity, thus initiating multiple small cavitations(Motosugi et al., 2005), and these multiple cavities coalesce with each other to form one large cavity. The current simulation calculates the most stable structure with the minimal potential energy composed of 40 cells and one cavity that should correspond to the eventual blastocyst morphology (see Discussion).

This modified Eq. (15) allows the computer simulation to successfully produce a structure recapitulating the mouse blastocyst(Fig. 2A,B). In fact, this additional elastic energy of the cell surface facing the blastocyst cavity has the effect of minimizing deviation from the relaxed cell surface area, and suggests that there should be some structure lining the blastocyst cavity with a tendency to keep its surface area closer to the relaxed condition. This result thus suggests the presence of the basement membrane at the basal side of the TE facing the blastocyst cavity, and/or possibly the difference in cell surface elastic properties or adhesive properties between the TE/primitive endoderm and the epiblast. Indeed, we find past(Nadijcka and Hillman, 1974; Thorsteinsdottir, 1992) and recent (Klaffky et al., 2001)studies clearly showing the presence of the basement membrane at the surface of the TE cells facing the blastocyst cavity, and the difference in its molecular components between the TE and the ICM. Thus, the basement membrane in this region may account for the altered surface properties required by the model and might be essential for maintaining the integrity of the cavity during blastocyst morphogenesis.

Under this simulation condition (Fig. 2A,B), however, the embryonic axis, i.e. the location of the inner cell cluster, is never fixed and remains variable during the process of calculation (data not shown), indicating that there is no certain orientation in which potential energy U in Eq. (15) becomes minimal. It is important to note, however, that there is no problem in forming the blastocyst structure itself, consistent with the fact that the ZP is not essential for the mouse embryo to develop to the blastocyst stage(Motosugi et al., 2005; Kurotaki et al., 2007). This simulation result rather predicts that, without the ZP, the embryonic axis of the mouse blastocyst can be oriented to any direction at an equal probability in terms of the potential energy (see below). This result leads us to include an additional element in the equation: (Ass. 5) the embryo is enclosed by a capsule (corresponding to the ZP in vivo) that constrains movement of the cells. This is expressed in the potential energy U as a restriction energy by a boundary, f(r_i), in which f(r_i) is a quadratic function, i.e. f(r_i)=1 expresses a quadratic surface (e.g. sphere or ellipsoid). Thus, Eq. (15) now takes the form:

In an initial attempt in which the surrounding capsule is assumed to be spherical, the axis of the embryo again remains variable(Fig. 2C). However, when the capsule is ellipsoidal (1.2: 1) with its long diameter vertical(Fig. 2D) or parallel to the z-axis (see Fig. 4),the position of the inner cell cluster, i.e. the axis of the blastocyst,always coincides eventually with the long axis of the ellipsoidal capsule. This indicates that the potential energy U of the cell aggregate becomes minimal when the position of the inner cell cluster is localized at one end of the long axis of the ellipsoidal capsule, thus eventually stabilizing the embryonic axis parallel to the long axis of the ellipsoidal ZP. This is independent of the initial position of cavity formation,regardless of whether it is initiated from a center (black circle in Fig. 2A; the outcome in Fig. 2B-D) or from a peripheral point (blue circle in Fig. 2A;the outcome in Fig. 2E,F).

We then decided to test experimentally these two predictions in two situations (when the ZP is absent or spherical, and when the ZP is ellipsoidal) by using time-lapse recording of the mouse preimplantation embryo(Fig. 3). The ratio of the longest to the shortest diameters of the ZP was on average 1.07 in 187 embryos analyzed at the 2-cell stage, and was essentially preserved until the mid-blastocyst stage (see also Motosugi et al., 2005). In 13 embryos in which the ratio was approximately 1.2, the Em-Ab boundary tended to be aligned with the shortest ZP axis(Fig. 3A,B, Embryo Group I; P=0.018, χ² analysis), thus confirming one of the two predictions of the computer simulation (i.e. when the ZP is ellipsoidal; see Fig. 2D-F and Fig. 4) as well as the mechanical constraint model for blastocyst morphogenesis proposed earlier(Alarcon and Marikawa, 2003; Kurotaki et al., 2007; Motosugi et al., 2005). The other prediction (when the ZP is absent or spherical, i.e. when there is no mechanical influence or bias by the ZP) was again experimentally verified by partially digesting the ZP to create a larger and spherical ZP that would exert essentially no mechanical or spatial constraints on the embryo(Fig. 3C and see Movie 1 in the supplementary material; showing the same experimental embryos). Note the spherical shape of the blastomeres of the 2-cell embryos(Fig. 3C, 47:10 post-hCG),which indicates the absence of mechanical constraints [see Motosugi et al.(Motosugi et al., 2005) for the ellipsoidal blastomeres of the normal 2-cell embryo]. In these experimental embryos with a spherical ZP, the tilt angle between the shortest ZP axis and the Em-Ab boundary (Fig. 3A) is randomly distributed(Fig. 3B, Embryo Group II; n=55, P=0.64, χ² analysis), which confirms the other prediction of the computer simulation when the ZP is spherical (see Fig. 2C).

The calculation process to find a stable state is exemplified in Fig. 4A (and see Movie 2 in the supplementary material). In this example, the aggregate is surrounded and restricted by an ellipsoidal (1.2:1) capsule, corresponding to an ellipsoidal ZP in vivo. A vertex (solid circle in Fig. 4A at t=0), replaced by a tetrahedron (see Fig. 1C), is enlarged until its volume reaches half the initial total volume (at t=500; see Fig. 1D). The embryonic axis(Em-Ab axis) keeps changing with respect to the outer capsule until the embryonic pole, the position of the ICM, is localized at one end of the long axis of the ellipsoidal capsule (t=2000). The structure is then`stabilized', remaining unchanged thereafter at least for a substantial time(t=2000-3500) equivalent to the time required to stabilize the structure (t=0-2000). Tracing of cells(Fig. 4A, marked in green, red,orange and blue; and four cells surrounding the orange cell at t=1000) further demonstrate that cells and the cavity gradually change their position relative to each other during the process; those cells marked either with a black circle or a black square at some point acquire a cell at their interface with the orange cell, and those marked with red circle and square disappear from the plane of observation (see Discussion). The enlarging cavity eventually becomes surrounded by an outer single-cell layer,with its position stabilized at one end of the long axis, while a cluster of cells forms on the opposite side (Fig. 4B,C and see Movie 3 in the supplementary material). Thus, the simulation recapitulates well the blastocyst morphology in vivo, with the outside cells corresponding to the TE and the inner population to the ICM.

Based on five simulations under various conditions with respect to the direction of the ellipsoidal ZP long axis or initial position of the cavity,the eventual cell numbers comprising the ICM and the TE are predicted to be on average 11.2 and 28.8, respectively. These numbers and their ratio are indeed verified by past and recent studies (Barlow et al., 1972; Dietrich and Hiiragi, 2007; Kelly et al.,1978), corresponding to those of the in vivo embryo. Intriguingly,despite the assumed equivalence in mechanical properties of the 40 cells in various simulations (Ass. 2), enlargement of the cavity enhances the difference in volume of the inner versus outer cells(Fig. 5), with the inner cell volume stabilizing at significantly lower levels (0.702±0.0088) than that of the outer cells (0.744±0.0136, average cell volume±s.d.; P<0.001, Student's t-test). This prediction based on the computer simulations is again confirmed by a recent study(Aiken et al., 2004)demonstrating that in vivo the average cell volume of the outer cells is significantly larger than that of the inner cells.

Several independent studies have recently suggested the absence of prepatterning in the mouse embryo until the blastocyst stage, and these led to the mechanical constraint model as a mechanism underlying Em-Ab axis formation in the mouse blastocyst (Alarcon and Marikawa, 2003; Kurotaki et al., 2007; Motosugi et al.,2005). The computer simulation in the present study first confirms that the mechanical constraint imposed on the embryo is indeed sufficient to orient the blastocyst axis. Furthermore, it illustrates the cellular mechanical properties minimally required for creating the blastocyst morphology, thereby allowing us to gain a deep insight into blastocyst morphogenesis from a mechanical point of view.

When the dynamics of the cell aggregate are restricted by the capsule of a deformed sphere, equivalent to the ZP in vivo, the inner cell cluster (ICM)consistently localizes to one end of the long axis of the ellipsoidal capsule. The simulated blastocyst contains the ICM formed by 11 or 12 cells, thus recapitulating well its composition in vivo. During the calculation process,most of the outside cells contribute to the TE, while the inside cells contribute to the ICM, except for one cell in this simulation that moves from inside to the outside TE (data not shown). However, a quantitative description of the contribution of cells to the ICM or TE awaits further studies with more-complex simulations (see below). The volume of the inner cells becomes significantly smaller than that of the outer cells during the simulation process, again similar to the difference observed in the blastocyst in vivo and despite the assumption of component cell equivalence in the simulations. Cellular polarization at the 8-cell stage and subsequent asymmetric division play a major role in the generation of asymmetry between inside and outside populations (Dietrich and Hiiragi,2007; Johnson and McConnell,2004; Ralston and Rossant,2008). The difference in cell size in those populations, for example, can be initiated by asymmetric divisions, and may be further enhanced by mechanical constraints, as seen in this study. Thus, the inside-outside asymmetry essential for blastocyst morphogenesis(Tarkowski and Wroblewska,1967) may emerge autonomously. This idea contrasts with the major role assigned to localized determinants in developmental mechanisms in non-mammalian `model' organisms. Whereas some reports(Niwa et al., 2005; Smith, 2005; Strumpf et al., 2005) claim a crucial role for the differential expression of Cdx2 and Oct4 (Pou5f1)molecules in specifying the fate of TE and ICM, recent studies, consistent with our present one, suggest that the initial difference might emerge as a result of cellular polarization (Ralston and Rossant, 2008), possibly includes stochastic mechanisms in the subsequent asymmetric divisions (Dietrich and Hiiragi, 2007), and is ultimately stabilized to accommodate intrinsic (gene expression pattern) and extrinsic (mechanical and structural integrity) cues.

Several issues remain to be addressed, largely owing to the limitations of current knowledge and simulation techniques. First, this simulation is not designed to recapitulate the developmental process of blastocyst morphogenesis because it lacks cell division and initial multiple cavities. This study allows us, however, to conclude that for an aggregate of 40 intrinsically equivalent cells with one cavity that has half the entire cell volume(conditions equivalent to the embryo at the mid- to late-blastocyst stage),the calculated structure has the minimal potential energy, regardless of the actual developmental process. Namely, the present computer simulation predicts that, even if the total cell number increases and multiple cavities are created during the developmental process, the eventual expanded blastocyst(composed of around 40 cells and most likely one cavity of half the entire cell volume) will have the calculated structure in terms of mechanical properties. Based solely on the mechanical properties, it is most likely that such an aggregate will eventually form a structure containing an inner cluster of cells localized to one end of its long axis, if the surrounding capsule is ellipsoidal. The process of the calculation (as shown in Fig. 4A) may thus appear less dynamic than development leading to blastocyst, as we observe only minor rearrangement of the cavity and the surrounding cells. The complete simulation to recapitulate blastocyst morphogenesis would certainly require integration of cell division and multiple cavity formation to provide dynamic aspects of the developmental process. Secondly, the TE:ICM cell volume ratio is 1.06 in the present simulation, whereas the actual cell volume ratio in the blastocyst in one study was reported as being 1.67(Aiken et al., 2004). This discrepancy might be partly due to less heterogeneity in the cellular population of the present simulation. Since embryonic cleavage in the mouse preimplantation stage is not precisely synchronous, the embryo at the blastocyst stage contains cells of various sizes and there may be dynamic changes in their relative position. Our simulation lacks this heterogeneous aspect, focusing on the eventual structure of the simplest population. Furthermore, the present simulations consider the cell membrane as a flat interface, without curvature, which may additionally account for the reduced difference in the TE:ICM cell volumes in the simulation. Thirdly, we assigned the surface energy coefficient of the outer surface and of the blastocyst cavity surface (σ_O and σ_I, respectively) at 0.7, and that of the interface energy between the neighboring cells at 1.0,based on the assumption that the surface tension of the single-membrane layer should be less than that of the double layer. More-complex simulations await measurement of the surface tension of the actual cell membrane in different regions of the embryo.

The present simulation also provides an insight into the mechanism by which the inside cells form a distinct cluster (ICM) instead of a continuous cell layer beneath the TE. This would indeed be the case only if the entire structure of the embryo were precisely radially symmetrical, and cavity formation were initiated at the very center of such a structure. However, in biological systems, such absolute precision is essentially impossible, and once there is the slightest deviation from the central location, the cavity never returns to the center but stabilizes at the periphery, forming a cluster of inside cells, the ICM. Taken together, the present findings using computer simulation reveal several unique features of early mammalian development, and will provide a conceptual basis for understanding the molecular mechanism of early embryonic patterning.

We thank R. Cassada, J.-E. Dietrich, B. Doak, M. Hoffman, H. Kirk and D. Solter for critical reading of the manuscript; T.H. is especially grateful to D. Solter for his generous support and continuous encouragement. This work was partially supported by a Japan JSPS Grant-in-Aid for Scientific Research (C)(to H.H.), and by the German Research Foundation, Priority Program 1109, and the Lalor Foundation (to T.H.).