PRESCRIBE: Predicting Single-Cell Responses with Bayesian Estimation

Dear Reviewer Z96M,

Thank you for your constructive feedback. We have addressed each of your points below.

Metrics and Their Descriptions

Q1/W1.1/W2.1. Gene Embedding Ablation Study

Thank you for the suggestion. To clarify, the ESM model operates on amino acid sequences to produce protein embeddings, not gene embeddings, so a direct comparison is not applicable.

However, to address the core of your question about the impact of the gene encoder, we conducted an ablation study comparing our model with three alternatives:

1. GEARS-GNN embeddings.
2. scFoundation embeddings.

Q2/W1.2: Choice of Distribution

We chose the Normal-Inverse-Wishart (NIW) distribution because our model, like the GEARS baseline, operates on preprocessed, normalized expression data, which is continuous. The NIW is the conjugate prior for a multivariate Gaussian likelihood, making it a natural fit for this data type.

In contrast, the Negative Binomial distribution is designed for modeling discrete count data. While adapting our evidence-based framework to a discrete likelihood is a good proposal for future research, it would require a non-trivial theoretical derivation.

Q3/W2.2/W2.3. Encoder Ablation Study

1.  Model Design: Our encoder in Equation. 7 models, both an additive effect (first term) and a non-linear interaction effect (second term). This architecture was adopted from prior work [2] and chosen for its effectiveness and simplicity. The encoder module is flexible, requiring only a deterministic mapping from inputs (c, x) to latent embeddings z (Appx.F. 6, Theorem 2).

2.  Experiment with a Set Transformer: Following your suggestion, we implemented a version of our model using a Set Transformer as the encoder. The results are presented below.

Our detailed adaptation is

$\sum f_{a1}(x_i) + f_{a2}(\sum f_{a1}(x_i)) \rightarrow f_{set}(CONCAT[f_{a1}(x_1),f_{a1}(x_2),...,f_{a12}(c)])$.

While this is a promising direction, the Set Transformer did not outperform our original model on this dataset, potentially due to the large data requirements for training a representative transformer.

Q4/W3.5: Comprehensive Non-Correlation Metrics (RMSE)

As requested, we now provide Root Mean Squared Error (RMSE) metrics (multiplied by 100).

W1: Originality and Significance

The three major contributions of our work are

1. Motivation Originality. Our paper tries to solve overconfidence in single-cell response prediction tasks, a crucial yet often overlooked limitation in previous research.
2. Theoretical Derivation Originality. The derivation and application of a multivariate evidence-based framework for single-cell perturbation prediction. This includes the novel multivariate loss function and the associated parameter update mechanisms.
3. Model Originality. The newly designed model is expected to achieve uncertainty-aware single-cell responses prediction and is expected to benefit by increasing the accuracy of predictive models.

W2.3: Control State Encoding

We apologize for the lack of clarity. The control profile $c_{i,j}$ for a given perturbation is the mean expression profile of all control cells in that experiment. We will state this explicitly in the revised manuscript.

W2.4: Decoder and Normalizing Flow Parameterization

The overall architecture is described in Lines 150-158, and hyperparameters are listed in Appendix B.

W3.1: Ablation Study on losses and other datasets

The $L_1$ term serves as the main regression target (like Mean Absolute Error for GEARS), so we did not ablate it.

The effect of the $L_2$ regularization term is studied in the hyperparameter search (Figure 5b), where setting $\lambda_1 = 0$ corresponds to its ablation.

The other ablations on the other two datasets are shown below:

W3.2: The paper reports correlation metrics as percentages, which is unconventional.

Regarding the calibration plot (Figure 5), the y-axis and x-axis show the percentile rather than PCC values, so it is not in the [-1, 1] interval.

W3.3: While PRESCRIBE seems to have better calibration than GEAR, it does not appear to be well-calibrated. Additionally, the calibration data points are few, with PRESCRIBE resulting in larger error bars.

(1) We think that one limiting factor is the calculation of the ground-truth E-distance itself, which requires a large number of cells to be genuinely accurate, a number often unavailable in real datasets. However, the consistent trend of performance improvement upon filtering (as shown in Table 2) demonstrates that our uncertainty metric has already been put into practical use.

