Predicting Time-Varying Flux and Balance in Metabolic Systems using Structured Neural ODE Processes

We thank the reviewer for the provided insightful comments and detailed feedback. Please find our clarifications below:

1. Little to no technical novelty.

The paper feels like the application of Neural ODE processes to DFBA, with minor tweaks to make it work better. Although the results are convincing, I am not sure that ICLR is a relevant venue for this work.

We would like to emphasize why the considered problem is impactful. Metabolic flux measurements, particularly at the single-cell level, are both expensive and technically challenging. Determining the compounds in a solution requires mass or NMR spectroscopy, followed by isotope subtyping or similar techniques to measure reaction fluxes. These methods rely on highly specialized instruments and experimental setups, resulting in very limited ground-truth data. Machine learning and data-driven methods are therefore highly suited to address these challenges. Please refer to our General Response 1 for further details.

While we acknowledge that several methodological components of our framework have appeared in earlier works individually and in different contexts, we argue that the combination of these components is novel in the context of our problem statement. Specifically, modeling flux-balance estimates as a random process with a biologically motivated approach is unique to our work. Please refer to our General Response 3 for additional context.

We thank the reviewer for the provided insightful comments and detailed feedback. Please find our clarifications below:

1. Can you please address the concern regarding the destructive measurement process of scRNA-seq?

It seems different cells at different time points are directly taken as a time series, which seems inaccurate. How do you address this challenge?

We address this challenge by formulating our task as a random process prediction task. A random process can be interpreted as either a distribution over trajectories or a distribution over values at each time step; we adopt the latter perspective. As the reviewer correctly noted, sampling from prior time steps does not yield precise values for the next step, which is why we predict a distribution instead of a fixed value. This distribution inherently captures the underlying uncertainty in the prediction. Our approach models random processes in a structured manner, focusing on distributions rather than trajectories. Please refer to our General Response 2.

2. The lack of structure might be the case for NODEP but not for other NODE models such as RNN-ODE (Rubanova et al.) or GRU-ODE (De Brouwer et al.).

While we agree that RNN-ODE and GRU-ODE introduce structure in their architectures, their outputs are trajectories rather than distributions. Our focus is on predicting distributions over time, not trajectories, as our problem is modeled as a random process. To the best of our knowledge, none of the existing neural process architectures incorporate structural modeling.

3. Most of the results are given in the form of training curves, which makes rapid evaluation of performance difficult.

The scales of some subplots are sometimes off, making it hard to distinguish between methods.

We chose to represent results as training curves to save space, given the comparison of three models across four metabolic pathways and five cases each. For irregularly sampled data, we have provided tables in the Appendix. Regarding the scale, we presented all plots on a log scale, but differences between metabolic pathways can lead to varying y-axis limits, which may appear inconsistent.

4. I am pretty sure that RNN-ODE and GRU-ODE could also be used for this problem.

I encourage the authors to expand their list of baselines, including RNNs or even one-step-ahead linear regressors.

While RNN-ODE and GRU-ODE could be adapted to this problem, these models predict values, not distributions. Estimating a distribution from these methods would require running each trajectory multiple times, which is computationally inefficient. Neural process models, by design, directly account for this need, which is why we focused on them instead. Please refer to our General Response 2.

5. I didn’t fully understand why the distribution of the latent is that important, given that you end up modulating it with a non-linear function for the final predictions. Could you please motivate this choice more rigorously?

In our initial experiments, aligning the latent distribution with the final output distribution led to a slight reduction in MSE, which motivated this choice. However, as the reviewer correctly pointed out, the theoretical justification for this choice remains inconclusive. We acknowledge the need for a more exhaustive investigation and plan to explore this in future studies.

We thank the reviewer for the provided insightful comments and detailed feedback. Please find our clarifications below:

1. The authors don't make a strong point for why the forecasting problem is specifically impactful.

The primary reason this problem is impactful is due to the high cost and difficulty of performing metabolic flux measurements, especially at the single-cell level. Determining flux requires specialized instruments and experimental setups, such as mass or NMR spectroscopy combined with isotope subtyping. As a result, ground truth data is scarce. Please refer to our General Response 1.

The value of the performance from the 80/20 split isn't limited to the 20% savings in cost but rather demonstrates that our model has captured the underlying dynamics, providing confidence in using the predictions for downstream analyses. While we considered other splits, the limited time span of our data (16 days) restricts the number of viable splits, as some context time points are required for training. In the Appendix subsection Effect of Varying Context Length, we explored these variations and selected the 80/20 split accordingly.

