Numerical study on spindle positioning using phase field method

A method of numerical simulation of cell division using phase fields is presented. The cell division plane is obtained as a result of the spindle position and orientation considered with the spatial distribution of the activated cortical force generators and the dividing cell shape. To exemplify the application of the proposed method, numerical simulations of the development of cysts and early embryos are performed. The numerical results demonstrate that the activated cortical force generators that are localized at the lateral cortices of the epithelial cells lead to the formation of a single central lumen. It is additionally shown that the linear distribution of the activated cortical force generators along the animal-vegetal axis of a spherical cell engenders a similar cell proliferation of the divided embryo generated by the 32 cell period in a sea cucumber.


Introduction
Precise positioning of the spindle is necessary for organizing multicellular tissues. Recently, it has been reported that the spindle position is determined by the spatial distribution of cortical force generators that are activated by cortical cues [1][2][3][4][5][6][7]. Activated cortical force generators (aCFGs) apply force to the spindle poles via astral microtubules (MTs), which spread radially from the poles. Controlling the spatial distribution of aCFGs leads to symmetric and asymmetric cell division in circular cells [7]. Moreover, the spindle positioning is known to be influenced by the cell shape, as well as by the spatial distribution of the aCFGs. For instance, the spindle axis has been observed to lie roughly parallel to the long axis of the cell [8].
For numerical study of the morphogenesis, a mathematical model that can systematically express the major developmental events, such as the cell size growth, cell death, cell differentiation, cell movement and cell division, is required. One of the authors has reported a mathematical model for a multicellular system based on the phase-field method [9], which can express these developmental events. However, even though the cortex of each cell can be easily described by the diffuse interface of the phase field of this model, the numerical simulation with the cell divisions in [9] is obtained without considering the spindle positioning but with using a plane of which orientation is randomly chosen, passing through the center of the dividing cell as the cell division plane. On the other hands, the remaining three authors proposed a mathematical model of cleavage based on the spindle pole movements [10].
In this study, we mathematically describe the mechanism of the spindle positioning [7,10] and incorporate it with the multicellular model [9]. Thus, using the proposed model, cell division considering the spindle positioning can be handled for any shape of cell. The basic concept underlying the spindle positioning in this paper was revealed in a previous study [7], in which the symmetric division in circular cells was computed by considering the spatial distribution of the aCFGs and the obtained spindle positions were compared with those that were experimentally determined.

M Akiyama et al
To show the applicability and validity of the presently proposed model, we present numerical results of the epithelial culture system, such as Madin-Darby canine kidney (MDCK) and Caco-2. These cell cultures have been extensively investigated because of their ability to reproduce tissue-like organization, such as sheets, cysts and tubules, when grown inside an extracellular matrix (ECM) [11][12][13]. During cyst development, each epithelial cell grows and divides, and fluid is secreted. The secretion of fluid induces the formation of a lumen enclosed by a monolayer of the epithelial cells. Each cell in the developed cyst or tubule has three types of membrane surfaces: an apical surface facing the lumen, a lateral surface adhering to adjacent cells and a basal surface adhering to the ECM. A tight junction also defines the boundary between the apical surface and the lateral surface. The establishment of such polarity within a single cell is a key feature of epithelial morphogenesis. In addition, Jaffe et al concluded that the apical surface facing a lumen is established between the two daughter cells as early as the first cell division [14] by the localization of the apical marker aPKC associating with the midbody prior to the occurrence of abscission, which is the last step of cytokinesis.
Several studies have focused on the components required for proper spindle orientation in cyst development [14][15][16]. One such component is the leucine-glycine-asparagine tripeptide in the N-terminal region (LGN).
LGN is known to be a component of the machinery for mitotic spindle positioning. It is interesting to note that the localization of LGN at the lateral cell cortex of a dividing cell has been observed during cyst development [16]. Moreover, the depletion of LGN in MDCK cells disrupts the spindle orientation and leads to inappropriate positioning of apical surfaces, thereby forming multi-lumen structures [15,16].
Although these cell cultures are much simpler than tissues, theoretically clarifying all development processes remains too complex. To date, cyst development has been numerically investigated using an in silico epithelial analogue [17] considering the hexagonal lattice sites. Each lattice size corresponds to the size of a cell. The state of each site changes according to the axiomatic operating principles that define the simulated cell actions, such as cell death and division, in response to the state of the surrounding sites. Dysregulation of the axiom was numerically shown to have a disruptive effect on cyst formation.
We assume that some intracellular materials, such as LGN, activate the CFGs. We thus studied how the cyst structure is affected by the aCFGs localized at the lateral cell cortex. The fluid secretion, as well as the changing shapes of cells and cell divisions, are considered in the model. From the numerical simulations performed both in two dimensions (2D) and three dimensions (3D), we found that the localization of the aCFGs on the lateral cortex is needed for the cyst formation with the single lumen.
Another example of the proposed model application is cell divisions in an early embryo. It has been reported that some materials localize in the cortex in the early embryo [18][19][20] and that disruptions of the localization can have a significant impact on the proper development of the embryo [21]. The detailed mech anisms of spatial segregation in the cortex are energetically investigated. For instance, Wang et al reported that cytoskeleton contractility plays important roles in the cortical clustering and the displacement of proteins [22]. However, for simplicity, we do not consider such detailed mechanisms; rather, we apply the fact that the localization occurs at the cortex. In this paper, we make an assumption that some localized materials nonuniformly activate the CFGs at the cortex. To express a blastomere, we consider a spherical region in which all cells are enclosed. The successive divisions occur in a certain period of time.

