The mechanism of embryonic polarity establishment in mammals has long been controversial. Whereas some claim prepatterning in the egg, we recently presented evidence that mouse embryonic polarity is not established until blastocyst and proposed the mechanical constraint model. Here we apply computer simulation to clarify the minimal cellular properties required for this morphology. The simulation is based on three assumptions: (1) behavior of cell aggregates is simulated by a 3D vertex dynamics model; (2) all cells have equivalent mechanical properties; (3) an inner cavity with equivalent surface properties is gradually enlarged. However, an initial attempt reveals a requirement for an additional assumption: (4) the surface of the cavity is firmer than intercellular surfaces, suggesting the presence of a basement membrane lining the blastocyst cavity, which is indeed confirmed by published data. The simulation thus successfully produces a structure recapitulating the mouse blastocyst. The axis of the blastocyst, however, remains variable,leading us to an additional assumption: (5) the aggregate is enclosed by a capsule, equivalent to the zona pellucida in vivo. Whereas a spherical capsule does not stabilize the blastocyst axis, an ellipsoidal capsule eventually orients the axis in accordance with its longest diameter. These predictions are experimentally verified by time-lapse recordings of mouse embryos. During simulation, equivalent cells form two distinct populations composed of smaller inner cells and larger outer cells. These results reveal a unique feature of early mammalian development: an asymmetry may emerge autonomously in an equivalent population with no need for a priori intrinsic differences.

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.

Three-dimensional vertex dynamics cell model

