A single-cell gene expression dynamics modeling method and system based on coupling condition flow matching

Abstract

This invention belongs to the fields of bioinformatics and machine learning, and relates to a method and system for modeling single-cell gene expression dynamics based on coupled conditional flow matching. The method first maps cell expression representations at each time point to a shared latent space; then, it calculates entropy-regularized optimal transport soft coupling between adjacent time points, and constructs a weighted flow matching supervision target based on this soft coupling to learn a time-varying velocity field; next, it learns the forward and backward velocity fields and uses their consistency constraints to improve the temporal coherence of soft coupling and local transitions; and it introduces distribution-level alignment and latent space dynamic regularization to globally anchor the long-range generation process, reducing distribution drift caused by accumulated integral errors. This invention can achieve interpolation and extrapolation prediction of single-cell distributions even without intermediate time points, and improves the temporal consistency of generated trajectories and the overall distribution fitting accuracy.

Description

Technical Field

[0001] This invention belongs to the interdisciplinary field of bioinformatics and machine learning, specifically involving a single-cell gene expression dynamics modeling method and system based on coupled conditional flow matching. It is applicable to the analysis and processing of single-cell RNA sequencing (scRNA-seq) time series data, realizing temporal interpolation, extrapolation prediction and cell development trajectory reconstruction of single-cell gene expression distribution. It belongs to the application of generative modeling and optimal transport technology in single-cell omics. Background Technology

[0002] Single-cell RNA sequencing technology overcomes the limitations of traditional batch sequencing, enabling the analysis of gene expression profiles at single-cell resolution. This reveals cellular heterogeneity and the dynamic transitions in cell states, making it a core tool for studying biological processes such as embryonic development, cell differentiation, and disease progression. However, in practical research, due to the destructive nature of single-cell sequencing, time-resolved scRNA-seq data can only be obtained at a few discrete time points. Furthermore, the cell sets at each time point are unpaired snapshot populations, resulting in a large number of unmeasured intermediate time points, making it difficult to directly reflect the continuous dynamic changes in cellular gene expression. Therefore, recovering continuous gene expression dynamics from sparse single-cell snapshot data through computational modeling methods, and predicting cell distribution at unmeasured time points, is a key research goal in the field of bioinformatics.

[0003] To address the aforementioned issues, existing technologies mainly fall into two research directions: one is snapshot alignment methods based on optimal transport (OT), which achieve temporal alignment at the distribution level by calculating transport coupling between cell populations at different time points, effectively reconstructing ancestor-descendant relationships among cells. However, these methods lack explicit temporal conditional dynamics models, only supporting short-distance temporal interpolation and unable to achieve long-term extrapolation prediction, and are difficult to capture cell-level local dynamic transition features. The other is models based on continuous-time generation dynamics, such as neural ODE and flow matching models, which learn temporal conditional vector fields and use ODE integrals to generate cell distributions at arbitrary time points, supporting interpolation and extrapolation tasks. However, these methods have many technical shortcomings when dealing with unpaired sparse snapshots: Poor discernibility of local transitions: The lack of paired snapshots leads to ambiguity in cell-level dependencies between adjacent time points. Existing methods often employ coarse-grained distribution matching or unidirectional constraints, resulting in insufficient determination of local cell transfer structures. In scenarios with large time intervals or trajectory branching / fusion, multiple solutions are easily generated, leading to poor temporal coherence of the generated results. Snapshot alignment and vector field learning are loosely coupled: the optimal transmission distribution alignment and the vector field learning of generative dynamics are mostly independent modules. OT is only used as a post-processing method and does not provide accurate cell-level supervision for vector field training. This leads to noise or instability in the training supervision of vector fields. Especially when the time interval is large, the model needs to infer cell pairing relationships and fitting dynamics at the same time, which further reduces the learning accuracy. Long-range generation suffers from distribution drift: continuous-time generation models achieve long-term prediction by repeatedly integrating the velocity field learned. Even small local regression errors accumulate with the integration process, causing the predicted cell population to gradually deviate from the true snapshot distribution manifold. Local supervision can only constrain the dynamics of adjacent time points and cannot explicitly punish the distribution mismatch at distant time points, making it difficult to suppress error accumulation.

[0004] Furthermore, existing flow-matching methods in single-cell dynamics modeling often rely on unidirectional alignment signals or Euclidean cost coupling, failing to fully utilize the bidirectional consistency information of time series data and lacking regularization constraints for latent space geometry, further exacerbating the drift problem in long-range generation. Therefore, there is an urgent need for a single-cell gene expression dynamics modeling method that can achieve strong coupling between snapshot alignment and vector field learning, strengthen local dynamic transition constraints, and anchor long-range generation distributions. This method would address the shortcomings of existing techniques in sparse unpaired scRNA-seq snapshot modeling, improve the prediction accuracy of temporal interpolation and extrapolation, and simultaneously ensure the temporal coherence of the generation trajectory. Summary of the Invention