(2) Actually, a large vertical error bar is a sign of good calibration. To conveniently visualize the data, we discretize continuous confidence scores and PCC values by binning them into percentile intervals, but variance within intervals is expected. Therefore, a better visualization might be to show 45-degree slanted error bars, which we are considering adopting.

We updated more standard calibration metrics using Expected Calibration Error (ECE), a standard metric for calibrated regression (Kuleshov et al., ICML 2018). Lower ECE values indicate better calibration.

The table below reports the ECE (multiplied by 100) for PRESCRIBE and baselines.

We will add the formal definition of ECE and a detailed discussion of these results to the revised manuscript.

W3.4: Pseudo E-distance Correlation*

Table 1(b) shows that the correlation between our predicted E-distance and the ground truth improves as the ground truth becomes more reliable (i.e., as $N$ increases). This suggests that the performance reported in Table 1(a) may be underestimated, as the number of available cells in the dataset limits the ground-truth E-distance calculation itself.

Thank you for the rebuttal and for providing additional results. The paper addresses some of my concerns but not all. I agree with reviewer UDmR that the clarity needs improvement, as the paper should be accessible to both ML and Biology audiences. Here are some of the outstanding issues:

• The rebuttal claims that lines 150-158 and the Appendix specify both the Flow $f_{\psi}$ and Decoder functions $f_{\beta}$. Unfortunately, only text descriptions of what those functions are provided, not actual specifications. Could you provide the actual specifications?
• The ablation study on the gene embeddings encoder is a good start. Note that ESM embeddings could also be used to extract gene embeddings since genes make proteins (see [1, 2]).
• Since the focus of the paper is on uncertainty prediction, it's still unclear why the proposed Normal Inverse Wishart approach is better than comparable alternatives, such as CellFlow [2], which models discrete count data instead of log-normalized counts. Besides the ECE metric, could you also provide uncertainty estimation beyond correlation metrics?
• Figure 4 and ECE metrics: Unfortunately, the proposed approach does not appear to be better calibrated than GEARS, although the ECE metric suggests some marginal improvement, which is unclear in terms of statistical significance.
• Table 1: The paper claims that pseudo-E-distance is correlated with E-distance; however, the correlation coefficient seems too small for it to be an effective measure of E-distance. Could you clarify?
• Also, in general, reporting correlation or ECE metrics as percentages is potentially misleading.

Given the above outstanding issues, I am leaning more towards maintaining my score.

References:

• [1] Rosen et al. 2023, "Universal Cell Embeddings: A Foundation Model for Cell Biology"
• [2] Klein et al. 2025, "CellFlow enables generative single-cell phenotype modeling with flow matching"

Category 2: Issues Requiring Further Discussion

We offer clarifications and provide multiple solutions for your consideration.

---

Q3.1 Comparison to CellFlow and Handling of Discrete Data

Clarification on CellFlow:

We thank you for mentioning CellFlow. It is high-quality research on generating post-perturbation profiles. However, we believe our framework and CellFlow serve different purposes and are not alternatives to each other.

1. Primary Goal: CellFlow is a generative model designed to generate single-cell phenotypes under complex perturbations and preserving heterogeneity. In contrast, PRESCRIBE is an uncertainty quantification framework that estimates the reliability of predictions for unseen perturbations by modeling their relationship to the training data.
2. Role of the Flow Module: In CellFlow, flow matching is used to improve the quality of the generative process. In our model, the flow module is used for density estimation in the latent space, which is central to our uncertainty quantification approach.
3. Type of Uncertainty: The uncertainty estimated by the two methods is complementary and not an alternative. CellFlow envisions its future stochastic extension as a valuable direction to model inherent randomness (data) uncertainty in a deeper way. Our framework, PRESCRIBE, is designed to capture wider uncertainty sources by quantifying both model uncertainty (density) and data uncertainty (entropy).

Clarification on Discrete Distributions:

We agree with your point about discrete count data. Our model is, in fact, compatible with this discrete distribution. Because the multivariate Poisson distribution also belongs to the natural exponential family, we can replace the current multivariate Gaussian output layer with a multivariate Poisson layer. This would involve swapping the Normal-Inverse-Wishart distribution for the appropriate conjugate prior. Our current implementation assumes a multivariate Gaussian distribution primarily to ensure a fair and direct comparison with our baselines.

Based on this, we propose the following solutions:

• Solution 1: We will explicitly mention the modeling of discrete count data as an important direction for our future work in the discussion section.

• Solution 2: We will expand the related work section to delineate better PRESCRIBE from generative frameworks like CellFlow and CPA, emphasizing that their primary goal is disentangling and modeling complex effects, whereas ours is quantifying prediction uncertainty.

Q3.2 Additional Uncertainty Metrics

We agree that a comprehensive evaluation of uncertainty is critical. However, for regression tasks, the set of established uncertainty metrics is less extensive than for classification. We currently report SPCC, ACC2, and ECE, and provide calibration plots, which are standard metrics referenced by others [2][3].

If these are insufficient, we would be grateful for any recommendations on other specific metrics you would find more informative for this task.

Q5. Underestimated Correlation and E-distance Estimation

We clarify that the low SPCC is caused by the poor quality of ground truth E, not by a flaw in our E-distance estimation. The original ground truth E was "noisy" because, to include all perturbation types, it had to rely on many sparse perturbations where true E cannot be reliably measured.

Table 1(b) validates this. When a more reliable ground truth E by calculated with high cell counts. While our estimated E-distance score remains the same, its correlation (SPCC) with this higher-quality ground truth E is significantly stronger. For clarity, we have copied both Table 1(a) and 1(b) below.

Impact of Ground Truth Quality (via Cell Count) on SPCC

Norman Dataset

Rpe1 Dataset

K562 Dataset

This demonstrates that our E-estimation is actually effective, and its performance was previously masked by a noisy ground truth.

---

We hope our responses and new results have addressed your concerns in Category 1. For the issues in Category 2, we have proposed several solutions and would greatly appreciate your feedback.

Sincerely,

The Authors

[1] Papamakarios et al. "Masked Autoregressive Flow for Density Estimation." NeurlPS (2018)

[2] Sun, E.D. et al. TISSUE: uncertainty-calibrated prediction of single-cell spatial transcriptomics improves downstream analyses. Nat Methods (2024)

[3] Rosen et al. "Kermut: Composite kernel regression for protein variant effects." NeurlPS (2024)

Dear Reviewer Z96M,

Thank you for explicitly listing the unresolved issues. We appreciate the opportunity to clarify these points.

We have structured our response into two categories based on your comments: Issues with Direct Solutions and Issues Requiring Deeper Discussion.

Category 1: Issues with Direct Solutions

---

Q1. Actual Specifications of the Flow and Decoder Modules

We provide the detailed specifications for the flow and decoder modules as requested.

1. Flow Module (src/nn/flow/mixflow.py)

• Overall Structure:

  • Number of Blocks: 10
  • Input/Output Dimension: $D$ (latent dimension)

• Composition of Each Block:

  1. Masked Autoregressive Flow (MAF)[1] Layer:
     • Network: 3 masked linear layers.
     • Dimensions: $[D, D \times \text{flow\_size}, D \times \text{flow\_size}, 2D]$
     • Activation: LeakyReLU
     • Scale Constraint: Enabled (constrain_scale: True)
  2. Batch Normalization (BatchNorm) Layer:
     • Dimension: $D$
     • Momentum: 0.5

All user-configurable hyperparameters are detailed in Appendix B. Other parameters use the default values from the source code.

2. Decoder Module (src/nn/output/decoder.py)

The decoder models the output distribution using a multivariate normal distribution with a conjugate prior.

• Output Class: MultiNormalOutput
• Parameters:
  • use_prior: True
  • output1: A linear layer mapping from latent dimension $D$ to output dimension $N$, predicting the mean $\boldsymbol{\mu}$.
  • output2: A linear layer mapping from latent dimension $D$ to $(N \times (N+1))/2$, predicting the elements of the covariance matrix $\boldsymbol{\Sigma}$.