In this model, a tissue composed of multiple cells is considered as a 3D-space tessellation consisting of polyhedra without gaps or overlaps. The interface and volume of the cells are expressed as x, y and z coordinates of the vertices(Fig. 1A). Spatial relationships between neighboring vertices are defined by surrounding cells(Honda et al., 2004). The vertices obey the equation of motion:
\[\ {\eta}{\ }d\mathbf{\mathit{r}}_{i}{/}dt=-{\nabla}_{i}U{\ }(i=1,{\ldots},n_{v}),\]
where ri is a 3D-positional vector of vertex i, and ▿i the nabla differential operator. The left side of Eq. (1) represents a viscous drag force proportional to the vertex velocity dri/dt with a positive constant η (an analog of the coefficient of viscosity). Vertices do not have mass (inertia), so that the motion of the vertices and cells is completely damped. The right side of Eq. (1) represents a potential force(driving force), i.e. minus the gradient of the potential U. The potential U includes various terms related to cell surface where ri is a 3D-positional vector of vertex i, and ▿i the nabla differential operator. The left side of Eq. (1) represents a viscous drag force proportional to the vertex velocity dri/dt with a positive constantη (an analog of the coefficient of viscosity). Vertices do not have mass(inertia), so that the motion of the vertices and cells is completely damped. The right side of Eq. (1) represents a potential force (driving force), i.e. minus the gradient of the potential U. The potential Uincludes various terms related to cell surface area, cell volume, and potential energy under restriction of the ZP, that are all expressed by vertex coordinates. Hence, U is a function of the vertex positions.
In the present study, the potential U contains terms of surface energy (US), elastic energy (UEV, UEI) and boundary energy due to restriction by the ZP(Ubound):
\[\ U=U_{\mathrm{S}}+U_{\mathrm{EV}}+U_{\mathrm{EI}}+U_{\mathrm{bound}}.\]
The potential US denotes the total surface energy of the cells:
\[\ U_{\mathrm{S}}={\sigma}_{\mathrm{S}}{\Sigma}^{n\mathrm{S}}_{k}S_{k}+{\sigma}_{\mathrm{O}}{\Sigma}^{n\mathrm{O}}_{k}O_{k}+{\sigma}_{\mathrm{I}}{\Sigma}^{n\mathrm{I}}_{k}I_{k}.\]
The first two terms in Eq. (3) represent the interface energy between neighboring cells and the surface energy between cells and the outside, where nS and nO are the numbers of polygons facing an adjacent cell and facing the outside, respectively. Sk is the surface area of a polygon k facing adjacent cells and σS is its interface energy per unit area. Ok is the surface area of a polygon k facing the outside and σO is its surface energy per unit area. The third term in Eq. (3) represents the interface energy between cells and the internal cavity (blastocyst cavity), where nI is the number of polygons facing the blastocyst cavity. Ik is the surface area of a polygon k facing the blastocyst cavity andσ I is its interface energy per unit area. The potential UEV contains two terms, the elastic energy of compression and expansion of the cells and that of the blastocyst cavity:
\[\ U_{\mathrm{EV}}={\kappa}_{V}{\Sigma}^{n}_{{\alpha}}(V_{{\alpha}}-V_{\mathrm{std}})^{2}+{\kappa}_{\mathrm{VI}}(V_{\mathrm{I}}-V_{\mathrm{Istd}})^{2},\]
where κV and κVI are the elastic constant of volumes of the cells and the blastocyst cavity, respectively, and nis the total number of cells. Vα and VI are the volumes of cell α and of the blastocyst cavity. Vstd and VIstd are the volumes of relaxed cells and of the relaxed cavity. Vstd is defined as the average volume of all cells. Eq. (4) imposes a constraint on the volumes of the cells and the blastocyst cavity. The potential UEI denotes the elastic energy of compression and expansion of the surface area of a polygon k facing the blastocyst cavity:
\[\ U_{\mathrm{EI}}={\kappa}_{\mathrm{I}}{\Sigma}^{n\mathrm{I}}_{k}(I_{k}-I_{\mathrm{std}})^{2},\]
where κI is the elastic constant of the surface of a polygon k facing the blastocyst cavity. Istd is the relaxed surface area of the cells facing the blastocyst cavity. Eq. (5)provides the cell surface constraint of a polygon k facing the blastocyst cavity. Restriction of the embryo by the ZP is described as an additional term Ubound:
\[\ U_{\mathrm{bound}}=w_{\mathrm{bound}}{\Sigma}^{n\mathrm{v}}_{{\iota}}1{/}[1+\mathrm{exp}(-{\lambda}f(\mathbf{\mathit{r}}_{i}))]\]
in which ri=(xi,yi, zi) and
\[\ f(\mathbf{\mathit{r}}_{i})=(x_{i}-x_{0})^{2}{/}a^{2}+(y_{i}-y_{0})^{2}{/}b^{2}+(z_{i}-z_{0})^{2}{/}c^{2}-1.\]
Ubound expresses the hindrance of outward motion of a vertex from the ellipsoid [the diameters of its main axes are a, b, cand its center is (x0, y0, z0)]. When a point riis outside or inside the ellipsoid, f(ri) is positive or negative,respectively. Vertex i moves without restriction when it is inside the ellipsoid (f(ri)<0),whereas it produces substantial potential energy when it crosses the ellipse boundary (f(ri)>0). Instead of a step function, i.e. {=1, f(ri)>0; =0, f(ri)<0}, we use a continuous analytic function, logistic function 1/[1+exp{-λ f(ri)}], where λ/4 is the slope of the logistic function at f(ri)=0, and wbound is the strength of the potential. Thus, Eq. (1)takes the form:

In order to reduce the parameter number without loss of generality, we introduce a new length unit R0 and rewrite Eq. (8) using new dimensionless quantities ri′,▿ i′(=∂/∂ri′), Sf′, Sm′, Vα′, Vstd′ and zα′ as follows:

Thus, Eq. (8) takes the form:

in which the new dimensionless quantities are defined as follows:

We take R0(=Vstd) =1, so that Vstd′=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′, wbound,λ′, a′, b′ and c′. Below, cell motions are measured in terms of the new length unit R0=Vstd 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).

Computer simulations

Vertex dynamics

