Recovering Time-Varying Networks From Single-Cell Data

We thank the reviewer for the thoughtful comments. We address the weaknesses and questions below.

1. This is indeed an important discussion. However, in this manuscript, we evaluate our method by finding support in already established benchmarks and known biological processes. Investigating novel edges or processes in detail would require further mechanistic studies, which are beyond the scope of this paper. Our goal is to demonstrate the model's ability to recover known regulatory interactions effectively.
2. We address this concern through gene set enrichment analysis on human datasets, showing that genes added between time points are enriched for relevant biological processes that align with the underlying biological response. This suggests that Marlene understands the biological context of the sample.
3. We mention Dictys in lines 59-60. Basically, while this is a dynamic method, it cannot be applied to our setting as it relies on a multimodal approach (using both scATAC-seq and scRNA-seq types of data which is less prevalent and more costly). The GRN inference in Dictys is achieved by leveraging TF binding sites obtained from scATAC-seq, so it has more information than scRNA-seq.
4. Perturbation analysis would indeed be a valuable addition. However, as noted in lines 513-515 in the manuscript, our goal was not to perform causal inference or perturbation studies, but to uncover the underlying GRNs given a biological dataset. We see this as a promising area for future work.

We hope these answer the reviewer's questions.

We thank the reviewer for the thoughtful and detailed feedback. We appreciate your acknowledgement of the methodological contributions in our paper. We hope to address the weaknesses below.

1. We appreciate the suggestion for direct experimental validation. In this manuscript, we evaluate our method by comparing it to well-established benchmarks and known biological processes. Investigating novel edges or processes in detail would require further mechanistic studies, which are beyond the scope of this paper. Our goal is to demonstrate the model's ability to recover known regulatory interactions effectively.
2. Scalability was not the primary focus of this manuscript. We discuss in the Limitations section that more efficient implementation (e.g., using FlashAttention) would lead to a better capability of the model to handle larger numbers of genes. We mainly view the current manuscript as a methodological contribution, but we would be happy to extend the method further to accommodate larger datasets.
3. Perturbation analysis is important, however, as noted in lines 513-515 in the manuscript, this is a different task. We were mainly concerned with uncovering the underlying GRNs from time series data which is already a significant goal. Causal inference and knockouts were not out goal and it would require other types of data, such as interventional data.
4. We thank the reviewer for highlighting the need to discuss the usefulness of time-varying GRNs. Since most biological processes are dynamic, time-varying GRNs are necessary to understand the evolution of regulatory interactions in the cell. For instance, understanding what TF-gene links are disrupted could help researchers discover drugs targets for specific TF-gene connections. We will expand our motivation section to add applications of time-varying GRNs.

Answers to the reviewer's questions can be found below.

1. We minimally tuned the parameters in this study (mainly learning rate), favoring significance of the results. While we do not use a hyperparameter tuning approach in this work, we believe that the results will improve with proper tuning.
2. Due to cost, most scRNA-seq datasets contain few time points, so we do not necessarily view this as a limitation of the method. However, other temporal units such as Structured State Space models (S4), can help with longer time series. Regarding the irregularity of time points, some of the datasets we used in the experiments contained irregular time points. For instance, the COVID-19 dataset contained time points from days 0, 2, 10, 28. Meanwhile, the time points in the aging dataset were separated by several years of age.
3. The number of cells of a single type varies greatly for each of the datasets (some only contain a few hundred, other thousands). Some types contained few cells, which is where meta-learning helps. We do not report the breakdown of the datasets by cell type, but we will add this information to the supplement.
4. Predicting cell type rather than exact gene expression is one way to mitigate this problem. The lung aging dataset is constructed by combining several datasets that come from different sources, thus it is quite prone to batch effects. However, as we show, Marlene performs better than most baselines for this dataset.
5. Connection strength could be correlated with the attention weights, however, it is unclear how meaningful the precise weights are for these matrices (interpretability of attention weights is limited). Therefore, we prefer to take the top x% of the edges and study them as a set, rather than looking at the weight. Directionality of the regulatory link is given by the attention (these matrices are not symmetric) and we analyze them as directed tuples in the results section.
6. In addition to static interactions, we also include gene set enrichment analysis of dynamic genes (i.e., genes that change between time points). This provides insight into the dynamic biological response and could, for example, help drug discovery efforts to target the right TFs.
7. Many of the baseline methods are static. Therefore, the lack of dynamic modeling might help explain some of the discrepancy. Also, the modeling scheme is different, as most of these baselines are not based on highly expressive nonlinear networks, but rely on simpler shallower approaches.

