Methods and systems for accurate and efficient simulation of multicellular morphology

Abstract

The application discloses a method and system for precisely and efficiently simulating multi-cell morphologies, which firstly gives initial states of a multi-cell system and phase field environment parameters; then semi-implicit discrete iteration is performed on a phase field evolution equation satisfied by the initial states to quickly solve morphologies (i.e. phase fields) of cells and cell groups. The application uses a diffusible field to describe cells, so that smooth, continuous and real cell characteristic cell morphology boundaries can be obtained; the application uses dense grids to accurately describe the morphology boundaries, and performs relative error control on cell volumes, so that higher reliability is achieved when calculating cell morphologies and intercellular interactions.

Description

TECHNICAL FIELD

[0001] The present application belongs to the technical field of computer simulation of life system, and relates to a simulation of multicellular morphology, in particular to a precise and efficient morphological simulation technology based on a phase field model for simulating real multicellular systems such as tissues, organs and embryos. BACKGROUND

[0002] Life systems are based on cells as the basic unit of activity, covering naturally occurring and artificially constructed and modified life systems. The mechanical and biochemical interactions of cell groups enable them to obtain functions beyond the single-cell level, realize the mutual combination of different cells in fate, function, and spatial architecture, and form a rich variety of multicellular life forms, such as animals, plants, and even some microorganisms (e.g., mushrooms and flagellates) that we see in daily life. For the most common multicellular animals, the interaction of cells endows specific spatial structures of tissues, organs, embryos, and even individuals; for example, the stomach cavity has a hollow structure, the intestinal tract has a tubular structure, and the lungs have a bifurcated structure. In addition to the systems that naturally exist in nature, many artificially constructed multicellular functional systems have emerged in recent frontier fields; for example, Xenobot robots that can move by themselves are constructed by Xenopus cells [S. Kriegman, D. Blackiston, M. Levin, J. Bongard. A scalable pipeline for designing reconfigurable organisms. Proc. Natl. Acad. Sci. U.S.A. 117 (4) (2020) 1853-1859], and synNotch machines that can self-organize into different shapes by encoding the specific adhesion of mouse cell surfaces [S. Toda, L. R. Blauch, S. K. Y. Tang, L. Morsut, W. A. Lim. Programming self-organizing multicellular structures with synthetic cell-cell signaling. Science 361 (6398) (2018) 156-162]. In addition, in the fields of medicine and developmental biology, organoids and embryoids obtained by culturing a small number of cells extracted from natural life systems in vitro also have similar morphological characteristics to real systems, and have been applied to drug screening and developmental mechanism research [J. Kim, B. K. Koo, J. A. Knoblich. Human organoids: Model systems for human biology and medicine. Nat. Rev. Mol. Cell Biol. 21 (2020) 571-584; N. Moris, K. Anlas, S. C. V. D. Brink, A. Alemany, J. S. Ghimire, T. Balayo, A. V. Oudenaarden, A. M. Arias. An in vitro model of early anteroposterior organization during human development. Nature 528(7812) (2020) 410-415; M. Rosner, M. Reithofer, D. Fink, M. Human embryo models and drug discovery. Int. J. Mol. Sci. 22(2) (2021) 637]. Similarly, tissue and organ engineering employs stem cells to re-culture in vitro into tissues and organs with certain functions and spatial structures, which are then transplanted back into patients to achieve a low-rejection therapeutic effect, such as treating skin burns [A. O. Lukomskyj, N. Rao, L. Yan, J. S. Pye, H. Li, B. Wang, J. J. Li. Stem cell-based tissue engineering for the treatment of burn wounds: A systematic review of preclinical studies. Stem Cell Rev. Rep. (2022)]

[0003] Morphologies of multicellular systems are ubiquitous in natural science research and real-world applications. However, the development of computer-aided simulation methods has many limitations, which restricts the scientific and industrial communities’ ability to compute, predict, and design morphologies of multicellular systems. The main cause of the backwardness of computational methods is the lack of experimentally verified model methods and the computing power to simultaneously compute cell populations at different scales, such as tissues, organs, and embryos.

[0004] There are several published computational models about multicellular population morphogenesis, the most widely used are coarse-grained model and vertex model. Coarse-grained model describes a cell with one particle, ignoring the complex and variable three-dimensional morphology of a cell [Z. Lv, J. Rosenbaum, S. Mohr, X. Zhang, D. Kong, Helen Preiss, S. Kruss, K. Alim, T. Aspelmeier, J. Grosshans. The emergent yo-yo movement of nuclei driven by cytoskeletal remodeling in pseudo-synchronous mitotic cycles. Curr. Biol. 30(13) (2020) 2564-2573; B. F. Nielsen, S. B. Nissen, K. Sneppen, J. Mathiesen, A. Trusina. Model to link cell shape and polarity with organogenesis. iScience 23 (2020) 100830]; vertex model reconstructs the cell boundary with Voronoi partition, describes a cell with a convex polyhedron (three-dimensional closest packing) including about 16 points [S. Alt, P. Ganguly, G. Salbreux. Vertex models: From cell mechanics to tissue morphogenesis. Phil. Trans. R. Soc. B 372 (2016) 20150520; K. Sato, D. Umetsu. A novel cell vertex model formulation that distinguishes the strength of contraction forces and adhesion at cell boundaries. Front. Phys. 9 (2021) 704878]. The aforementioned two models greatly simplify the cell morphology, limit the cell dimension, so the computational cost is very low, but at the same time, the ability to accurately describe the cell morphology is weakened. On the other hand, phase field model describes the cell as a diffusible fluid with a three-dimensional grid, giving the cell boundary the ability to smoothly deform.Recently, the phase field model has been demonstrated by simulating the embryonic development of nematodes, proving its ability to accurately simulate the morphological and mechanical interactions of real cells and cell populations. However, its three-dimensional mesh calculation and partial differential equation description still result in huge computational costs, and it is currently limited to two-dimensional or small-cell calculations, which restricts its widespread application [J.Jiang,K.Garikipati,S.Rudraraju.A diffuse interface framework for modeling the evolution of multi-cell aggregates as a softpacking problem driven by the growth and division of cells.Bull.Math.Biol.81(8)(2019)3282-3300;X.Kuang,G.Guan,MKWong,LYChan,Z.Zhao,C.Tang,L.Zhang.Computable early Caenorhabditis elegans embryo with a phase field model.PLoS Comput.Biol.18(1)(2022)e1009755]. Therefore, developing a precise and efficient phase-field model and method for simulating multicellular morphology would enable the scientific and industrial communities to obtain precise and efficient morphological calculation capabilities for naturally existing and artificially constructed and modified multicellular systems, providing computer-based simulation predictions and practical guidance for multicellular system research. Summary of the Invention