Cell boundary surface areas (Sk, Ok,Ik), cell volumes (Vα), and the volume of the blastocyst cavity (VI) are calculated by the vertex coordinates, providing the potential (U) and its partial derivatives with respect to xi, yi, and zi (▿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 δ.

Initial structure of the cell aggregate and its development in computer simulations

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(Istd) is made to increase to half the initial volume of the cell aggregate (Istd=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.

Parameter values used in the simulations

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 σOI, κV, κVII and wbound 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 κVV=12) in this study allows both flexibility and rigidity. The coefficients of surface energy are chosen to be smaller(σoI=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 (Istd) 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 Istd 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 Istd=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 wbound andλ allow a wide range of values. When wbound is too large (e.g. wbound=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 wbound=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(Vstd) is 1. The step size h for the Runge-Kutta method is 0.005 to enable smooth changes in the process.

Fig. 1.

Model and elements of the computer simulation. (A) Polyhedral cells and polygonal surfaces are defined by vertices. (B) Reconnection of neighboring vertices. (C) Formation of a tetrahedral intercellular space at a vertex, corresponding to the blastocyst cavity. (D) Process of expanding the blastocyst cavity in the computer simulation. Volume of the intercellular space (VIstd) is forced to increase linearly until t=500.

Fig. 1.

Model and elements of the computer simulation. (A) Polyhedral cells and polygonal surfaces are defined by vertices. (B) Reconnection of neighboring vertices. (C) Formation of a tetrahedral intercellular space at a vertex, corresponding to the blastocyst cavity. (D) Process of expanding the blastocyst cavity in the computer simulation. Volume of the intercellular space (VIstd) is forced to increase linearly until t=500.

Computation

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.

Partial digestion of the zona pellucida of the mouse embryo

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% CO2 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.

Accordingly, the potential energy of the aggregate, U, consists of three surface energy terms and two elastic energy terms. These terms are expressed by the positional vectors ri(i=1,..., nv) of nv vertices of the polyhedra. The vertex obeys the equation of motion:
\[\ d\mathbf{\mathit{r}}_{i}{/}dt=-{\nabla}_{i}U,\]
which is expressed using dimensionless variables, and▿ i is a nabla differential operator with respect to ri. The vertices move according to Eq.(13) so that the potential energy U becomes smaller(Honda et al., 2004). In addition to the equation of motion, our model involves an elementary process of reconnection of neighboring vertices(Honda et al., 2004). When the distance of the two neighboring vertices is shorter than a critical length,δ, their relationship changes, including reconnection(Fig. 1B).
The potential energy U is:
\[\ U={\Sigma}_{k}S_{k}+{\sigma}_{\mathrm{O}}{\Sigma}_{k}O_{k}+{\sigma}_{\mathrm{I}}{\Sigma}_{k}I_{k}+{\kappa}_{\mathrm{V}}{\Sigma}_{{\alpha}}(V_{{\alpha}}-1)^{2}+{\kappa}_{\mathrm{VI}}(V_{\mathrm{I}}-V_{\mathrm{Istd}})^{2}.\]
The first three terms are the surface energy of the interfaces between neighboring cells (Sk), an outer surface energy of the interfaces between cells and the externals (Ok), and the inner surface energy of the interfaces between a blastocyst cavity and cells(Ik), respectively. The next two terms are the elastic energies accompanied by a change in volume of the individual cells(Vα; the volume of the relaxed cell is assumed to be 1) and of the blastocyst cavity (VI; the volume of the relaxed blastocyst cavity is VIstd). Parameters,σ O, σIOI), κV andκ VI are weights applied to the respective terms.

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 (VIstd),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(VIstd=20), reaching the expanded blastocyst stage. The simulation includes an additional phase, after the increase in the volume(VIstd) 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).

Fig. 2.

Mouse blastocyst morphology under various conditions. The initial cell aggregate (A) and the simulated blastocyst at t=2000 (B-F)are shown in cross-sectional views of xz-(left column) and yz-(right column)planes. (A) The initial cell aggregate, in which either a central (black circle) or peripheral (blue circle) vertex is replaced by a tetrahedron(t=0), followed by enlargement under various conditions specified in B-F. (B) An aggregate without zona pellucida (ZP) and with cavitation initiated from a central vertex. (C) An aggregate enclosed with a spherical ZP with cavitation initiated from a central point. (D) An aggregate enclosed with the ellipsoidal (long y-axis) ZP and cavitation from the center. (E)The ellipsoidal (long z-axis) ZP with cavitation initiated from a peripheral vertex. (F) The ellipsoidal (long y-axis) ZP with cavitation from a peripheral point.

Fig. 2.

Mouse blastocyst morphology under various conditions. The initial cell aggregate (A) and the simulated blastocyst at t=2000 (B-F)are shown in cross-sectional views of xz-(left column) and yz-(right column)planes. (A) The initial cell aggregate, in which either a central (black circle) or peripheral (blue circle) vertex is replaced by a tetrahedron(t=0), followed by enlargement under various conditions specified in B-F. (B) An aggregate without zona pellucida (ZP) and with cavitation initiated from a central vertex. (C) An aggregate enclosed with a spherical ZP with cavitation initiated from a central point. (D) An aggregate enclosed with the ellipsoidal (long y-axis) ZP and cavitation from the center. (E)The ellipsoidal (long z-axis) ZP with cavitation initiated from a peripheral vertex. (F) The ellipsoidal (long y-axis) ZP with cavitation from a peripheral point.

A computer simulation is carried out based on the above three assumptions,i.e. this study is initiated from the simplest assumption. However, this results in failure, with the computer calculation disturbing tessellation and freezing during the simulation process. Several further attempts led us to additional complexity, by modifying (Ass. 3) and taking into consideration a difference in mechanical properties between the cell and the cavity: (Ass. 4)there is some component lining the cavity surface to make it firmer and more stabilized than the interface between cells. In practice, the additional elastic energy term of the cell surface facing the blastocyst cavity,κ Ik(Ik-Istd)2, is included in Eq. (14),in which Ik is the surface area of the cells facing the blastocyst cavity, Istd is the surface area in relaxed condition, and κI is an elastic constant of the surface area. Thus, the potential energy U takes the form:
\[\ U={\Sigma}_{k}S_{k}+{\sigma}_{\mathrm{O}}{\Sigma}_{k}O_{k}+{\sigma}_{\mathrm{I}}{\Sigma}_{k}I_{k}+{\kappa}_{\mathrm{V}}{\Sigma}_{{\alpha}}(V_{{\alpha}}-1)^{2}+{\kappa}_{\mathrm{VI}}(V_{\mathrm{I}}-V_{\mathrm{Istd}})^{2}+{\kappa}_{\mathrm{I}}{\Sigma}_{k}(I_{k}-I_{\mathrm{std}})^{2}.\]

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(ri), in which f(ri) is a quadratic function, i.e. f(ri)=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).