Multicellular model based on the phase-field method
For numerical study of morphogenesis, we use the phase-field method for multicellular systems [9]. For self-containedness, a brief explanation of the method is provided in this subsection. Note that all model equations presented in this paper are written with dimensionless parameters.
In this model [9], a vector variable u(r, t) = (u 1 (r, t), · · · , u M (r, t)) is introduced, where M is the total number of cells in the system, r denotes the position, and t represents the time. The component u m (r, t) (m = 1, · · · , M ) describes the shape of the mth cell, such that u m (r, t) u c in the region with the mth cell (<u c in the region not taken up by the mth cell), with a constant u c ∈ (0, 1).
The time evolution of the mth cell shape is described by the following equation: where the coefficients τ u , ε u , α, β, γ, η, and V m are positive constants. The variable ψ(r, t) is defined as ψ(r, t) = m h(u m (r, t)), where h(u) = u 2 (3 − 2u). The volume of the mth cell, defined as v(u m ) = Ω h(u m )dr, is assumed to be conserved, where Ω denotes the area of the system. It should be noted that (1) is not explicitly written with the variables u m (m = m and m = 1, · · · , M), which denote the other cell shapes; rather, it is written with the variable ψ. Thus, one can easily design a program that does not consume computer memory and is appropriate for parallel computing.
If f u is a constant f , (1) is known as the Allen-Cahn equation in the field of material science and can be easily solved in one dimension as This means that the front moves and expands the region of u m = 1 (u m = 0) if f > 0 (f < 0). Therefore, the volume of the mth cell v(u m ) approaches the target volume V m by the first term on the right side of (2). The second term on the right side of (2) describes the condition in which the cell does not overlap other cells. This is because the term −β(ψ − h(u m )) is negative in regions where other cells exist. The remaining terms of (2) express the effect of adhesion. The third term on the right side of (2) acts similar to negative diffusion in the region where the cell adhesion occurs. The fourth term on the right side of (2) is introduced to prevent divergence due to the third term of (2) with the condition whereby γ > η.
The values of parameters can be selected depending on the characteristics of the shapes of the target cells. For example, if the target cells can easily change their volumes, the parameter α should be set as a small value. Moreover, the value of ε u concurrently controls the thickness of cell membrane denoted by the diffuse interface of u m and the cell deformability from the circle in 2D or the sphere in 3D. Thus, when the cells have the thick membrane and easily deform their shape from the circle or the sphere, one needs to set V m large as well as ε u .