• Prior: NormalWishartPrior

3. Prior Specification (src/distributions/multinormal.py)

• Class: NormalWishartPrior
• Parameters:
  • evidence: 1
  • N: Output dimension
  • $\boldsymbol{\chi_1}$: The prior's mean vector.
  • $\boldsymbol{\chi_2}$: The prior's second moment matrix (Prior Covariance + $\boldsymbol{\chi_1}\boldsymbol{\chi_1}^T$).

---

Q2. Ablation Study of ESM2 Embeddings

We supplement the ESM2 embeddings ablation study here. (Metrics are in their original value without multiplying by 100)

Future users can easily switch between different gene embedding backbones (e.g., scGPT, scFoundation, GNN, UCE) by modifying the backbone argument in the configuration. This information will be added to the supplementary section.

---

Q4. Statistical Significance

We performed one-sided paired t-tests comparing our model against baselines.

P-values from One-Sided Paired t-test (Our Model vs. Baseline)

---

Q6. Results Format

Thank you for the suggestion. In the revised manuscript, we will report all metrics (PCC, ACC, ECE, MSE, etc.) as their original decimal values, formatted to four decimal places, without multiplying by 100.

Dear Reviewer tFuR,

Thank you for your insightful comments and valuable suggestions. Our point-by-point responses are detailed below.

Metrics and Their Descriptions

---

Q1. Scalability: Given the use of multivariate NatPN with normalizing flows, how does PRESCRIBE scale with increasing gene dimension or perturbation complexity? What are the memory/runtime bottlenecks?

We benchmarked PRESCRIBE's runtime while varying the model's latent dimension, $N$. The results below show that the average batch time scales in a near-linear fashion.

---

Q2. All benchmark datasets seem to share a common structure and origin. Did you test generalization on held-out cell types or across studies to evaluate transfer performance?

To supplement our standard benchmarks, we performed a new cross-cell-type generalization experiment.

Experimental Design: We used the dataset from Nadig et al. (2025), which contains CRISPR screens in Jurkat and HepG2 cells. We created a transfer task where the model was trained on all perturbations in HepG2 and a subset of perturbations in Jurkat. It was then tested on the held-out perturbations in the Jurkat cell line. The model uses a shared encoder with cell-type-specific flow and decoder modules.

Results & Interpretation: As the table below shows, PRESCRIBE's uncertainty metric remains well-calibrated (SPCC > 0.38) in this transfer task. And by filtering out the 5-10% of predictions with the highest uncertainty, the performance on the remaining predictions consistently improves.

---

Q3 / W2. In the filtering-based accuracy gains, what are the characteristics of the discarded predictions?

Low confidence metric our model predicted induced from two uncertainty sources:

For seen or familiar perturbations, low confidence means either high prediction variance (like GEARS) (learnt by entropy) or low effect magnitude (supervised by ranking loss). (Fig.1 & Table 1).

For unseen perturbations,  we expected the flow module to generate low evidence or likelihood so that it constructs a common feature for those epistemic uncertainties (Fig.3(a)) rather than an arbitrary prediction.

W1. The paper is technically very dense... This may limit accessibility...

Thank you for this feedback. We will revise the introduction and methods to improve accessibility. We will clarify the core contributions for different audiences:

• For the Machine Learning Audience: They can ignore E-distance and prior settings to directly use a multivariate nature posterior network, which mainly focuses on estimating variance and epistemic uncertainties.
• For the Computational Biology Audience: We will highlight the practical utility, thus we consider the magnitude of the effect by using the E-distance concept.

W3. While Spearman correlation is used... no standard calibration metrics are utilized.

We agree and have now benchmarked our model using Expected Calibration Error (ECE), a standard metric for calibrated regression (Kuleshov et al., ICML 2018). Lower ECE values indicate better calibration.

The table below reports the ECE (multiplied by 100) for PRESCRIBE and baselines.

We will add the formal definition of ECE and a detailed discussion of these results to the revised manuscript.