Fig. 3.

Experimental verification of the predictions based on the computer simulation. (A) Scheme illustrating the tilt angle (dark green)between the shortest ZP axis at the 2-cell stage and the Em-Ab boundary in the blastocyst. (B) The proportion of mouse embryos with tilt angles within a certain range, in relation to the ratio of the longest to the shortest diameter of the ZP. Embryos in Group I have the longest ZP diameter, 20%longer than that of the shortest. Embryos in Group II have the spherical ZP that results from its partial digestion and enlargement. Numbers within the bars indicate the number of embryos with tilt angles within the range specified on the right. (C) Sequential DIC images of embryos with enlarged and spherical ZP developing from the 2-cell to the blastocyst stage. Solid and dashed colored lines indicate, for each embryo, the shortest ZP axis at the 2-cell stage and the Em-Ab boundaries at the mid-blastocyst stage,respectively. The tilt angle for each embryo is measured by superimposing these two lines. White arrowheads indicate the ZP. In each frame, time is given in hours:minutes after human chorionic gonadotropin (hCG) injection. Scale bar: 50 μm.

Fig. 3.

Experimental verification of the predictions based on the computer simulation. (A) Scheme illustrating the tilt angle (dark green)between the shortest ZP axis at the 2-cell stage and the Em-Ab boundary in the blastocyst. (B) The proportion of mouse embryos with tilt angles within a certain range, in relation to the ratio of the longest to the shortest diameter of the ZP. Embryos in Group I have the longest ZP diameter, 20%longer than that of the shortest. Embryos in Group II have the spherical ZP that results from its partial digestion and enlargement. Numbers within the bars indicate the number of embryos with tilt angles within the range specified on the right. (C) Sequential DIC images of embryos with enlarged and spherical ZP developing from the 2-cell to the blastocyst stage. Solid and dashed colored lines indicate, for each embryo, the shortest ZP axis at the 2-cell stage and the Em-Ab boundaries at the mid-blastocyst stage,respectively. The tilt angle for each embryo is measured by superimposing these two lines. White arrowheads indicate the ZP. In each frame, time is given in hours:minutes after human chorionic gonadotropin (hCG) injection. Scale bar: 50 μm.

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, χ2 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, χ2 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.

