Modeling Cell Dynamics and Interactions with Unbalanced Mean Field Schrödinger Bridge

Summary and Strengths

• Thank you. We sincerely appreciate your recognition of the strengths in our paper. In this study, we address the challenge of jointly modeling cell-cell interactions and unbalanced stochastic biological effects using multiple single-cell snapshots data. This introduces a new framework in the integration of interpretable cellular interactions with complex cell dynamics within the optimal transport framework, offering both great theoretical and practical potential for exploring the spatiotemporal dynamics and interactions of single cells.

• Training complexity and hyperparameter selection

• Thank you. Here we adopted the two-phase pre-training strategy detailed in Algorithm 2. This strategy effectively decouples the optimization problem: We first initialize the velocity ($v_\theta$) and growth ($g_\theta$) networks and the interaction potential ($\Psi_\theta$), providing a robust starting point for stable end-to-end training. This significantly reduces sensitivity to challenges associated with PINN losses. We will expand the discussion of this strategy in the manuscript.
• Regarding hyperparameter selection, we found that lowering the mass loss weighting during pre-training improves distribution matching and performance. Additionally, extending training epochs for larger datasets enhances results. As detailed in Appendix C.2 and summarized in Table 14, we maintained most key hyperparameters consistent across diverse datasets, from synthetic to real-world scRNA-seq data, except for the scale parameter ($d_\text{cutoff}$), which is computed directly from each dataset. We also found that increasing the number of training epochs has no adverse effect, allowing us to use a default of (500, 100, 500) epoches. We will include the more detailed discussion of the hyperparameter selection in the paper.
• To support this, we successfully applied these hyperparameters (Epoches {500, 100, 500}, with other parameters same in the table) to two additional datasets—a 100D EB dataset and a zebrafish spatiotemporal (ST) dataset—with results wil be presented below.
• This demonstrates that CytoBridge does not require extensive, dataset-specific tuning and that the presented hyperparameters serve as a reliable default configuration, enhancing the method's practical applicability.

The paper does not report wall-clock or GPU times. Long training cycles could limit practical adoption, especially given the emergence of simulation-free single-cell methods that train much faster (e.g., OT-CFM [1], OT-MFM [2]).

• Thank you. We fully agree that flow matching's simulation-free approach enables much faster computation. To our knowledge, jointly learning growth, interaction, and velocity from snapshot data in a simulation-free manner is challenging but represents a highly promising future research direction. Below, we present the computational time and memory usage for each method on the 100D EB dataset. Our results show that CytoBridge requires a training time of 11 minutes, comparable to other simulation-based methods, with a peak GPU memory usage of 6.3 GB. On the largest pancreatic $\beta$-cell differentiation data with 51,274 cells across 8 time points, CytoBridge requires a training time of 23 minites, with a peak GPU memory of 3.1 GB. Given the dataset's dimensionality and scale, we consider this computational cost to be acceptable. We will list the training time and memory usage in Appendix C.

• Omitted baselines – While MIOFlow (Huguet et al., 2022) is included, newer simulation-free and manifold-aware baselines such as Metric Flow Matching [2] are missing. Because Metric Flow Matching is currently state of the art on the EB dataset, it should be discussed and compared.

• Thank you for pointing this out. Metric Flow Matching is a very important and outstanding contribution. This is a critical comparison that will greatly strengthen our paper. We have now benchmarked CytoBridge against Metric Flow Matching and found that our method achieves improved performance on the key datasets. We show the $\mathcal{W}_1$ distance on the 100D EB dataset as below, and we will add the results of Metric Flow Matching to all relevant tables in the revised version.
• We also applied CytoBridge to a zebrafish spatiotemporal transcriptomics dataset (Dev. Cell, 2022, PMID: 35512701), using 5.25 hpf and 10 hpf as input. We first compared the $\mathcal{W}_1$ distances in the physical space (denoted as 'Space') and gene expression space (denoted as 'Gene') between CytoBridge and other methods. As shown below, CytoBridge achieves improved performance over other methods. To enhance interpretability, we identified top 200 genes influenced most significantly by cell-cell interactions. Subsequent enrichment analysis of these genes revealed key pathways associated with the biological processes. Through interpretable analysis, we identified pathways related to somite development which are critical to zebrafish embryonic development. These findings align with known biological processes, underscoring the framework’s applicability to spatially resolved data.
• We also applied CytoBridge to the mouse hematopoiesis dataset, identifying pathways associated with positive regulation of leukocyte activation and cell-cell adhesion, key processes in hematopoiesis, interactions, and differentiation. Moreover, we calculated the gradients of growth network with respect to genes to identify key genes that contribute to growth dynamics, including Meis1, Nfkb1, and Rap1b. Meis1 regulates hematopoietic stem cell proliferation and self-renewal, Nfkb1 promotes cell cycle entry and survival as a signaling hub, and Rap1b drives cell division via pathways like MAPK, consistent with established biological knowledge. This demonstrates CytoBridge’s ability to uncover biologically relevant mechanisms.
• So, CytoBridge’s preliminary success on spatial data and its robust biological interpretability across diverse datasets highlight its potential. We believe it will inspire important follow up research in the spatiotemporal dynamics research. We will clarify these points in the revised manuscript.