"How good is good enough?" is indeed an important question. Addressing it requires downstream analysis of predictions, such as comparing ground truth metabolic-flux/balance values or biomass objective values derived from predictions. As noted in the Conclusions section, this would require time-varying scRNA-seq data aligned with ground truth metabolic data. Unfortunately, such datasets are not publicly available due to the challenges of data collection, so we left this for future work. Meanwhile, we evaluate models based on MSE and prioritize those with lower values. Our model's ability to generalize the dynamics enables practitioners to perform gene-knockout studies or predict biomass yield without solving the dFBA optimization problem, which often requires domain expertise to formulate.

2. I did not fully understand the inputs and outputs of the method.

We train separate models for each task: gene expression, metabolic flux, metabolic balance, metabolic flux with gene knockouts, and metabolic balance with gene knockouts. Each model adapts to regular or irregular sampling. During training, the encoder takes in context points, calculates the latent variables, samples the latent $l_0$ and $d$, and passes them to the decoder. The decoder predicts latent values for target points via a neural ODE, from which it derives distributions for target timepoints. During testing, we evaluate all available points, including unseen timepoints during training.

3. Could the problem be set up as a regular prediction task with simpler architectures?

This is not a regular prediction task because we cannot track gene expression of individual cells over time; measurements destroy the cells. Consequently, flux and balance estimates cannot be determined for single cells across time, and the resulting trajectories are random processes with uncertainty about future values. Modeling this as a random process, which predicts a distribution instead of a point, is a natural fit.

To achieve similar results with LSTMs or GRUs, multiple runs per trajectory would be required to estimate a distribution at each timestep, which is computationally inefficient. Furthermore, the chemical kinetics underlying these systems are commonly formulated as ODEs, motivating the use of neural-ODE models.

4. Baselines are from NODEP; what about trying simpler methods or time-series foundation models?

Although simpler methods could be explored, this is not purely a time-series problem but a random process problem, where the predictions are distributions rather than points. Moreover, chemical kinetics are inherently described by ODEs, which aligns well with NODEP-based approaches. To estimate distributions using simpler time-series methods, multiple runs would be needed, making them less efficient. Please refer to our General Response 2.

5. Unfortunately, the authors do not show any interpretation of the learned systems.

We chose a latent architecture for two reasons: (a) scalability, given the large number of genes/metabolites, and (b) suitability for irregularly sampled data, as shown in prior works like [1]. Since precise formulations of metabolic pathways via chemical kinetics are still ongoing, we did not focus on interpreting latent mappings but agree that exploring interpretability would be valuable.

[1] Rubanova, Y., Chen, R. T. Q., & Duvenaud, D. (2019). Latent ODEs for Irregularly-Sampled Time Series.

6. Misc

We appreciate the suggestion for ablation studies, and will incorporate in our revisions. We chose not to explore comparisons on other tasks because it diverges from our primary motivation of predicting metabolic flux and balance. However, we understand the importance of broader comparisons and will consider this in subsequent work.

We thanks the reviewer for their further comments. Please find our responses below:

1. This limitation seems to be mostly addressed by scFEA in your current work

To an extent you're correct, however, in the real world settings, scRNA-seq data isn't usually considered for metabolic analysis, instead what we usually have are urine or blood samples from which we get the metabolomic profile. This is aligned with our methodology where our flux model is trained directly on flux time samples.

We use the scFEA framework to get estimated flux and balance data because the real world data is difficult to procure, and thus in our context the scFEA should be considered as a proxy for the real world flux-balance data-- something we assume in the section Problem Formulation.

Although, we did perform experiments with single-cell gene-expression data for precisely for the reason that it's a real world data and is highly correlated with the amount of metabolic pathway molecules-- discussed in the section 4.2 Gene-expression.

Ideally, we would've preferred working directly with real world flux-balance data. However we were unable to find such open sourced datasets, and as a result we resorted to scFEA estimates. We discuss this in more detail in the section Conclusions, along with open questions. And we're trying our level best to procure real world data but till then, we think it's best to use scFEA as a proxy for real world data to make methodological progress.

2. I still believe this would lend itself very nicely to a simple predictive architecture such as LSTM or GRU (ignoring the rest of the observed distribution).