Model for spindle positioning
Now let us consider the situation in which the cell division of the mth cell begins at t. To calculate the spindle position using the cell shapes obtained by (1), we use the fact that the cortex of the mth cell can be expressed as a function of u m , such as H(u m ) = u m (1 − u m ), since u m has an interface with a thickness of the order of ε u . In other words, H(u m ) has a positive value of order 1 in the region where the cortex of the mth cell exists, whereas H(u m ) becomes almost zero in the other region.
The spindle poles, pulled by the MTs radiating from them, are assumed to be connected to each other by a spring, as shown in figure 1(a). Therefore, the MTs bound to the aCFGs at the mth cell cortex in the region denoted by Ω 1 (Ω 2 ) in figure 1(b) can exert a force on the spindle pole p 1 ( p 2 ).
We assume that the time evolution of the spindle positioning is much faster than that of the shape change. In other words, the shapes of the cells and the spatial distributions of the aCFGs are assumed to not change during spindle positioning. This assumption means that only the final positions of the poles are used in the results of spindle positioning presented in this paper. To avoid confusion, we use another expres-sion, t s , to denote the time in the calculations of spindle positioning. That is, during the calculation of the spindle positioning, the time t s changes instead of t . Once the positions of the poles reach a steady state, the mth cell is divided into two daughter cells. Then, the calcul ations of the cell shapes, (1), are resumed, in which time t changes instead of t s .
In many organisms, oscillations have been observed during the process of spindle repositioning [23][24][25]. It has been theoretically demonstrated that the cooperative attachment and detachment of aCFGs to the MTs lead to spontaneous oscillations [26]. It is also known that the spindle positioning switches direction from centering to posteriorly displacing in the Caenorhabditis elegans embryo [27]. In addition, dynamical models of reaction network have been reported for studying the mechanism of proper spindle positioning [28,29]. However, since this paper focuses on how the plane of cell division is controlled by the spatial distribution of aCFG and the shape of the dividing cell, we do not consider such detailed dynamics.
The positions of the poles p 1 and p 2 at t s are denoted by r 1 (t s ) and r 2 We assume that the magnitude of pulling force from each MT depends linearly on its length [30] as well as the aCFGs density [7]. Thus, the force applied to the poles p i (i = 1, 2) can be given as follows [7]: where Ω i denotes the region shown in figure 1(b).
The density of the aCFGs at the cell cortex of the m th cell is written as follows: where the function H(u m ) = u m (1 − u m ) denotes the cell cortex of the mth cell, as mentioned in the first paragraph of this subsection. The spatial distribution of the aCFGs is decided by the function G. For instance, by setting G to be uniform throughout the dividing cell, the situation in which the cortex aCFGs are uniformly distributed at the cortex of the dividing cell can be addressed. The time evolution of each pole is determined by the summation of the forces exerted by the MTs and by the spring connecting the poles. It can thus be presented as follows: where µ, s and σ are positive constants. Equations (5) and (6) are solved until the positions of the poles reach a steady state, r * 1 and r * 2 .
The cell division of the mth cell can be expressed by setting two new phase fields, u m1 (r, t) and u m2 (r, t), from u m (r, t). The shapes of the two daughter cells can be expressed as: where the step function χ(r; r * 1 , r * 2 ) becomes one (zero) in the region Ω 1 (Ω 2 ) with the steady spindle position r * 1 and r * 2 . In the simulation, the step function χ(r; r 1 , r 2 ) in (7) and (8) is approximated by the following function: where ε c is a positive constant. The function g(r; r 1 , r 2 ) is defined as: where the line (plane) in 2D (3D) at which g(r; r 1 , r 2 ) = 0 represents the line (plane) bisecting the one connecting r 1 and r 2 at the center, and g(r; r 1 , r 2 ) becomes positive (negative) in the region Ω 1 (Ω 2 ).