Minor typo

• Thank you. We apologize for the typo. We will thoroughly review the manuscript to address any additional errors and will update all related presentation to ensure clarity.

Thank you again!

Summary and Strenghts

• Thank you. We sincerely appreciate your recognition of the strengths in our paper.
• scRNA-seq offers only a static snapshot of gene expression, lacking the temporal information needed to understand cell state transitions. Our work, similar to recent advances like WaddingtonOT (Cell, 2019) and Moscot (Nature, 2025), focuses on temporally resolved scRNA-seq to capture gene expression profiles across multiple time points.
• We address the challenge of jointly modeling cell-cell interactions and unbalanced stochastic biological effects using multiple single-cell snapshots data. This is an early and essential step in the integration of interpretable cellular interactions with complex cell dynamics within the optimal transport framework, offering both great theoretical and practical potential for exploring the spatiotemporal dynamics and interactions of single cells.

Compared to Works in Trajectory Inference

• Thank you for your valuable feedback. We are pleased to clarify the distinctions between various methods in trajectory inference and their applicable contexts. In single-cell transcriptomics, methods vary depending on the data type:
• Single Time-Point Snapshot Data: In this context, trajectory inference methods are broadly categorized into two groups:
  • Pseudotime-Based Methods (e.g., Monocle, Slingshot, PAGA, Diffusion Map) infer trajectories by ordering cells along a pseudotime axis but often require specifying a starting point, limiting their flexibility.
  • RNA Velocity-Based Methods (e.g., scVelo, VeloAE, VeloVI, DeepVelo) leverage spliced and unspliced counts to estimate dynamics but are constrained by the need for such data. Both approaches infer dynamics from existing data without generating new cell states.
• Multi-Time-Point Snapshot Data: Our work addresses a distinct scenario involving single-cell sequencing data across multiple time points. Due to the destructive nature of technology, this can be framed as a generative modeling task, where dynamics are inferred from distributions at different time points and could be interpolated at unseen time points. Once learned, these dynamics enable the generation of intermediate cell states. Methods like optimal transport, which have gained significant attention in computational systems biology and machine learning, are well-suited for this task.
• Evaluation and Comparability: Our generative approach, unlike pseudotime or RNA velocity methods, evaluates performance using metrics for generated distributions, which traditional methods cannot produce. This fundamental difference in data and modeling objectives, along with our hold-out experiments, renders direct comparisons with traditional trajectory inference methods infeasible. Very few studies, to our knowledge, compare these method types in the same context. A recent review (PMID: 40422408) discusses modeling approaches for different data types. The table below summarizes these distinctions:

• We greatly appreciate your suggestion to consider the TICCI method. TICCI is a highly significant and outstanding work. Moreover, TICCI uses CellChat to compute intercellular communication patterns, and this also can be integrated into CytoBridge’s function $k$, enhancing its compatibility with biological priors and practical value. Although direct comparisons are inappropriate due to differing objectives and evaluation metrics, we are very happy to incorporate TICCI as a reference in our study. We have applied TICCI to our datasets to validate our biological intepretation and will discuss it as a relevant work in the revised manuscript. Additionally, we will expand the discussion to provide a more thorough comparison between these methods.

Limitated improvement, Biological Significance and Interpretation