We appreciate the suggestion and agree that LSTM/GRU-based approaches could serve as a valuable baseline. However, adapting these methods for our specific setting would require additional effort beyond a straightforward prediction task using LSTMs or GRUs, as these methods do not inherently consider the underlying dynamics or provide distributional outputs. That said, we do recognize the potential value in exploring this direction and will seriously consider it for future work. Thank you for highlighting this opportunity!

Thank you for your response, clarifying some aspects of the proposed method.

1. Why is it impactful? Cost and complexity

I understand that collecting the right type of data for the metabolic flux and balance prediction can be complex and costly. However, this limitation seems to be mostly addressed by scFEA in your current work, reducing the problem to collecting time-series scRNA-Seq data.

3. Regular prediction task with simpler architectures. Distribution-level problem

I agree that measurements are destructive so that you don't have paired data available. But if I understand correctly, you currently measure prediction quality exclusively at the level of the mean squared error on the mean of the distribution. I am fine with that, but as a prediction target, I still believe this would lend itself very nicely to a simple predictive architecture such as LSTM or GRU (ignoring the rest of the observed distribution).

I would like to thank the authors again for the clarifications provided. Sadly, overall, I maintain my evaluation. To change my score, the task itself would have to better motivated and explored in-depth to be convincing, and the technical novelty seems limited. In particular, the baseline suggestions raised by me and other reviewers remain unaddressed in the current version of the paper.

We thank the reviewer for the provided insightful comments and detailed feedback. Please find our clarifications below:

1. Limited advancement in the individual methodological components used.

While we agree that our framework uses several methodological components that could already be found in earlier works on an individual basis in different contexts, we believe that the sum total of those components is still novel with respect to our problem statement and biological motivation of modeling the flux-balance estimates as a random process. Please refer to our General Response 3.

2. Does the model only take time as input? My understanding is that $F$ is the decoder. Is using just time sufficient?

We agree that we could have been clearer here. The reviewer is correct that $F$ is the trained decoder, taking $t$ and sampled latent $l_{0}, d$ as input. During inference, the input should be some prior context trajectory $(t_{i}, y_{i})_{i \in I_C}$ to determine the latent and $t$ for which we want to get the samples. We will clarify this in the revised version of our paper.

3. What do the $Y$ variables mean? $G, M_{f}, M_{b}, \tilde{M_{f}}, \tilde{M_{b}}$ are not clear.

The $Y$ variables refer to the different types of data we consider in our modeling. Note that $Y$ is a probability distribution in our context. The different types of data are:

• $G$: Gene-expression.
• $M_{f}$: Metabolic-flux.
• $M_{b}$: Metabolic-balance.
• $\tilde{M_{f}}$: Metabolic-flux for the gene-knockout case.
• $\tilde{M_{b}}$: Metabolic-balance for the gene-knockout case.

We train a separate model for each type of data above. For each kind of data, we have samples from a particular time, which are understood as being generated by the underlying probability distributions.

4. The metrics are limited, why didn’t you try Wasserstein distance and/or maximum mean discrepancy (MMD) distance?

We agree that several other metrics could have been reported as well. However, as mentioned in the subsection of Experimental Settings under the heading Evaluation Metric, we can interpret the MSE loss as the squared distance between the means of distributions, which satisfies us with the legitimacy of our metric. Additionally, another reason was the use of test-MSE as a metric in prior neural-processes papers, making it an apples-to-apples comparison with our task.

5. What is the setup for predicting metabolic flux and metabolic balance with gene knockouts?

Is this for test gene knockouts that are seen during training, or gene knockouts unseen by the model during training? This is not entirely clear.

During test time, we consider gene-knockout configurations that aren't present during training. We wanted to see whether the framework generalizes to unseen knockout configurations and unseen timesteps, and we observe that it does. This is very useful from a practitioner's point of view because a sufficiently well-trained model can generate samples for newer gene-knockout configurations, saving the time needed to perform the gene-knockout experiments.

6. Equations throughout the text don't have labels -- it would be helpful to numerically label all equations and reference them in the text.

We thought the equations could be read within the context of the text around them. However, we agree that labeling them would make it easier for the reader. We will revise them accordingly.

7. The learning objective is not very clear.

Our learning objective is an ELBO loss, with a slight modification where we use the context points to calculate the prior of the latent variables $L_{0}(l_{0}|\mathcal{C})$ and $D(d|\mathcal{C})$. This is similar to the loss used in the Neural-ODE Process paper.