[0005] The purpose of this invention is to address the lack of accurate and efficient tools for calculating multicellular morphology in the existing technology. Based on an optimized and validated phase-field model, this invention provides a method and system for accurately and efficiently simulating multicellular morphology. This system is used to achieve dynamic simulation of the morphology of real multicellular systems at the scale of tissues, organs, and embryos, and can accurately and efficiently calculate the mechanical interactions and morphological movements of more than 100 cells.

[0006] To achieve the above objectives, the present invention adopts the following technical solutions.

[0007] In this invention, a cell is described by a diffusible three-dimensional continuous scalar field, with values ​​of 1 inside the cell and 0 outside the cell; this is achieved by identifying 0-1 transition regions or isosurfaces (φ). C This allows you to obtain smooth shape boundaries for cells.

[0008] In the present application, the biophysical parameters for characterizing cells are set to include the boundary and viscosity coefficient of the environment where the cells are located, the intercellular attractive force of each cell, the surface tension / rigidity, the motion noise (endogenous noise + exogenous noise), the volume, etc. Among them, the environmental boundary can be set to be of any shape or an infinite space; regarding the intercellular attractive force of each cell, a specific value can be assigned to each specific pair of cells, and there can also be asymmetric attractive forces between cells; regarding the volume, it can be set to a constant value or a function that changes over time to simulate life processes such as apoptosis, growth, etc. that are accompanied by changes in volume.

[0009] The present application first provides a method for accurately and efficiently simulating the morphology of a multicellular system composed of a plurality of non-dividing cells, which comprises the following steps:

[0010] S1 setting the initial state of a plurality of non-dividing cells and the phase field environment parameters;

[0011] S2 iteratively solving the following phase field evolution equation according to the initial state of the cells and the phase field parameters to obtain the morphology (i.e. phase field) of the plurality of non-dividing cells and cell groups;

[0012]

[0013] In the formula, φ i represents the phase field of the i-th cell; F ten represents the surface tension of the cell, which makes the cell surface area as small as possible to tend to a sphere, thereby reducing the deformation amount caused by external forces, and thus also representing the rigidity of the cell; F atr represents the attractive force between the i-th cell and the remaining cells; F rep represents the repulsive force between the i-th cell and the remaining cells and the constraint boundary; F vol represents the force for controlling the volume of the cell; t represents time; τ represents the viscosity coefficient of the environment; ξ i (t) represents the motion noise assigned to the i-th cell, and κ represents the noise intensity quantization constant.

[0014] Therefore, by the above method, the morphology of each cell in the multicellular system under the set environmental conditions can be obtained, and the biophysical parameters related to the multicellular system are provided for generating ideal morphology and motion patterns of multicellular robots, thereby efficiently and accurately obtaining an ideal multicellular system.

[0015] In the above step S1, the initial state of the cells includes the initial boundary shape and position of the cells. The phase field environment parameters include the constraint boundary (φ e ), the surface tension quantization constant (β) of the cell, the attractive force quantization constant (σ i,j ) between cells, the viscosity coefficient quantization constant (τ) of the environment, and the volume function (Vi (t))、motion noise intensity quantification constant (K) and the like.

[0016] In the above step S2, the cell mechanics under the phase field framework includes cell surface tension (or rigidity) (formula (2)) and attraction force between cells (formula (3)):

[0017]

[0018]

[0019] In the formula, N represents the total number of cells in the multicellular system; Δ represents the Laplace operator, represents the gradient operator; β represents the cell surface tension quantification constant (or rigidity); σ i,j represents the attraction force quantification constant between cell i and cell j; W (φ i ) = φ i 2 (φ i -1) 2 is a double potential well function, which separates the phase field into two states of 0 (outside the cell) and 1 (inside the cell), and its derivative is W' (φ i ) = 2φ i (φ i -1) (2φ i -1); c represents a quantification constant for controlling the width of the phase field boundary.

[0020] The phase field of the cell group needs to be mutually repulsive to avoid overlapping; similarly, the phase field of the cell and the constraint boundary cannot overlap (formula (4)):

[0021]

[0022] In the formula, g represents the repulsion force quantification constant for preventing overlapping between the phase fields of the cells, g e represents the repulsion force quantification constant for preventing overlapping between the phase field of the cell and the phase field of the space outside the constraint boundary; φ e is a given constraint boundary phase field, which is 0 inside the boundary and 1 outside the boundary.

[0023] In the present application, the cell volume control adopts relative error precision control (formula (5)) to avoid the phenomenon of phase field disappearance caused by too small grid:

[0024]

[0025] In the formula, M represents the volume constraint strength quantification constant; represents the unit vector perpendicular to the boundary of the cell phase field and pointing to the inside of the cell; r represents the position vector of any point in space; V i(t) represents the ideal volume given to the i-th cell at time t.

[0026] The present application adopts numerical method (here refers to second-order semi-implicit numerical format) to discretize the phase field evolution equation (formula (1)), so that the numerical calculation still ensures the calculation stability when using lower time resolution (time step), thereby realizing the calculation acceleration.

[0027] Therefore, the evolution equation (formula (1)) is rewritten as:

[0028]

[0029]

[0030] The semi-implicit evolution equation is obtained through the second-order semi-implicit discretization processing, as shown in formula (8):

[0031]

[0032] In the above formula (8), the upper right subscript m represents the calculation step number (m is the m-th step, and m+1 is the m+1-th step); the so-called semi-implicit refers to the implicit processing of the linear term βΔφ i in formula (6) as The nonlinear term F i is processed through the explicit Adams-Bashforth method as [Yang X. Error analysis of stabilized semi-implicit method of Allen-Cahn equation. Discrete Continuous Dyn. Ser. B 11 (4) (2009) 1057-1070]. The above semi-implicit evolution equation has higher stability, allowing the use of larger time step. The stability term is an additional introduced dissipation term, which balances the calculation instability caused by the explicit processing of the nonlinear term and improves the stability of the semi-implicit evolution equation; wherein S is a positive value parameter (i.e. stability term coefficient) for maintaining the calculation stability. In summary, the second-order time discretization can successfully obtain higher numerical precision. According to the above steps, iterative calculation is carried out until the multi-cell system reaches a stable state or the morphological characteristics change significantly, i.e. the cell morphology of each cell is obtained, and the cell population morphology is obtained by all cell morphologies. Here, the steady state can be defined as the root mean square speed of all cell centroids being less than a set threshold.

[0033] The present application further provides a system for accurately and efficiently simulating the morphology of a multi-cell system composed of a plurality of non-dividing cells, which comprises:

[0034] The first initial parameter generation module is used to provide the initial state of several non-dividing cells and phase field environment parameters.