• Thank you. We address your comments as follows:
• Comprehensive Benchmarking: In our study, we compared CytoBridge against nine representative methods (including two newly added methods), including neural ODE-based, flow matching-based, multi-marginal, and unbalanced approaches. To our knowledge, these encompass the state-of-the-art in the field. CytoBridge demonstrates clear performance improvements across multiple datasets, underscoring its effectiveness.
• Novel Contribution: And we believe our work is the very first to provide a framework that simultaneously learns explicit cell-cell interaction dynamics alongside unbalanced and stochastic effects, a critical biological process difficult be modeled by existing methods. This enables direct inference of cell-cell interactions from snapshot data, offering novel biological insights that constitute the primary significance of our contribution.
• Strong Potential: As noted by Reviewer 6tJf, CytoBridge holds substantial potential for spatiotemporal transcriptomics. We have incorporated an additional spatial transcriptomics dataset and conducted further analyses, demonstrating CytoBridge’s strong biological relevance across both spatial and non-spatial data from three key perspectives:
  • Detection of Cell-Cell Interactions Across Datasets: CytoBridge’s interaction term enables the identification of the presence or absence of cell-cell interactions in a dataset. In our study, we analyzed four biological datasets: mouse hematopoiesis, EB, pancreatic differentiation, and EMT datasets. The first three datasets involve organoid-like systems cultured in vitro, where cell-cell interactions are plausible due to their three-dimensional, multicellular organization mimicking tissue-like environments. Conversely, the EMT dataset, derived from an in vitro cell line experiment, consists of isolated cells with minimal physical or signaling interactions. Our algorithm successfully identified weak interaction patterns in the EMT dataset (Appendix B.5), aligning with biological expectations, while detecting stronger interactions in the organoid-like datasets, consistent with their multicellular nature.
  • Interpretable Analysis of Different Contributions: To enhance interpretability, we identified top 200 genes influenced most significantly by cell-cell interactions. Subsequent enrichment analysis of these genes revealed key pathways associated with the biological processes. Applying this to the mouse hematopoiesis dataset, we identified pathways closely linked to hematopoiesis, interactions, and differentiation (Please see our reply to Reviewer idfk). Moreover, we calculated the gradients of growth network with respect to genes to identify key genes that contribute to growth dynamics, including Meis1, Nfkb1, and Rap1b. Meis1 regulates hematopoietic stem cell proliferation and self-renewal, Nfkb1 promotes cell cycle entry and survival as a signaling hub, and Rap1b drives cell division via pathways like MAPK, consistent with established biological knowledge. This demonstrates CytoBridge’s ability to uncover biologically relevant mechanisms.
  • Extension to Spatiotemporal Transcriptomics: We applied CytoBridge to a zebrafish spatiotemporal transcriptomics dataset (Dev. Cell, 2022, PMID: 35512701), using 5.25 hpf and 10 hpf as input. We find that CytoBridge is better over other methods. Through interpretable analysis, we identified pathways related to somite development which are critical to zebrafish embryonic development (Please see our reply to Reviewer 6tJf). These findings align with known biological processes, showing the framework’s applicability to spatially resolved data.
• So, CytoBridge’s success on spatial data and its robust biological interpretability across diverse datasets highlight its potential. We believe it will inspire important follow up research in the spatiotemporal dynamics research. We will clarify these points in the revised manuscript.

Biological Motivation and Practical Impact.

• Thank you. As outlined in the related works section of our manuscript, simultaneously modeling cell-cell interactions and complex dynamics (e.g., unbalanced stochastic effects) within the optimal transport framework is well-motivated for temporally resolved snapshot data. We have demonstrated the effectiveness of our algorithm across five biological datasets,including one spatiotemporal datasets, highlighting its practical impact and biological relevance.

scale and complexity

• Thank you. We clarify that CytoBridge was evaluated on large-scale, real-world datasets, such as the pancreas dataset with 51,274 cells across 8 time points, and also a spatiotemoral transcriptomics achieving state-of-the-art performance. To address complexity, we further added experiments with 100 dimensional space, which is comparable to the most of methods in this domain, confirming robust performance in high-dimensional settings (results included in the revised manuscript). We comprehensively evaluated computational time and memory usage, confirming they are acceptable (Please see our response to Reviewer rSSy).
• Thus, we believe CytoBridge is both methodologically innovative and capable of handling the scale and complexity of current single-cell datasets.

