Recovering Time-Varying Networks From Single-Cell Data

## Abstract

Abstract Gene regulation is a dynamic process that underlies all aspects of human development, disease response, and other key biological processes. The reconstruction of temporal gene regulatory networks has conventionally relied on regression analysis, graphical models, or other types of relevance networks. With the large increase in time series single-cell data, new approaches are needed to address the unique scale and nature of this data for reconstructing such networks. Here, we develop a deep neural network, Marlene, to infer dynamic graphs from time series single-cell gene expression data. Marlene constructs directed gene networks using a self-attention mechanism where the weights evolve over time using recurrent units. By employing meta learning, the model is able to recover accurate temporal networks even for rare cell types. In addition, Marlene can identify gene interactions relevant to specific biological responses, including COVID-19 immune response, fibrosis, and aging.

## 1 Introduction

Biological systems are dynamic, changing over time in response to various stimuli and events. To construct accurate models of biological activity during development, disease progression, treatment response, and other biological processes, it is essential to track their evolution over time (Bar-Joseph et al., ). Studying the regulation of these dynamic processes is key for understanding the underlying mechanisms that drive the response and for identifying potential interventions that can serve as cures for diseases (Silverman et al., ).

Much of the research in this area is focused on the reconstruction of regulatory networks (Karlebach & Shamir, ; Badia-I-Mompel et al., ). These networks comprise a subset of proteins known as transcription factors (TFs), which regulate the activity of all other genes and proteins within the cell. However, these gene regulatory networks (GRNs) are not static. Instead, both the active nodes (proteins) and the edges (genes) change over time (Gasch et al., ; Luscombe et al., ). To reconstruct such networks, researchers often integrate static data—such as the type of nodes in the network—with dynamic data, such as time series measurements of node activity (gene expression profiles). Early work in this area employed microarrays and ChIP-chip data
(Lee et al., ; Wang et al., ; Blais & Dynlacht, ; Gilchrist et al., ) followed by time series
next-generation RNA-seq data (Wang et al., ), and
most recently, single-cell RNA-seq (scRNA-seq) data
(Ding et al., ; Nguyen et al., ; Matsumoto et al., ).

Several computational methods have been proposed over the last two decades to reconstruct such dynamic GRNs
(Bar-Joseph et al., ; Yosef et al., ; Schulz et al., ; Yan et al., ).
Some of these methods utilized
time-varying graphical models including Hidden Markov models, Markov random fields, and Dynamic Bayesian Networks
(Ahmed & Xing, ; Song et al., ; Dondelinger et al., ; Zhu & Wang, ). Other approaches attempted to use regression or to learn temporal precision matrices using extensions of the graphical lasso algorithm (Wang et al., ; Hallac et al., ).

While such models successfully reconstructed some processes (Ahmed & Xing, ; Kim et al., ), they are less suitable for the more recent types of data, most notably scRNA-seq time series data. First, the much larger size of the data presents a challenge
for traditional graphical models. In addition, prior methods do not directly account for the fact that multiple cells are profiled for each time point. Finally, prior methods do not utilize the ability to learn much larger models such as neural networks which have been shown to improve many other learning tasks (Ching et al., ; Angermueller et al., ).

Very recently, a few methods have been proposed for using deep learning to recover static GRNs (Shrivastava et al., ; Shu et al., ). However, they cannot be directly used to capture dynamic GRNs (i.e., enforcing learning between time points). Two recent approaches, Dictys and CellOracle (Wang et al., ; Kamimoto et al., ) can be used to infer dynamic GRNs, however, these methods rely on a different type of data (ATAC-seq) which is harder to obtain and less prevalent.

Beyond the realm of biology, the inference of dynamic graphs using neural networks has garnered significant attention. This problem has found applications in diverse domains, including information retrieval, molecular graphs, and traffic forecasting (Zhu et al., ; Zhang et al., ). While there are similarities between these problems and the dynamic GRN problem, there are also significant differences that make it hard to extend these methods for time series scRNA-seq. The problem of inferring temporal graphs is usually defined by recovering a series of graph adjacency matrices $\mathbf{A}_{t}\in\mathbb{R}^{n\times n}$ where $n$ is the number of nodes and each node is a $k$-dimensional feature vector. However, when dealing with scRNA-seq data, the problem becomes: given a gene expression matrix $\mathbf{X}_{t}\in\mathbb{R}^{c\times g}$, where $c$ is the number of cells (samples) and $g$ is the number of genes (features), we are interested in recovering gene networks $\mathbf{A}_{t}\in\mathbb{R}^{g\times g}$, i.e., graphs of features rather than nodes (cells).

In this paper, we present a novel deep learning framework that effectively addresses the challenges discussed above for reconstructing dynamic GRNs. Our contribution is three-fold. First, we demonstrate that existing deep learning methods for temporal graph structure learning can be adapted for scRNA-seq data analysis. To achieve this, we perform a gene featurization step by leveraging set-like architectures such as DeepSets or Set Transformers (Zaheer et al., ; Lee et al., ).
Second, we construct dynamic graphs by applying a self-attention mechanism (Bahdanau et al., ) to these gene feature vectors. To model dynamics, we draw inspiration from EvolveGCN where a gated recurrent unit (GRU) evolves the weights of a graph neural network (Pareja et al., ). However, unlike EvolveGCN, our approach uses a GRU to evolve the weights of key and value projection matrices in the self-attention module. This allows for the construction of dynamic graphs that capture regulatory interactions over time. Lastly, GRNs are highly dependent on cell functions, hence, separate GRNs need to be learned for each cell type. A single scRNA-seq dataset may combine cells of multiple types, some of which are rare cell populations. To this end, we employ a model-agnostic meta-learning (MAML) (Finn et al., ) training procedure by treating each cell type as a “task” to be learned. With this approach, the model quickly adapts to tasks with few samples, enabling the reconstruction of dynamic graphs even for rare cell types.