Dear Reviewer UDmR,

Thank you for your thoughtful review and constructive comments. We have addressed each of your points below.

Metrics and Their Descriptions

---

Responses to Questions

Q1. How is the PCC calculated in Fig. 1? Is it across gene LFCs for a given perturbation?

Yes. This is a standard evaluation metric in the field. We use log-fold changes (LFCs) for two reasons:

1.  Biological Relevance: Accurately predicting the change induced by a perturbation is more useful since the downstream study on the perturbation effect is conducted by comparison.
2.  Scoring Robustness: Evaluating PCC directly on post-perturbation expression can lead to artificially inflated scores due to the high sparsity in single-cell data.

Q2. Why should perturbation embeddings be additive in Equation 7? It seems like it should be a set transformer or something similar.

1.  Model Design: Our encoder in Equation. 7 models both an additive effect (first term) and a non-linear interaction effect (second term). This architecture was adopted from prior work [1] and chosen for its effectiveness and simplicity. The encoder module is flexible, requiring only a deterministic mapping from inputs (c, x) to latent embeddings z (Appx.F. 6, Theorem 2).

2.  Experiment with a Set Transformer: Following your suggestion, we implemented a version of our model using a Set Transformer as the encoder. The results are presented below.

Our detailed adaptation is

$\sum f_{a1}(x_i) + f_{a2}(\sum f_{a1}(x_i)) \rightarrow f_{set}(CONCAT[f_{a1}(x_1),f_{a1}(x_2),...,f_{a12}(c)])$.

While this is a promising direction, the Set Transformer did not outperform our original model on this dataset, potentially due to the large data requirements for training a representative transformer.

Q3. What is the intuition for setting $N_H=N$?

Our choice of $N_H = N$ is based on an empirical evaluation of four heuristics proposed in the NatPN paper. We tested these strategies on the Norman dataset to select the optimal setting for our model.

As shown, the 'exp-half' ($N_H = e^{N/2}$) and 'constant' ($N_H = 1$) settings led to training collapse. The 'normal' setting ($N_H=e^{log(\sqrt{4\pi})N}$) underperformed the 'exp' setting ($N_H = N$) (Please refer Table 2 in the paper). We will include these analyses in the appendix.

Q4. It is unclear what you are showing in Table 1(b).

In Table 1(b), $N$ represents the number of cells sampled from the experimental data to calculate the "ground truth" E-distance. The purpose of this table is to show that the correlation between our predicted E-distance and the ground truth improves as the ground truth becomes more reliable (i.e., as $N$ increases). This suggests that the performance reported in Table 1(a) may be underestimated, as the ground-truth E-distance calculation itself is limited by the number of available cells in the dataset.

Q5. How does an untrained model (PRESCRIBE-Null) have such good predictive performance as is shown in Table 2?

The PRESCRIBE-Null model, while its network weights are random, makes predictions by leveraging the prior distribution, which is conditioned on the control cell profile. However, its performance is substantially lower than a more meaningful baseline like 'AverageKnown', which highlights the value of our fully trained model.

Q6. Is Ours-MLP trained with the E-distance ranking loss? If not, it seems like it should be for a fair ablation.

Thank you for this valuable point. For a fair comparison, we re-trained the 'Ours-MLP' baseline using the same E-distance ranking loss as our whole model. The updated results are below.

---

Responses to Weaknesses

W1. The primary weakness of this paper is its clarity... the methods section... inundates the reader in equations...

We sincerely apologize for this and will revise the manuscript to improve clarity.

(1)&(2). In the original paper, we have compiled a comprehensive notation table in Appendix A. We recognize that overloading the symbol $N$ created unnecessary confusion (Q4), and we apologize for this oversight. In the revision, we will implement the following changes: The symbol $N$ will be used exclusively to represent the latent dimension. A new, distinct symbol (e.g., $N_{cell}$) will be introduced for the number of sample cells used in the E-distance calculation.

(3). We will reframe the method to clarify our method: to create a unified confidence score that handles multiple sources of uncertainty.