Thank you again!

Thank you for your continued feedback on our manuscript. We appreciate the time and effort you have spent in reviewing our work. We have carefully considered your comments. Below, we respond to your concerns, providing clarifications and evidence where necessary, and respectfully seek further guidance on points that remain unclear.

As far as I can tell results compared to DeepRUOT in all real datasets do not demonstrate significant improvement. However CytoBridge is the slowest and least memory efficient approach.

• Thank you for your feedback on our CytoBridge method. We respectfully note that, contrary to your comment, our results on real datasets, particularly zebrafish spatial transcriptomics data, show significant improvement over DeepRUOT. As shown below, CytoBridge achieves a Space score of 0.035, a [calculate: (0.261-0.035)/0.261 ≈ 86.6%] 86.6% improvement over DeepRUOT’s 0.261, while maintaining a competitive Gene score (4.712 vs. 4.745). This substantial gain likely stems from our incorporation of biological priors, aligning with Reviewer 6tJf’s observation of our method’s strong potential in spatial transcriptomics. Additionally, addressing Reviewer rSSy’s concerns, our computational cost remains comparable to baselines, with acceptable or no significant increase in runtime or memory usage. |Methods (ST data)|Space|Gene| |-------|-----|----| |MMFM|0.247|5.208| |SF2M|0.265|5.423| |MetaFM|0.268|5.413| |uAM|0.177|6.499| |UOT-FM|0.227|5.173| |MetricFM|0.273|5.366| |MIOflow|0.263|4.720| |Tigon|0.352|4.979| |DeepRUOT|0.261|4.745| |Ours|0.035|4.712|

Biological insight

• Thank you for your feedback on the biological significance of our work. In response to your earlier comment requesting evaluation “with respect to known biology,” we have demonstrated CytoBridge’s strengths in (1) detecting cell-cell interactions across datasets, (2) providing interpretable analysis of distinct contributions, and (3) extending to spatiotemporal transcriptomics. These results align with known biological insights. To ensure we meet your expectations, could you kindly provide specific guidance on what constitutes “novel understanding” in this context or suggest particular analyses we could perform?

PCA space and origincal space

Thank you for your comment regarding PCA space and the original space. We are happy to clarify that, following standard preprocessing in the field, we first select 2,000 highly expressed genes in the original gene space and apply PCA to reduce dimensionality, mitigating dataset noise. PCA defines a linear mapping via a loading matrix, which is reversible. Thus, any gradient vector computed in the PCA space can be transformed back to the original high-dimensional gene space by multiplying with the transpose of the loading matrix. This ensures our model’s efficient computation in low-dimensional space while preserving full biological interpretability in the original gene space.

Scability. I am afraid 50K is a very very small number of cells in the single cell world.

• Thank you for your comment on scalability. We clarify that, in the context of temporally resolved single-cell data using dynamics-based trajectory inference, our analysis of 50K cells is consistent with or exceeds the scale of datasets used in recent works published at NeurIPS, ICLR, and ICML in 2024 and 2025, where the largest methods handled similar or smaller cell counts. These represent outstanding and representative contributions in the field. We are concerned that scaling to 3M cells may extend beyond the scope of our current study and the typical benchmarks in this field. Could you specify which methods or datasets you believe we should benchmark against to address this concern? | Method | Max cells | Highest Dim | Publication | |----------------------------|-----------|-------------|----------------| | MetricFM | 33,701| 100dim| NeurIPS 2024 | | Topological SB| 18,203 | 2dim| ICLR 2025 spotlight | | Wasserstein Lagrangian Flows| 33,701| 100dim | ICML 2024 | | UOT-FM | 20,519 | 50dim| ICLR 2024| | JKO-net|18,203 | 5 dim | NeurIPS 2024 oral | DeepRUOT | 10,998 | 10dim| ICLR 2025 oral|

Thank you for your valuable feedback and we are looking forward to hearing from you!

Sincerely,

The Authors

I am afraid the plan will only partially address my concerns as again you are building on explaining the potential of the framework rather than providing an application which provides novel or concrete biological insight.