[0035] The first cell morphology update module is used to iteratively solve the following phase field evolution equation based on the initial state of the cells and phase field environment parameters to obtain the morphology of several non-dividing cells and cell populations.

[0036]

[0037] In the formula, φ i F represents the phase field of the i-th cell; ten Indicates cell surface tension; F atr F represents the attraction between the i-th cell and the remaining cells; rep F represents the repulsive force between the i-th cell and the other cells and the constraint boundary; vol The force controlling cell volume is represented by t; time is represented by τ; the quantization constant of the environmental viscosity coefficient is represented by ξ. i (t) represents the motion noise assigned to the i-th cell, and κ represents the noise intensity quantization constant.

[0038] This invention further addresses the case of multicellular systems containing divisible cells by providing a method for accurately and efficiently simulating multicellular morphology. First, it provides the initial cell state and phase field environmental parameters of a system containing at least one divisible cell. Then, using the phase field of divisible cells in the current iteration as the parent cell phase field, it divides the parent cell phase field using a given plane to obtain the initial phase fields of the two daughter cells corresponding to the parent cell in the current iteration. Next, based on the initial phase fields of each daughter cell generated in the current iteration and the cell phase fields of undivided cells, iteratively solving the phase field evolution equations satisfied by each cell yields the cell morphology and movement of each cell in the current iteration, thus realizing the morphological evolution of the entire multicellular system.

[0039] Cell division is an iterative process, and the cells that divide in the previous iteration can serve as the parent cells for the next iteration. Let there be Q cells after step L; of which P cells (1 ≤ P ≤ Q) divide in the (L+1)th iteration. Then, the cell morphology evolution steps of the multicellular system in the (L+1)th iteration are as follows:

[0040] S1′ takes any divisible cell p as the parent cell and obtains the initial phase field of two daughter cells by dividing in a given plane; repeat this operation until P cells have completed division, p = 1, 2, ..., P;

[0041] S2' splits the initial phase field of all daughter cells (numbering 2P) and the cell phase field of the non-splitting cells (numbering Q-P) according to the L+1th iteration process, and iteratively solves the phase field evolution equation satisfied by each cell to obtain the cell morphology and movement of each cell in the L+1th iteration process, that is, to realize the morphological evolution of the entire multicellular system.

[0042] In the above step S1', the initial phase fields of the two daughter cells are as follows:

[0043]

[0044]

[0045] In the formula, represents the phase field of the pth parent cell, and represents the initial phase field of the two newly born daughter cells in the L+1th iteration; Ω represents a given calculation grid, and r represents the position vector of any point in space, represents the phase field center (center of mass) of the parent cell; ε represents a quantitative constant for adjusting the gap between the phase fields of the two daughter cells; n represents a unit vector along the cell division direction; g is determined by minimizing the function , wherein V p,D1 and V p,D2 are the volumes of the two daughter cells, which can be obtained by experiment or arbitrarily set according to the purpose of morphological evolution of the multicellular system.

[0046] In the above step S2', each cell satisfies the phase field evolution equation (1). The semi-implicit evolution equation (formula (8)) is obtained by using the semi-implicit discrete processing method given above, and then the initial phase fields of all daughter cells (numbering 2P) and the cell phase field of the non-splitting cells (numbering Q-P) are obtained according to the L+1th iteration process splitting, and the semi-implicit evolution equation (8) is iteratively solved until the multicellular system evolution meets the set requirements.

[0047] The present application further provides a system for accurately and efficiently simulating multicellular morphology, which comprises:

[0048] A second initial parameter generation module for generating an initial cell state and phase field environment parameters of a multicellular system containing at least one splittable cell;

[0049] A cell splitting module for splitting the parent cell phase field in the phase field of all cells in the current iteration process to obtain the initial phase field of the daughter cells corresponding to the parent cell in the current iteration process by giving a plane;

[0050] A second cell morphology updating module is configured to perform iterative solution of the phase field evolution equation satisfied by each cell according to the initial phase field of each offspring cell generated in the current iteration process and the cell phase field of the non-dividing cell, to obtain the cell morphology and movement of each cell in the current iteration process, and to realize morphology evolution of the entire multicellular system.

[0051] Compared with the prior art, the present application has the following beneficial effects:

[0052] (1) The present application firstly gives the initial state of the multicellular system and the phase field environment parameters; then performs iterative solution of the phase field evolution equation satisfied by each cell to obtain the cell morphology (i.e. phase field); and uses the diffusible field to describe the cell, so that the smooth, continuous and real cell morphology boundary can be obtained; compared with other methods of simplifying the cell into a point (coarse-grained model) or a convex polyhedron (vertex model), the present application uses dense grids to accurately describe the morphology boundary and controls the relative error of the cell volume, so that the present application has higher reliability in calculating the cell morphology and intercellular interaction.

[0053] (2) The present application can calculate the morphology dynamics of a natural living multicellular system; for example, using the nematode embryo development system as a case and precision calibration basis, the present application can reproduce the cell morphology and movement, reproduce the embryo structure, include the conserved cell connection map observed in experiments, and inversely deduce the cell mechanics distribution of the real system through comparison of the experimental structure.

[0054] (3) The present application can calculate the morphology dynamics of an artificial living multicellular system; for example, using the synNotch artificial system as a case, the present application can reproduce the rich morphology dynamics characteristics exhibited by the system under different viscosity settings, different cell proportions and initial states; so as to provide biophysical parameters with guiding significance for artificial multicellular biological engineering and assist in carrying out artificial multicellular biological engineering.

[0055] (4) The present application has excellent precision and high efficiency, can accurately capture the cell morphology characteristics (such as cell contact relationship and cell contact area), and thus assist in judging the mechanical state of the cell; in addition to the calculation precision, the present application needs very short time to process a cell system with a scale of about 100 cells, and can realize the scaling of calculation to carry out parameter scanning and other virtual experiments.

[0056] (5) The present application can be applied to other artificial multicellular systems (such as organoids, embryo-like embryos, etc.), has wide market demand and application prospect in assisting the construction, modification, design of artificial multicellular systems and in organ tissue engineering, drug screening, etc., and is suitable for popularization and use in the field. BRIEF DESCRIPTION OF DRAWINGS

[0057] Figure 1 Morphology output vs. experimental observation for synNotch artificial system simulation (Example 1).

[0058] Figure 2 Flowchart for nematode embryonic development system simulation (Example 3).

[0059] Figure 3 Cell division sequence for nematode embryonic development system simulation (Example 3).

[0060] Figure 4 Cell morphology output for nematode embryonic development system simulation (Example 3).

[0061] Figure 5 Simulation-experiment effect comparison chart for phase field model and method for precise and efficient simulation of multicellular morphology. DETAILED DESCRIPTION