[0005] To achieve the above-mentioned objectives, the present invention adopts the following technical solution: A single-cell gene expression dynamics modeling method based on coupled conditional flow matching, the method comprising the following steps: Acquire single-cell gene expression snapshot datasets at multiple time points, and encode the single-cell gene expression vectors in the single-cell gene expression snapshot datasets at each time point into a shared latent space representation; For any two sets of shared latent space representations of adjacent time points, the pairing weight matrix is ​​obtained by calculating the entropy-regularized optimal transmission soft coupling based on the cost matrix. Based on the pairing weight matrix, bridge samples and target velocities at intermediate time points are constructed, and the time-conditional velocity field is regressed and trained using the bridge samples and target velocities (in order to learn the latent space velocity field that changes with time). The forward and backward velocity fields of the temporal conditional velocity field are jointly learned, and the local transition discernibility under unpaired supervision is improved by the consistency constraint of the forward and backward velocity fields. To address long-range generation, we introduce distribution-level alignment loss and latent space dynamic regularization loss to anchor the predicted distribution generated by the time-conditional velocity field integral (i.e., the predicted distribution of single-cell gene expression) to the observed snapshot distribution (i.e., the distribution of single-cell gene expression in the single-cell gene expression snapshot dataset) and suppress distribution drift. During the inference phase, the shared latent space representation of a single cell at the starting time point is integrated using the time-conditional velocity field after training to obtain the shared latent space representation of a single cell at the target time point, and then decoded to obtain the single-cell gene expression prediction result at the target time point.

[0006] Furthermore, the shared latent space representation is obtained by a variational autoencoder (VAE), specifically including: inputting a high-dimensional single-cell gene expression vector into the VAE encoder, encoding it into a Gaussian latent variable distribution, and sampling it using a reparameterization technique to obtain a low-dimensional latent vector. The VAE decoder then reconstructs the low-dimensional latent vector back into the gene expression space. By jointly training the VAE with reconstruction loss and KL divergence regularization, the mapping from the high-dimensional gene expression space to the low-dimensional shared latent space is realized, solving the problems of high dimensionality and sparsity in direct modeling in the gene space.

[0007] Furthermore, the cost matrix adopts a phased construction strategy: in the training warm-up phase, the Euclidean cost matrix is ​​constructed using the Euclidean distance of the latent vectors at adjacent time points to achieve initial stable training of the model; after the warm-up phase, a bidirectional fused cost matrix is ​​constructed using the one-step prediction error that fuses the forward velocity field and the backward velocity field, so that the coupled construction is more in line with the real dynamic transfer law of cells.

[0008] Furthermore, the bidirectional fusion cost matrix is ​​obtained as follows: based on the currently learned forward velocity field pairs at time points. One-step Euler prediction is performed on the latent vector at (first time point) to obtain the time point. The predicted latent vector at the second time point; based on the currently learned backward velocity field at the time point. Perform a one-step Euler prediction on the latent vectors to obtain the time points. The predicted latent vector; the squared Euclidean distance of the forward prediction error and the squared Euclidean distance of the backward prediction error are weighted equally and summed as the pair of adjacent time points. The bidirectional fusion cost is composed of the bidirectional fusion costs of each time point pair, forming a bidirectional fusion cost matrix. The corresponding cost matrix elements are the bidirectional fusion costs of the forward prediction error and the backward prediction error for the same time point pair.

[0009] Furthermore, the entropy-regularized optimal transport soft coupling is solved using the Sinkhorn iterative algorithm. Taking the uniform distribution of the two groups of disordered cell batches at two adjacent time points as marginal constraints, the inner product of the cost matrix and the coupling matrix is ​​minimized, and an entropy regularization term is added. The resulting pairing weight matrix is ​​a non-negative matrix, and its element values ​​represent time points. and The probability of soft pairing between disordered batches of cells is used to transform the alignment at the distribution level into soft pairing supervision at the cell level.

[0010] Furthermore, the bridging samples and target velocities are constructed as follows: endpoint pairs are sampled according to the probability distribution of the pairing weight matrix. ,in for The latent vector, for Latent vectors; then randomly sample coefficients , indicating that uniform random sampling is performed on the interval [0, 1] to calculate the bridging time. and bridging latent vector To form bridging samples Divide the difference of the endpoint latent vectors by the time interval. Construct the constant target velocity corresponding to the bridging sample. .

[0011] Furthermore, the regression training specifically involves: using bridging samples... Input the forward velocity field, output the predicted velocity, and calculate the squared Euclidean distance between the predicted velocity and the target velocity as the forward flow matching loss; pair the endpoints... The backward bridging samples and backward target velocities are interchanged and constructed in the same way. The backward velocity field is input to calculate the backward flow matching loss. The forward flow matching loss and the backward flow matching loss are summed and used as the overall training target for local flow matching, so as to realize the regression training of the time-conditional velocity field.

[0012] Furthermore, to reduce the computational complexity of local supervision, a Top-K truncation approximation is applied to the pairing weight matrix: for each row of the pairing weight matrix, only the K largest pairing weight values ​​are retained, and the remaining weights are reset to zero. Endpoint pairs are sampled based on the truncated weight matrix, and the flow matching loss is calculated, reducing the computational complexity from O(B) to O(K). 2 The computational efficiency is reduced to O(BK) (where B is the batch size), which improves computational efficiency while ensuring training accuracy.

[0013] Furthermore, the distribution-level alignment loss specifically involves: starting from the latent vector at the initial time point, generating predicted latent vectors for each training time point through the integration of the forward velocity field's ODE (Ordinary Differential Equation), which are then mapped by the decoder to a predicted cell set in the gene expression space; calculating the entropy-regularized OT (Optimal Transport) cost between the predicted cell set and the corresponding time point observation snapshot cell set to measure the distribution difference between the two; and averaging the entropy-regularized OT cost over all training time points to form the distribution-level alignment loss, constraining the predicted distribution to be anchored to the observation snapshot distribution.

