MorphoSim: an efficient and scalable phase-field framework for accurately simulating multicellular morphologies

Abstract

The phase field model can accurately simulate the evolution of microstructures with complex morphologies, and it has been widely used for cell modeling in the last two decades. However, compared to other cellular models such as the coarse-grained model and the vertex model, its high computational cost caused by three-dimensional spatial discretization hampered its application and scalability, especially for multicellular organisms. Recently, we built a phase field model coupled with in vivo imaging data to accurately reconstruct the embryonic morphogenesis of Caenorhabditis elegans from 1- to 8-cell stages. In this work, we propose an improved phase field model by using the stabilized numerical scheme and modified volume constriction. Then we present a scalable phase-field framework, MorphoSim, which is 100 times more efficient than the previous one and can simulate over 100 mechanically interacting cells. Finally, we demonstrate how MorphoSim can be successfully applied to reproduce the assembly, self-repairing, and dissociation of a synthetic artificial multicellular system - the synNotch system.

Fig. 1. The phase-field framework for C. elegans early embryogenesis.

a Flow chart. b Comparison between ab initio simulation and in vivo ground truth of C. elegans embryonic morphologies from 2- to 8-cell stages. The embryonic structures at 2-, 3-, and 4-cell stages are two-dimensional and shown by fluorescence images; the ones at 6-, 7-, and 8-cell stages are three-dimensional and shown by segmented images. The cells’ identities and corresponding colors and existing stages (represented by the total cell number of the stage) are denoted on the right. All the subfigures, except the schematics for the initial state setting in (a), are adapted from ref. with granted permission.

Fig. 2. Framework improvement by the addition of stabilization term.

a Computing time compared to the baseline at different stages in the first-order scheme. b Computing time compared to the baseline at different stages in the second-order scheme. In both (a) and (b), the baseline is set as the computing time when $\delta t=0.3$ and no stabilization term is added. In b, the dashed line is the curve $t_{\mathrm{c}}^{\prime }=0.3\times \delta t^{-1}$ and the inset is the same figure presented in a dual logarithmic coordinate system. c Embryonic morphologies from 1- to 8-cell stages in the second-order scheme with $\delta t=2.0$. The cells’ identities and corresponding colors are denoted on the right; ${\sigma }_{\mathrm{S}}$ and ${\sigma }_{\mathrm{W}}$ denote the relatively strong (i.e., ${\sigma }_{\mathrm{S}}=0.5$) and weak (i.e., ${\sigma }_{\mathrm{W}}=0$) cell–cell attraction inferred at the 4-cell stage respectively (Supplementary Fig. 2a, b and Supplementary Note 1).

Fig. 3. Cell disappearance that happens when the cell size is too small.

a Gradual disappearance of the P4 cell during the 24-cell stage; the P4 cell is colored dark red and pointed by an arrowhead. b Phase field distribution in 1D (x axis) and 2D (xy plane) of a free cell with a radius of 4 μm (upper two rows) and 5 μm (lower two rows). The single-cell simulation is conducted with open boundaries and illustrated with an in silico interval of 20 from left to right. In the 1D distribution, the boundary of a cell ($0.07\le \varphi <0.93$) is painted with light gray while the interior ($\varphi \ge 0.93$) and exterior ($\varphi <0.07$) of a cell are painted with dark gray and white respectively. c A scanning on parameter c with different cell radii R_min; yellow and blue indicate if the cell disappears or not in a single-cell simulation.

Fig. 4. The graphical user interface of MorphoSim.

a The interface and instruction of inputs required. b, c The simulation inputs for 8-cell C. elegans embryogenesis with relatively strong (i.e., ${\sigma }_{\mathrm{ABpl},\mathrm{E}}={\sigma }_{\mathrm{S}}=0.5$) and weak (i.e., ${\sigma }_{\mathrm{ABpl},\mathrm{E}}={\sigma }_{\mathrm{W}}=0.5$) adhesion in ABpl-E contact respectively. (d) The initial state (in silico time = 0) of the 8-cell embryo. e, f The final state (in silico time = 15,000) of the 8-cell embryo with relatively strong (i.e., ${\sigma }_{\mathrm{ABpl},\mathrm{E}}={\sigma }_{\mathrm{S}}=0.5$) and weak (i.e., ${\sigma }_{\mathrm{ABpl},\mathrm{E}}={\sigma }_{\mathrm{W}}=0.5$) adhesion in ABpl-E contact respectively.

Fig. 5. The 1- to 102-cell C. elegans embryogenesis.

a The cell lineage tree reconstituted from in vivo experiment. The cell divisions from the same founder cell (i.e., AB, MS, E, C, D, or P4) and in the same generation are regarded as synchronous; the cell cycle length of each division group is approximated with the shortest one among the cells obtained from experimental measurement. The in vivo time axis is placed on the left and the cell division order is labeled near each cell division group. b The embryonic morphologies generated by phase-field simulation, using the cell division order and axis and volume segregation ratio measured experimentally as input. At each stage defined by the cell division order in (a), the in silico time (i.e., simulation time) from 1- to 8-cell stages is introduced in Supplementary Note 1, while the one after the 8-cell stage is set to reach the equilibrium ($\overline{v}<1\times 10^{-4}$).

Fig. 6. In silico simulation (left panel) and in vivo experiment (right panel) of the synNotch systems.

In each row, the in silico structures are shown from time points 0 to 50,000 in a step of 10,000; the computational domain is outlined with a black cube. In the first row, since in silico time = 40,000 (the fifth column), a part of red cells, which contact at least one green cell, are painted blue considering that the red cells (Type 1) receive contact-dependent signaling from the green cells (Type 2) and undergo subsequent differentiation into the blue cells (Type 3), in accordance to the experimental description. In both the fourth and fifth rows, the first in silico structure is adopted from the last one in the third row; for the former one, the cells with z < 0 or without contact with the cell aggregate are removed and shown by gray shadow. For the in silico structures, the cell types and corresponding colors and value assignments on cell–cell adhesion are listed in Table 2. For the in vivo structures, the scale bars represent 100 μm in reality; the fluorescence colors are in line with the ones used for the in silico structures, except that the red cells in the fourth and fifth rows also correspond to the blue ones in silico as they are the same in the cell–cell adhesion program. The images of in vivo experiment are from ref. , reprinted with permission from AAAS.