[0062] The technical solutions of the embodiments of the present application will be described clearly and completely in combination with the accompanying drawings. Obviously, the described embodiments are only some of the embodiments of the present application, not all the embodiments. Based on the embodiments in the present application, all other embodiments obtained by those of ordinary skill in the art without creative labor belong to the present application.

[0063] Example 1

[0064] This embodiment is directed to a multicellular system composed of several non-dividing cells, and the method for precise and efficient simulation of multicellular morphology provided by the present application is explained in detail. Here, the cells cannot divide, but can only deform and move under the action of random motion noise (Gaussian white noise is used here) and mechanical interaction.

[0065] At the same time, this embodiment takes the synNotch artificial system as an example to explain in detail the method for precise and efficient simulation of multicellular morphology provided by the present application in the aspect of artificial system.

[0066] SynNotch is one of the most advanced synthetic multicellular machines, which controls cell adhesion and cell-cell interaction by encoding surface adhesion proteins in mouse cells [S. Toda, L. R. Blauch, S. K. Y. Tang, L. Morsut, W. A. Lim. Programming self-organizing multicellular structures with synthetic cell-cell signaling. Science 361(6398) (2018) 156-162]. By encoding different types of adhesion in two types of cells, the cell population can obtain separate mosaic structures and inner and outer layered structures; the manufactured multicellular machine can self-repair when encountering external mechanical damage, and automatically degrade when the adhesion code is removed.

[0067] The method for accurately and efficiently simulating multicellular morphology provided by the embodiment comprises the following steps:

[0068] S1 Given a number of non-dividing cell initial states and phase field environment parameters.

[0069] The cell initial state is: first, uniformly scatter points (≥120 cells) in three-dimensional space (xyz rectangular coordinate system) as the initial position of the cell population phase field, and then give each initial point a spherical phase field with a radius of R as the initial cell morphology. Here, the uniformly scattered points are within the artificially set constraint boundary.

[0070] The set phase field environment parameters include: constraint boundary (φ e ), phase field boundary thickness quantization constant (c), cell surface tension quantization constant (β), constraint boundary and cell repulsion quantization constant (g e ), cell and cell repulsion quantization constant (g), cell and cell attraction quantization constant (σ i,j ), environmental viscosity coefficient quantization constant (τ), volume constraint strength quantization constant (M), cell volume function (V i (t)), noise intensity quantization constant of Gaussian white noise (κ), etc., as shown in Table 1. Cell and cell attraction quantization constant (σ i,j) represents the cell adhesion between different cells, which is related to the target cell, according to the real experimental configuration (from [S. Toda, L. R. Blauch, S. K. Y. Tang, L. Morsut, W. A. Lim. Programming self-organizing multicellular structures with synthetic cell-cell signaling. Science 361 (6398) (2018) 156-162]) and the specific settings are shown in Table 2.

[0071] The embodiment further sets the spatial resolution (spatial step size δl) and the temporal resolution (temporal step size δt).

[0072] Table 1 Phase field environment parameters for simulating synNotch artificial system

[0073]

[0074]

[0075] S2 iteratively solves the corresponding phase field evolution equation according to the initial state of the cell and the phase field environment parameters to obtain the morphology (i.e. phase field) of a number of non-dividing cells and cell groups.

[0076] In this embodiment, the constraint boundary (φ e ) is set as a cubic space consistent with the experimental observation scale, and the edge length of the cubic space is 64 microns. The multicellular system is set to be unable to divide, and the simulation of the split tiling phenomenon (three types of cells), the split tiling phenomenon (two types of cells), and the inner-outer layered phenomenon (two types of cells) is realized through different types of adhesion coding.

[0077] Since each cell satisfies the phase field evolution equation (1), the noise used in this embodiment is Gaussian white noise, δt represents the time step; u i (t) represents the random noise of the centroid velocity of the i-th cell, which is obtained by randomly assigning values satisfying the standard normal distribution to x, y, and z directions respectively. The semi-implicit evolution equation (formula (8)) is obtained by using the semi-implicit discrete processing method given above, the cell adhesion is set according to the split tiling phenomenon (three types of cells, including the initial and final two cases), the split tiling phenomenon (two types of cells), and the inner-outer layered phenomenon (two types of cells) in Table 2, and then the semi-implicit evolution equation (8) is iteratively solved, and the operation is performed according to the virtual time of 50000. The operation results are shown in Figure 1 Table 2. The cell type division and time consumption are shown in Table 2.

[0078] Table 2 Cell type partition and time consumption record of synNotch artificial system

[0079]

[0080] Note: 3 types of cells (blue) in the segregation mosaic phenomenon (three types of cells) are defined as type 1 cells (red) in contact with type 2 cells (green), which are colored from virtual time = 40000.

[0081] Note: The calculation unit uses a GPU (NVIDIA Tesla P100)

[0082] From Figure 1 As can be seen from the first to third rows, all three phenomena can be simulated and reproduced, and are highly consistent with experimental observations.

[0083] For the simulation of the self-repairing phenomenon (two types of cells), the cell type partition and time consumption record are shown in Table 2. At the beginning of the simulation, the simulation end state of the inner-outer stratification phenomenon (two types of cells) is used as the initial cell population structure, and a small amount of cells free outside the central cell population and the centroid position below z = 0 (i.e., the midpoint of the longitudinal coordinate) are deleted to simulate a cell population subjected to mechanical cutting, and then the semi-implicit evolution equation (8) is iteratively solved according to the virtual time of 50000, and the calculation result is shown in Figure 1 The fourth row. The cell type partition and time consumption record are shown in Table 2. It can be seen that the simulation of the self-repairing phenomenon can be realized.

[0084] For the simulation of the automatic degradation phenomenon (two types of cells), the cell type partition and time consumption record are shown in Table 2. At the beginning of the simulation, the simulation end state of the inner-outer stratification phenomenon (two types of cells) is used as the initial cell population structure, and the intercellular adhesion is set to the lowest adhesion value of the system 0.3, as shown in Table 2; then the semi-implicit evolution equation (8) is iteratively solved according to the virtual time of 50000, and the calculation result is shown in Figure 1 The fifth row. The cell type partition and time consumption record are shown in Table 2. It can be seen that the originally stratified cell population gradually returns to a scattered and random state, which is highly consistent with experimental observations.

[0085] It can also be seen from the present embodiment that when calculating a noisy system of 120 cells or less, it takes no more than 13.5 hours to run from the initial state to a significant change in shape on the GPU and reproduce experimental observations. It shows that the method provided by the present application for accurately and efficiently simulating the morphology of multiple cells has high computational efficiency.

[0086] Example 2