[0014] Furthermore, the latent space dynamic regularization loss specifically involves: at each training time point, obtaining the predicted latent vector set obtained by velocity field integration, and the observed latent vector set obtained by the VAE encoder directly encoding the observed cells at that time point; calculating and minimizing the entropy regularization OT difference between the two sets of latent vector sets in the latent space to form the latent space dynamic regularization loss, suppressing the generated latent vectors from deviating from the latent space geometry learned by the VAE, and ensuring the consistency between latent dynamics and representation learning.

[0015] Furthermore, the overall training of the method employs a two-stage process: Phase 1: VAE characterization pre-training. Single-cell gene expression data from all time points are mixed and the VAE is trained until convergence, resulting in a stable encoder and decoder, realizing the mapping of high-dimensional gene expression to low-dimensional latent space. The second phase involves joint dynamics learning. In each training iteration, adjacent time points are sampled to batch calculate OT soft coupling and optimize local flow matching loss and VAE reconstruction loss; simultaneously, according to a preset period... The distribution-level alignment loss and latent space dynamic regularization loss are calculated and backpropagated to achieve joint optimization of local dynamics learning and global distribution anchoring, balancing training efficiency and long-term prediction stability.

[0016] A single-cell gene expression dynamics modeling system based on coupled conditional flow matching, the system comprising: The representation learning module is used to acquire single-cell gene expression snapshot datasets at multiple time points and encode the single-cell gene expression vectors in the single-cell gene expression snapshot datasets at each time point into a shared latent space representation. The soft coupling building module is used to represent two sets of shared latent spaces for any pair of adjacent time points, calculate entropy-regularized optimal transmission soft coupling based on the cost matrix, and obtain the pairing weight matrix. The local flow matching training module is used to construct bridging samples and target velocities at intermediate time points based on the pairing weight matrix, and to perform regression training on the time-conditional velocity field using the bridging samples and target velocities. The bidirectional consistency constraint module is used to jointly learn the forward and backward velocity fields of the temporal conditional velocity field, and improve the local transition discernibility under unpaired supervision by the consistency constraints of the forward and backward velocity fields. The global anchoring module is used to introduce distribution-level alignment loss and latent space dynamic regularization loss into long-range generation, so that the predicted distribution generated by the time-conditional velocity field integral is anchored to the observation snapshot distribution and the distribution drift is suppressed. The inference generation module is used to integrate the shared latent space representation of a single cell at the starting time point using the time-conditional velocity field after training during the inference phase, to obtain the shared latent space representation of a single cell at the target time point, and decode it to obtain the single-cell gene expression prediction result at the target time point.

[0017] Furthermore, the functions of each of the above modules are as follows: Representation learning module: It is used to encode high-dimensional single-cell gene expression vectors into low-dimensional shared latent space representations and decode the latent space representations back into the gene expression space. This module has a built-in VAE encoder and decoder, and pre-training is achieved through reconstruction loss and KL divergence regularization to output the latent vector representations of cells at each time point. The soft coupling building module is used to calculate entropy-regularized OT soft coupling between latent space representations at adjacent time points and output the paired weight matrix. This module has a built-in cost matrix construction unit and a Sinkhorn iterative solution unit, which supports the phased construction of Euclidean cost matrix and bidirectional fused cost matrix. The paired weight matrix of soft coupling is obtained through Sinkhorn iteration. Local flow matching training module: used to construct bridging samples and target velocities based on the pairing weight matrix and train the time-conditional velocity field; this module has built-in bridging sample construction unit, target velocity calculation unit and flow matching loss calculation unit, supports Top-K truncation approximation, and realizes regression training of velocity field; Bidirectional Consistency Constraint Module: Used to jointly train the forward and backward velocity fields and apply consistency constraints; this module optimizes the forward and backward flow matching loss simultaneously to ensure the temporal consistency of the bidirectional velocity fields and improve the recognizability of local transitions. Global Anchoring Module: Used to calculate the distribution-level alignment loss and latent space dynamic regularization loss and constrain the long-range generation process; This module has built-in distribution-level alignment loss calculation unit and latent space dynamic regularization loss calculation unit, and participates in the joint optimization of the model according to a preset period to achieve global anchoring of long-range generation and suppress distribution drift; The inference generation module is used to integrate and decode the latent vector at the starting time point to generate single-cell expression predictions at the target time point. This module has a built-in ODE integration unit. Starting from the latent vector at the starting time point, it obtains the latent vector at the target time point through ODE integration of the forward velocity field. Then, it is mapped to the gene expression vector by the decoder of the representation learning module and outputs the final prediction result.