Fig. 4.

Computer simulation of mammalian blastocyst morphology. Drawing (top left) shows cross-sections by xz- and yz-planes of the simulated aggregate.(A) An example of the calculation process to a stable state, viewed from cross-sections of the yz-plane. Numbers indicate the time point of the simulation (t). A vertex (a black circle in the sample at t=0) is replaced by a tetrahedron (see Fig. 1C) and is enlarged until its volume reaches half the initial total volume (at t=500; see Fig. 1D). The blastocyst axis keeps changing during t=500-2000 (see text) until it is localized and stabilized at one end of the long axis of the ellipsoidal ZP(t=2000), when it no longer migrates (t=2000-3500). Several cells are marked (green, red, orange and blue) to illustrate their movement. In addition, four cells surrounding the orange cell at t=1000 are marked by black and red squares and circles. Some of these cells are not visible at certain time points because they are moving in 3D. (B) The simulated blastocyst at t=2000 in cross-sectional views of xz- and yz-planes. (C) A stereoscopic view of the blastocyst at t=2000, in which the ZP and some of the TE cells are removed for internal view.

Fig. 4.

Computer simulation of mammalian blastocyst morphology. Drawing (top left) shows cross-sections by xz- and yz-planes of the simulated aggregate.(A) An example of the calculation process to a stable state, viewed from cross-sections of the yz-plane. Numbers indicate the time point of the simulation (t). A vertex (a black circle in the sample at t=0) is replaced by a tetrahedron (see Fig. 1C) and is enlarged until its volume reaches half the initial total volume (at t=500; see Fig. 1D). The blastocyst axis keeps changing during t=500-2000 (see text) until it is localized and stabilized at one end of the long axis of the ellipsoidal ZP(t=2000), when it no longer migrates (t=2000-3500). Several cells are marked (green, red, orange and blue) to illustrate their movement. In addition, four cells surrounding the orange cell at t=1000 are marked by black and red squares and circles. Some of these cells are not visible at certain time points because they are moving in 3D. (B) The simulated blastocyst at t=2000 in cross-sectional views of xz- and yz-planes. (C) A stereoscopic view of the blastocyst at t=2000, in which the ZP and some of the TE cells are removed for internal view.

Fig. 5.

Volume of inner cells and outer cells of the blastocyst after simulations. The data are based on a total of five simulations under various conditions (differing in direction of the long axis of the ellipsoidal ZP, and in the initial position of the cavity). Average volumes (±s.d.)of the inner cells (black bars) and of the outer cells (gray bars) are 0.702±0.0088 (n=56) and 0.744±0.0136 (n=144),respectively.

Fig. 5.

Volume of inner cells and outer cells of the blastocyst after simulations. The data are based on a total of five simulations under various conditions (differing in direction of the long axis of the ellipsoidal ZP, and in the initial position of the cavity). Average volumes (±s.d.)of the inner cells (black bars) and of the outer cells (gray bars) are 0.702±0.0088 (n=56) and 0.744±0.0136 (n=144),respectively.

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.).