[0087] The embodiment is directed to a multicellular system composed of a plurality of non-dividing cells, and provides a system for accurately and efficiently simulating multicellular morphology, which comprises:

[0088] A first initial parameter generation module is configured to give initial states of a plurality of non-dividing cells and phase field environment parameters.

[0089] A first cell morphology update module is configured to iteratively solve the phase field evolution equation (1) satisfied by the initial states of the cells and the phase field environment parameters to obtain the morphology of the plurality of non-dividing cells and cell groups.

[0090] The system for accurately and efficiently simulating multicellular morphology can be set in a computer with computing power. The first initial parameter generation module is configured to give initial states of a plurality of non-dividing cells and phase field environment parameters. The first cell morphology update module is configured to iteratively solve the phase field evolution equation corresponding to the multicellular system to obtain the morphology of the plurality of non-dividing cells and cell groups according to the initial states of the cells and the phase field environment parameters. The specific operation can refer to steps S1 and S2 in the method for accurately and efficiently simulating multicellular morphology given in Embodiment 1.

[0091] Embodiment 3

[0092] The embodiment is directed to a case where the multicellular system contains dividing cells, and provides a method for accurately and efficiently simulating multicellular morphology. First, the initial cell state and phase field environment parameters of the multicellular system containing at least one dividing cell are given. Then, the phase field of the dividing cell in the phase field of all cells in the current iteration process is taken as the mother cell phase field, and the mother cell phase field is divided by a given plane to obtain the initial phase field of the two daughter cells corresponding to the mother cell in the current iteration process. Then, the initial phase field of each daughter cell generated in the current iteration process and the cell phase field of the non-dividing cell are used to iteratively solve the phase field evolution equation satisfied by each cell to obtain the cell morphology of each cell in the current iteration process, i.e., to realize the morphology evolution of the entire multicellular system.

[0093] The embodiment takes the embryonic development of C. elegans (i.e., the cell division process) as an example to illustrate the method for accurately and efficiently simulating multicellular morphology in natural systems.

[0094] The embodiment provides a method for accurately and efficiently simulating the embryonic development of C. elegans. The cell division sequence (including cell identity and division time), cell division direction, and volume allocation ratio are obtained by experiment shooting. Then, the cell morphology and cell movement are obtained by the method for accurately and efficiently simulating multicellular morphology. The flowchart is shown in Figure 2 .

[0095] The cell division sequence data used in this embodiment comes from [G.Guan,MKWong,VWSHo,X.An,LYChan,B.Tian,Z.Li,LHTang,Z.Zhao,C.Tang.System-level quantification and phenotyping of early embryonic morphogenesis of Caenorhabditis elegans.bioRxiv(2019)776062]. In this literature, green fluorescent protein (GFP) is used to label the cell nucleus, and three-dimensional time-lapse fluorescence imaging is used to trace cell lineage, thereby obtaining the cell division sequence (including cell identity and division time). See Figure 2 As shown.

[0096] Data on cell division direction and volume fraction were obtained from [J. Cao, G. Guan, VWS Ho, MK Wong, LY Chan, C. Tang, Z. Zhao, H. Yan. Establishment of a morphological atlas of the Caenorhabditis elegans embryo using deep-learning-based 4D segmentation. Nat. Commun. 11(2020)6254]. In this paper, red fluorescent protein (mCherry) was used to label the cell membrane to segment cells and obtain volume fraction; at the same time, green fluorescent protein (GFP) was used to label the cell nucleus to obtain the cell division direction.

[0097] In this embodiment, as Figure 3 As shown, the cell lineage tree represents 24 cell divisions before the 102-cell stage, serving as the phase division for phase-field simulation; except for the 6th and 7th cell stages, which are calculated up to the quasi-steady state (defined as: in, For the root mean square rate of the centroids of all cells, t q The calculation time for the quasi-steady-state point) and the 8-cell phase is a virtual time of 15000. All phases are calculated to the steady state.

[0098] The method for accurately and efficiently simulating multicellular morphology provided in this embodiment includes the following steps:

[0099] Step 1: Given the initial state of the initial cells and the phase field environment parameters.

[0100] In this embodiment, the initial morphology of the initial cell P0 is given, including the initial boundary shape and position of the cell, such as...Figure 2

[0101] In this embodiment, the set phase field environment parameters include: constraint boundary (φ e ), phase field boundary thickness quantization constant (c), cell surface tension quantization constant (β), constraint boundary and cell repulsion quantization constant (g e ), cell and cell repulsion quantization constant (g), cell and cell attraction quantization constant (σ i,j ), environment viscosity coefficient quantization constant (τ), volume constraint strength quantization constant (M), cell volume function (V i (t)), noise intensity quantization constant of Gaussian white noise (κ), etc., as shown in Table 3. Cell and cell attraction quantization constant (σ i,j ) represents the cell adhesion between different cells, which is related to the target cell, and is configured according to the real experiment (see from the literature [K. Yamamoto, A. Kimura. An asymmetric attraction model for the diversity and robustness of cell arrangement in nematodes. Development 144(23) (2017) 4437-4449; P. Dutta, D. Odedra, C. Pohl. Planar asymmetries in the C. elegans embryo emerge by differential retention of a PAR at cell-cell contacts. Front. Cell Dev. Biol. 7 (2019) 209]) and the following analysis.

[0102] The constraint boundary (φ e ) is set as an ellipsoidal eggshell consistent with experimental observation, and the semi-axes of the ellipsoid in the x / forward and backward axis, y / left and right axis, and z / ventral and dorsal axis are 28.28 microns, 12.87 microns, and 18.81 microns, respectively. The place 10.21 microns away from the origin on both sides of the y / left and right axis is truncated by a section parallel to the xy plane to simulate the external extrusion of the embryo in the real fluorescence shooting process.

[0103] The embodiment further sets the spatial resolution (spatial step size δl) and the time resolution (time step size δt).

[0104] Table 3 Phase field environment parameters of the simulated nematode embryo development system

[0105]

[0106] ​Note: (1) The virtual time length m is judged at steady state t The number of calculation steps for cyclic solution of formula (8) is Where FLOOR is the floor function; every time steps, it is determined whether the cell movement reaches steady state;

[0107] (2) The root mean square rate of the multicellular system at steady state Where r c,i represents the phase center (centroid) of the i-th cell.

[0108] Step two, cell division

[0109] In this step, the phase field of all cells in the current iteration process that can be divided is the mother cell phase field, which is divided by a given plane to obtain the initial phase field of the two daughter cells corresponding to the mother cell in the current iteration process; then, according to the initial phase field of each daughter cell generated in the current iteration process and the cell phase field of the non-dividing cell, the phase field evolution equation satisfied by each cell is iteratively solved to obtain the cell morphology and movement of each cell in the current iteration process, that is, the morphological evolution of the entire multicellular system is realized.