Spindle positioning
To elucidate the effects of the spatial distributions of the aCFGs and the cell shapes on the spindle positioning, the function G and the cell shape u 1 are set to be different in the panels of figures 2(a) and (b). The other parameters are set to be the same. The shapes of the cells expressed by the phase field u 1 were prepared in advance by solving (1). The spindle is initially placed at the center of each cell, denoted by the × symbol in figure 2. Equations (5) and (6) are solved until the spindle poles reach their stationary positions. We have checked that the initial positions of the spindle poles do not affect their stationary positions as shown in figures S1 and S2, see supplemental information (stacks.iop.org/PhysBio/16/016005/mmedia). Figure 2(a) shows the final position of the spindle in a circular cell. The function G is set as a function of the angle from the y-axis, θ, as a + (θ/π) 2 , a + |θ/π| and a + |θ/π|, from left to right, respectively, where a is a constant. The definition of θ (−π θ π) is found in the left panel of figure 2(a). From figure 2(a), one can easily observe that the spindle is situated on the lower part of cell, where G is large, and is perpend icular to the y-axis regardless of the form of the function G.
The panels of figure 2(b) show results for the spindle pole position in cells with various shapes. The cell shapes are made in advance by solving (1). In order to make the various cell shapes, a new term proportional to h(q(r)) is added to (2), where q(r) is a new variable which is independent on time. This new term represents that the cells escape from (to) the region where q = 1 (q = 0). The regions where q = 1 are colored by gray in figure 2(b). The function G is set to be a constant everywhere in the simulation box of the results shown in figure 2(b). That is, the aCFGs are set to be uniformly distributed at the cell cortex. In such cases, the spindle is placed parallel to the longitudinal direction of each cell, since the forces applied to the spindle poles depends linearly on the distance from the cell cortex. Despite the cell shape, the spindle is located around the cell center because of the balance between the forces applied to the poles.
In addition, we numerically verified how the spindle positions are affected by the magnitudes of the parameters that determine the spindle character, i.e. σ and s . It is found that the orientations of spindles and the midpoints between the positions of the two poles are not affected by changing the values of s and σ, whereas the spindle becomes longer for a large value of s and a small value of σ. That is, the plane of cell division is not affected by changing the values of s and σ because this plane is determined by the line perpendicularly bisecting the line connecting the poles.
Incidentally, substituting (4) to (3), one can find that components of the force F i = (F (i) x , F (i) y ) (i = 1, 2) are decoupled. Thus, the force can be easily estimated by analytical expression in some case. For example, the situation of figure 2(a) is one of the typical cases, in which the spindle is placed perpendicular to y-axis within a circular cell. Figure 3 shows the graphs of F The × symbol indicates the center of the cell. The parameters for computing u 1 are set as follows: a simulation box of Ω = 6 × 6 in size, a spatial grid size of δ = 5.00 × 10 −2 , and a time increment of dt = 1.00 × 10 −2 , τ u = 1.00, ε 2 u = 1.00 × 10 −3 , α = 1.00 × 10, β = γ = η = 0.00 and V 1 = 2.00. For calculation of the pole positions, the parameters are set as follows: the time increment is dt s = 1.00 × 10 −1 , a = 2.00 × 10 −1 , µ = 1.00, σ = 1.00 × 10 −1 and s = 0.00.  (3) and (4) with the phase field and by (b) analytical estimations. The dotted lines, the thick solid lines, the thin solid lines indicate the y-components of the forces, F y , in the cases of G(θ) = a + (θ/π) 2 , a + |θ/π| and a + |θ/π|, respectively. ure 3(a) are performed with the same settings with those for figure 2(a). In the analytical estimation, the diffuse interface of the phase field is approximated by the sharp interface. The details of the analytical estimation are described in the supplemental information. In spite of the rough estimation, the graphs in figure 3(b) capture the features of those in figure 3(a).