• For seen or familiar perturbations, E-distance is a powerful metric because it considers both prediction variance (like GEARS) and the predicted effect magnitude. (Fig.1 & Table 1), So identifying this case is our primary goal.

• For unseen perturbations, our NatPN-based approach is designed to map these cases to a low E-distance, which constructs a common feature for them (Fig.3(a)), which is better than arbitrary prediction.

  To demonstrate this works, we supplemented an experiment that shows our model's advantage at identifying these cases (PCC<0.2 & E>0). Analysis conducted on 600+ test points across three datasets.

  % of Data Filtered by Metric	5%	10%	15%	20%
  Ours Discovery Rate (%)	17.1	22.9	25.7	31.4
  GEARS Discovery Rate (%)	8.2	12.2	16.3	24.5

W2. It does seem to be the case that uncertainty estimates from their method are more calibrated... But it's not the case that they are calibrated particularly well either...

We think that one limiting factor is the calculation of the ground-truth E-distance itself, which requires a large number of cells to be genuinely accurate, a number often unavailable in real datasets. However, the consistent trend of performance improvement upon filtering (as shown in Table 2) demonstrates that our uncertainty metric has already been put into practical use.

W3. It is good to see though that there predictive performance... is quite strong. In fact, the paper seems to undersell how strong the predictive performance is.

Our primary focus was on validating our uncertainty quantification framework, so we emphasized the consistent trend of performance improvement upon filtering.

W4. Since GEARS uncertainty estimates are so bad, I think two simple additional baselines should be included...

Thank you for suggesting these excellent baselines. We have run both experiments.

(a) GEARS Predicted Effect Magnitude as Uncertainty: We calculated the E-distance between GEARS's predictions and the control state to use as an uncertainty score.

The results show this is indeed a strong baseline. We think in perturbation models (different from just a regression task). E-distance is a more comprehensive metric than variance alone. It is not a parallel concept to variance; instead, it's an inclusive one, capturing both variance (via self-distance) and effect magnitude (via pairwise-distance). The following correlations between these components and GEARS's original variance score support this:

(b) Embedding Similarity as Uncertainty: We also evaluated using the cosine similarity between gene embeddings of test and training perturbations as a confidence score. We used the Rpe1 and K562 datasets, as defining similarity for the double knockouts in the Norman dataset is ambiguous.

This straightforward similarity-based approach does not yield a well-calibrated uncertainty metric, highlighting the need for a more sophisticated method.

---

Responses to Limitations

Limitations are not discussed in the conclusion at this point.

You can find limitations in Appendix J, where we discuss the model's reliance on the quality of perturbation embeddings and the challenges that remain in fully disentangling epistemic uncertainty.

[1] SEASON COMBINATORIAL INTERVENTION PREDICTIONS WITH SALT & PEPER. ICLR (2024)

Dear Reviewer HczP,

Thank you for all your helpful comments. The following are our responses, following the order of questions, weaknesses, and limitations.

Q1/W1.Why did you focus on extending NatPN with energy distances?

Our motivation was to construct trustworthy perturbation models, using NatPN as it is efficient, flexible, and can capture both model and data uncertainty. We used the E-distance because previous perturbation-related research[1] suggests that the perturbation E-distance can reflect data quality, and we further find that it can efficiently discover unsatisfactory results, so we try to combine them.

We admit that the link with pairwise distance with evidence causes some confusion, as we want to incorporate them together. Actually, for broad ML users, they don't need to consider E-distance and L3; they can just use our extended multivariate NatPN. For perturbation models, we prioritize the utility, so we introduce pairwise distance and E-distance.

Q2.The results of Table 3 are confusing - why random..

This table aims to ablate the impact of filtering. Performance can be significantly enhanced by filtering out the few worst predictions, such as those with negative PCCs in Fig. 1(a). Therefore, as the number of filtered items increases, there is a higher probability of removing the worst predictions, which in turn increases performance. Thus, the stable trend is more significant.

Q3.Why is it important to keep the evidence in the range [N, 2N]?