[0110] Cell division is an iteration process, and the cells divided in the last iteration process can be used as the mother cells in the next iteration process; let there be Q cells after step L; among them, there are P cells (1≤P≤Q) that divide in the iteration process of step L+1, then the cell morphological evolution steps of the multicellular system in the iteration process of step L+1 are as follows:

[0111] S1' takes any cell p that can be divided as a mother cell, and divides it by a given plane to obtain the initial phase field of two daughter cells; repeat this operation until P cells complete division, p = 1, 2, …, P;

[0112] S2' according to the initial phase field of all daughter cells (the number is 2P) obtained by division in the iteration process of step L+1 and the cell phase field of non-dividing cells (the number is Q-P), iteratively solve the phase field evolution equation satisfied by each cell to obtain the cell morphology and movement of each cell in the iteration process of step L+1, that is, realize the morphological evolution of the entire multicellular system.

[0113] In the above step S1', the initial phase fields of the two daughter cells are as follows:

[0114]

[0115]

[0116] In the formula, represents the phase field of the p-th mother cell, and Let represent the initial phase field of the two newborn daughter cells in the (L+1)th iteration; Ω represents the given computational grid, and r represents the position vector of any point in space. ε represents the phase field center (centroid) of the maternal cell; ε represents the quantization constant that adjusts the gap between the phase fields of the two daughter cells; n represents the unit vector along the cell division direction; b is obtained by minimizing the function Determined, among which, V p,D1 and V p,D2 The volume of two daughter cells was obtained experimentally in this embodiment.

[0117] In step S2′ above, each cell satisfies the phase field evolution equation (1). The semi-implicit evolution equation (Equation (8)) is obtained by using the semi-implicit discrete processing method given above. Then, based on the L+1 step iteration process, the initial phase field of all daughter cells (number of 2P) and the cell phase field of undivided cells (number of QP) are obtained by splitting. The semi-implicit evolution equation (8) is solved iteratively until the evolution of the multicellular system meets the set requirements.

[0118] During the simulation of the above cell division, it was found that key developmental processes include the establishment of the anterior and posterior body axis and the ventral-dorsolateral body axis in the 4-cell stage of the nematode embryo: such as Figure 5 As shown, the angle between the centroid of ABa cells and P2 cells and the x-axis is less than 10°; the angle between the centroid of ABP cells and EMS cells and the z-axis is less than 10°.

[0119] The establishment of cell contact relationships during the 4-cell phase includes the establishment of 5 cell contact pairs: ABa-ABp, ABa-EMS, ABP-P2, ABP-EMS, and EMS-P2.

[0120] Cell adhesion at the 4-cell stage: when the intercellular attraction σ i,j When all values ​​are set to 0, the calculated intercellular contact area is smaller than the experimental contact area; in this case, the quantitative requirement is the relative deviation of all five contact areas. The average value is less than 0, where h refers to the h-th contact surface. By uniformly increasing the intercellular attraction, the cell contact areas obtained from simulation and experiment can be made more similar; at this point, the quantitative requirement is that the average deviation of all 5 contact areas decreases after increasing the attraction.

[0121]

[0122] 4. The viscosity between specific cells (EMS) and specific cells (P2) during cell phase 4 is lower than the viscosity of other cell contact surfaces (ABa-ABp, ABa-EMS, ABP-P2, ABP-EMS): When the global cell viscosity (attraction) is set to a value that minimizes the overall difference in contact area among all cells (σ...S ) the deviation of EMS-P2 contact area is larger than that of other 4 contacts; at this time, the quantitative requirement is the average deviation of the first 4 contacts less than 20%, the deviation of EMS-P2 more than 40%. By adjusting the viscosity (σ W ) of EMS-P2, the deviation of this contact area can be made smaller, and the deviations of the final 5 pairs of contact areas are all less than 20%.

[0123] The above key developmental phenomena include the specific nematode embryo morphology and cell contact map at the 6, 7, and 8 cell stages. Among them, the end time of the 6 and 7 cell stages (the division time of the next cell) is set as the time point when the system kinetic energy reaches the quasi-steady state (defined as: wherein, is the root mean square speed of all cell centers, t q is the time point of quasi-steady state), and the virtual time (step size x step number) of the end of the 8 cell stage is 15000. In the simulation of the 6, 7, and 8 cell stages, the viscosity between sister cells (two cells derived from the same mother cell) is set to be weaker σ W , and the viscosity between non-sister cells (two cells derived from different mother cells) is set to be stronger σ S . In the simulation of the 8 cell stage, the cell pair ABpl-E is set to σ S , which cannot reproduce the stable real embryo structure, and is set to σ W , which can, as shown in Figure 5 .

[0124] Through the above analysis, the present embodiment can completely reconstruct the conserved cell contact map at the 2, 3, 4, 6, 7, and 8 cell stages. Specifically, the conserved contact of cell pairs is established by experimentally observing multiple embryos (all embryo samples exist), non-conserved contact (part of the embryo samples exist), and conserved non-contact (all embryo samples do not exist). In the multi-cell system obtained by simulation, all conserved contacts must be reproduced in the simulation structure, and all conserved non-contacts are not allowed to be established in the simulation, indicating that the method provided by the present application simulates the multi-cell morphology consistent with the experiment.

[0125] The included 4-cell stage EMS-P2 cells have low adhesion to the contact surface of the 8-cell stage ABpl-E cell pair, which has been confirmed by biological experiments [K. Yamamoto, A. Kimura. An asymmetric attraction model for the diversity and robustness of cell arrangement in nematodes. Development 144 (23) (2017) 4437-4449; P. Dutta, D. Odedra, C. Pohl. Planar asymmetries in the C. elegans embryo emerge by differential retention of aPARs at cell-cell contacts. Front. Cell Dev. Biol. 7 (2019) 209], indicating that the method provided by the present application can capture the real cell mechanical interaction, such as Figure 5 as shown.

[0126] Using a GPU (NVIDIA Tesla P100), the noise-free system was calculated within 102 cells according to the above steps, from the initial state to the steady state. The simulation results are shown in Figure 4 . The stable calculation of 25 stages can be achieved (see Table 4), and the time-consuming of calculation to the steady state is not more than 1 hour, and the total time-consuming is not more than 6.5 hours, indicating the high efficiency of the method provided by the present application.

[0127] Table 4: Simulation stage and time-consuming record of nematode embryo development system

[0128]

[0129]

[0130] Note: The calculation unit uses GPU (NVIDIA Tesla P100)