Cyst formation
As an application of the proposed model, the development of a cyst enclosing a fluid-filled lumen is investigated in this subsection. During development, each epithelial cell in ECM grows and divides with the secretion of fluid. We use the facts that the apical surface facing a lumen is established between the two daughter cells as early as the first cell division [14], and that the depletion of LGN localized at the lateral cortex disrupts the spindle orientation and leads to inappropriate positioning of apical surfaces, thereby forming multi-lumen structures [16].
To incorporate the effect of the localized aCFGs on the lateral surface of a dividing cell, we assume the form of the function G as follows: where ρ 0 and ρ e are positive constants and e η (r, t) is given as: (r, t)). (12) The value e η is negative in the region where the cells adhere, and it has a zero value in the other region. In fact, the functional derivative of the integral over r of (12) is identical to the third term of the right side of (2). Substituting (11) into (4), and solving (5) and (6), the plane of division of the cell can be obtained considering the aCFGs localized on the lateral surfaces.
A new variable s(r, t) is introduced to represent the lumen location, such that the region where a lumen exists is indicated by s s l (< s l in the region without a lumen), with a constant s l ∈ (0, 1). The time evolution equation for s is given as follows: where τ s , ε s , ξ and β s are positive constants. With a large value of the constant ξ, liquid is prone to secretion. The second term of (14) express that the lumen does not overlap the cells. Similarly, to express that the cells do not overlap the lumen, a new term, −β s h(s) , is also added to the right side of (2) as: (15) Figure 4 shows the numerical results in 2D. We start the simulations with the data of a circle cell (v(u 1 ) = v ini ), placed at the center of the simulation box. The cell division of the mth cell is assumed to start when the mth cell volume v(u m ) becomes larger than the positive constant V d (< V m ). The three upper panels in figures 4(a) and (b) show the shapes of the cells upon cell division, as well as the stationary positions of the spindle poles within the dividing cells. At the next time step of each of these panels, two daughter cells are generated by (7) and (8). At the same time, to express the apical cortex formation at the midbody of a dividing cell [14], we set a region, s = 1.00, for a circle (sphere) of radius r s in 2D (3D). The center of this circle (sphere) is the midpoint between the stationary positions of the spindle poles obtained by solving (5) and (6). The lowest panels show the structures with 15 cells, which are obtained by performing the simulations until t = 260.00.
We set ρ e as 5.00 and 1.00 × 10 −3 for the simulations shown in figures 4(a) and (b), respectively, where the other parameters are the same for both cases. In figure 4(a), a single lumen enclosed by a monolayer is formed, whereas figure 4(b) shows the time evolution of a cyst containing multiple lumens. In these simulations, the tight junctions are not considered. Therefore, if further simulations are performed without further cell division, the region where the value of s is not zero expands outside the cyst by breaking the cell adhesion. As discussed in the previous subsection, the plane of cell division is not affected by changing the values of parameters related to spindle positioning, s and σ.
From the simulations shown in figure 4, we found that the localized aCFGs cause two actions needed for cyst formation with a single lumen. One is cell division with a spindle perpendicular to the apico-basal axis for generating the monolayer of epithelial cells. If the value of ρ e is not adequately large, the spindles tend to be oriented along the long axis of the cell, as shown in the panels at t = 85.12 and t = 88.10 in figure 4(b). In other words, if the value of ρ e is sufficiently large, the spindle is located perpendicular to the apicobasal axis during development of the cyst structures, and thus a monolayer of cells surrounding a single lumen is formed. The other is asymmetric positioning of the spindle for generating a centrally localized lumen. Comparing the top panels of figures 4(a) and (b) reveals that the centers of the spindles for the case of a large value of ρ e is located nearer to the cyst center than for the case of a small value of ρ e . Thus, the small regions of s = 1.00 created at the centers of the spindles during cell divisions are incorporated into the central lumen if ρ e is large. Conversely, the small regions are easily isolated if the value of ρ e is small. These results are consistent with the experimental findings whereby depleting LGN caused misorientation of spindles and the formation of multiple lumens [16]. Figure 5 shows a snapshot of a 3D simulation at t = 250.00. The gray surfaces in the upper panel in figure 5 are contour plots of u m = 2.00 × 10 −1 (m = 1, · · · , 40). From left to right, the lower panels of figure 5 show cross-sectional views at x = 3, x = 5, x = 7 and x = 9. The fluid-filled region described by the order parameter s is not shown, although its time evolution was also calculated. This region is enclosed by a monolayer of cells. As in the 2D simulations, we found that a cyst enclosing a single lumen is formed in 3D if the localization of the aCFGs at the lateral surface of the dividing cell is adequately strong.