We apply our meta learning approach for inferring temporal gene regulatory networks (Marlene) to three publicly available scRNA-seq datasets. The first is a time series SARS-CoV-2 mRNA vaccination dataset of human peripheral blood mononuclear cells (PBMCs) (Zhang et al., ). The second dataset is a human lung aging atlas from the Human Cell Atlas Project (Regev et al., ; Sikkema et al., ). The third dataset is from a study of lung fibrosis using a mouse lung injury model (Strunz et al., ). All three datasets incorporate several time points, thus enabling a longitudinal analysis of the relevant biological responses through the inference of dynamic, cell type-specific GRNs. As we show, our method is able to reconstruct accurate networks for these datasets, significantly improving upon prior methods proposed for this task.

## 2 Methods

### 2.1 Problem Setup

Consider a gene expression matrix $\mathbf{X}\in\mathbb{R}^{c\times g}$ where $c$ is the number of cells and $g$ is the number of genes. In the human genome, $g$ varies from $25,000$ to $30,000$, while the number of cells could be between a couple thousand to a few million. In the setting of dynamic graphs, we assume the existence of a time point for each row (cell), leading to a time series $\widetilde{\mathbf{X}}:=\{\mathbf{X}_{1},\dots,\mathbf{X}_{T}\}$ with $\mathbf{X}_{t}\in\mathbb{R}^{c_{t}\times g}$. Here, the number of cells $c_{t}$ may vary with $t$. We are interested in recovering a series of directed graphs $\widetilde{\mathcal{G}}:=\{\mathcal{G}_{1},\dots,\mathcal{G}_{T}\}$ where each $\mathcal{G}_{t}=\{\mathcal{N},\mathcal{E}_{t}\}$. The set of nodes is the set of genes, i.e., $\mathcal{N}=[g]$, and we assume this set is static over time. The dynamic edge sets $\mathcal{E}_{t}=\{(u,v,w)\}_{u,v\in\mathcal{N},w\in\mathbb{R}}$ denote directed weighted links between genes, where the source $u$ is the gene that regulates the expression of its target $v$, and $w$ is the strength of this relationship. The source nodes are called transcription factor genes (TFs).

We can alternatively characterize each graph $\mathcal{G}_{t}$ by the corresponding adjacency matrix $\mathbf{A}_{t}\in\mathbb{R}^{g\times g}$. Denote $\widetilde{\mathbf{A}}:=\{\mathbf{A}_{1},\dots,\mathbf{A}_{T}\}$. Since TFs control the expression of their target genes, the underlying GRNs should, in principle, allow the recovery of the full expression profile for a cell. In other words, $\widetilde{\mathbf{X}}=f(\widetilde{\mathbf{X}}^{\text{TF}},\widetilde{\mathbf
{A}})$ where $\widetilde{\mathbf{X}}^{\text{TF}}$ denotes the expression of all TFs. The function $f$ is unknown as it involves intricate interactions among genes, including combinatorial effects. For instance, certain scenarios exist where TFs cooperate with each-other to activate a gene, while in other instances, the activation requires some TFs to be active and others to be inactive (Wise & Bar-Joseph, ).

In existing deep learning literature, $f$ is sometimes modeled using autoencoders (Seninge et al., ; Shu et al., ; Wang et al., ). However, reconstructing the full gene expression vector is challenging as the data is extremely sparse and conventional reconstruction losses, such as mean squared error, tend to emphasize overall averages.
Since GRNs are dependent on cell function, we hypothesize that simplifying the problem by predicting cell types may improve the accurate recovery of GRNs. In other words, given a temporal batch of cells of the same type $\widetilde{x}^{\text{TF}}\in\mathbb{R}^{T\times\text{batch size}\times\text{n.
TFs}}$, we consider the classification problem $y=f(\widetilde{x}^{\text{TF}},\widetilde{\mathbf{A}})$ where $y$ is the known cell type label for the batch. Finally, the task of learning $\widetilde{\mathbf{A}}$ given a batch $\widetilde{x}$ becomes

$$ $\operatorname*{arg\,min}_{\widetilde{\mathbf{A}}}\text{CrossEntropyLoss}(y,f( \widetilde{x}^{\text{TF}},h(\widetilde{x}))),\qquad\widetilde{\mathbf{A}}:=h( \widetilde{x})$ (1) $$

for choices of functions $f$ and $h$.

### 2.2 Architecture of Marlene

Figure: Figure 1: Overview of Marlene. Marlene takes as input gene expression data in the form of a cell-by-gene matrix. It then performs gene featurization via the pooling by multihead attention (PMA) mechanism which returns a gene feature matrix. This matrix is then inputted into a self-attention module to obtain a gene network in the form of an adjacency matrix. The weights of the self-attention module evolve from one time point to the next via a gated recurrent unit (GRU). The expression of transcription factors and the recovered graph are used to reconstruct the full gene expression vector. Finally, the reconstructed matrix is used to predict the cell type for the batch. The network is trained in a model-agnostic meta-learning fashion where each cell type is treated as a “task” to be learned, thus enabling the model to quickly adapt to cell types with low representation.
Refer to caption: x1.png

In this work, we propose a neural network architecture called Marlene that effectively learns dynamic GRNs (Figure 1). Marlene consists of three main steps. The first two steps address the choice of $h$ in (1), while the last addresses the choice of $f$.

In the first step, we apply a gene featurization step by treating a batch of cells as a set of elements. The DeepSet architecture introduces a pooling operator that allows the neural network to be invariant to the order of input samples, effectively treating the input as a set (Zaheer et al., ). Similarly, the Set Transformer architecture is designed to process sets of data via attention-based operators that are permutation invariant (Lee et al., ). Specifically, the pooling by multihead attention (PMA) aggregation scheme introduced in Set Transformers outputs a matrix of $k$ vectors $\mathbf{H}\in\mathbb{R}^{k\times g}$ for an arbitrary input set $\mathbf{X}\in\mathbb{R}^{c\times g}$. Each of the $k$ output vectors has a specific meaning such as statistics of the input data. By applying the PMA operator to a temporal batch of cells, we obtain a gene-by-feature matrix $\mathbf{G}=\mathbf{H}^{\top}\in\mathbb{R}^{g\times k}$ that encodes information about the cells from each time point. PMA consists of a multihead attention block (MAB) where the input $\mathbf{X}$ consists of key vectors, and the query is a learnable set of $k$ vectors $\mathbf{S}$. In this work, we use a shared PMA layer for all time points assuming that the specific key statistical properties are invariant (though their value obviously changes for different time points). Given $\widetilde{x}=[x_{1},\dots,x_{T}]$, we have