As per my previous comments: The the significance of the growth relevant genes compared to other hits is unclear and recovering generic pathways, genes, does not provide novel understanding over these settings.

On a similar note the "Detection of Cell-Cell Interactions Across Datasets" (referring to the magnitude justification lacking actual detection of interactions presented in the authors' response: "CytoBridge’s interaction term enables the identification of the presence or absence of cell-cell interactions in a dataset. In our study, we analyzed four biological datasets: mouse hematopoiesis, EB, pancreatic differentiation, and EMT datasets. The first three datasets involve organoid-like systems cultured in vitro, where cell-cell interactions are plausible due to their three-dimensional, multicellular organization mimicking tissue-like environments. Conversely, the EMT dataset, derived from an in vitro cell line experiment, consists of isolated cells with minimal physical or signaling interactions. Our algorithm successfully identified weak interaction patterns in the EMT dataset (Appendix B.5), aligning with biological expectations, while detecting stronger interactions in the organoid-like datasets, consistent with their multicellular nature.")
is insuffucient as a justification for the method.

Summary and Strengths

• Thank you. We sincerely appreciate your recognition of the strong potential of our method. We address the challenge of jointly modeling cell-cell interactions and unbalanced stochastic biological effects using multiple single-cell snapshots data. This makes our work a new method in the integration of interpretable cellular interactions with complex cell dynamics within the optimal transport framework, offering both great theoretical and practical potential for exploring the spatiotemporal dynamics and interactions of single cells.

Weaknesses 1 and 2: Biological significance and Interactions

• Thank you for your insightful comments. We fully agree that applying our CytoBridge framework to spatiotemporal transcriptomics data holds strong potential, and we appreciate the opportunity to address this point. Below, we outline our perspective on this application and its challenges, and explain our focus on scRNA-seq data analysis in the current work:
  • Challenges in ST: Extending CytoBridge to ST data introduces complexities, particularly in aligning temporal and spatial scales across time-series slices due to technical and biological variability. Non-rigid or non-linear spatial alignment are needed and pose great difficulties.
  • Limited Methods for Joint Modeling: To our knowledge, very few studies simultaneously model spatiotemporal batch effects, cell-cell interactions and complex dynamics (e.g., growth and stochastic) in this context, making the application of our framework to spatial transcriptomics directly very challenging. Applying CytoBridge to scRNA-seq data simplifies complexity, prioritizing cellular interaction modeling. We focus on this in the current work, designating detailed exploration of spatiotemporal data as future work.
• Meanwhile, to show CytoBridge’s potential, we are very happy to present an example applying it to a zebrafish ST dataset (Please see below) (Dev. Cell, 2022, PMID: 35512701). And we believe CytoBridge hold strong biological significance from three key perspectives, even in the non-spatial data:
  • Detection of Cell-Cell Interactions Across Datasets: Even in non-spatial dataset, CytoBridge’s ability to model interaction term enables the identification of the presence or absence of cell-cell interactions in a dataset. In our study, we analyzed four biological datasets: mouse hematopoiesis, EB, pancreatic differentiation, and EMT datasets. The first three datasets involve organoid-like systems cultured in vitro, where cell-cell interactions are plausible due to their three-dimensional, multicellular organization mimicking tissue-like environments. Conversely, the EMT dataset, derived from an in vitro cell line experiment, consists of isolated cells with minimal physical or signaling interactions. Our algorithm successfully identified weak interaction patterns in the EMT dataset (Appendix B.5), aligning with biological expectations, while detecting stronger interactions in the organoid-like datasets, consistent with their multicellular nature.
  • Interpretable Analysis of Interaction Contributions: To enhance interpretability, we identified top 200 genes influenced most significantly by cell-cell interactions. Subsequent enrichment analysis of these genes revealed key pathways associated with the biological processes. On mouse hematopoiesis dataset, we identified pathways closely linked to hematopoiesis, interactions, and differentiation (Please see the reply to Reviewer idfK). This demonstrates CytoBridge’s ability to uncover biologically relevant mechanisms.
  • Extension to Spatiotemporal Transcriptomics: We applied CytoBridge to a zebrafish spatiotemporal transcriptomics dataset, using 5.25 hpf and 10 hpf as input. We first compared the $\mathcal{W}_1$ distances in the physical space (denoted as 'Space') and gene expression space (denoted as 'Gene') between CytoBridge and other methods. CytoBridge achieves better performance over other methods. Through interpretable analysis, we identified pathways related to somite development which are critical to zebrafish embryonic development. These findings align with known biological processes, showing the framework’s applicability to spatially resolved data.

• So while more systematic ST dataset application requires further development especially in spatiotemporal alignment, CytoBridge’s preliminary success on spatial data and its robust biological interpretability across diverse datasets highlight its potential. We believe it will inspire important follow up research in the spatiotemporal dynamics research. We will clarify these points in the revised manuscript.

The significance of the RBF kernel.

• Thank you. The RBF kernel transforms the one-dimensional cell distance $d_{ij}$ into a high-dimensional feature vector. This allows our model to learn interactions between different cells by encoding interaction distances across multiple scales, rather than enforcing interactions only between similar cells. We will add the detailed discussion in the revised version.

Typo, citation, reference, presentation, evaluation metrics.

• Thank you. We apologize for the typographical error and for confusion. We will thoroughly review the manuscript to address any additional errors. We will update all related presentation to ensure clarity.

Algorithm stability and complexity.

• Thank you. We fully agree that flow matching is a very important area to expore, to our knowledge, simulation free jointly learning growth, interaction and velocity from snapshots data is challenging but can be a very intersting future direction. CytoBridge is robust and does not require extensive, dataset-specific tuning, enhancing the method's practical applicability. We have provided a detailed response regarding this to Reviewer rSSy. Please refer to there. Thank you for your understanding :).

Including UOT-FM

• Thank you. The UOT-FM is an outstanding and very important work. We have now benchmarked UOT-FM. We show the $\mathcal{W}_1$ distance on the mouse dataset as below. And also on 100D EB dataset (Please see our reply to Reviewer rSSy). We will add the results of UOT-FM to all relevant tables in the paper. The results show that CytoBridge is better on the key datasets.

Waddington landscape.

• Thank you. The low-lying regions on this landscape correspond to areas of high probability density. These regions represent the stable states or attractors of the system. We apologize for the missing axis labels. We will clearly label the axes to make the visualization easier to understand.

Description of the plots.

• Thank you. In the mouse dataset, the regions with high learned growth rates correspond to the hematopoietic stem cell (undifferentiated cell type) populations. This is biologically consistent with the lineage tracing barcode results (Nat. Mach. Intell., 2023, PMID: 38274364). We will add more discussion in the revised version.
• We agree that stream plots can be challenging to interpret. Regarding the visualization of scores, we will highlight the cell types with higher stability. Additionally, we will highlight Figure 9’s analysis in the main text, which correlates the learned interaction force with cell transition velocity, directly illustrating whether interactions promote or inhibit differentiation, providing a more intuitive biological interpretation.

Questions: Random batch method

• Thank you. We have conducted an ablation study to analyze its impact, and we will include the full results in the revised Appendix. Generally, increasing $p$ slightly improves performance on evaluation metrics by providing a more accurate approximation of the interaction term. However, gains diminish at higher $p$ values, while computational costs increase. Thus, we selected $p=16$ as the default, balancing robust performance with computational efficiency (Please see our reply to Reviewer idfk, thanks:)).