The reviewer has praised the work, and their main concerns were with additional experimentation or extensions of the work. Given this and the comments, we feel like the score of 3 is too low for this work. We would appreciate it if the reviewer could increase their score.

We thank the reviewer for their comments. The points raised regarding the interpretation of the adjacency matrix are valid, and we address them below.

First, we acknowledge that attention weights do not inherently capture regulatory interactions. The limitations of interpreting attention weights have also been detailed in papers like "Attention is not Explanation" [1]. However, this is why we believe the softmax + reconstruction step add a level of biological plausibility. The softmax plays a double role: 1) it enforces a probabilistic interpretation. Each gene's regulatory influence is weighted across all TFs, i.e., the TF with most influence would be assigned a higher weight. 2) Most importantly, the softmax does not allow negative weights. Indeed, if the weights are close to 0, then they contribute minimally in the reconstruction step, otherwise, the larger the magnitude, the more effect this particular link will have in the reconstruction step. Having negative weights on top of the complex nonlinearities in the model, would make the interpretation of such weights very complicated. In short, we believe that the attention can be interpreted as an adjacency matrix because it is directly being treated as one through a message passing step during reconstruction.

We agree with the reviewer that interpreting the exact values of these weights is limited, which is why we only take the top x% of the edges as the most influential and discard the specific weight. Given this, the computation of AUROC and AUPR values is also not very informative or meaningful. For this reason, we prefer direct overlap tests, like Fisher's exact test. If, instead, the suggestion is to compute these metrics by increasing the number of selected edges, this is also problematic as the set of ground truth edges is very small. Hence, increasing the number of selected edges would simply increase false positives.

Reducing the number of rows in A (i.e., deleting genes) would impact the performance of gene enrichment analysis. We showed in the results section that the genes selected by the method are relevant to the response being studied (and most of these gene sets are small), suggesting that the relevant genes are indeed picked by the method.

The surrogate task of predicting cell type was chosen in order to align the regulatory interactions with specific cell types which is the goal of the method.

Finally, regarding the question on few time points, these are typical for scRNA-seq data. Most datasets do not extend beyond a few time points, which is one of the ongoing challenges in the field. Regarding the suggestion to run the method independently for each time point, we note that, for the very first time point in each series, there is no influence from the future. Despite this, the overlap with the ground truth databases is still very significant as shown in the figures. Hence, treating Marlene as a static method should still recover relevant edges. It is unclear, however, if the model will still capture the dynamics of the system (i.e., the consecutive edge sets may not be related to each other).

We hope these address the reviewer concerns. Please let us know if you have additional questions.

[1] https://arxiv.org/abs/1902.10186

We thank the reviewer for praising the work. We hope to address their concerns below.

1. We apologize for the lack of clarity on this point. The connection to EvolveGCN is the dynamic adaptation of weights. In most neural networks, weights are updated via some gradient-based algorithm. In EvolveGCN and Marlene, some of the network weights are predicted by the model itself. In the case of Marlene, the GRU module predicts the weights of the linear layers in the self-attention module for the next time point. These weights are then multiplied with the input x to output keys, queries. We hope this clarifies the use of the GRU.
2. This is a good point. As is common in deep learning literature, neural networks are sensitive to initialization and this is also true for Marlene. However, we found that many TFs and genes were consistently selected across runs.
3. Good catch; this is not explained in the paper. During evaluation, we obtain $A_t$ for multiple batches and then average the outputs to obtain a single $A_t$ for each $t$. We have now clarified this in the main text.

We thank the reviewer for suggesting the additional BEELINE datasets. We looked at BEELINE datasets from GEO # GSE75748. However, these datasets are derived from a single cell type, hESCs. Currently, Marlene is designed to predict the cell type which we claim helps with the discovery of cell type specific GRNs. In the case of hESCs, there is only one cell type, so this data does not leverage the benefits of meta learning or Marlene for cell specific GRNs.

We hope this clarifies some of the reviewer's comments.