Early cleavage of a sea cucumber
In this subsection, we consider the cell divisions within a spherical region of constant volume V s = 4πR 3 s /3. We regard this spherical region as a blastomere. The variable s(r, t), which was introduced in the previous subsection for describing the lumen, is initially set to one (zero) if the distance from the center of the simulation box is larger (smaller) than R s with an interface with a thickness of the order of ε s and does not change throughout the simulation. Therefore, the cells whose shapes are automatically computed by solving (1) with (15) are unable to penetrate into the outside region of the sphere of s.
The cell division occurs at the same time in all cells at an interval of T. We set the objective size of the daughter cell, m, which is denoted as V m in (15), to be 0.9 times its volume immediately after the cell division. Figure 6 shows the results where G is a linear function of z , which connects G = A at the top of the sphere of s, z = L/2 + R s , and G = B at the bottom of the sphere of s, z = L/2 − R s , where L denotes a length of a side of the simulation box. We consider the case B > A > 0. We take the z -axis along one side of the simulation box as shown in figure 6. In other words, the function G is written as G = (A − B)(z − L/2 + R s )/(2R s ) + B. We start the simulations with the data u 1 , which is one (zero) in (out) of a sphere, of which the volume v ini is chosen to be smaller than V s . Before the first cell division, the volume of the cell changes to approach its objective size, V 1 . The cell division occurs at an interval of T = 30.00. The snapshots in figure 6 from the top left to the bottom right are the results at t = 27.00, 57.00, 87.00, 117.00, 147.00 and 177.00, respectively, whereas the cell divisions occur five times at t = 30.00, 60.00, 90.00, 120.00 and 150.00.
the first cell division is parallel to the z -axis. The two daughter cells produced by the first cell division take hemispherical forms. In each daughter cell, the direction of the spindle is perpendicular to the z -axis and parallel to the longitudinal direction of the semicircle, which is the cross-section of each daughter cell with a plane parallel to the x-y plane. Therefore, the plane of the second cell division is orthogonal to that of the first cell division. The remaining cell divisions are also analogously explained. If the differences between A and B are not large enough, the spindles fail to align properly as shown in figure S3.

Discussion and conclusion
We proposed a mathematical model that can serve as a powerful tool for numerical studies of morphogenesis. The model combines the mechanism of the spindle positioning [7,10] with the multicellular model [9]. Using the proposed method, one can perform numerical simulations of multicellular systems by taking cell divisions into account. This paper focuses on the final positions of poles by using the spatial distributions of the aCFGs and the cell shapes, even though it has been reported that spindle oscillations and switch directions have been observed during the process of spindle positioning [23][24][25][26][27] and that spindle positioning is regulated by the reaction network [28,29].
Numerically, we can reproduce the symmetric and asymmetric cell division, as reported in [7]. If the aCFGs are uniformly distributed at the cortex, the plane computed by using the proposed model is placed at the cell center and is perpendicular to the longitudinal axis of the cell. These numerical findings are consistent with the experimental observations [8].
To demonstrate the applicability and validity of the proposed model, we presented the numerical results of two developmental systems: cyst development, and early embryo development. Both systems are very important for studying morphogenesis. It has been reported that the localization at the cell cortex plays important roles in both systems. We numerically demonstrated that a cyst with a single liquid-filled lumen can be formed with the aCFGs localized at the lateral surfaces. The localization of the aCFGs orients the spindle perpendicular to the apico-basal axis, and asymmetric positioning of the spindle generates a centrally localized apical surface, which is analogous to experimental observations [14,16]. The result of the cell division enclosed in the spherical region is impressive. When the aCFGs are linearly distributed, we found the cell divisions to be similar to those in the early embryo of the sea cucumber. Further results with different forms of the function of G will be reported elsewhere.
The proposed model can express developmental events, such as cell size growth, cell death, cell differentiation, cell movement, and cell division, by considering the spindle positioning and shape of the dividing cell. Thus, we believe that the various systems can be investigated by using the proposed model.