Computation of the interaction

• Thank you. Forces need be calculated at each intermediate time step of the ODE solver, using the current, simulated particle positions. Relying only on fixed, real data points can not capture the dynamic influence of these forces between observed time points.

Initial weights.

• Thank you. Mathematically, we use the empirical measure $\mu_0^N=\frac{1}{N}\sum_{i=1}^N\delta_{x_i}$ to approximate the true distribution $\mu_0$, hence the uniform weights. This convergence to the true distribution is guaranteed when $N \rightarrow \infty$. We'll add this explanation to the paper. Thank you for pointing this out!

Target weights.

• Thank you. As detailed in Appendix A.4, this weight, derived from the number of closest real data points, encourages the growth network to assign higher weights to particles moving into denser state space regions. We'll elaborate more on the point in the revised version.

Algorithm - line 5.

• Yes, correct! In pre-training, we optimize the local mass constraint ($M_t$) for more flexible, fine-grained guidance. Line 12's notation means we sample $(x_0, x_1)$ pairs from the optimal transport plan between consecutive time distributions to build Brownian bridges for score matching. We'll add the explanations in the revised version.

Thank you again!

Dear authors,

Thank you very much for your further elaboration. I will outline a couple of additional points below.

• Thanks for posting additional publications. I find these papers relevant, but many of them have an aspect in common: they use ligand-receptor information (I mentioned this concept in my answer above as well). To me, this kind of biological evidence is important as a before make interactions in dissociated data somewhat identifiable. I do see that using this information may go beyond the methodological scope of the paper, but in future revisions, I believe this should at least be part of the ground truth formulation.