[0018] Compared with the prior art, the beneficial effects of the present invention are as follows: 1) Achieving strong coupling between snapshot alignment and vector field learning: The distribution alignment of entropy regularized optimal transmission is transformed into a cell-level soft-weighted flow matching supervision target, replacing the traditional coarse-grained distribution matching, providing accurate cell-level supervision for velocity field training, solving the problem of loose coupling between the two, and improving the accuracy of local dynamics learning; 2) Improve the discernibility of local transitions under unpaired supervision: By simultaneously learning the forward and backward velocity fields and applying symmetric supervision / consistency constraints, the bidirectional information of the time series is used to constrain the local cell transfer structure, which solves the problem of local transition ambiguity under unpaired snapshots and significantly improves the temporal coherence of the generated trajectory. 3) Effectively suppress distribution drift in long-term generation: The dual global anchoring of distribution-level alignment constraints and latent space dynamic regularization is introduced. On the one hand, it constrains the predicted distribution to fit the observation snapshot distribution, and on the other hand, it ensures the consistency between latent dynamics and VAE latent space geometry. It suppresses the accumulation of integral error from both distribution and latent space dimensions, and greatly improves the robustness of long-term extrapolation prediction. 4) Balancing modeling accuracy and computational efficiency: Stable initialization of the model is achieved through phased cost matrix construction, the computational complexity of local flow matching is reduced by Top-K truncation approximation, and the training strategy of "local iterative optimization + global periodic update" balances training efficiency and global constraint effect, so that the model has good engineering applicability while ensuring high accuracy. 5) Strong generalization ability and wide applicability: This invention can process time-series scRNA-seq data from different biological scenarios, including rapid embryonic development and long-term cell reprogramming. It also supports two types of tasks: time interpolation (completing intermediate unmeasured time points) and time extrapolation (predicting future time points). The generated cell distribution has a high degree of matching with the real distribution and can accurately reconstruct the global trajectory structure of cell development, providing biologists with a reliable computational tool for analyzing the dynamic changes of cells. Attached Figure Description

[0019] Figure 1 This is a schematic diagram of the overall process of the single-cell gene expression dynamics modeling method based on coupled conditional flow matching of the present invention.

[0020] Figure 2This is a schematic diagram of the module composition of the single-cell gene expression dynamics modeling system based on coupled conditional flow matching of the present invention. Detailed Implementation

[0021] The technical solutions of the embodiments of the present invention will be clearly and completely described below with reference to the accompanying drawings. Obviously, the described embodiments are only some embodiments of the present invention, and not all embodiments. Based on the embodiments of the present invention, all other embodiments obtained by those skilled in the art without creative effort are within the scope of protection of the present invention.

[0022] The terms used in this invention are defined as follows: Single-cell gene expression snapshot dataset: refers to a collection of single-cell gene expression data obtained using scRNA-seq technology at multiple discrete physical time points. Each time point corresponds to an unordered cell population with no pairing relationships between cells, denoted as . ,in For physical time points, For this time point Gene expression vector of each cell This represents the number of cells at that time point.

[0023] Shared latent space representation: refers to the low-dimensional latent vector obtained by mapping high-dimensional gene expression vectors through a VAE. Cells at all time points are mapped to the same latent space, denoted as . This enables dimensionality reduction and feature extraction of high-dimensional data.

[0024] Entropy-regularized optimal soft-coupling: refers to the optimal transmission coupling matrix solved under entropy regularization constraints, whose elements are non-negative soft pairing probabilities rather than hard 0-1 pairings, which can better adapt to the heterogeneity and sampling noise of single-cell data.

[0025] Time-conditional velocity field: refers to a velocity field defined by time. and latent vector As input, it is a function that outputs the rate of change of cell state in the latent space, divided into forward velocity field and output velocity field. (From past to future) and backward velocity field (From the future to the past) is the core of describing the continuous dynamics of gene expression in single cells.

[0026] Bridging samples: These are intermediate time points and corresponding latent vectors constructed through linear interpolation between latent vector endpoint pairs at adjacent time points. They are used to simulate continuous temporal dynamics and provide continuous training samples for flow matching.

[0027] The single-cell gene expression dynamics modeling method based on coupled conditional flow matching of the present invention, such as Figure 1 As shown, it includes the following steps: Acquire single-cell gene expression snapshot datasets at multiple time points, and encode the single-cell gene expression vectors in the single-cell gene expression snapshot datasets at each time point into a shared latent space representation; For any two sets of shared latent space representations of adjacent time points, the pairing weight matrix is ​​obtained by calculating the entropy-regularized optimal transmission soft coupling based on the cost matrix. Based on the pairing weight matrix, bridge samples and target velocities at intermediate time points are constructed, and the time-conditional velocity field is regressed and trained using the bridge samples and target velocities. The forward and backward velocity fields of the temporal conditional velocity field are jointly learned, and the local transition discernibility under unpaired supervision is improved by the consistency constraint of the forward and backward velocity fields. To introduce distribution-level alignment loss and latent space dynamic regularization loss into long-range generation, the predicted distribution generated by time-conditional velocity field integral is anchored to the observation snapshot distribution and distribution drift is suppressed. During the inference phase, the shared latent space representation of a single cell at the starting time point is integrated using the time-conditional velocity field after training to obtain the shared latent space representation of a single cell at the target time point, and then decoded to obtain the single-cell gene expression prediction result at the target time point.

[0028] In one embodiment, the shared latent space is represented by a variational autoencoder, which encodes gene expression vectors into a Gaussian latent variable distribution and samples latent vectors, and then the decoder reconstructs the latent vectors back into the gene expression space.

[0029] In one embodiment, the cost matrix adopts a phased construction strategy: during the training warm-up phase, it is constructed using the Euclidean distance of latent vectors at adjacent time points; after the warm-up phase, it is constructed using the one-step prediction error of fusing the forward and backward velocity fields to form a bidirectional fused cost matrix.

