Cell development process dynamic modeling method and device based on time sequence single cell transcriptome data and medium

Abstract

The invention provides a cell development process dynamic modeling method and device based on time sequence single cell transcriptome data and a medium, and relates to the crossing field of bioinformatics and computational biology. The method comprises the following steps: constructing a Shenchang differential equation learning framework; adjusting parameters of the single cell development state change model based on the Shenxuan differential equation learning framework so as to construct a population cell development state change model; obtaining a cell specific gene regulation network and a population cell gene regulation network based on the population cell development state change model so as to predict occurrence opportunity of cell lineage differentiation and a molecular decision mechanism of cell differentiation; therefore, the problems of incomplete modeling mechanism, insufficient noise processing and lack of energy principle in the existing cell development process are solved.

Description

TECHNICAL FIELD

[0001] The present application relates to the field of bioinformatics and computational biology, in particular to a cell development process dynamic modeling method and device based on time-series single-cell transcriptome data and a medium. BACKGROUND

[0002] In recent years, the rapid development of single-cell RNA sequencing (scRNA-seq) technology enables researchers to analyze the heterogeneity of tissues and organs at single-cell resolution. Especially in developmental biology, by collecting samples at multiple time points, time-series single-cell transcriptome data is obtained, which provides an unprecedented data basis for studying cell fate determination, lineage differentiation and tissue formation. However, due to the destructive sampling characteristics of single-cell sequencing, i.e., each cell can only be measured once and its subsequent state cannot be tracked, traditional static data analysis methods are difficult to reveal the continuous evolution dynamics of cell state.

[0003] Therefore, the academic community has proposed the concepts of "cell trajectory inference" or "pseudotime ordering", aiming to construct the evolution path between cells according to gene expression similarity. Early methods such as Monocle, Slingshot, PAGA, etc. mainly rely on geometric connection in graph structure or dimension reduction space for path inference, although they can reveal the topological structure, they lack modeling of the physical mechanism of the development process, and are difficult to explain dynamic behaviors such as bifurcation decision and speed change.

[0004] In addition, some studies attempt to introduce dynamic models, such as methods based on RNA velocity (scVelo, VELOVI), which estimate the instantaneous gene expression change direction by the ratio of precursor mRNA (messenger ribonucleic acid) to mature mRNA, and preliminarily realize the modeling of the velocity field. In addition, there are also works trying to apply optimal transport theory to time-series cell distribution matching (such as Waddington-OT, CellOT) to infer population evolution paths. However, existing methods still generally face problems such as incomplete modeling mechanism, insufficient noise processing, and lack of energy principle. SUMMARY

[0005] To solve the problems of incomplete modeling mechanism, insufficient noise processing and lack of energy principle in existing cell development process modeling, a dynamic modeling method of cell development process based on time series single cell transcriptome data is proposed; the evolution of cell population in gene expression space is regarded as a controlled random dynamic process, and under the guidance of the principle of energy minimization, the end-to-end learning of high-dimensional trajectory is realized through neural differential equation.

[0006] The first aspect of the present application provides a dynamic modeling method of cell development process based on time series single cell transcriptome data, characterized by comprising:

[0007] constructing a single cell development state change model based on time series single cell transcriptome data;

[0008] constructing a neural differential equation learning framework;

[0009] adjusting the parameters of the single cell development state change model based on the neural differential equation learning framework to construct a population cell development state change model;

[0010] obtaining cell-specific gene regulation network and population cell gene regulation network based on the population cell development state change model to predict the occurrence time of cell lineage differentiation and the molecular decision mechanism of cell differentiation.

[0011] In one possible design, the single cell development state change model is constructed based on time series single cell transcriptome data, comprising:

[0012] Based on time series single cell transcriptome data, a stochastic differential equation (SDE) is used to model the cell development process to obtain the single cell development state change model:

[0013] ,

[0014] wherein, represents the M-dimensional gene expression vector of the cell at time t, the drift term represents the driving term of the gene regulation network, which describes the deterministic influence of gene interaction on cell state, and f is the velocity field; the diffusion term , is the noise intensity matrix, which reflects the technical noise and biological random fluctuation, , represents the standard Wiener process, which simulates the random disturbance in the development process; represents an M-row S-column real number, represents an S-row real number.

[0015] In one possible design, the constructing the neural ordinary differential equation learning framework includes:

[0016] constructing a population cell probability density distribution evolution simplified equation for representing population cell paths;

[0017] defining a population cell total energy cost discrete equation from an initial distribution to a terminal distribution based on the population cell probability density;

[0018] defining a loss function of the neural ordinary differential equation learning framework based on the population cell total energy cost discrete equation.

[0019] In one possible design, the constructing the population cell probability density distribution evolution simplified equation includes:

[0020] defining an evolution of a population cell probability density function following an extended Fokker-Planck equation, constructing a population cell probability density distribution evolution initial equation based on the population cell probability density;

[0021] defining an effective velocity field, simplifying the population cell probability density distribution evolution initial equation based on the effective velocity field to obtain the population cell probability density distribution evolution simplified equation.

[0022] In one example, the constructing the population cell probability density distribution evolution initial equation based on the population cell probability density includes introducing a source term representing an effect of a local proliferation rate on a population density, the population cell probability density distribution evolution initial equation being:

[0023] ,

[0024] wherein, is a population probability density, is a diffusion coefficient matrix, , is an M by M real number, denotes a divergence;

[0025] defining an effective velocity field , simplifying the population cell probability density distribution evolution initial equation based on the effective velocity field to obtain the population cell probability density distribution evolution simplified equation, .

[0026] In one possible design, the defining the population cell total energy cost discrete equation from an initial distribution to a terminal distribution includes:

[0027] The definition of the cell development trajectory follows the principle of energy optimization, and the total energy cost equation of the population cells from the initial distribution to the terminal distribution is defined based on the set of all possible continuous evolution paths of the population cells from the initial distribution to the terminal distribution.

[0028] The initial high-dimensional space particle point set is randomly sampled, and the continuous total energy cost equation is converted into the total energy cost discrete equation according to the evolution of the flow simulation probability density of the particles.