Aiken, C. E., Swoboda, P. P., Skepper, J. N. and Johnson, M. H. (
2004
). The direct measurement of embryogenic volume and nucleo-cytoplasmic ratio during mouse pre-implantation development.
Reproduction
128
,
527
-535.
Alarcon, V. B. and Marikawa, Y. (
2003
). Deviation of the blastocyst axis from the first cleavage plane does not affect the quality of mouse postimplantation development.
Biol. Reprod.
69
,
1208
-1212.
Alarcon, V. B. and Marikawa, Y. (
2005
). Unbiased contribution of the first two blastomeres to mouse blastocyst development.
Mol. Reprod. Dev.
72
,
354
-361.
Aziz, M. and Alexandre, H. (
1991
). The origin of the nascent blastocoele in preimplantation mouse embryos: ultrastructural cytochemistry and effect of chloroquine.
Roux's Arch. Dev. Biol.
200
,
77
-85.
Barlow, P., Owen, D. A. and Graham, C. (
1972
). DNA synthesis in the preimplantation mouse embryo.
J. Embryol. Exp. Morphol.
27
,
431
-445.
Beddington, R. S. and Robertson, E. J. (
1999
). Axis development and early asymmetry in mammals.
Cell
96
,
195
-209.
Calarco, P. G. and Brown, E. H. (
1969
). An ultrastructural and cytological study of preimplantation development of the mouse.
J. Exp. Zool.
171
,
253
-283.
Chisholm, J. C., Johnson, M. H., Warren, P. D., Fleming, T. P. and Pickering, S. J. (
1985
). Developmental variability within and between mouse expanding blastocysts and their ICMs.
J. Embryol. Exp. Morphol.
86
,
311
-336.
Chroscicka, A., Komorowski, S. and Maleszewski, M.(
2004
). Both blastomeres of the mouse 2-cell embryo contribute to the embryonic portion of the blastocyst.
Mol. Reprod. Dev.
68
,
308
-312.
Dietrich, J. E. and Hiiragi, T. (
2007
). Stochastic patterning in the mouse pre-implantation embryo.
Development
134
,
4219
-4231.
Gardner, R. L. (
1997
). The early blastocyst is bilaterally symmetrical and its axis of symmetry is aligned with the animal-vegetal axis of the zygote in the mouse.
Development
124
,
289
-301.
Gardner, R. L. (
2001
). Specification of embryonic axes begins before cleavage in normal mouse development.
Development
128
,
839
-847.
Gueth-Hallonet, C., Antony, C., Aghion, J., Santa-Maria, A.,Lajoie-Mazenc, I., Wright, M. and Maro, B. (
1993
). gamma-Tubulin is present in acentriolar MTOCs during early mouse development.
J. Cell Sci.
105
,
157
-166.
Hiiragi, T. and Solter, D. (
2004
). First cleavage plane of the mouse egg is not predetermined but defined by the topology of the two apposing pronuclei.
Nature
430
,
360
-364.
Hiiragi, T., Alarcon, V. B., Fujimori, T., Louvet-Vallee, S.,Maleszewski, M., Marikawa, Y., Maro, B. and Solter, D.(
2006
). Where do we stand now? Mouse early embryo patterning meeting in Freiburg Germany (2005).
Int. J. Dev. Biol.
50
,
581
-588.
Honda, H., Tanemura, M. and Nagai, T. (
2004
). A three-dimensional vertex dynamics cell model of space-filling polyhedra simulating cell behavior in a cell aggregate.
J. Theor. Biol.
226
,
439
-453.
Johnson, M. H. and McConnell, J. M. (
2004
). Lineage allocation and cell polarity during mouse embryogenesis.
Semin. Cell Dev. Biol.
15
,
583
-597.
Kelly, S. J., Mulnard, J. G. and Graham, C. F.(
1978
). Cell division and cell allocation in early mouse development.
J. Embryol. Exp. Morphol.
48
,
37
-51.
Klaffky, E., Williams, R., Yao, C. C., Ziober, B., Kramer, R. and Sutherland, A. (
2001
). Trophoblast-specific expression and function of the integrin alpha 7 subunit in the peri-implantation mouse embryo.
Dev. Biol.
239
,
161
-175.
Kurotaki, Y., Hatta, K., Nakao, K., Nabeshima, Y. and Fujimori,T. (
2007
). Blastocyst axis is specified independently of early cell lineage but aligns with the ZP shape.
Science
316
,
719
-723.
Louvet-Vallee, S., Vinot, S. and Maro, B.(
2005
). Mitotic spindles and cleavage planes are oriented randomly in the two-cell mouse embryo.
Curr. Biol.
15
,
464
-469.
Motosugi, N., Bauer, T., Polanski, Z., Solter, D. and Hiiragi,T. (
2005
). Polarity of the mouse embryo is established at blastocyst and is not prepatterned.
Genes Dev.
19
,
1081
-1092.
Nadijcka, M. and Hillman, N. (
1974
). Ultrastructural studies of the mouse blastocyst substages.
J. Embryol. Exp. Morphol.
32
,
675
-695.
Niwa, H., Toyooka, Y., Shimosato, D., Strumpf, D., Takahashi,K., Yagi, R. and Rossant, J. (
2005
). Interaction between Oct3/4 and Cdx2 determines trophectoderm differentiation.
Cell
123
,
917
-929.
Ohno, Y. and Isoda, K. (
1977
).
Suuchi Keisan Handbook
. Tokyo: Ohmsha.
Piotrowska, K., Wianny, F., Pedersen, R. A. and Zernicka-Goetz,M. (
2001
). Blastomeres arising from the first cleavage division have distinguishable fates in normal mouse development.
Development
128
,
3739
-3748.
Piotrowska, K. and Zernicka-Goetz, M. (
2001
). Role for sperm in spatial patterning of the early mouse embryo.
Nature
409
,
517
-521.
Piotrowska-Nitsche, K., Perea-Gomez, A., Haraguchi, S. and Zernicka-Goetz, M. (
2005
). Four-cell stage mouse blastomeres have different developmental properties.
Development
132
,
479
-490.
Ralston, A. and Rossant, J. (
2008
). Cdx2 acts downstream of cell polarization to cell-autonomously promote trophectoderm fate in the early mouse embryo.
Dev. Biol.
313
,
614
-629.
Rossant, J. and Tam, P. P. (
2004
). Emerging asymmetry and embryonic patterning in early mouse development.
Dev. Cell
7
,
155
-164.
Smith, A. (
2005
). The battlefield of pluripotency.
Cell
123
,
757
-760.
Smith, R. and McLaren, A. (
1977
). Factors affecting the time of formation of the mouse blastocoele.
J. Embryol. Exp. Morphol.
41
,
79
-92.
Strumpf, D., Mao, C. A., Yamanaka, Y., Ralston, A.,Chawengsaksophak, K., Beck, F. and Rossant, J. (
2005
). Cdx2 is required for correct cell fate specification and differentiation of trophectoderm in the mouse blastocyst.
Development
132
,
2093
-2102.
Tarkowski, A. K. and Wroblewska, J. (
1967
). Development of blastomeres of mouse eggs isolated at the 4- and 8-cell stage.
J. Embryol. Exp. Morphol.
18
,
155
-180.
Thorsteinsdottir, S. (
1992
). Basement membrane and fibronectin matrix are distinct entities in the developing mouse blastocyst.
Anat. Rec.
232
,
141
-149.
Torres-Padilla, M. E., Parfitt, D. E., Kouzarides, T. and Zernicka-Goetz, M. (
2007
). Histone arginine methylation regulates pluripotency in the early mouse embryo.
Nature
445
,
214
-218.
Wiley, L. M. and Eglitis, M. A. (
1981
). Cell surface and cytoskeletal elements: cavitation in the mouse preimplantation embryo.
Dev. Biol.
86
,
493
-501.