• With respect to potentials, I think I did not express myself clearly enough, and I apologise for that. I appreciate synthetic experiments with known potentials. However, eventually, I do not think they reflect intricate biological potential. Since a ground truth potential is not available unless ligand-receptor informed, I find this experiment a bit difficult to factor into the downstream story.

Having said this, I commend the authors for all the detail they provided in the rebuttal. The amount of work was great, and I hope my criticism is taken as an encouragement for the next iteration of revisions. I confirm the score increase I discussed in my previous answer and remain open minded to a better assessment upon interaction with the other reviewers.

Thank you, and best wishes for the next steps.

Thank you for your valuable feedback. We agree that incorporating ligand-receptor information is crucial for enhancing the biological relevance of our model. We believe our framework offers a good foundation for integrating such information, and we plan to explore this direction in future work to further refine our approach. Thank you once again for your insightful comments, which have greatly enriched our manuscript.

Summary and Strengths

• Thank you. We sincerely appreciate your recognition of the strengths in our paper. In this paper, we address the challenge of jointly modeling cell-cell interactions and unbalanced stochastic biological effects using multiple single-cell snapshots data. This makes our work an early and essential step in the integration of interpretable cellular interactions with complex cell dynamics within the optimal transport framework, offering both great theoretical and practical potential for exploring the spatiotemporal dynamics and interactions of single cells.

W1: Hyperparameter of baselines

• Thank you. We fully agree that a rigorous and fair evaluation is essential. To address this, we have made every effort to enhance the baseline comparisons by thoroughly tuning key baselines. Specifically, we have implemented the following improvements：
  • We have increased the network sizes of the baseline models to ensure their parameter counts are comparable to CytoBridge.
  • We have extended the training epochs for these models beyond default settings to better suit different datasets.
  • We have incorporated a learning rate scheduler into the training process to optimize performance.
• These adjustments have enhanced the performance of most baselines. Nevertheless, our algorithm achieves improved performance over tuned baselines, indicating the robust contribution of CytoBridge. Here we present our tuned results ($\mathcal{W}_1$) on the mouse hematopoiesis data as below, and we will update all relevant tables in the manuscript with the results from these tuned baselines.

W2: Training Complexity

• Thank you. We fully agree that directly optimizing four neural networks could be challenging. To address this, here we adopted the two-phase pre-training strategy detailed in Algorithm 2. This strategy effectively decouples the optimization problem: We first initialize the velocity ($v_\theta$) and growth ($g_\theta$) networks and the interaction potential ($\Psi_\theta$), providing a robust starting point for stable end-to-end training. This significantly reduces sensitivity to challenges associated with PINN losses. We will expand the discussion of this strategy in the manuscript.
• Regarding hyperparameter selection, we found that lowering the mass loss weighting during pre-training improves distribution matching and performance. Additionally, extending training epochs for larger datasets enhances results. As detailed in Appendix C.2 and summarized in Table 14, we maintained most key hyperparameters consistent across diverse datasets, from synthetic to real-world scRNA-seq data, except for the scale parameter ($d_\text{cutoff}$), which can be computed directly from each dataset. We also found that increasing the number of training epochs has no adverse effect, allowing users to use a default of (500, 100, 500) epoches. We will include the more detailed discussion of the hyperparameter selection in the paper.
• To support this, we successfully applied these hyperparameters (Epoches {500, 100, 500}, with other parameters same in the table) to two additional datasets—a 100D EB dataset and a zebrafish spatiotemporal (ST) dataset—with results wil be presented below.

Q1 : Evaluate the integral during training?