[0029] In one example, the definition of the cell development trajectory follows the principle of energy optimization, and the unbalance optimal transport theory framework is introduced. The total energy cost equation of the population cells from the initial distribution to the terminal distribution is defined based on the set of all possible continuous evolution paths of the population cells from the initial distribution to the terminal distribution. The total energy cost equation of the population cells from the initial distribution to the terminal distribution is:

[0030] ,

[0031] wherein, is the initial distribution, is the terminal distribution, is the energy cost, T is the time corresponding to the terminal distribution, is an M-row real number, is the kinetic energy cost required for driving, is the mass change cost of the proliferation activity, is a set weight coefficient;

[0032] The total energy cost discrete equation is:

[0033] .

[0034] In one possible design, the definition of the loss function of the neural ordinary differential equation learning framework based on the total energy cost discrete equation includes:

[0035] The reconstruction loss of the time sequence probability density distribution is defined based on the total energy cost of the population cells from the initial distribution to the terminal distribution to define the loss function of the neural ordinary differential equation learning framework.

[0036] The reconstruction loss of the time sequence probability density distribution includes adjacent time reconstruction error and long-term time reconstruction error.

[0037] In one example, the loss function of the neural ordinary differential equation learning framework is:

[0038] ,

[0039] wherein, a reconstruction loss for the time-series probability density distribution, a reconstruction loss weight coefficient for the set time-series probability density distribution;

[0040] ,

[0041] wherein, denotes the probability density distribution corresponding to the real gene expression state x at time denotes the probability density distribution corresponding to the real gene expression state x at time denotes the probability density distribution corresponding to the gene expression state predicted by the neural ordinary differential equation learning framework at time denotes the probability density distribution corresponding to the gene expression state predicted by the neural ordinary differential equation learning framework at time denotes the probability density distribution corresponding to the real gene expression state x at time denotes the probability density distribution corresponding to the gene expression state predicted by the neural ordinary differential equation learning framework at time denotes the probability density distribution corresponding to the gene expression state denotes the loss corresponding to the real gene expression state x at time denotes the loss corresponding to the real gene expression state x at time denotes the loss corresponding to the real gene expression state x at time denotes the loss corresponding to the real gene expression state x at time

[0042] To enhance robustness, reconstruction errors at adjacent time points and long-term time points are calculated to capture short-term and long-term consistency of trajectory-driven probability density distribution. By backpropagation, all parameters are jointly optimized to achieve the best approximation of real developmental trajectories.

[0043] In one possible design, the method of adjusting parameters of the single-cell developmental state change model based on the neural ordinary differential equation learning framework to construct a population-cell developmental state change model comprises:

[0044] Based on the drift term and the diffusion term of the single-cell developmental state change model, a velocity field of population state change is obtained, a proliferation field is introduced, the velocity field and the proliferation field are input into the neural ordinary differential equation learning framework, the velocity field and the proliferation field of population state change are adjusted based on the loss function of the neural ordinary differential equation learning framework, and the velocity field and the proliferation field of population state change after feature extraction are used to construct the population-cell developmental state change model.

[0045] In one example, based on the drift term and the diffusion term of the single-cell developmental state change model, a velocity field is obtained, a proliferation field The above field is input as a neural ordinary differential equation learning framework, and the speed field and the proliferation field of the population state change are optimized based on the neural ordinary differential equation learning framework loss function to obtain the optimized speed field and the optimized proliferation field ; based on the optimized speed field and the proliferation field of the population state change, a population cell development state change model is constructed

[0046] ;

[0047] The population cell development state change model is simplified as .

[0048] The initial cell state is converted into a predicted trajectory by integration, thereby operationalizing the theoretical dynamics while maintaining interpretability through the underlying physical structure. Therefore, the neural speed field can serve as both a computational engine for high-dimensional reasoning and a differentiable representation of the single-cell development process.

[0049] In one possible design, obtaining a cell-specific gene regulatory network and a population cell gene regulatory network based on the population cell development state change model comprises:

[0050] Based on the adjusted speed field, the Jacobian matrix under the single cell state is analyzed by partial derivative decomposition;

[0051] The time-resolved cell-specific gene regulatory network is obtained according to the obtained Jacobian matrix;

[0052] According to the obtained Jacobian matrix, the population level gene regulatory network is captured by aggregating the single cell network;

[0053] Based on the cell-specific gene regulatory network and the population cell gene regulatory network, the timing of cell lineage differentiation and the molecular decision mechanism of cell differentiation are predicted.

[0054] In one example, the Jacobian matrix of the cell is analyzed by partial derivative decomposition:

[0055] ,

[0056] wherein, quantifies the instantaneous regulatory effect of the gene on the velocity component of the cell , is the velocity component of the cell , is the velocity component of the cell , is the velocity component of the cell , is the velocity component of the cell represents the Jacobian matrix of the cell , represents the optimized cell 's effective velocity field;

[0057] By , the time-resolved cell-specific gene regulatory network (GRN) is obtained, wherein, represents the specific gene regulatory network of cell k, represents the absolute value of the interaction strength, represents the activation / inhibition of the sign, is a sign function, represents the Hadamard product, and the final directed GRN is given by the Hadamard product ;

[0058] By aggregating single-cell networks to capture population-level GRN, wherein, represents the population-level GRN, and N represents the total number of cells, the directionality of cell specificity is retained in

[0059] The second aspect of the application provides a cell development process dynamic modeling device based on time-series single-cell transcriptome data, comprising a storage device, a processor, and a computer program stored in the storage device and executable on the processor, characterized in that the processor executes the computer program to realize the steps of the cell development process dynamic modeling method based on time-series single-cell transcriptome data according to the first aspect of the application.

[0060] The third aspect of the application provides a computer storage medium comprising a computer program, which, when executed by a processor, realizes the steps of the cell development process dynamic modeling method based on time-series single-cell transcriptome data according to the first aspect of the application.

[0061] Advantages:

[0062] The application provides a cell development process dynamic modeling method based on time-series single-cell transcriptome data, which solves the problems of incomplete existing cell development process modeling mechanism, insufficient noise processing, and lack of energy principle, and has significant effects compared with the prior art:

[0063] (1) significantly improves the accuracy and biological rationality of cell development trajectory reconstruction;