When the t-distribution degenerates into a Gaussian, the model behaves like GEARS, which only mimics variance. This means it tends to interpolate and underestimate epistemic uncertainty in missing or sparse areas of the perturbation space, e.g., for new perturbation inputs.

Q4.Why the form of the encoder mentioned in section 3.3? Is it appropriate to assume additive effects?

1.  Model Design: The encoder module actually is flexible, requiring only a deterministic mapping from inputs (c, x) to latent embeddings z (Appx.F. 6, Theorem 2), so we take a simple and efficient one.

2. We do not find that the E-distance has a significant relationship with the interaction type. Here are the mean E-distances for several common interaction types analyzed on Norman.

W2.There is an application to only three datasets of one data and perturbation type, when ..

Actually, our baselines, like scGPT and GEARS, all use these three typical datasets, so we followed the same convention. As for transfer performance, we consider two scenarios.

The first is perturbation transfer, such as predicting the effects of new genetic or drug perturbations; the adaptation is simple, replacing $x_i$ with the new perturbation embeddings.

The second is cell type or basal state transfer, which is trickier since different cell lines have very different systems. We conducted a trial where we held out some perturbations in one cell type and used another known cell type to predict them. The design is that two cell types share the same encoder but have different flow and decoder modules. We used the dataset from Nadig et al. (2025), which contains CRISPR screens in Jurkat and HepG2 cells as an example.

Results & Interpretation: As the table below shows, PRESCRIBE's uncertainty metric remains well-calibrated (SPCC > 0.38) in this transfer task. And by filtering out the 5-10% of predictions with the highest uncertainty, the performance on the remaining predictions consistently improves.

W3.Following on the previous point, claims of PRESCRIBE being well-calibrated are limited to the narrow set of empirical results included in the paper. Results like Figure 3b seem cherry-picked.

(1) We think that one limiting factor is the calculation of the ground-truth E-distance itself, which requires a large number of cells to be genuinely accurate, a number often unavailable in real datasets. However, the consistent trend of performance improvement upon filtering (as shown in Table 2) demonstrates that our uncertainty metric has already been put into practical use.

(2) Additionally, we updated more standard calibration metrics using Expected Calibration Error (ECE), a standard metric for calibrated regression (Kuleshov et al., ICML 2018). Lower ECE values indicate better calibration.

The table below reports the ECE (multiplied by 100) for PRESCRIBE and baselines.

We will add the formal definition of ECE and a detailed discussion of these results to the revised manuscript.

W4.Model choices that could have been insightful..

We provide ablation studies below.

W5.There is some confusion about the $\lambda$ s in equation 9 ...

(1) We do not add a $\lambda$ coefficient for $L_1$ because it is the central optimization function, and we wanted to reduce the hyperparameter tuning burden; other parameters are all smaller than 1.

(2) We think the coefficient does not affect this core design: when $\partial L_1/\partial p_v$ approaches zero, $\lambda_3 \partial L_4/\partial p_v \neq 0$ when $\lambda_3 \neq 0$.

Limitation The recent paper [1] explores the performance issues of many of the baseline models and datasets used in this paper.

After reading this paper, we agree that the reference frame is essential in perturbation models. We think your concern may be that the E-distance is only correlated with the old metrics. Therefore, we conducted a quick preliminary study where we used $\mu_{all}$ to recalculate the GEARS metric and obtained the following results:

We find that although the new E-distance is more correlated with the new PCC, it does not dominate the strong E-distance signal with PCC(new), as the old E-distance is still highly correlated with the new PCC.

And since the paper was published very recently, we do not have time to update all our research, but would like to include it in the future.

[1] Peidli, S., Green, T.D., Shen, C. et al. scPerturb: harmonized single-cell perturbation data. Nat Methods.(2024)

Dear Reviewer HczP:

We've updated our general response with a more detailed and clear explanation of our core idea and modules, based on your feedback. We would be grateful if you could review our responses one last time before the discussion period ends and let us know if we've addressed your concerns. Your guidance has been extremely helpful.

Best regards,

Authors

Dear Reviewer HczP,