• Yes, the integral is evaluated during each training step. This is performed numerically using a standard ODE solver. Computed over ODE paths, this step incurs minimal overhead, ensuring computational efficiency without impeding the training process.

Q2 : Interpretation of the Cell Interaction

• Thank you. We fully acknowledge the significance of providing a clear biological interpretation of the cell-cell interactions modeled by CytoBridge. Below, we present the biological interpretability of CytoBridge from three key perspectives:
  • Detection of Cell-Cell Interactions Across Datasets: CytoBridge’s interaction term enables the identification of the presence or absence of cell-cell interactions in a dataset. In our study, we analyzed four biological datasets: mouse hematopoiesis, EB, pancreatic differentiation, and EMT datasets. The first three datasets involve organoid-like systems cultured in vitro, where cell-cell interactions are plausible due to their three-dimensional, multicellular organization mimicking tissue-like environments. Conversely, the EMT dataset, derived from an in vitro cell line experiment, consists of isolated cells with minimal physical or signaling interactions. Our algorithm successfully identified weak interaction patterns in the EMT dataset (Appendix B.5), aligning with biological expectations, while detecting stronger interactions in the organoid-like datasets, consistent with their multicellular nature.
  • Interpretable Analysis of Interaction Contributions: To enhance interpretability, we identified top 200 genes influenced most significantly by cell-cell interactions. Subsequent enrichment analysis of these genes revealed key pathways associated with the biological processes. Applying this to the mouse hematopoiesis dataset, we identified pathways related to positive regulation of leukocyte activation and cell-cell adhesion which are closely linked to hematopoiesis, interactions, and differentiation. This demonstrates CytoBridge’s ability to uncover biologically relevant mechanisms.
  • Extension to Spatiotemporal Transcriptomics: Our framework is readily extensible to ST data. We applied CytoBridge to a zebrafish spatiotemporal ST dataset (Dev. Cell, 2022, PMID: 35512701). Further details are provided in our response to Reviewer 6tJf.
• So the ability to detect interaction patterns, identify key contributing genes, and uncover relevant biological pathways highlights its significant potential. We will expand the discussion of these points in the revised manuscript to ensure clarity.

CytoBridge scale to higher dimensions?

• Thank you. Our method maintains scalability comparable to these existing neural ODE-based frameworks. In TrajectoryNet, MIOflow, Neural LSB and DeepRUOT, these methods have been tested up to 100 dimensions.
• Following your suggestion, we have now evaluated CytoBridge on the EB data, extending the input from 50 to 100 PCs. The results ($\mathcal{W}_1$, time), reported below, demonstrate that our algorithm remains effective at this dimensionality.

• This finding provides empirical evidence for the robustness of our method in scaling to more complex, higher-dimensional scenarios. However, at dimensions around 1000, our approach encounters challenges similar to those faced by other methods, highlighting the need for simulation-free training strategies, such as flow matching, to effectively handle ultra-high-dimensional settings.

Choice of p in RBM

• Thank you. We have conducted an ablation study to analyze its impact. We report the $\mathcal{W}_1$ distances as below and will include the full results in the revised Appendix. Generally, increasing $p$ slightly improves performance on evaluation metrics by providing a more accurate approximation of the interaction term. However, gains diminish at higher $p$ values, while computational costs increase. Thus, we selected $p=16$ as the default, balancing robust performance with computational efficiency.

RBM inference vs standard inference?

• Thank you. We have conducted an ablation study to assess the RBM impact on model performance. Results show that the full model (without RBM) and the RBM model ($p=16$) yield comparable performance, with the full model slightly better at certain time points. However, RBM significantly enhances computational efficiency, reducing the interaction term’s complexity from $\mathcal{O}(N^2)$ (e.g., 10^6 for 1000 particles) to $\mathcal{O}(pN)$. The tables demonstrate substantial reductions in memory and inference time with RBM. Given its comparable accuracy and scalability for large-scale single-cell datasets, RBM is a justified and essential component of our framework. Results will be detailed in the revised Appendix C.2.

Def 3.2. Figure 1 is never referenced in the text.

• Thank you. We apologize for the typographical error. We will thoroughly review the manuscript to address any additional errors.

Limitations

• Thank you. We will include the limitations and future work in the main text.

Thank you again!