[0064] The application models the evolution of cell state as a controlled stochastic differential equation, while introducing three types of physical field functions, namely velocity field, noise field and proliferation field, to fully characterize the three core mechanisms of gene regulation driving force, expression fluctuation and cell proliferation effect. Compared with traditional methods (such as Monocle, RNA velocity) that rely only on geometric structure or local velocity estimation, the present scheme can more realistically reflect the behavior of cells in complex dynamic processes such as bifurcation, convergence and acceleration, significantly improving the biological consistency and prediction ability of trajectory inference.

[0065] (2) Efficiently differentiable solution of high-dimensional nonlinear systems is achieved;

[0066] In view of the problem that the dimension of gene expression space is high and the traditional PDE (partial differential equation) is difficult to solve, the application adopts a neural ODE (ordinary differential equation) framework for meshless modeling, parameterizes the physical field function by using a deep neural network, and realizes end-to-end training through automatic differentiation. This method avoids the "curse of dimensionality" of constructing a grid or performing density propagation in a high-dimensional space, greatly improves the computational efficiency and model scalability, and is suitable for real single-cell data with thousands to tens of thousands of gene dimensions.

[0067] (3) The development path is given a clear physical meaning and energy optimality;

[0068] The application introduces the non-equilibrium optimal transport theory into single-cell trajectory modeling, and defines an energy functional containing kinetic energy, noise dissipation and proliferation cost as an optimization objective. This energy minimization principle ensures that the inferred development path is the most energy-efficient evolution path, not only enhancing the stability of the model, but also making the results more consistent with the universal law of biological systems pursuing energy efficiency, improving the interpretability and robustness of the modeling results.

[0069] (4) Good modularity and scalability, suitable for various development scenarios;

[0070] The modules of the application, such as the dimension reduction method, the density estimation method and the neural network structure, support multiple alternative implementations, and users can flexibly choose according to the data characteristics. In addition, the model naturally supports joint modeling of multiple batches and multiple condition data, and can be used to compare the differences in development trajectories under different genotypes, drug treatments or disease states, and has wide application prospects. BRIEF DESCRIPTION OF DRAWINGS

[0071] Figure 1 A cell development process dynamic modeling method flowchart based on time-series single-cell transcriptome data provided by the embodiments of the application;

[0072] Figure 2A flowchart for optimizing a single cell development state change model parameter based on a neural ordinary differential equation learning framework is provided for the embodiments of the present application.

[0073] Figure 3 A schematic diagram of a cell development process dynamic modeling device based on time series single cell transcriptome data is provided for the embodiments of the present application.

[0074] Figure 4 A schematic diagram of a computer storage medium is provided for the embodiments of the present application. DETAILED DESCRIPTION

[0075] The technical solutions in the embodiments of the present application will be described clearly and completely below with reference to the drawings in the embodiments of the present application. Obviously, the described embodiments are only part of the embodiments of the present application, rather than 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 work belong to the scope of protection of the present application. It should be noted that the embodiments in the present application and the features in the embodiments can be combined with each other without conflict.