[0030] In one embodiment, the bidirectional fusion cost matrix is ​​obtained as follows: a one-step Euler prediction is performed on the latent vector at the first time point based on the forward velocity field to obtain the predicted latent vector at the second time point; a one-step Euler prediction is performed on the latent vector at the second time point based on the backward velocity field to obtain the predicted latent vector at the first time point; the weighted sum of the forward prediction error and the backward prediction error is used as the bidirectional fusion cost; and the bidirectional fusion cost of each time point pair constitutes the bidirectional fusion cost matrix.

[0031] In one embodiment, the entropy-regular optimal transport soft coupling is solved by Sinkhorn iteration, and the pairing weight matrix represents the soft pairing distribution between disordered cell batches at two time points.

[0032] In one embodiment, the bridging sample and the target velocity are constructed as follows: endpoint pairs are sampled according to the pairing weight matrix, and then the bridging time and bridging latent vector are obtained by randomly sampling coefficients, and a constant target velocity is constructed by dividing the difference between the endpoint latent vectors by the time interval.

[0033] In one embodiment, the regression training includes: regressing the target velocity from the forward velocity field to obtain the forward flow matching loss, and constructing the backward flow matching loss by swapping the forward endpoint and the backward endpoint, with the sum of the forward flow matching loss and the backward flow matching loss serving as the local flow matching training target.

[0034] In one embodiment, a Top-K truncation approximation is applied to the pairing weight matrix, which includes: retaining only the K largest pairing weights in each row of the pairing weight matrix and calculating the flow matching loss.

[0035] In one embodiment, the distribution-level alignment loss includes: starting from the latent vector at the starting time point, generating the predicted latent vector at each training time point through the ODE integral of the forward velocity field, and mapping it to the predicted cell set in the gene expression space by the decoder; calculating the entropy-regularized OT cost between the predicted cell set and the cell set observed at the corresponding time point; and averaging the entropy-regularized OT cost at all training time points to form the distribution-level alignment loss.

[0036] In one embodiment, the latent space dynamic regularization loss includes: at each training time point, minimizing the OT difference in the latent space between the latent vector set obtained by the velocity field integral and the latent vector set obtained by the encoder encoding the observed cells at each time point, so as to suppress the generated latent vectors from deviating from the geometry of the latent space.

[0037] In one embodiment, the training process of the method includes two stages: first, VAE representation pre-training is performed; then, in each iteration, adjacent time points are sampled to calculate OT soft coupling and optimize local flow matching loss, and the distribution-level alignment loss and latent space dynamic regularization loss are calculated and backpropagated according to a preset period.

[0038] In one embodiment, taking a zebrafish embryonic development time-series scRNA-seq dataset as an example, the specific implementation steps of the single-cell gene expression dynamics modeling method based on coupled conditional flow matching of the present invention are described in detail. The dataset contains single-cell gene expression snapshots at 12 physical time points, with 2000 highly variable genes and a total of 38731 cells. The goal is to achieve interpolation prediction for unmeasured time points and extrapolation prediction for future time points. This embodiment includes the following specific steps: Step 1: Data preprocessing and VAE representation pre-training.

[0039] Data preprocessing: The library size of the raw scRNA-seq data was normalized and then transformed by log (1+x) to screen out 2000 highly variable genes (HVGs) and obtain a standardized gene expression matrix with an expression vector dimension of 2000 for each cell.

[0040] VAE Model Construction: A variational autoencoder is constructed. The encoder uses a 3-layer fully connected network with an input dimension of 2000, hidden layer dimensions of 512 and 256, and outputs the mean and variance vectors of the latent space. The latent space dimension is set to 50. The decoder uses a 3-layer fully connected network with an input dimension of 50, hidden layer dimensions of 256 and 512, and output dimension of 2000 to reconstruct the gene expression vector.

[0041] VAE pre-training: Standardized representation data from all 12 time points were mixed and randomly divided into training and validation sets. The Adam optimizer was used with a learning rate of 2×10⁻⁶. -3 The training was conducted for 50 epochs with a batch size of 256. The loss function was jointly optimized using reconstruction loss (MSE) and KL divergence regularization. in, For encoder parameters, For decoder parameters, For batch size, The weight of the KL divergence regularization term is set to 0.01. The latent vector is obtained by reparameterizing the Gaussian latent variable distribution output by the encoder. The reconstructed gene expression vector output by the decoder. Let be the gene expression vector of the b-th single cell in the batch. The encoder is given an approximate posterior distribution, where For latent vectors, It is a d-dimensional standard Gaussian prior distribution.

[0042] Latent vector encoding: The standardized expression vectors at each time point are input into the trained VAE encoder to obtain the latent vector representation of each cell, with a dimension of 50, thus completing the mapping from the high-dimensional gene expression space to the low-dimensional shared latent space.

[0043] Step 2: Constructing entropy-normalized optimal transmission soft coupling.