$$ $\widetilde{\mathbf{G}}=\text{PMA}(\widetilde{x})^{\top}:=\text{MAB}(\mathbf{S} ,\widetilde{x})^{\top}:=(\mathbf{M}+\text{rFF}(\mathbf{M}))^{\top}$ (2) $$

where $\mathbf{M}=\mathbf{S}+\text{Multihead}(\mathbf{S},\widetilde{x},\widetilde{x})$ and rFF is a row-wise feedforward layer. For completeness, these operations are defined in the appendix.

Note that if cells of different types are mixed in the same batch, the statistics derived by the PMA step may not capture cell type-specific information. Consequently, in a single batch, we only include cells of one type.

In the second step, Marlene learns temporal adjacency matrices using a self-attention mechanism. To model dynamics, we draw inspiration from EvolveGCN which performs model adaptation using a GRU (Pareja et al., ). Unlike EvolveGCN, which uses the GRU to update the weights of a graph convolution layer, we use a GRU to evolve the key and query projection weights of the self-attention module. Since most time series we deal with contain very few time points, a GRU should suffice and not suffer from vanishing gradient problems. Similar to EvolveGCN, we apply a summarization step via top $k$ pooling to reduce the gene feature matrix to a square matrix for the GRU (appendix).

More precisely, we initialize self-attention weights $\mathbf{W}^{Q}_{0},\mathbf{W}^{K}_{0}\in\mathbb{R}^{k\times k}$ and two recurrent units $\text{GRU}_{Q},\text{GRU}_{K}$. Given the time sequence of gene feature matrices $\widetilde{\mathbf{G}}$ obtained from the previous step, temporal adjacency matrices are constructed in the following recurrent fashion for all $t\in[T]$:

$$ $\displaystyle\mathbf{Z}_{t}=$ $\displaystyle\text{TopK}(\mathbf{G}_{t})$ (3) $\displaystyle\mathbf{W}^{Q}_{t}=\text{GRU}_{Q}(\mathbf{Z}_{t},\mathbf{W}^{Q}_{ t-1}),$ $\displaystyle\quad\mathbf{W}^{K}_{t}=\text{GRU}_{K}(\mathbf{Z}_{t},\mathbf{W}^ {K}_{t-1})$ (4) $\displaystyle\mathbf{Q}_{t}=\mathbf{G}_{t}\mathbf{W}^{Q}_{t},$ $\displaystyle\quad\mathbf{K}_{t}=\mathbf{G}_{t}\mathbf{W}^{K}_{t}$ (5) $\displaystyle\mathbf{A}_{t}=\text{softmax}$ $\displaystyle\left(\frac{\mathbf{Q}_{t}\mathbf{K}_{t}^{\top}}{\sqrt{k}}\right).$ (6) $$

Here, $\mathbf{W}^{Q}_{t}$ and $\mathbf{W}^{K}_{t}$ serve as hidden states for the respective GRUs. The GRUs dynamically adapt self-attention weights, influencing which TFs specific genes should attend to in subsequent time steps. Consequently, the evolution of these weights is constrained. We also restrict the columns of $\mathbf{A}_{t}$ (i.e., sources) to $p$ known TFs in the TRRUST database (Han et al., ) which greatly reduces the number of parameters to be learned. Therefore, in our implementation $\mathbf{A}_{t}\in\mathbb{R}^{g\times p}$.

Next, we perform a gene expression reconstruction step based on the expression of TFs and the inferred adjacency matrices. This is followed by any number of fully connected layers with nonlinear activation functions $\sigma$. Finally, we sum across output vectors to obtain a logit vector with the same dimension as the number of cell types in the data:

$$ $\tilde{y}=\text{Pool}(\text{Linear}(\dots\sigma(\text{Linear}(\widetilde{x}^{ \text{TF}}\widetilde{\mathbf{A}}^{\top})))).$ (7) $$

Network depth can be introduced at all three levels by stacking MAB layers during gene featurization, stacking GRUs, or stacking linear layers at the end.

### 2.3 Meta learning for rare cell types

ScRNA-seq datasets often originate from biological samples that exhibit cellular heterogeneity, potentially containing multiple distinct cell types. Some of these cell subpopulations are rare and are represented by a small number of cells in the sample (Jindal et al., ). Since we are concerned with the discovery of cell type-specific temporal GRNs, learning such large graphs for these rare cell types may not be feasible and lead to overfitting. Since many interactions are shared across cell types (Chasman & Roy, ), we employ the model-agnostic meta-learning framework (MAML) (Finn et al., ). MAML is specifically designed to enable neural networks to adapt to novel tasks with limited training samples (i.e., few shot learning). By treating each cell type as a “task”, the MAML training paradigm facilitates the recovery of dynamic graphs for rare cell types. We begin by adapting model parameters through multiple optimization steps using a batch of support examples (cells). These adapted parameters are then evaluated on a separate set of query cells, followed by a meta-update.

During the adaptation step, we perform gradient descent, while for the meta-update, we employ the Adam optimizer (Kingma & Ba, ). During training, gradient clipping proves crucial to prevent overfitting of the MAML adaptation step to the cell type under consideration.

## 3 Experiments

**Table 1: Time series scRNA-seq datasets used in this study.**
|  | Number of | Metadata |  |  |  |  |
| --- | --- | --- | --- | --- | --- | --- |
| Dataset | Cells | Genes^1 | TFs | Cell Types | Time Points | Sample |
| SARS-CoV-2 | $113,271$ | $1899$ | 556 | 7 | d0, d2, d10, d28 | PBMCs (human) |
| HLCA | $27,953$^2 | $2433$ | 674 | 11 | Ages $<35,35-50,\geq 50$ | lung (human) |
| Fibrosis | $22,758$ | $1217$ | 433 | 6 | PBS, d3, d7, d10, d14, d21, d28 | lung (mouse) |