[0076] The existing population evolution modeling method based on optimal transport regards the cell population at different time points as a probability distribution, and uses the Wasserstein distance (also known as the earth mover's distance (EMD)) to find the minimum cost cell state transition scheme, so as to infer the development path. Its advantage lies in modeling from the perspective of population evolution, which naturally adapts to the destructive sampling characteristics of single cell data. Although various cell trajectory inference methods have been proposed, there are still the following key technical problems in practical application, which limit their accuracy and interpretability in complex development systems:

[0077] Traditional optimal transport (OT) methods usually assume that the system is in thermodynamic equilibrium, has no source and no sink, and cannot describe the mass change caused by cell proliferation, and lack explicit modeling of the internal noise structure;

[0078] Lack of unified dynamic mechanism modeling: existing methods are mostly based on geometric structure or local velocity estimation, and fail to regard cell development as a complete stochastic dynamic system driven by gene regulation, disturbed by noise and regulated by proliferation, resulting in insufficient explanation of nonlinear phenomena such as bifurcation and convergence.

[0079] Ignoring the energy optimization principle in population evolution: cell state transition is not an arbitrary path, but should follow a certain "most economical" evolution strategy. Existing methods lack a clear energy minimization loss function as a trajectory optimization criterion, resulting in a deviation of the inferred path from the true biophysical law.

[0080] Model solving difficulty in high-dimensional space: The dimension of gene expression space is extremely high (usually >1000), and traditional PDE solving methods based on grid or density propagation face the "curse of dimensionality", making it difficult to achieve efficient and stable numerical simulation.

[0081] Insufficient modeling of noise and proliferation effect: Most methods assume constant noise or ignore its spatial heterogeneity, and do not explicitly model the effect of cell proliferation on population density, resulting in significant bias in trajectory estimation in stem cell expansion, tumor growth, etc.

[0082] Poor model scalability and flexibility: Existing frameworks often rely on specific assumptions (such as linear dynamics, equilibrium transport), making it difficult to adapt to diverse developmental patterns, and lack modular design, limiting the ability to integrate with other prior knowledge.

[0083] In summary, the existing technology has not yet provided a unified modeling framework that can reflect the intrinsic physical mechanism of cell development, efficiently solve in high-dimensional space, and have good interpretability.

[0084] Therefore, the embodiments of the present application provide a dynamic modeling method of cell development process based on time series single cell transcriptome data to solve the above problems. In the following, the embodiments of the present application will be described in detail with reference to the accompanying drawings Figures 1-2 The scheme of the embodiments of the present application will be described in detail.

[0085] Reference Figure 1 The embodiments of the present application provide a dynamic modeling method of cell development process based on time series single cell transcriptome data, comprising:

[0086] Constructing a single cell development state change model based on time series single cell transcriptome data;

[0087] Constructing a neural ordinary differential equation learning framework;

[0088] Adjusting the single cell development state change model parameters based on the neural ordinary differential equation learning framework to construct a population cell development state change model;

[0089] Obtaining cell-specific gene regulation network and population cell gene regulation network based on the population cell development state change model to predict the timing of cell lineage differentiation and the molecular decision mechanism of cell differentiation.

[0090] Time-series single-cell transcriptome data: Refers to a dataset obtained by measuring gene expression in individual cells at multiple consecutive time points using single-cell RNA sequencing technology. These data capture the dynamic changes in gene expression during development, differentiation, or response to external stimuli, and can be used to infer trajectories of cell state transitions, identify key time-point events (such as cell fate decisions), and typically include a gene expression matrix, time labels, and cell metadata. Time-series designs allow for the analysis of continuous evolution of cell states rather than static snapshots, revealing molecular dynamics during development or disease progression.

[0091] Single-cell developmental state change model: Refers to a mathematical model or computational framework that describes the evolution of an individual cell's internal state (such as gene expression profile, metabolic state, or phenotype) over time during development. This model is typically based on differential equations, stochastic processes, or machine learning methods (such as hidden Markov models) that simulate the transition of cells from a starting state to a terminal state. The goal is to quantify the rate, direction, and driving factors of cell state changes, for example, by fitting time-series single-cell transcriptome data to infer cell-specific kinetic parameters.

[0092] Neural ordinary differential equation learning framework: Refers to a machine learning framework that combines neural networks and ordinary differential equations (ODEs), where neural networks are used to parameterize ODEs functions to learn the continuous-time evolution of dynamic systems. In this framework, the weights of the neural network are optimized through backpropagation or adjoint methods to fit observed data (such as time-series single-cell data) to infer the differential equation representation of the system state (such as cell gene expression). This framework is suitable for modeling high-dimensional, nonlinear biological processes such as cell developmental trajectories and can handle irregularly time-sampled data.

[0093] Population cell developmental state change model: Refers to a computational model that extends from the single-cell model to describe the statistical evolution of the overall state of an entire cell population (such as a tissue or culture) during development. The goal is to simulate population-level behavior such as lineage differentiation proportions, population synchrony, or instability, and integrate single-cell data to infer population parameters such as differentiation rates or intercellular variability.

[0094] Cell-specific gene regulatory network: Refers to a gene regulatory network (GRN) constructed for an individual cell, which represents the regulatory relationships between genes within that cell in a graph structure (such as the activation or inhibition interactions between transcription factors and target genes). Network nodes represent genes, and edges represent regulatory strengths, which are typically inferred from single-cell data (such as based on co-expression, information theory, or machine learning methods). This network emphasizes that different cells have different gene regulatory relationships and can be used to explain the specific behavior of individual cells during development or disease.

[0095] Population cell gene regulatory network: refers to the gene regulatory network constructed at the cell population level, which infers the consensus or average regulatory relationship by integrating multiple single-cell data (such as clustering or average expression profile). This network captures the universal regulatory patterns in the population, ignoring the heterogeneity of cell regulatory relationships.

[0096] Timing of lineage differentiation and molecular decision-making mechanism of cell differentiation: refers to the time point (timing) at which the lineage differentiation event (such as the differentiation of pluripotent stem cells into specific cell types) occurs during cell development, and the molecular cascade events (decision-making mechanism) that control the differentiation. This includes identifying key signaling pathways, transcription factor dynamics, epigenetic modifications, and gene expression switches, which together determine cell fate selection. By analyzing time-series single-cell data or perturbation experiments, the differentiation timing (such as through pseudo-time analysis) can be quantified and the decision-making rules (such as bistable networks or feedback loops) can be revealed, thereby understanding developmental errors or disease mechanisms.

[0097] Specifically, the present scheme first receives "time-series single-cell transcriptome data input" in the data preparation and preprocessing stage, and performs necessary preprocessing and dimensionality reduction. In the preferred embodiment of the present application, the method of autoencoder is used for dimensionality reduction processing to extract key features and reduce data dimension. Subsequently, kernel density estimation (KDE) method is used for population probability density estimation to construct the distribution model of cell state at different time points, providing reliable data support for subsequent construction of random dynamics model integrating driving, noise and proliferation.

[0098] Next, enter the core modeling stage. In order to solve the above random dynamics model, a meshless high-dimensional solution framework based on NeuralODE is designed. This framework uses deep neural network learning algorithm to simulate the dynamic changes in cell development process through Neural ODE. In order to ensure the accuracy and biophysical significance of the model, a set of loss functions are defined as optimization objectives. Through this principle, the learning process of neural network can be effectively constrained to ensure that the output results conform to the biological laws.

[0099] Finally, after obtaining the modeling results of each part, further analysis and interpretation are carried out to realize the interpretable dynamic trajectory modeling and bifurcation prediction. This stage not only includes the visualization display of cell development trajectory, but also involves the prediction of key gene dynamics and lineage bifurcation points, thereby providing a powerful tool for in-depth understanding of cell development mechanisms.

[0100] The present application models the cell development as a stochastic dynamic system controlled by intrinsic regulatory forces and extrinsic noise by introducing stochastic differential equation, thus naturally explaining the bifurcation behavior on the development path - when the velocity field appears unstable equilibrium point or saddle point, the system will evolve along different stable manifold, corresponding to different fate determination. The present scheme is designed based on the following principles:

[0101] Cell development is the result of the interaction of deterministic process driven by gene regulatory network and random noise, and stochastic differential equation can reflect both characteristics;

[0102] The evolution of the probability density of the population cells can effectively avoid the technical limitations of single-cell tracking, and the Fokker-Planck equation is used to realize the macroscopic description of the population behavior;

[0103] When the driving term of the gene regulatory network meets certain conditions (such as the existence of multi-stable state), the probability density distribution will show bifurcation structure, which corresponds to the key node of cell fate differentiation.

[0104] By introducing the noise term and the proliferation term, the present scheme solves the problem that the existing method cannot quantify the influence of random disturbance on cell fate and ignores the contribution of proliferation to population evolution, and at the same time, through the modeling of bifurcation phenomenon, the precise capture of the cell fate differentiation process is realized.

[0105] In one embodiment, based on the time-series single-cell transcriptome data, a single-cell development state change model is obtained by modeling the cell development process using a stochastic differential equation:

[0106] ,

[0107] wherein, represents the M-dimensional gene expression vector of the cell at time t, the drift term represents the driving term of the gene regulatory network, which describes the deterministic influence of gene interaction on cell state, and f is the velocity field; the diffusion term , is the noise intensity matrix, which reflects the technical noise and biological random fluctuation; , represents a standard Wiener process, which simulates the random disturbance in the development process; represents an M-row S-column real number, represents an S-row real number; the noise term is designed as a variable related to the cell state, which can describe the source of cell heterogeneity, which follows the known biological principle, for example, the noise is high in the metastable state before bifurcation occurs.

[0108] In one embodiment, the neural stochastic differential equation learning framework is constructed, comprising:

[0109] constructing an evolution simplified equation of the probability density distribution of the population cells, for representing a population cell path;

[0110] defining a total energy cost discrete equation of the population cells from an initial distribution to a terminal distribution based on the probability density of the population cells;

[0111] defining the neural ordinary differential equation learning framework loss function based on the total energy cost discrete equation, for training the neural ordinary differential equation learning framework.

[0112] In an embodiment, the constructing an evolution simplified equation of the probability density distribution of the population cells comprises:

[0113] defining that the evolution of the population cell probability density function follows an extended Fokker-Planck equation, constructing an evolution initial equation of the probability density distribution of the population cells based on the probability density of the population cells;

[0114] defining an effective velocity field, and simplifying the evolution initial equation of the probability density distribution of the population cells based on the effective velocity field to obtain the evolution simplified equation of the probability density distribution of the population cells.

[0115] Specifically, since single-cell RNA sequencing has destructive sampling characteristics, it is impossible to track the same cell longitudinally, and therefore it is impossible to directly observe individual paths To this end, the present scheme replaces the tracking of individual cells by constructing an evolution equation of the probability density distribution of the population cells.

[0116] The constructing an evolution initial equation of the probability density distribution of the population cells based on the population cell probability density comprises: considering that there are asymmetric division and proliferation behaviors of cells in the development process, introducing a source term , representing the influence of the local proliferation rate on the population density, performing population probability density estimation by kernel density estimation, and the evolution initial equation of the probability density distribution of the population cells is:

[0117] ,

[0118] wherein, is the population probability density, is a diffusion coefficient matrix, , is an M-row and M-column real number, represents a divergence; represents a cell division hotspot, represents an apoptosis region; from a biological point of view, strictly explains the local change of cell density caused by the cell cycle; the model uniformly depicts three biological mechanisms: gene regulation drive (through the velocity field ); molecular noise influence (through diffusion coefficient matrix ); cell proliferation effect (through source term );

[0119] defining an effective velocity field , , the effective velocity field incorporates the combined influence of deterministic forces and noise-induced drift, where the logarithmic gradient term measures the density-driven "statistical pressure" induced by randomness;

[0120] Based on the effective velocity field, the initial equation of evolution of the probability density distribution of the population of cells is simplified to obtain a simplified equation of evolution of the probability density distribution of the population of cells, .

[0121] In an embodiment, the defined discrete equation of total energy cost of the population of cells from the initial distribution to the terminal distribution includes:

[0122] The definition of the cell development trajectory follows the principle of energy optimization, and the total energy cost equation of the population of cells from the initial distribution to the terminal distribution is defined based on the set of all possible continuous evolution paths of the population of cells from the initial distribution to the terminal distribution;

[0123] Randomly sampling a set of particle points in a high-dimensional space at an initial time, and simulating the evolution of the probability density according to the flow of the particles, the continuous total energy cost equation is converted into the total energy cost discrete equation.

[0124] Specifically, given a time-series single-cell data set , the probability density of the cell population at each time point , the present application defines that the cell development trajectory should follow the principle of energy optimization, that is, the system minimizes the weighted sum of total kinetic energy, noise dissipation and proliferation cost during the evolution process.

[0125] The non-equilibrium optimal transport theory framework is introduced, and the total energy cost equation of the population of cells from the initial distribution to the terminal distribution is defined based on the set of all possible continuous evolution paths of the population of cells from the initial distribution to the terminal distribution, and the total energy cost equation of the population of cells from the initial distribution to the terminal distribution is:

[0126] ,

[0127] wherein, is the initial distribution, is the terminal distribution, is the energy cost, T is the time corresponding to the terminal distribution, is an M-row real number, is the kinetic energy cost required for driving, a quality change cost of proliferation activity, a weight coefficient of setting;

[0128] The total energy cost equation of population cells from the initial distribution to the final state distribution comprehensively considers the kinetic energy cost required to drive , the quality change cost of proliferation activity and the weight coefficient balance the contributions of gene expression dynamics and proliferation; obviously, this is the non-equilibrium optimal transport cost defined in the high-dimensional space (single-cell state usually contains thousands of genes), and its direct calculation is not feasible. Therefore, it is intended to overcome this problem by trajectory-based Monte Carlo integration, by randomly sampling a set of particle points in the high-dimensional space at the initial time , and using the flow of particles to simulate the evolution of the probability density, converting the continuous energy cost into a discrete form:

[0129] The total energy cost discrete equation is:

[0130] ;

[0131] Embed this solver in a deep neural network framework, parameterize , and by neural networks; combine the Monte Carlo integration approximation of the energy functional with the Fokker-Planck constraint to realize the end-to-end optimization of energy.

[0132] Cell development, as a natural process, follows the principle of energy consumption minimization. The non-equilibrium optimal transport model describes this feature through the energy functional, and allows the existence of source terms (proliferation) and non-conservative forces (non-gradient velocity field), which is more consistent with the real development biology scenario;

[0133] Monte Carlo estimation avoids grid division in high-dimensional space by random sampling, significantly reducing computational complexity;

[0134] Deep neural networks use the nonlinear fitting ability of neural networks to achieve flexible parameterization of high-dimensional velocity fields.

[0135] This scheme solves the problem of solving the optimal transport model in high-dimensional space, and at the same time, through the introduction of the energy optimization target, makes the cell trajectory inference more consistent with the internal rules of the biological process, and improves the physical consistency of the trajectory inference.

[0136] In one embodiment, the loss function of the neural ordinary differential equation learning framework is defined based on the total energy cost discrete equation, comprising:

[0137] The loss function of the neural ordinary differential equation learning framework is defined based on a reconstruction loss of a time-series probability density distribution based on a total energy cost of population cells from an initial distribution to a final distribution;

[0138] The reconstruction loss of the time-series probability density distribution comprises adjacent time reconstruction errors and long-term time reconstruction errors.

[0139] Specifically, to solve the model solving problem in a high-dimensional gene space, the neural ordinary differential equation learning framework is constructed, and a deep neural network is used to parameterize a velocity field and a proliferation field . The velocity field , the proliferation field , the drift term , and the diffusion term are obtained by modeling; the velocity field , the proliferation field , the drift term , and the diffusion term are taken as inputs, and the initial cell state is converted into a predicted trajectory by integration, so as to operationalize the theoretical dynamics while maintaining interpretability through the physical basis structure of

[0140] . Therefore, the neural velocity field can serve as both a computational engine for high-dimensional reasoning and a differentiable representation of single-cell developmental processes.

[0141] ,

[0142] wherein the reconstruction loss of the time-series probability density distribution is

[0143] ,

[0144] wherein denotes a probability density distribution corresponding to a real gene expression state x at time , denotes a probability density distribution corresponding to a gene expression state predicted by the neural ordinary differential equation learning framework at time , denotes a probability density distribution corresponding to a real gene expression state x at time , denotes a probability density distribution corresponding to a gene expression state predicted by the neural ordinary differential equation learning framework at time ; denotes a loss corresponding to a real gene expression state x at time , Indicates time The loss corresponding to the true gene expression state x;

[0145] To enhance robustness, calculate adjacent time points. and long-term time point The reconstruction error is minimized, and the short-term and long-term consistency of the trajectory-driven probability density distribution is captured. All parameters are jointly optimized through backpropagation to achieve the best approximation of the true developmental trajectory. End-to-end differentiable modeling is achieved using Neural ODE, avoiding the grid limitations of traditional PDE solutions; physical constraints and energy regularization are combined to prevent overfitting and improve generalization ability.

[0146] In one implementation, the step of optimizing the parameters of a single-cell developmental state change model based on a neural constant differential equation learning framework to construct a population cell developmental state change model includes:

[0147] The velocity field of population state change is obtained based on the drift and diffusion terms of the single cell developmental state change model. A proliferation field is introduced, and the velocity field and the proliferation field are input into the neural network constant differential equation learning framework. The velocity field and the proliferation field of population state change are characterized based on the loss function of the neural network constant differential equation learning framework. Based on the characteristic velocity field and proliferation field of population state change, a population cell developmental state change model is constructed.

[0148] Specifically, based on the drift term of the model of changes in the developmental state of a single cell. and diffusion terms Obtain the velocity field of the group state change Introducing a breeding field The above velocity field and breeding farm As input, by passing the initial cell state x through Integrals are transformed into predicted trajectories, thus operationalizing theoretical dynamics, while simultaneously... The physical underlying structure remains interpretable. Therefore, the neural velocity field can serve as both a computational engine for high-dimensional reasoning and a differentiable representation of single-cell development processes.

[0149] refer to Figure 2, the training process of the meshless high-dimensional solving framework based on the Neural ODE includes: obtaining time series observation values of population density by using kernel density estimation, obtaining time series prediction values of population density based on the Neural ODE learning framework and the time series observation values of population density, performing developmental trajectory prediction based on the time series observation values of population density and the time series prediction values of population density, obtaining an optimal model describing the cell development process, and realizing interpretable dynamic trajectory modeling and bifurcation prediction; constructing a loss function of the Neural ODE learning framework based on reconstruction loss and energy loss, and optimizing 、 and obtaining 、 and , constructing a population cell development state change model:

[0150] ;

[0151] The population cell development state change model is simplified as: .

[0152] In an embodiment, a method for obtaining a cell-specific gene regulatory network and a population cell gene regulatory network based on the population cell development state change model includes:

[0153] Based on the adjusted velocity field, the Jacobian matrix in the state of a single cell is analyzed by partial derivative decomposition;

[0154] The time-resolved cell-specific gene regulatory network is obtained according to the obtained Jacobian matrix;

[0155] The population-level gene regulatory network is captured by aggregating the single-cell network according to the obtained Jacobian matrix;

[0156] The cell lineage differentiation occurrence time and the molecular decision mechanism of cell differentiation are predicted based on the cell-specific gene regulatory network and the population cell gene regulatory network.

[0157] Specifically, the built population cell development state change model can not only generate continuous cell trajectories, but also identify local stability, fixed points and bifurcation structures by analyzing the Jacobian matrix of the velocity field, so as to predict the occurrence time and molecular decision mechanism of lineage differentiation. For example, when the velocity field appears a saddle point or a Hopf bifurcation, the system will enter a multi-stable evolution stage, corresponding to a selection window of cell fate. Once the trained velocity field model , it can be used to reconstruct the population-level and single-cell-resolved gene regulatory network (GRN).

[0158] The Jacobian matrix in the state of a single cell is analyzed by partial derivative decomposition:

[0159] ,

[0160] where, quantifies the instantaneous regulatory effect of gene on the velocity component of cell , is the first component of the velocity field, corresponding to the rate of change of gene expression; represents the Jacobian matrix of cell , represents the effective velocity field of the optimized cell ; By , the time-resolved cell-specific gene regulatory network is obtained, where,

[0161] represents the specific gene regulatory network of cell k, the absolute value of represents the interaction strength, the sign of represents activation / inhibition, is the sign function, represents the Hadamard product, and the final directed GRN is given by the Hadamard product ;

[0162] By aggregating single-cell networks to capture population-level GRNs, where, represents the population-level GRN, and N represents the total number of cells, the cell-specific directionality is retained.

[0163] Traditional gene regulatory network reconstruction methods (such as SCENIC, GENIE3) rely on static correlations or co-expression patterns, making it difficult to capture time-varying regulatory logic during development. The principle of this scheme is:

[0164] The gradient of the velocity field directly reflects the dependence of gene expression rate on other genes, revealing the instantaneous regulatory direction and strength between genes, and is more dynamic.

[0165] The velocity field is learned from the entire developmental trajectory, and the reconstructed GRN naturally has spatiotemporal resolution.

[0166] Solves the problem of missing cell heterogeneity and dynamics in regulatory network reconstruction, and realizes high-resolution cell-specific time-varying network construction through velocity field decomposition, providing a precise tool for analyzing underlying regulatory mechanisms.

[0167] The embodiment of the present application also provides a cell development process dynamic modeling device based on time-series single-cell transcriptome data, referring to Figure 3 , comprising a storage device, a processor, and a computer program stored in the storage device and executable on the processor, wherein the processor executes the computer program to implement the steps of the cell development process dynamic modeling method based on time-series single-cell transcriptome data provided by the embodiment of the present application. Wherein, the same or corresponding contents in the scheme shown in Figures 1-2 may refer to the detailed description of the above Figures 1-2 , which will not be repeated here.

[0168] The embodiment of the present application also provides a computer storage medium, referring to Figure 4 , comprising a computer program, which is executed by a processor to implement the steps of the cell development process dynamic modeling method based on time-series single-cell transcriptome data provided by the embodiment of the present application. Wherein, the same or corresponding contents in the scheme shown in Figures 1-2 may refer to the detailed description of the above Figures 1-2 , which will not be repeated here.

[0169] Embodiment:

[0170] The cell development process dynamic modeling method based on time-series single-cell transcriptome data provided by the present application has high flexibility and modularization characteristics. On the premise of ensuring the consistency of the overall modeling logic, various components can be implemented by using multiple technical paths. In order to fully expand the technical boundaries of the present application, improve the applicability, and enhance the stability of the patent claims, the following provides several alternative technical implementation schemes for the key links. These alternative schemes can all achieve the basic purpose of the present application, that is, modeling the cell development process from time-series single-cell transcriptome data and analyzing biological driving factors such as velocity field, proliferation field, and noise field.

[0171] In the preferred embodiment of the present application, the original high-dimensional gene expression data (usually with dimensions of thousands to tens of thousands of genes) is first processed by dimension reduction to alleviate the curse of dimensionality, improve the calculation efficiency, and reduce noise interference. Although principal component analysis as a typical linear dimension reduction method is preferentially adopted, the person skilled in the art can select any one or more of the following alternative schemes according to the specific application scenario:

[0172] Alternative schemes of data dimension reduction method:

[0173] Nonlinear manifold learning method: including t-distributed stochastic neighbor embedding, uniform manifold approximation and projection, which can be used to preserve the local cell neighborhood structure, and is particularly suitable for development systems with complex bifurcation or circular trajectory;

[0174] Autoencoder and its variants: learn nonlinear low-dimensional latent space representation through deep neural network, support end-to-end training, and can be jointly optimized with subsequent Neural ODE module;

[0175] Diffusion map: construct Markov transition graph based on cell-cell similarity, extract low-dimensional coordinates aligned with developmental kinetics, and naturally adapt to the geometric structure of the Fokker-Planck equation;

[0176] Latent variable model: such as scVI, DESTINY, and other generative models based on variational inference, which can model technical noise and batch effects while reducing dimensionality.

[0177] The low-dimensional latent space obtained by the above dimensionality reduction methods can be used as the input space for subsequent SDE modeling and Neural ODE solving. At this time, the velocity field and the proliferation field are defined on the latent space, and their physical meaning can be interpreted by mapping them back to the original gene space through the decoder.

[0178] Alternative solutions for probability density calculation methods:

[0179] This application relies on the empirical probability density of the cell population at each time point as a supervisory signal or physical constraint condition. In the preferred embodiment, the kernel density estimation method is used to construct the continuous density distribution. However, according to the data sparsity, dimensionality, and computational resource constraints, the following alternatives can also be applicable:

[0180] k-Nearest Neighbor-based Density Estimation: Estimate local density using the inverse relationship between the density of each cell in its neighborhood, suitable for irregular distribution and does not require assumptions about kernel function form;

[0181] Graph-based density estimation: define graph Laplacian on cell similarity graph (such as k-nearest neighbor graph), simulate density evolution through heat kernel diffusion, and naturally compatible with optimal transport framework;

[0182] Mixed model method: use Gaussian mixture model or Dirichlet process mixture model to cluster cell states and fit global density, suitable for data with obvious subpopulation structure;

[0183] Deep generative model estimation: use normalized flow, generative adversarial network, or diffusion model to learn complex density distribution, support efficient sampling and density evaluation;

[0184] Histogram and gridding method: divide the grid in the low-dimensional latent space and count the cell frequency, combine with smoothing filtering (such as Gaussian filtering) to obtain discrete density representation, suitable for visualization and rapid verification.

[0185] It is worth noting that the results output by all density estimation methods can be used to construct data fitting terms in the loss function or as boundary / initial conditions for the Fokker-Planck equation. The choice of different methods does not affect the overall architecture of this application, but only involves the trade-off between accuracy, stability, and computational overhead.

[0186] Alternative solutions for neural network structure:

[0187] This application uses deep neural network to parameterize the drift term, diffusion term, and proliferation field. In the preferred embodiment, a multi-layer perceptron is used as the basic structure. However, depending on the characteristics of the input data, spatio-temporal dependencies, and the availability of prior knowledge, the following neural network architectures can also be used as effective alternatives:

[0188] Residual network: Introducing skip connections to alleviate the problem of gradient vanishing, improving the training stability of deep networks, suitable for learning complex nonlinear fields;

[0189] Graph neural network: If an adjacency graph has been constructed between cells (such as based on expression similarity), local interactions can be modeled using graph convolutional networks or graph attention networks to reflect the influence of the cell microenvironment;

[0190] Conditional neural network structure: Input time as a conditional variable (such as through time embedding or FiLM mechanism) to explicitly model the time-varying characteristics of the field;

[0191] Shared weight structure: For physical quantities such as velocity field and proliferation field that have potential correlations, a shared bottom feature extraction network can be designed to improve parameter utilization efficiency;

[0192] Physical information neural network: Directly calculate the Fokker-Planck equation residual at the network output and add it to the loss function as a soft constraint to enhance physical consistency;

[0193] Transformer architecture: When dealing with large-scale cell collections, use self-attention mechanisms to capture long-range cell dependencies, suitable for modeling the co-evolution of tissue scales.

[0194] In addition, the activation function of the neural network can also be flexibly selected, including nonlinear functions such as ReLU, Swish, and Tanh, or custom activation forms designed according to biological priors (such as sigmoid response functions to simulate gene regulation threshold effects). These alternative network structures can all achieve parameterized modeling of field variables, and the choice should consider the temporal characteristics of the data, cell-cell relationships, and computational resources, with the core being to ensure that the network can accurately fit the mapping relationship between field variables and gene expression, time.

[0195] Alternative solutions for optimal transport solving method:

[0196] In the process of solving the energy functional of the non-equilibrium optimal transport model, in addition to the meshless high-dimensional solver based on Monte Carlo estimation, the entropy regularization optimal transport method can also be used: by introducing an entropy regularization term into the energy functional, the original problem is transformed into a form that can be efficiently solved by the Sinkhorn algorithm, reducing the computational complexity, especially suitable for high-dimensional space scenarios with high requirements for solving speed. This method sacrifices a small amount of accuracy for the improvement of computational efficiency, and has obvious advantages in large-scale single-cell data processing.

[0197] Alternative of loss function:

[0198] In defining the loss function of the neural network, in addition to using the combination of mean square error and KL divergence, the Wasserstein distance can be used to replace the KL divergence as the probability density distribution matching loss. The Wasserstein distance is more sensitive to local differences in distribution, and can more robustly measure the difference between two probability distributions, especially suitable for scenarios where the data distribution has multi-peak or long-tail characteristics, reducing the model bias caused by outliers, and improving the stability of trajectory estimation.

[0199] Although the present application is described herein with reference to particular embodiments, it will be understood that these examples are merely illustrative of the principles and applications of the present application. It will therefore be understood that numerous modifications can be made to the illustrative embodiments, and that other arrangements can be devised without departing from the spirit and scope of the application as defined by the appended claims. It will be understood that different dependent claims and features described herein can be combined with each other in ways other than those described in the original claims. It will also be understood that features described in connection with individual embodiments can be used in other described embodiments.

Claims

1. A method for dynamic modeling of cell development processes based on time-series single-cell transcriptome data, characterized in that, include: A model of single-cell developmental state changes was constructed based on time-series single-cell transcriptome data; Constructing a learning framework for divine ordinary differential equations; The parameters of the single-cell developmental state change model are adjusted based on the aforementioned neural constant differential equation learning framework to construct a population cell developmental state change model. Based on the population cell development state change model, cell-specific gene regulatory networks and population cell gene regulatory networks are obtained to predict the timing of cell lineage differentiation and the molecular decision-making mechanism of cell differentiation.

2. The method for dynamic modeling of cell development processes based on time-series single-cell transcriptome data according to claim 1, characterized in that, The model for single-cell developmental state changes based on time-series single-cell transcriptome data includes: Based on time-series single-cell transcriptome data, a model of single-cell developmental state changes was obtained by using stochastic differential equations to model the cell development process.

3. The method for dynamic modeling of cell development processes based on time-series single-cell transcriptome data according to claim 1, characterized in that, The framework for constructing a learning framework for divine ordinary differential equations includes: A simplified equation for the evolution of the probability density distribution of the population cells is constructed to represent the population cell paths; Based on the probability density of the population cells, a discrete equation for the total energy cost of the population cells from the initial distribution to the final state distribution is defined. The loss function of the learning framework of the divine constant differential equation is defined based on the discrete equation of the total energy cost.

4. The method for dynamic modeling of cell development processes based on time-series single-cell transcriptome data according to claim 3, characterized in that, The simplified equation for the evolution of the probability density distribution of the constructed population of cells includes: The evolution of the population cell probability density function is defined to follow the extended Falk-Planck equation, and the initial equation for the evolution of the population cell probability density distribution is constructed based on the population cell probability density. Define an effective velocity field, and simplify the initial equation for the evolution of the probability density distribution of the population cells based on the effective velocity field to obtain the simplified equation for the evolution of the probability density distribution of the population cells.

5. The method for dynamic modeling of cell development processes based on time-series single-cell transcriptome data according to claim 3, characterized in that, The discrete equation defining the total energy cost of the population cells from the initial distribution to the final distribution includes: The cell development trajectory is defined according to the principle of energy optimization. The total energy cost equation for the distribution of cells from the initial distribution to the final distribution is defined based on the set of all possible continuous evolutionary paths of the population cells from the initial distribution to the final distribution. The initial high-dimensional space particle point set is randomly sampled, and the evolution of the probability density is simulated based on the particle flow, transforming the continuous total energy cost equation into the discrete total energy cost equation.

6. The method for dynamic modeling of cell development processes based on time-series single-cell transcriptome data according to claim 3, characterized in that, The loss function defined based on the discrete equation of the total energy cost in the learning framework of the divine ordinary differential equation includes: Based on the total energy cost of the population cells from the initial distribution to the final distribution, the reconstruction loss of the temporal probability density distribution is introduced, and the loss function of the aforementioned neural ordinary differential equation learning framework is defined. The reconstruction loss of the time-series probability density distribution includes reconstruction errors at adjacent times and reconstruction errors at longer times.

7. The method for dynamic modeling of cell development processes based on time-series single-cell transcriptome data according to claim 6, characterized in that, The step of adjusting the parameters of the single-cell developmental state change model based on the neural constant differential equation learning framework to construct a population cell developmental state change model includes: The velocity field of population state change is obtained based on the drift and diffusion terms of the single cell developmental state change model. A proliferation field is introduced, and the velocity field and the proliferation field are input into the neural network constant differential equation learning framework. The velocity field and the proliferation field of population state change are adjusted based on the loss function of the neural network constant differential equation learning framework. The population cell developmental state change model is constructed based on the adjusted velocity field and proliferation field.

8. The method for dynamic modeling of cell development processes based on time-series single-cell transcriptome data according to claim 7, characterized in that, The process of obtaining cell-specific gene regulatory networks and population cell gene regulatory networks based on the population cell developmental state change model includes: Based on the adjusted velocity field, the Jacobian matrix in a single cell state is analyzed by partial derivative decomposition. The time-resolved cell-specific gene regulatory network was obtained based on the Jacobian matrix. Based on the obtained Jacobian matrix, the gene regulatory network of population cells is captured by aggregating single-cell networks; Predicting the timing of cell lineage differentiation and the molecular decision-making mechanism of cell differentiation based on cell-specific gene regulatory networks and population cell gene regulatory networks.

9. A device for dynamic modeling of cell development processes based on time-series single-cell transcriptome data, comprising a storage device, a processor, and a computer program stored in the storage device and executable on the processor, characterized in that, The processor executes the computer program to implement the steps of the method for dynamic modeling of cell development processes based on time-series single-cell transcriptome data as described in any one of claims 1 to 8.

10. A computer storage medium comprising a computer program, characterized in that, When the computer program is executed by the processor, it implements the steps of the method for dynamic modeling of cell development processes based on time-series single-cell transcriptome data as described in any one of claims 1 to 8.

Images