We understand that chasing down your reply is not our job, and we do not intend to add any pressure on your busy schedule. However, as we are getting closer to the end of the discussion phase, we would really appreciate it if you could be so kind as to let us know if we have properly addressed your comments and questions in the rebuttal and if anything can be further clarified.

Many thanks in advance!

Authors

Dear Reviewer YtU5,

Thank you for all your helpful comments. The following are our responses.

Metrics and Their Descriptions

Q1. The author(s) initially present experimental results for all the datasets and then focus the rest of the exploration on the Norman dataset. What is the reason for focusing on the Norman dataset for the later analyses? Are the results robust across all three datasets? Similar experimental results for the other two datasets (e.g., in the Appendix) would be helpful here and could improve the clarity and significance of this work. If this is not feasible, an explanation for why this analysis is specific to the Norman dataset is necessary for transparency.

(1) We conducted experiments across all datasets, with the exception of the N-increasing analysis, generalization analysis, and ablation study. The following are the reasons and supplement experiments.

(2) We supplement the N-increasing analysis for all datasets here.

(3) For the generalization analysis, we focused on the Norman dataset because it involves double knock-outs, which allows us to mimic different degrees of generalization. The other datasets involve single knock-outs, for which the degree of generalization is harder to define. Therefore, for those datasets, we provide the average training and test set confidence scores instead.

(4) For the ablation study, constrained by time, we only provide an ablation on the loss function here.

(5) The hyperparameters are different from those for the Norman dataset, but the tuning pipeline remains the same.

Q2/W1.2.In the presentation of the method in 3.1, the author(s) fix the proportion between, and (Eq. 4). Could the authors provide the theoretical or empirical justification for this specific modeling choice? A brief ablation study showing that the model is robust to this choice, or that this choice is optimal, would strengthen the clarity of the model design choices here.

(1) This modeling choice is based on previous research[1].

(2) The Problem (Ambiguity): The original loss function ($\mathcal{L}_{1}$) does not depend on the individual parameters that control uncertainty ($\kappa$ and $\nu$). Instead, it depends on a combination of them. This means that infinitely many different combinations of these parameters can produce the exact same loss value, making it impossible for the model to learn the correct, unique uncertainties. This issue is called "degeneration".

(3) The Solution (Fixing the Ratio): To solve this, the authors propose to link, or "couple," the two parameters with a constant ratio, $r$, such that $\nu_{i}=2\kappa_{i}$. This removes the ambiguity and eliminates the need for a separate, flawed regularization term that was used in prior work.

Q3.How robust is the model to misspecification of its Normal-Inverse-Wishart distribution? It may be nice to know how strong the regularization of the term is in handling any model misspecification here, since it encourages the pseudo E-distance derived from the Normal-Inverse-Wishart prior to match the training data E-distance. A discussion of its robustness or a simple experiment on synthetic data with a known non-Gaussian ground truth (but still reasonable for single cell settings) to demonstrate that not only the predictions but also the uncertainty scores are robust would clarify the limitations and strengths of this paper.

This is useful advice, and we conducted the following experiments in response.

In single-cell preprocessing, data is usually first log-normalized. Although raw count data is very sparse, it sometimes contains a few very large counts. Thus, we first normalize to increase differences, then take the logarithm to reduce the long tail. In this way, single-cell data can be treated as having a somewhat Gaussian distribution. However, as you suggested, in a non-Gaussian setting, we did not perform log-normalization. Instead, we used min-max linear normalization, which yields a very asymmetric distribution as our input. Other operations remained the same. We still use the Norman dataset as an example, and the following are the results.

We conducted two trials: the first uses a non-normal input and a non-normal E-distance, while the second uses a non-normal input and a normal E-distance.

The results show that our model cannot handle non-Gaussian distributions. We also note that using an E-distance derived from a normal distribution assumption deteriorates model performance in this setting. We think this may be because the E-distance is a distance calculation suitable for Gaussian distributions. This inspires the idea that the E-distance could be replaced with another distribution distance metric, but that is outside the scope of this paper. However, we sincerely thank you for your advice and will include this limitation regarding the distribution assumption in our paper.

[1] Multivariate Deep Evidential Regression.