To validate our approach, we use three public scRNA-seq datasets (Table 1): a human SARS-CoV-2 mRNA vaccination dataset, a lung aging atlas (The Human Lung Cell Atlas—HLCA), and a mouse lung fibrosis dataset (Zhang et al., ; Sikkema et al., ; Strunz et al., ). To assess the quality of the inferred networks, we draw upon two databases of TF-gene regulatory interactions, which have been curated from the scientific literature—TRRUST and RegNetwork (Han et al., ; Liu et al., ). For the human genome, TRRUST contains 8427 unique validated regulatory edges, while RegNetwork contains $150,405$. Note that certain edges lack a corresponding TF or gene in the expression data, so the numbers used for the analysis are smaller. We used only the genes that were present in TRRUST for all three datasets. For the mouse lung dataset we used the corresponding mouse networks for both databases. To match the number of links in these databases, we selected the top 2% of edges for all methods. For Marlene, this was done by sparsifying the self-attention matrix to retain only the top scoring edges. The significance of the overlap was carried out via Fisher’s exact test (Fisher, ). All $p$-values were corrected for multiple testing using the Benjamini-Hochberg procedure (Benjamini & Hochberg, ).

We compare Marlene against several popular static gene regulatory network inference methods included in the BEELINE benchmark (Pratapa et al., ) and beyond, such as PIDC, GENIE3, GRNBoost2, SCODE, and DeepSEM (Chan et al., ; Huynh-Thu et al., ; Moerman et al., ; Matsumoto et al., ; Shu et al., ), which are applied independently to each time point. DeepSEM is a deep generative model based on structural equation modeling. We also compare against time-varying graphical lasso (TVGL)(^1^11We used the implementation of [https://github.com/fdtomasi/regain](https://github.com/fdtomasi/regain)), a method that models temporal precision matrices (Hallac et al., ), and to a deep neural network that utilizes the S4 module (GraphS4mer) (Tang et al., ; Gu et al., ).

We train Marlene using a batch size of 16 cells and also set the number of seeds in the PMA layer to 16. For MAML, we use 5 inner steps. The model with the lowest loss is selected for GRN inference. For the meta-update, we use a decaying learning rate starting with ${1}\mathrm{e}{-4}$, while for the inner step we use ${1}\mathrm{e}{-3}$ for both datasets. Experiments were performed using an NVIDIA RTX 3060 and took only a few minutes per run.

### 3.1 Case study 1: SARS-CoV-2 vaccination

The SARS-CoV-2 vaccination dataset consists of peripheral blood mononuclear cells (PBMCs) from six healthy donors at four time points (days 0, 2, 10, and 28) (Zhang et al., ). Day 0 samples were obtained before vaccination. We removed the “Other” cell type group and kept the remaining seven. These include B cells, dendritic cells (DC), monocytes (Mono), natural killer cells (NK), and various types of T cells.

#### 3.1.1 Marlene recovers accurate gene regulatory networks

Analysis results using Marlene and prior methods is presented in Figure 2. As can be seen, Marlene outperformed competing methods in the dynamic GRN inference task for 5 of the cell types, yielding statistically significant results across time points. Specifically, for B cells, Marlene successfully identified more than 800 regulatory links within the RegNetwork framework at each time point (FDR $\leq{1}\mathrm{e}{-67}$), surpassing the performance of the second-best method, SCODE, which detected 579 links (day 2, FDR $\leq{1}\mathrm{e}{-15}$). Analogous findings were observed for natural killer cells, where Marlene identified over 600 RegNetwork links at each time point (FDR $\leq{1}\mathrm{e}{-18}$). In comparison, the second-ranking method, SCODE, showed a significant overlap for only one time point (day 28). Upon examining the TRRUST database, we observed less pronounced differences in the results. Nonetheless, Marlene obtained higher overlap for 5 out of 7 cell types followed by GENIE3, which performed well for monocytes and the “Other T” category.

Figure: Figure 2: Overlap analysis of the SARS-CoV-2 vaccination dataset. Showing $-\log_{10}(\text{FDR})$ values from a Fisher’s exact test measuring the overlap between predicted TF-gene interactions in reconstructed networks and two TF-gene interaction databases (TRRUST, RegNetwork). Best performing method is starred.
Refer to caption: x2.png

#### 3.1.2 Marlene recovers realistic dynamic transitions

The analysis so far has primarily focused on individual time points. Next, we turned our attention to assessing the quality of graph transitions between consecutive time points. Specifically, we examined whether the learned graphs demonstrated smooth transitions over time. To evaluate this, we computed the intersection-over-union (IoU) score for edges between time points $t$ and $t+1$ (Figure 3a). Notably, our findings revealed that for most cell types, Marlene exhibited the lowest IoU score during the initial period (days $0\to 2$), followed by higher scores during days $2\to 10$, and $10\to 28$. This pattern aligns with our expectations, as variations in gene expression are likely to be most pronounced during the early post-vaccination period (days $0\to 2$). From the methods we compare against, GENIE3, GRNBoost2, SCODE, PIDC, and DeepSEM exhibited significantly lower IoU scores across all temporal transitions, likely due to their lack of dynamic modeling and the fact that they were ran independently per time point. TVGL, on the other hand, showed high IoU scores which remained close to constant over time. Finally, Graphs4mer displayed a reduction in IoU scores over time, which is unlikely given the expected immediate immune response.

Figure: Figure 3: Temporal analysis of the predicted gene regulatory networks for the SARS-CoV-2 vaccine dataset. (a) Intersection-over-union (IoU) scores between consecutive graphs. (b) For each method, top 3 MSigDB terms enriched for genes that were regulated at day 2 but not day 0.
Refer to caption: x3.png

Next, we sought to assess the quality of the TF-gene regulatory links added between time points. For brevity, we focused specifically on the initial temporal transition (days $0\to 2$), as this period is likely to witness a more significant biological response. For each cell type, we took note of all the genes that were regulated by some TF at day 2 but not at day 0. Using this set of genes $z$, we performed gene set enrichment analysis (GSEA) (Subramanian et al., ; Fang et al., ) using the molecular signatures database (MSigDB) (Liberzon et al., ). Through permutation tests, GSEA assigns an enrichment score (ES) to $z$ reflecting its overrepresentation within the MSigDB gene set collection. We found that for many cell types, genes added by Marlene at day 2 greatly overlapped with COVID-19 and SARS-CoV-2 related gene sets. For instance “Interferon Gamma Response”, which was identified as a SARS-CoV-2 antiviral response (Hilligan et al., ), was significantly enriched in dendritic cells (15 genes, FDR $={1}\mathrm{e}{-6}$). Similarly, “TNF-alpha Signaling via NF-kB”—a pathway involved in the immune response and inflammation (Hayden & Ghosh, )—was enriched in several cell types, as well as processes such as “Apoptosis” (cell death) and “p53 Pathway” (inhibits replication of infected cells) (Elmore, ; Harris & Levine, ). Other methods, while being enriched for relevant terms, showed a smaller gene overlap for these types (Figure 3b) or were not consistent across cell types (e.g., DeepSEM, SCODE).

Overall, these results suggest that Marlene is able to capture both known TF-gene links, but also genes that are relevant to the response being studied.

### 3.2 Case study 2: Aging and senescence in the lung

The Human Lung Cell Atlas (HLCA) is a large data integration effort by the Human Cell Atlas Project (Sikkema et al., ; Regev et al., ). This data combines scRNA-seq samples from 107 individuals spanning an age range of 10 to 76 years, making it particularly attractive for studying aging and senescence (a form of aging characterized by the absence of cell division) (van Deursen, ; SenNet Consortium, ).

Figure: Figure 4: Results on the HLCA dataset. (a) FDR corrected $p$-values of Fisher exact tests reflecting the number of links that overlap with TRRUST and RegNetwork databases. (b) Top 3 Jensen Diseases terms enriched for genes added between the first and second age group.
Refer to caption: x4.png

We split the atlas into three age groups at 35 and 50 years old, thus forming a pseudotime series of length 3. We removed smokers from the dataset as these will likely confound the results. To accommodate the data in the GPU, we randomly selected cells from 11 cell types, including type II pneumocytes, endothelial cells, and monocytes.

Similar to the vaccination dataset, we begin the analysis by evaluating the set of regulatory links using the TRRUST and RegNetwork databases. For this dataset we find that Marlene and SCODE are the top two performing methods (Figure 4a). For some of the cell types, Marlene achieves significant results, recovering more than 1000 RegNetwork links (classical monocytes, FDR $={1}\mathrm{e}{-76}$). Even for cell types with fewer cells, such as non-classical monocytes (with only 138 cells for the second age group), Marlene still recovered more than 800 known TF-gene links for each transition (FDR $\leq{1}\mathrm{e}{-27}$). SCODE performed well for some cell types such as CD1c-positive myeloid dendritic cells and CD4-positive, alpha-beta T cells. For all other methods, the overlap was smaller (Appendix Figure 6a). Note that while SCODE is comparable for the static network (single time point) inference task, it does not utilize dynamic information.

We next examined the ability of different methods to capture the dynamics of the biological processes. For this, we looked at graph transitions. IoU scores show that only the temporal methods (Marlene, TVGL, Graphs4mer) capture the smooth temporal transition between time points, while other methods, including SCODE, achieve low IoU scores (Appendix Figure 6b). We performed GSEA using Jensen Diseases gene set to see if genes added by Marlene in these transitions were enriched for any age-related diseases (Pletscher-Frankild et al., ; Grissa et al., ). We found that Marlene added genes are enriched for several diseases such as arthritis, lung disease, and coronary artery disease. Other dynamic baselines were also enriched for relevant terms, but contained fewer marker genes (Figure 4b).

Finally, we also investigated whether the genes regulated at different age groups were enriched for senescence. Cellular senescence refers to a permanent arrest of cell division triggered by the accumulation of DNA damage (Suryadevara et al., ). The absence of cell division can detrimentally impact tissue regeneration and repair, thereby contributing to various age-related diseases. Here, we use the SenMayo gene set which contains 125 genes reported to be enriched for senescence (Saul et al., ). Only 81 of these genes overlapped with our data. We found that for 4 cell types, there was an increase in SenMayo gene regulation at the oldest age group (age $>50$), suggesting that senescent cells accumulate with age as hypothesized (Figure 5).

Figure: Figure 5: Enrichment for senescence using the SenMayo set. For 4 cell types, there was statistically significant enrichment for the oldest age group. We only used the top 200 regulated genes.
Refer to caption: x5.png

### 3.3 Case study 3: Fibrosis in a mouse lung injury model

Next, we evaluated whether Marlene could perform effectively across different species by analyzing a dataset from a mouse model of lung injury induced by the chemotherapeutic agent bleomycin (Strunz et al., ). The dataset included seven time points: one pre-treatment and six post-treatment intervals. After filtering out cell types with low representation and genes with low counts, we retained six cell types, including B cells, T cells, and macrophages.

In this analysis, Marlene outperformed competing methods in four of the six cell types when benchmarked against the RegNetwork database, specifically in alveolar epithelial cells, dendritic cells, endothelial cells, and macrophages. For T cells, TVGL showed slightly better performance. When evaluated against the TRRUST database, SCODE performed well in four cell types, while Marlene surpassed it in the remaining two. The differing results between two databases may reflect their incomplete coverage, highlighting the need for further refinement.

Finally, all static baselines, including SCODE, showed low IoU scores across time points, indicating their inability to capture temporal evolution. In contrast, Marlene, showed increasing IoU scores over time, suggesting ongoing lung regeneration following treatment which slowly stabilizes. Figures illustrating these findings are provided in the appendix.

## 4 Discussion

Gene regulation is a dynamic process that underlies all biological systems and responses. Understanding which TFs regulate which genes, and when this regulation occurs, provides insights into these dynamic processes which can ultimately lead to better treatment options.

To improve on current methods for reconstructing time varying regulatory networks, we use the expressive capabilities of deep neural networks to model the dynamic regulation of genes. Specifically, we focused on inferring dynamic networks from scRNA-seq data.

Our proposed method, Marlene, constructs dynamic graphs from time series data. Marlene begins with a set pooling operator based on PMA to extract gene features. These gene features are subsequently used to construct dynamic graphs via a self-attention mechanism. The weights of the self-attention block are updated through the use of GRUs. Additionally, by employing MAML, we help Marlene uncover graphs even for rare cell types However, Marlene optimizes the prediction of cell type label rather than gene expression. As such, Marlene is not currently equipped to determine the impact of perturbations including gene knockouts or overexpression experiments. Exploring the integration of causal inference capabilities into Marlene represents a promising direction for future research.

We demonstrated the effectiveness of Marlene in recovering dynamic GRNs using three datasets: a SARS-CoV-2 vaccination dataset, a lung aging atlas, and a mouse dataset of fibrosis. In all three datasets, Marlene successfully identified many validated TF-gene links from the TRRUST and RegNetwork databases across various cell types. It also accurately modeled the temporal dynamics of these connections. Some prior methods ignored the temporal aspect, leading to very little similarity between consecutive time networks. Other methods integrated all time points together, leading to very similar networks for each time point. In contrast, Marlene accurately recovered the variation dynamics, which is often characterized by strong rewiring following treatment that later stabilizes. In addition, Marlene identified many relevant edges. For instance, in the lung aging data, several dynamic edges were enriched for age-related diseases, such as arthritis. Meanwhile, in the SARS-CoV-2 data, these dynamic links were enriched for immune response processes. Prior methods captured some known edges, however, the overall results were less significant. By providing better models to explain disease and vaccine response, researchers can zoom in on the specific mechanisms targeted which in turn can lead to better treatments.

## 5 Limitations

While successful, Marlene has a few limitations. The datasets we used in this study, while typical for scRNA-seq time series, consisted of only a few time points. For longer sequences, the GRU operation may suffer from vanishing gradient problems (Pascanu et al., ). In such scenarios, the S4 module may be preferred as it has been shown to model long sequences better than traditional GRUs (Gu et al., ). In addition, using a large number of genes for training, results in quadratic growth in memory consumption due to the need to store adjacency matrices. This led us to restrict the set of genes for each of the two studies. A more efficient implementation or alternative approaches such as FlashAttention (Dao et al., ) can lead to better ability to utilize all genes profiled.

## Appendix A Appendix

### A.1 Set Transformer Operations

We redefine the Multihead and rFF operations from Set Transformers (Lee et al., ) to the ones used for Marlene here.

First, we define the Attention operation. Let $Q\in\mathbb{R}^{k\times g}$ be the query matrix of $k$ elements and $g$ dimensions. The Attention operation used for MAB is

$$ $\text{Attention}(Q,K,V)=\text{softmax}\left(\frac{QK^{\top}}{\sqrt{g}}\right)V$ (8) $$

where the key and value matrices are $K,V\in\mathbb{R}^{c\times g}$. Next, the Multihead attention operation with $h$ heads (Vaswani et al., ) is given by

$$ $\text{Multihead}(Q,K,V)=\text{concat}(O_{1},\dots,O_{h})W^{O}$ (9) $$

where $O_{j}=\text{Attention}(QW^{Q}_{j},KW^{K}_{j},VW^{V}_{j})$ for weight matrices $W^{Q}_{j},W^{K}_{j},W^{V}_{j}\in\mathbb{R}^{g\times g/h}$ and $W^{O}\in\mathbb{R}^{g\times g}$ (these matrices are not to be confused with self-attention weights used consequently for Marlene). In our implementation, $k$ is the number of seeds or output vectors used for the PMA layer. This is a hyperparameter that corresponds to the number of “statistic” vectors we expect to learn from data. Finally, rFF is a feedfoward layer such as a linear layer.

### A.2 EvolveGCN Operations

Here, we introduce the GRU and topK pooling operations used in the second step of Marlene.

The topK pooling operation is needed to summarize nodes into $k$ representative ones (Cangea et al., ; Pareja et al., ). Here $k$ is the same as the number of seeds used for PMA. Given an input $\mathbf{G}\in\mathbb{R}^{g\times k}$ and a learnable vector $q$, the TopK operation performs the following steps:

$$ $\displaystyle\rho$ $\displaystyle=\frac{\mathbf{G}q}{\|q\|}$ $\displaystyle i$ $\displaystyle=\text{Top-k-indices}(\rho)$ $\displaystyle\mathbf{Z}$ $\displaystyle=[\mathbf{G}\odot\tanh{(\rho)}]_{i}.$ $$

At time step $t$, given a pooled matrix $\mathbf{Z}_{t}$ and hidden state $\mathbf{W}_{t-1}$ (i.e., self-attention weights $\mathbf{W}^{Q}_{t-1}$ or $\mathbf{W}^{K}_{t-1}$), the standard GRU operation is:

$$ $\displaystyle r_{t}$ $\displaystyle=\sigma(M_{ir}\mathbf{Z}_{t}+b_{ir}+M_{hr}\mathbf{W}_{t-1}+b_{hr})$ $\displaystyle z_{t}$ $\displaystyle=\sigma(M_{iz}\mathbf{Z}_{t}+b_{iz}+M_{hz}\mathbf{W}_{t-1}+b_{hz})$ $\displaystyle n_{t}$ $\displaystyle=\text{tanh}(M_{in}\mathbf{Z}_{t}+b_{in}+r_{t}\odot(M_{hn}\mathbf {W}_{t-1}+b_{hn}))$ $\displaystyle\mathbf{W}_{t}$ $\displaystyle=(1-z_{t})\odot n_{t}+z_{t}\odot\mathbf{W}_{t-1}$ $$

where $\sigma$ is the sigmoid function and $\odot$ is the Hadamard product. See also Paszke et al. ().

### A.3 Supplementary Figure for the HLCA dataset

Figure: Figure 6: (a) FDR corrected $p$-values of Fisher exact tests reflecting the number of links that overlap with the two TF-gene databases. (b) IoU scores across time reflecting the overlap between consecutive graphs.
Refer to caption: x6.png

### A.4 Analysis of the fibrosis dataset

Figure: Figure 7: (a) FDR corrected $p$-values of Fisher exact tests reflecting the number of links that overlap with the two mouse databases. (b) IoU scores across time reflecting the overlap between consecutive graphs.
Refer to caption: x7.png

## Peer Reviews

confidence: 4
contribution: 2
initial_rating: 3
presentation: 3
questions: The model is aimed at time series, but the experimental data very few data points. If the method was applied to each time point separately (as just the first time point, eliminating the recurrent part), would the performance be affected substantially?
soundness: 2
strengths: The paper focuses on an important, though well-studied task of regulatory gene network inference. It focuses on time-series, single-cell data, an increasing available, more detailed view of gene expression. 

The approach goes beyond simple application of existing deep learning models by using an architecture in which the projection matrices for calculating attention are evolve as part of an RNN.
summary: The paper aims to predict a series of graphs describing gene expression regulation by transcription factors from a series of single-cell gene expression data. The method uses attention-based architecture within each time point, with a recurrent component linking the time points.
weaknesses: The architecture relies on interpreting the attention matrix A generated within the model as the adjacency matrix of the regulatory network, with the model itself being trained on a surrogate task of predicting cell type (y). The assumption that A, used in this surrogate task, will capture direct regulatory interactions is not very well justified in the manuscript. Would using a matrix with different dimensionality (e.g. having #columns the same, to match # of TFs, but with different number of rows), result in similar performance on the surrogate task? If yes, what beyond shape leads to the interpretation of A as the adjacency matrix? For example, is there a reason to apply softmax to rows of the matrix?

The experimental results are missing key details relating to the performance on regulatory network inference. Key statistics related to network-wide performance in discovering edges are not reported: it would be helpful to see AUROC and AUPR values.

---

confidence: 4
contribution: 4
initial_rating: 8
presentation: 3
questions: 1. BEELINE does have a few dataset that has multiple time points (For example, hESC). I wonder how Marlene perform on those datasets in BEELINE. You can use the non-chipseq ground truth, which is very similar with the method I use. It would be great if you can use their metric (EPR/AUPRR etc)
soundness: 4
strengths: This is a very solid piece of research. The motivation is well explained and attractive. The idea is novel and has potential to be practically useful. The description of the method is very clear and the logic flows very well. The experiment part is comprehensive. 

In terms of the method itself, using PMA or Set transformer to convert the expression to gene feature matrix eliminate the axis of cells so downstream analysis could focus on the attentions on genes. The use of GRU in the next step is not so intuitive but seems to have literature support. The task of predicting cell labels is also very clever in this case because single cell data is noisy and complete reconstruction is more prone to error. Also, control the sparsities of the adjacency matrix by using the top k edges is also very inspiring.
summary: This paper introduces an attention and graph based approach called Marlene, to infer dynamic GRN from time series (longitudinal) single cell RNA sequencing data. The author provided an in-depth overview of the GRN inference field and explained the need of developing novel methods to capture the time-varying gene networks. 

Marlene starts with a pooling by multihead attention (PMA) layer which transform gene expression profiles at multiple time points to gene feature matrix at multiple time points. Then self attention are applied on each of these gene embeddings and the adjacency matrix is calculated based on the learned transformation. The weights of the self attention module are connected with GRU units so that information could be shared across time points. Finally, the extracted GRN are used to regenerate the expressed data, and the generated features are used to predict the cell type label. The entire model is trained to predict cell type. 

In term of experiment, the author managed to squeeze results from 3 experiments into this paper. These 3 cases are SARS-CoV-2 Vaccination, Aging and Lung, and fibrosis a mouse lung injure model. All 3 experiments are solid and supporting the claim.
weaknesses: 1. I would like to see a more clear explanation on Equation 3-6. What exactly are the rational of using GRU here beyond it was used in EvolveGCN? Do we have any physical meaning on this operation on Equation 3-6? 
2. Algorithm stability is a key metric in BEELINE. Could you comment on the stability of Marlene?
3. The output of Equation 6 is the adjacency matrix at timepoint t. Then, at least for the evaluation you have performed in this study, you must have transformed multiple At into one At. Could you explain in detail how you did that?

---

confidence: 4
contribution: 2
initial_rating: 3
presentation: 4
questions: 1	How sensitive is Marlene to the choice of hyperparameters, such as the number of attention heads, hidden units, or depth of the neural network layers?
Explanation: The performance of deep learning models often depends on selecting hyperparameters. The paper does not provide a comprehensive analysis of how different hyperparameter settings affect Marlene's ability to recover accurate GRNs. Understanding the model's sensitivity to these choices is important for assessing its robustness and guiding practical applications.

	2	Can Marlene effectively handle datasets with a large number of time points or with irregular time intervals between samples?
Explanation: The study's datasets contain a relatively small number of time points (3-7). It is unclear how well the model would scale to longer time series, which are common in many biological processes. Additionally, the paper does not discuss how Marlene would handle irregularly sampled data, where the time intervals between consecutive points vary. Addressing these scenarios is crucial for the model's general applicability.

	3	How does the model's performance change when dealing with datasets of different sizes, both in terms of the number of cells and the number of genes?
Explanation: The paper reports results on three specific datasets but does not provide a systematic analysis of how the model's performance scales with the size of the input data. Understanding Marlene's data efficiency and ability to handle datasets of varying sizes is important for assessing its practical utility and guiding data collection efforts.

	4	How does Marlene handle technical noise and batch effects that are common in scRNA-seq data?
Explanation: Single-cell RNA sequencing data often contains technical noise and batch effects that can confound the analysis of biological variation. The paper does not explicitly discuss how Marlene deals with these issues or if any preprocessing steps (e.g., normalization, batch correction) were applied to the input data. Clarifying the model's robustness to these factors is vital for its reliable application to diverse datasets.

	5	Can the model provide insights into the strength and directionality of the inferred regulatory interactions?
Explanation: While Marlene outputs weighted directed graphs, the paper focuses primarily on evaluating the presence or absence of edges against existing interaction databases. It does not delve into how well the model captures the strength and directionality of the regulatory relationships. Providing a more detailed analysis of these aspects could enhance the interpretability and biological relevance of the inferred networks.

6    What are the downstream applications of time-varying GRNs? Currently, the authors evaluate according to the recovery of (static) curated interactions. How can the method be evaluated in a way that evaluates its utility for real downstream applications?

7   The paper shows an extremely large difference between Marlene and alternative methods. This is surprising because all methods perform largely the same tasks with only minor methodological tweaks. What can explain the large differences between models?
soundness: 2
strengths: Strengths:
* Marlene effectively leverages recent advances in deep learning, such as self-attention mechanisms and recurrent units, to model dynamic GRNs from scRNA-seq data. The approach uses set-based architectures to handle multiple cells per time point.
* Employing meta-learning (MAML) enables Marlene to reconstruct accurate networks even for rare cell types by treating cell types as tasks. This enhances the model's ability to handle heterogeneous cell populations.
* The model demonstrates strong empirical results on three diverse datasets, outperforming several static and temporal baselines in recovering known regulatory interactions and identifying relevant biological processes.
summary: This paper proposes a novel deep learning framework called Marlene for reconstructing dynamic gene regulatory networks (GRNs) from time-series single-cell RNA sequencing (scRNA-seq) data. The authors aim to address the challenges of modeling the temporal evolution of GRNs using the increasingly available scRNA-seq data, which existing network inference methods are not well-equipped to handle.

Problem Significance: 
Studying the dynamic regulation of biological processes is crucial for understanding the mechanisms driving responses such as development, disease progression, and treatment outcomes. Reconstructing accurate temporal models of GRNs can provide key insights to identify potential interventions and treatments. However, the increasing scale of scRNA-seq data and the presence of multiple cells profiled per time point pose challenges for traditional GRN inference approaches. The development of Marlene tackles an important open problem in the field.

Model Input/Output: 
The input to Marlene is a time series of gene expression matrices {X^1, ..., X^T}, where each matrix X^t has dimensions (cells x genes) for a given time point t. The number of cells may vary per time point.
The output is a series of directed weighted graphs {G^1, ..., G^T} representing the GRN at each time point. The nodes are genes (assumed to be the same across time) and the edges capture regulatory relationships between transcription factors (TFs) and target genes. The graphs are characterized by adjacency matrices {A^1, ..., A^T}.
For evaluation, Marlene's inferred networks are benchmarked against curated TF-gene interaction databases (TRRUST and RegNetwork) using Fisher's exact test to assess overlap significance. The dynamics captured between time points are analyzed by computing intersection-over-union (IoU) scores.
weaknesses: Weaknesses:
* The paper primarily evaluates using overlap analysis with existing incomplete databases of static interactions. More direct experimental validation of novel predicted regulatory links would strengthen the findings. 
* Potential limitations in scaling Marlene to a very large number of genes are not thoroughly discussed. In experiments, the quadratic memory usage from adjacency matrices led to gene filtering.
* The model currently lacks the ability to predict the effects of perturbations like TF knockouts, which could enhance its utility for causal inference and treatment design.
* The usefulness of time-varying GRNs is not discussed.

Overall, the paper makes a valuable methodological contribution to a significant problem and demonstrates promising results. However, there is room for further validation and discussion of scalability limitations. More concrete insights into Marlene's novel interactions and their biological significance would also improve the impact. Overall, this is a solid paper that advances the field of GRN inference from single-cell data.

---

confidence: 4
contribution: 3
initial_rating: 3
presentation: 2
questions: 1. While the paper presents a technical approach, it could benefit from a deeper discussion on the biological implications of the findings and how they align with or differ from current scientific understanding.
2. In the experimental phase, mainly presents results based on metrics, demonstrating that it outperforms other methods. However, for solving a specific GRN problem, we are more concerned with whether the inferred GRN can identify some key genes or transcription factors (TFs) that can be further analyzed downstream. For a biological problem, merely comparing metrics does not adequately demonstrate the model’s performance.
3. The paper emphasizes using GRU to capture the temporal dynamics of the GRN. However, there are already many studies on dynamic GRNs, such as Dictys: dynamic gene regulatory network dissects developmental continuum with single-cell multiomics. I believe it would be more convincing to compare with these dynamic methods.
4. In GRN inference problems, perturbation experiments for some key genes are also an important downstream analysis. Including some perturbation experiments could help validate the accuracy of the inferred GRN.
soundness: 3
strengths: 1. propose a novel deep learning framework that employs self-attention mechanisms and GRUs to model the dynamics of gene regulatory networks.
2. The method is tested across three diverse datasets, which demonstrate the generalization ability of the model.
summary: The paper titled "Recovering Time-Varying Networks from Single-Cell Data" introduces Marlene, a deep neural network designed to infer dynamic gene regulatory networks from time series single-cell gene expression data. The authors propose a model that employs self-attention mechanisms and recurrent units to construct directed gene networks, with weights that evolve over time. A significant contribution is the integration of meta-learning to enable accurate recovery of temporal networks, even for rare cell types. The study demonstrates Marlene's effectiveness in identifying gene interactions relevant to specific biological responses across multiple datasets, showcasing the potential of deep learning in unraveling the dynamics of gene regulation and its implications for understanding complex biological processes.
weaknesses: 1. While the paper presents a technical approach, it could benefit from a deeper discussion on the biological implications of the findings and how they align with or differ from current scientific understanding.
2. In the experimental phase, mainly presents results based on metrics, demonstrating that it outperforms other methods. However, for solving a specific GRN problem, we are more concerned with whether the inferred GRN can identify some key genes or transcription factors (TFs) that can be further analyzed downstream. For a biological problem, merely comparing metrics does not adequately demonstrate the model’s performance.
3. The paper emphasizes using GRU to capture the temporal dynamics of the GRN. However, there are already many studies on dynamic GRNs, such as Dictys: dynamic gene regulatory network dissects developmental continuum with single-cell multiomics. I believe it would be more convincing to compare with these dynamic methods.
4. In GRN inference problems, perturbation experiments for some key genes are also an important downstream analysis. Including some perturbation experiments could help validate the accuracy of the inferred GRN.