[0131] By controlling the spatial resolution (spatial step size δl, for example, δl≤0.5 microns) and the time resolution (time step size δt, for example, δt≤1.5), the calculation accuracy is ensured, and under the experimental conditions of C. elegans embryo (cell division sequence, cell division direction, volume distribution ratio), there is no calculation overflow, and the 1-8 cell stage embryo cell morphology observed in the experiment (such as Figure 5 ) can be reproduced, and some key development processes of the real embryo can be deduced.

[0132] The method for accurately and efficiently simulating the morphology of a multicellular system can conform to the real world.

[0133] Embodiment 4

[0134] The embodiment provides a system for accurately and efficiently simulating the morphology of a multicellular system, which comprises the following:

[0135] The second initial parameter generation module is configured to input the initial cell state of the multicellular system comprising at least one cell that can split and the phase field environment parameters.

[0136] The cell splitting module is configured to split the mother cell phase field into the initial phase fields of two daughter cells corresponding to the mother cell in the current iteration process by using the given splitting plane.

[0137] The second cell morphology updating module is configured to obtain the cell morphology and movement of each cell in the current iteration process by iteratively solving the phase field evolution equation satisfied by each cell according to the initial phase fields of the daughter cells generated in the current iteration process and the cell phase field of the non-splitting cell, that is, realize the morphology evolution of the entire multicellular system.

[0138] The system for accurately and efficiently simulating the morphology of a multicellular system can be arranged in a computer with computing capability, the initial cell state of the multicellular system comprising at least one cell that can split and the phase field environment parameters are input by using the second initial parameter generation module, the initial phase fields of two daughter cells corresponding to the mother cell in the current iteration process are obtained by splitting the mother cell phase field by using the cell splitting module, and the cell morphology and movement of each cell in the current iteration process are obtained by iteratively solving the phase field evolution equation satisfied by each cell according to the initial phase fields of the daughter cells generated in the current iteration process and the cell phase field of the non-splitting cell by using the second cell morphology updating module. The specific operation can refer to steps one and two in the method for accurately and efficiently simulating the morphology of a multicellular system in embodiment 3.

[0139] Those skilled in the art will understand that the embodiments described herein are intended to help the reader understand the principles of the present application and should be understood as not limiting the scope of protection of the present application to such specific statements and embodiments. Those skilled in the art can make various other specific modifications and combinations according to the technical inspiration disclosed in the present application without departing from the essence of the present application, and these modifications and combinations still fall within the scope of protection of the present application.

Claims

1. A method for accurate and efficient simulation of multicellular morphology, characterized in that, For a multicellular system consisting of several non-dividing cells, the steps are as follows: S1: Given the initial state of several non-dividing cells and the phase field environment parameters; S2: According to the initial state of the cells and the phase field environment parameters, the following phase field evolution equation is iteratively solved to obtain the morphology of several non-dividing cells and cell groups; (1)； wherein represents the phase field of the i-th cell; represents the cell surface tension; represents the attractive force between the i-th cell and the rest of the cells; represents the repulsive force between the i-th cell and the rest of the cells and the confining boundary; represents the force controlling the cell volume; represents time; represents the environmental viscosity quantification constant; represents the motion noise given to the i-th cell, represents the noise intensity quantification constant; The cell volume control satisfies the following equation: (5)； wherein represents a volume-constrained strength quantification constant; represents a unit vector normal to the cell phase boundary pointing towards the cell interior; r represents the position vector of an arbitrary point in space; represents the ideal volume of the i-th cell at time t; The phase field evolution equation is discretized by semi-implicit method to obtain the semi-implicit evolution equation, as shown in equation (8): (8)； (7)； wherein, φi denotes the phase field of the i-th cell; the upper right superscript m denotes the number of calculation steps; denotes the environmental viscosity quantization constant; denotes the time resolution; denotes the cell surface tension quantization constant; S denotes the stabilization term coefficient; denotes the quantization constant that controls the boundary width of the phase field; denotes the repulsive force quantization constant that prevents the overlap between the phase field of the cell and the phase field of the space outside the constraint boundary; denotes the repulsive force quantization constant that prevents the overlap between the phase fields of the cells; is the phase field of the given constraint boundary; denotes the attractive force quantization constant between the cell i and the cell j; denotes the volume constraint strength quantization constant; denotes the ideal volume of the i-th cell at time t; denotes the motion noise given to the i-th cell, denotes the noise strength quantization constant.

2. The method for precise and efficient simulation of multicellular morphology of claim 1, wherein, The initial state of the cell includes the initial boundary shape, position of the cell; The phase field environment parameters include the constraint boundary, the cell surface tension constant, the cell-cell attraction force constant, the environmental viscosity constant, the cell volume function, and the motion noise intensity constant.

3. The method for precise and efficient simulation of multicellular morphology of claim 1, wherein, In the above step S2, the cell mechanics under the phase field framework includes the cell surface tension and the attraction force between cells as follows: (2)； (3)； wherein represents the total number of cells in the multicellular system; represents the Laplacian operator, represents the gradient operator; represents the cell surface tension constant; represents the attractive force constant between cell i and cell j; is a double-well function that separates the order parameter into two states, 0 and 1, and its derivative is 0 represents the extracellular space and 1 represents the intracellular space; represents a quantization constant that controls the width of the order parameter interface.

4. The method for precise and efficient simulation of multicellular morphology of claim 1, wherein, The phase field of the cell group and the phase field of the cell and the constraint boundary satisfy the following equation: (4)； wherein represents a repulsive force per unit constant that prevents overlap between cell phase fields; represents a repulsive force per unit constant that prevents overlap between cell phase fields and phase fields outside the confinement boundary; is 0 inside the boundary and 1 outside the boundary for a given confinement boundary phase field.

5. A system for accurate and efficient simulation of multicellular morphology, characterized in that, For a multicellular system consisting of several non-dividing cells, including: A first initial parameter generation module for giving the initial state of several non-dividing cells and the phase field environment parameters; A first cell morphology update module for iteratively solving the following phase field evolution equation according to the initial state of the cells and the phase field environment parameters to obtain the morphology of several non-dividing cells and cell groups; (1)； wherein represents the phase field of the i-th cell; represents the cell surface tension; represents the attractive force between the i-th cell and the rest of the cells; represents the repulsive force between the i-th cell and the rest of the cells and the confining boundary; represents the force controlling the cell volume; represents time; represents the environmental viscosity quantification constant; represents the motion noise given to the i-th cell, represents the noise intensity quantification constant; The cell volume control satisfies the following equation: (5)； wherein represents a volume-constrained strength quantification constant; represents a unit vector normal to the cell phase boundary pointing towards the cell interior; r represents the position vector of an arbitrary point in space; represents the ideal volume of the i-th cell at time t; The phase field evolution equation is discretized by semi-implicit method to obtain the semi-implicit evolution equation, as shown in equation (8): (8)； (7)； wherein, φi denotes the phase field of the i-th cell; the upper right superscript m denotes the number of calculation steps; μ denotes an environmental viscosity quantization constant; Δt denotes a time resolution; σ denotes a cell surface tension quantization constant; S denotes a stabilization term coefficient; β denotes a quantization constant that controls the width of the phase field boundary; represents a repulsive force quantization constant that prevents the overlap between the phase field of the cell and the phase field of the space outside the constraint boundary; represents a repulsive force quantization constant that prevents the overlap between the phase fields of the cells; φb denotes a given constraint boundary phase field; σij denotes an attractive force quantization constant between the cell i and the cell j; represents a volume constraint strength quantization constant; Videalsi(t) denotes a given ideal volume of the i-th cell at time t; represents a motion noise given to the i-th cell, represents a noise strength quantization constant.