[0044] For adjacent time points (like and , and (etc.), construct the cost matrix and solve for entropy-regularized optimal transmission soft coupling respectively: Cost matrix construction: The first 30 epochs of training are a warm-up phase, using a Euclidean cost matrix with elements of... ,in for The One latent vector, for The One latent vector; after the warm-up phase, construct a bidirectional fusion cost matrix, first through the forward velocity field. right A one-step Euler prediction yields Through the backward velocity field right A one-step Euler prediction yields ,in The time interval is represented by the elements of the bidirectional fusion cost matrix: Sinkhorn iterative solution for soft coupling: using uniform distribution as the marginal constraint and setting entropy regularization coefficients. (C is the cost matrix, (The median) is used to solve the entropy-normalized optimal transport problem using the Sinkhorn iterative algorithm: in, For uniform marginal distribution, The set of coupling matrices that satisfy the marginal constraints, for One of the coupling matrices, Coupling matrix elements, The batch size is determined, and the process converges after 50 iterations, yielding the paired weight matrix (i.e., the entropy-normalized optimal transmission coupling matrix). Its elements for No. individual cells and No. The soft pairing probability of each cell.

[0045] In this context, "coupling" in "the set of coupling matrices satisfying marginal constraints" refers to the soft pairing relationship between two groups of disordered cell batches at two adjacent time points; in the latent space representation, it is the coupling between one set of latent vectors and another set of latent vectors. More specifically, it refers to the coupling at time point t. a mini-batch and time point t b The soft pairing matrix is ​​formed by the pairwise matching weights of each cell / latent vector between mini-batches.

[0046] Step 3: Velocity field training for local flow matching.

[0047] Velocity field model construction: The forward and backward velocity fields use the same residual network structure, with time as the input. Sinusoidal temporal embedding (64 dimensions) and latent vectors The concatenated vector has a dimension of 50, the hidden layer consists of 3 residual blocks, each containing layer normalization, SiLU activation and a linear layer, and the output dimension is 50, consistent with the dimension of the latent space. Construction of bridging samples and target velocities: based on paired weight matrix probability distribution sampling endpoint pair Random sampling Calculate bridging time and bridging latent vector Construct target speed ; Top-K truncation approximation: Set K=32, retain the 32 largest weight values ​​in each row of the paired weight matrix, and set the rest to zero. Sampling is based on the truncated matrix to reduce computational complexity. Flow matching loss calculation: Input the forward velocity field to obtain the predicted velocity. Calculate the forward flow matching loss By swapping endpoints, construct a backward bridging sample and a backward target velocity, and input the backward velocity field to calculate the backward flow matching loss. The total loss for local flow matching is .

[0048] Step 4: Consistency constraints for the two-way velocity field.

[0049] In each training iteration, the parameters of the forward and backward velocity fields are simultaneously optimized. The sum of the forward and backward flow matching losses is used as the local training objective. The velocity field parameters are updated through backpropagation to ensure that the prediction results of the bidirectional velocity fields satisfy temporal consistency: that is, the forward velocity field changes from... arrive The prediction and the backward velocity field from arrive The predictions are matched with each other. Through this symmetric supervision constraint, the local cell transfer structure under unpaired snapshots is clarified, and the identifiability is improved.

[0050] Step 5: Global anchoring constraints generated over a long period.

[0051] Long-range generation refers to the process of generating multiple subsequent time points, especially single-cell distributions at target time points far from the initial observation point, by continuously integrating the learned temporal condition velocity field from the initial observation point.

[0052] Set global update cycle =5, meaning that the global anchoring loss is calculated and propagated every 5 local iterations: Calculation of distribution-level alignment loss: from the starting time point Starting from the latent vectors, the predicted latent vectors for all 12 training time points are generated through RK4 integration of the forward velocity field (step size 0.1). These latent vectors are then mapped to the predicted cell set in the gene expression space by the VAE decoder. ;calculate With observation snapshot cell set The entropy regularization OT cost is averaged across all time points to obtain the distribution-level alignment loss. ; Latent space dynamic regularization loss calculation: For each training time point t, obtain the set of predicted latent vectors generated by integration. The set of latent observation vectors directly encoded by the encoder ,in The learnable parameters of the encoder are represented. The entropy-regularized OT costs of the two parameters in the latent space are calculated and summed to obtain the latent space dynamic regularization loss. ; Joint loss optimization: The VAE reconstruction loss, local flow matching loss, distribution-level alignment loss, and latent space dynamic regularization loss are weighted and summed to form the overall training loss. The weighting coefficient is set as follows: , , The Adam optimizer and OneCycleLR learning rate scheduling are used for optimization, with a maximum learning rate of 2×10. -3 The gradient clipping was set to 1.0, and the training lasted for a total of 300 epochs.

[0053] Step 6: Prediction generation during the inference phase.

[0054] Interpolation prediction: for unmeasured intermediate time points (e.g.) , , ), from the starting time point Starting from the latent vector, the latent vector set at the target time point is obtained by integrating the forward velocity field through RK4 to the target time point; Extrapolation prediction: For future points in time (e.g.) , Using the same ODE integration method, integrate from the latent vector at the most recent observation time point to the target time point; Decoding output: Input the latent vector set of the target time point into the VAE decoder, map it back to the gene expression space, and obtain the final single-cell gene expression prediction result, realizing the distribution prediction and trajectory reconstruction of unmeasured time points.

[0055] In one embodiment, corresponding to the above method, a single-cell gene expression dynamics modeling system based on coupled conditional flow matching is provided, the module composition of which is as follows: Figure 2 As shown, the specific implementation of each module is as follows: Representation learning module: Includes VAE encoder and VAE decoder, both of which are composed of fully connected networks and activation functions. After pre-training, it realizes bidirectional mapping between gene expression vectors and latent vectors, and outputs the latent vector representation of cells at each time point; The soft coupling building block includes a cost matrix construction unit and a Sinkhorn iterative solution unit. The cost matrix construction unit supports the phased generation of Euclidean cost and bidirectional fusion cost. The Sinkhorn iterative solution unit realizes the solution of entropy regular optimal transmission soft coupling and outputs a paired weight matrix. Local flow matching training module: includes a bridging sample construction unit, a target velocity calculation unit, and a flow matching loss calculation unit. The bridging sample construction unit generates bridging samples using linear interpolation. The target velocity calculation unit calculates a constant target velocity. The flow matching loss calculation unit supports Top-K truncation approximation and outputs the local flow matching loss. Bidirectional consistency constraint module: It works in conjunction with the local flow matching training module to synchronously receive the forward and backward flow matching losses, achieve symmetric supervision through summation, and backpropagate the gradients to the forward and backward velocity fields to achieve consistent training of the bidirectional velocity fields. Global anchoring module: includes a distribution-level alignment loss calculation unit and a latent space dynamic regularization loss calculation unit. The distribution-level alignment loss calculation unit measures the difference between the predicted and observed distributions in the gene space, and the latent space dynamic regularization loss calculation unit measures the difference between the predicted and observed latent vectors in the latent space. It outputs the global anchoring loss according to a preset period. The inference generation module contains an ODE integration unit and uses the RK4 integration algorithm to integrate from the latent vector at the starting time point to the target time point. The obtained latent vector is then input into the decoder of the representation learning module 1, and the single-cell gene expression prediction result at the target time point is output.

[0056] The modules communicate with each other via data interfaces to transfer latent vectors, weight matrices, loss values, and other data. The entire system is based on the PyTorch / TensorFlow framework and can be deployed on servers, workstations, and other computer devices, supporting batch processing of scRNA-seq time series data. For the specific working process of each module, please refer to the corresponding processes in the aforementioned method embodiments.

[0057] Another embodiment of the present invention provides a computer device, including a memory and a processor. The memory stores a computer program, which, when executed by the processor, implements the single-cell gene expression dynamics modeling method based on coupled conditional flow matching described in any of the preceding claims. The computer device may be a desktop computer, laptop computer, server, cloud computing platform, or other device with data processing capabilities. The processor may be a computing chip such as a CPU, GPU, or TPU, and the memory may be a storage medium such as RAM, ROM, solid-state drive, or hard disk drive.

[0058] Another embodiment of the present invention provides a computer-readable storage medium storing a computer program, which, when executed by a processor, implements the single-cell gene expression dynamics modeling method based on coupled conditional flow matching described in any of the preceding claims. The computer-readable storage medium can be a volatile or non-volatile storage medium, including various media capable of storing program code such as a USB flash drive, portable hard drive, read-only memory (ROM), random access memory (RAM), magnetic disk, or optical disk.

[0059] The above methods and systems were used to conduct experimental validation on three types of scRNA-seq datasets: zebrafish embryos, fruit flies, and cell reprogramming. Compared with existing mainstream methods (WOT, scNODE, PRESCIENT, VGFM), the results were evaluated using 2-Wasserstein distance (distribution matching degree, lower is better) and average pairwise matching. Using distance (point-level similarity, lower is better) as the evaluation metric, experimental results show that: Interpolation prediction: This invention applies 2-Wasserstein distance to three types of datasets. The distances were all significantly lower than those of the comparative methods and superior to the OT-based method WOT, indicating that the present invention has higher distribution prediction accuracy when supplementing the intermediate unmeasured time points; Extrapolation prediction: The OT-based method cannot perform extrapolation. The prediction error of this invention is much lower than that of continuous dynamic models such as scNODE, PRESCIENT, and VGFM. The degree of error accumulation is significantly reduced, indicating that the global anchoring constraint effectively suppresses the distribution drift of long-range generation. Trajectory Reconstruction: The PAGA connectivity graph and IM distance are used to evaluate the topological consistency of the trajectory. After the missing time points are filled in, the IM distance of the present invention is significantly lower than that of the comparison method, which can better preserve the global connectivity of cell development and generate a stronger temporal coherence of the trajectory.

[0060] Meanwhile, ablation experiments show that removing any one of the bidirectional velocity field, global alignment constraint, or latent space dynamic regularization will cause the model performance to decrease to varying degrees, proving that each core design of this invention is key to performance improvement.

[0061] The single-cell gene expression dynamics modeling method and system based on coupled conditional flow matching of this invention can be directly applied to scRNA-seq time series data analysis in the field of bioinformatics, providing a reliable computational modeling tool for the study of biological processes such as embryonic development, cell differentiation, tumor microenvironment evolution, and cell reprogramming. Furthermore, the method of this invention is based on deep learning and optimal transfer technology, can be deployed on conventional computer equipment, has high computational efficiency and strong robustness, and possesses good industrial applicability and promotional value.

[0062] The above description is merely a preferred embodiment of the present invention and is not intended to limit the scope of protection of the present invention. Any modifications, equivalent substitutions, improvements, etc., made within the spirit and principles of the present invention should be included within the scope of protection of the present invention.

Claims

1. A method for modeling single-cell gene expression dynamics based on coupled conditional flow matching, characterized in that, Includes the following steps: Acquire single-cell gene expression snapshot datasets at multiple time points, and encode the single-cell gene expression vectors in the single-cell gene expression snapshot datasets at each time point into a shared latent space representation; For any two sets of shared latent space representations of adjacent time points, the pairing weight matrix is ​​obtained by calculating the entropy-regularized optimal transmission soft coupling based on the cost matrix. Based on the pairing weight matrix, bridge samples and target velocities at intermediate time points are constructed, and the time-conditional velocity field is regressed and trained using the bridge samples and target velocities. The forward and backward velocity fields of the temporal conditional velocity field are jointly learned, and the local transition discernibility under unpaired supervision is improved by the consistency constraint of the forward and backward velocity fields. To introduce distribution-level alignment loss and latent space dynamic regularization loss into long-range generation, the predicted distribution generated by time-conditional velocity field integral is anchored to the observation snapshot distribution and distribution drift is suppressed. During the inference phase, the shared latent space representation of a single cell at the starting time point is integrated using the time-conditional velocity field after training to obtain the shared latent space representation of a single cell at the target time point, and then decoded to obtain the single-cell gene expression prediction result at the target time point.

2. The method according to claim 1, characterized in that, The shared latent space is represented by a variational autoencoder, which encodes gene expression vectors into a Gaussian latent variable distribution and samples latent vectors, and then the decoder reconstructs the latent vectors back into the gene expression space.

3. The method according to claim 1, characterized in that, The cost matrix adopts a phased construction strategy: in the training warm-up phase, it is constructed using the Euclidean distance of latent vectors at adjacent time points; after the warm-up phase, it is constructed using the one-step prediction error of fusing the forward and backward velocity fields to form a bidirectional fusion cost matrix. The bidirectional fusion cost matrix is ​​obtained as follows: a one-step Eulerian prediction is performed on the latent vector at the first time point based on the forward velocity field to obtain the predicted latent vector at the second time point; a one-step Eulerian prediction is performed on the latent vector at the second time point based on the backward velocity field to obtain the predicted latent vector at the first time point; the weighted sum of the forward prediction error and the backward prediction error is used as the bidirectional fusion cost; and the bidirectional fusion cost of each time point pair constitutes the bidirectional fusion cost matrix.

4. The method according to claim 1, characterized in that, The entropy-regular optimal transmission soft coupling is solved by Sinkhorn iteration. The pairing weight matrix represents the soft pairing distribution between disordered cell batches at two time points. The bridging sample and target velocity are constructed as follows: the endpoint pairs are sampled according to the pairing weight matrix, and the coefficients are randomly sampled to obtain the bridging time and bridging latent vector. The constant target velocity is constructed by dividing the difference between the endpoint latent vectors by the time interval. The regression training includes: regressing the target velocity from the forward velocity field to obtain the forward flow matching loss, and constructing the backward flow matching loss by exchanging endpoints. The sum of the forward flow matching loss and the backward flow matching loss is used as the local flow matching training target.

5. The method according to claim 1, characterized in that, The Top-K truncation approximation is applied to the pairing weight matrix. The Top-K truncation approximation includes: retaining only the K largest pairing weights in each row of the pairing weight matrix and calculating the flow matching loss.

6. The method according to claim 1, characterized in that, The distribution-level alignment loss includes: starting from the latent vector at the initial time point, generating predicted latent vectors for each training time point through ODE integration of the forward velocity field, and mapping them to a predicted cell set in the gene expression space via the decoder; calculating the entropy regularization OT cost between the predicted cell set and the observed snapshot cell set at the corresponding time point; averaging the entropy regularization OT cost over all training time points to form the distribution-level alignment loss; the latent space dynamic regularization loss includes: at each training time point, minimizing the OT difference in the latent space between the latent vector set obtained by velocity field integration and the latent vector set obtained by the encoder encoding the observed cells at each time point, so as to suppress the generated latent vectors from deviating from the geometry of the latent space.

7. The method according to claim 1, characterized in that, The training process of the method includes two stages: first, VAE representation pre-training is performed; then, in each iteration, adjacent time points are sampled to calculate OT soft coupling and optimize local flow matching loss, and the distribution-level alignment loss and latent space dynamic regularization loss are calculated and backpropagated according to a preset period.

8. A single-cell gene expression dynamics modeling system based on coupled conditional flow matching, characterized in that, include: The representation learning module is used to acquire single-cell gene expression snapshot datasets at multiple time points and encode the single-cell gene expression vectors in the single-cell gene expression snapshot datasets at each time point into a shared latent space representation. The soft coupling building module is used to represent two sets of shared latent spaces for any pair of adjacent time points, calculate entropy-regularized optimal transmission soft coupling based on the cost matrix, and obtain the pairing weight matrix. The local flow matching training module is used to construct bridging samples and target velocities at intermediate time points based on the pairing weight matrix, and to perform regression training on the time-conditional velocity field using the bridging samples and target velocities. The bidirectional consistency constraint module is used to jointly learn the forward and backward velocity fields of the temporal conditional velocity field, and improve the local transition discernibility under unpaired supervision by the consistency constraints of the forward and backward velocity fields. The global anchoring module is used to introduce distribution-level alignment loss and latent space dynamic regularization loss into long-range generation, so that the predicted distribution generated by the time-conditional velocity field integral is anchored to the observation snapshot distribution and the distribution drift is suppressed. The inference generation module is used to integrate the shared latent space representation of a single cell at the starting time point using the time-conditional velocity field after training during the inference phase, to obtain the shared latent space representation of a single cell at the target time point, and decode it to obtain the single-cell gene expression prediction result at the target time point.

9. A computer device, characterized in that, It includes a memory and a processor, the memory storing a computer program that, when executed by the processor, implements the method of any one of claims 1-7.

10. A computer-readable storage medium, characterized in that, The computer-readable storage medium stores a computer program that, when executed by a processor, implements the method according to any one of claims 1-7.

Images