6. A method for accurate and efficient simulation of multicellular morphology, characterized in that, For the case of containing dividing cells in the multicellular system, first, the initial cell state and the phase field environment parameters of the system containing at least one dividing cell are given; Then, the initial phase field of the two daughter cells corresponding to the mother cell in the current iteration process is obtained by dividing the mother cell phase field with a given plane; Then, according to the initial phase field of each daughter cell generated in the current iteration process and the cell phase field of the non-dividing cell, the phase field evolution equation satisfied by each cell is iteratively solved to obtain the cell morphology and motion of each cell in the current iteration process, that is, the morphological evolution of the entire multicellular system is realized; The phase field evolution equation is represented as: (1)； wherein represents the phase field of the i-th cell; represents the cell surface tension; represents the attractive force between the i-th cell and the rest of the cells; represents the repulsive force between the i-th cell and the rest of the cells and the confining boundary; represents the force controlling the cell volume; represents time; represents the environmental viscosity quantification constant; represents the motion noise given to the i-th cell, represents the noise intensity quantification constant; The cell volume control satisfies the following equation: (5)； wherein represents a volume-constrained strength quantification constant; represents a unit vector normal to the cell phase boundary pointing towards the cell interior; r represents the position vector of an arbitrary point in space; represents the ideal volume of the i-th cell at time t; The phase field evolution equation is discretized by semi-implicit method to obtain the semi-implicit evolution equation, as shown in equation (8): (8)； (7)； wherein, φi represents the phase field of the i-th cell; the upper right superscript m represents the number of calculation steps; represents an environmental viscosity quantization constant; represents a time resolution; represents a cell surface tension quantization constant; S represents a stabilization term coefficient; represents a quantization constant that controls the boundary width of the phase field; represents a repulsive force quantization constant that prevents overlap between the phase field of the cell and the phase field of the space outside the constraint boundary; represents a repulsive force quantization constant that prevents overlap between the phase fields of the cells; is a given constraint boundary phase field; represents an attractive force quantization constant between cell i and cell j; represents a volume constraint strength quantization constant; represents a given ideal volume of the i-th cell at time t; represents a motion noise given to the i-th cell, represents a noise strength quantization constant.

7. The method for accurately and efficiently simulating the morphology of a multicellular system according to claim 6, characterized in that, The cell division is an iterative process, the cell after the last iteration process can be used as the mother cell of the next iteration process; there are Q cells after the Lth step; wherein P cells divide in the L+1th iteration process, The cell morphological evolution step of the multicellular system in the L+1th iteration process is as follows: S1′: Take any one of the dividing cells p as the mother cell, and divide it by a given plane to obtain the initial phase field of the two daughter cells; The initial phase field of the two daughter cells is as follows: (9)； (10)； In the formula, This represents the phase field of the p-th maternal cell. and Indicates the first The initial phase field of two newly generated daughter cells in the step iteration; Represents a given computational grid. Let represent the position vector of any point in space. Represents the phase field center of the maternal cell; This represents the quantization constant that regulates the gap between the phase fields of two daughter cells; Represents a unit vector along the direction of cell division; By minimizing the function It is confirmed that, among them, and The volume of two daughter cells; This operation is repeated until P cells complete division, ; S2′: According to the initial phase field of all daughter cells obtained by dividing in the L+1th iteration process and the cell phase field of the non-dividing cell, the phase field evolution equation satisfied by each cell is iteratively solved to obtain the cell morphology and motion of each cell in the L+1th iteration process, that is, the morphological evolution of the entire multicellular system is realized.

8. A system for accurate and efficient simulation of multicellular morphology, characterized in that, For the case of containing dividing cells in the cell system, including: A second initial parameter generation module is configured to generate initial cell states and phase field environment parameters of a multicellular system containing at least one splittable cell; A cell splitting module is configured to split a mother cell phase field of a splittable cell in all cell phase fields in a current iteration process to obtain initial phase fields of two daughter cells corresponding to the mother cell in the current iteration process by a given splitting plane; A second cell morphology update module is configured to generate initial phase fields of each daughter cell and cell phase fields of non-splittable cells in the current iteration process, and to solve a phase field evolution equation satisfied by each cell to obtain cell morphologies and movements of each cell in the current iteration process, i.e., to realize morphological evolution of the entire multicellular system. The phase field evolution equation is expressed as: (1)； wherein represents the phase field of the i-th cell; represents the cell surface tension; represents the attractive force between the i-th cell and the remaining cells; represents the repulsive force between the i-th cell and the remaining cells and the confining boundary; represents the force controlling the cell volume; represents time; represents the environmental viscosity quantification constant; represents the motion noise given to the i-th cell, represents the noise intensity quantification constant; The cell volume control satisfies the following equation: (5)； wherein represents a volume-constrained strength quantification constant; represents a unit vector normal to the cell phase boundary pointing towards the cell interior; r represents the position vector of an arbitrary point in space; represents the ideal volume of the i-th cell at time t; The phase field evolution equation is discretized by a semi-implicit method to obtain a semi-implicit evolution equation, as shown in equation (8): (8)； (7)； wherein, φi represents the phase field of the i-th cell; the upper right superscript m represents the number of calculation steps; represents an environmental viscosity quantization constant; represents a time resolution; represents a cell surface tension quantization constant; S represents a stabilization term coefficient; represents a quantization constant that controls the boundary width of the phase field; represents a repulsive force quantization constant that prevents overlap between the phase field of the cell and the phase field of the space outside the constraint boundary; represents a repulsive force quantization constant that prevents overlap between the phase fields of the cells; is a given constraint boundary phase field; represents an attractive force quantization constant between cell i and cell j; represents a volume constraint strength quantization constant; represents a given ideal volume of the i-th cell at time t; represents a motion noise given to the i-th cell, represents a noise strength quantization constant.

Images