Pruning neural network models for gene regulatory dynamics using data and domain knowledge

More biological background: Thank you for the suggestion. With the additional space, we will add further introductory background including a simple figure of how proteins and genes are interacting, and how such an interaction is represented in the prior matrices, and further explanations on the biological mechanism of gene regulation and its impact on disease, in particular cancer, where these regulatory networks are perturbed leading to the particular diseases, which makes GRN analysis so relevant. If you have further specific requests on what needs elucidation, we would appreciate your feedback.

Can knowledge be updated by the model: Yes, definitely. While the prior knowledge is a great starting point, DASH sparsifies based on both the prior knowledge and the patterns it learns from the data itself. If there are new patterns in the data that are not present in the prior knowledge, then the final model will have learned updated relationships on top of the prior domain knowledge.

Simulating effect of activating or deactivating genes: Yes, this is something the PHOENIX model itself (i.e. the base model on which we benchmark the different sparsification strategies) is able to do. We can intervene on an input and set it to zero and analyze the effect on the predicted dynamics.

Calculating $\Omega$: We are sorry for the confusion, during a revision of the paper we changed the notation inconsistently. The prior should be $C\in \mathbb{R}^{k\times k},$ where $k$ is the number of genes. Similarly $\boldsymbol{{\Omega^{(t)}_1}^{\intercal} \Omega_1^{(t)}} \in \mathbb{R}^{k\times k}$. So the product of the score matrix reflecting the first layer should align with the prior input-input relationship matrix $\boldsymbol{C}$. Similarly, $\boldsymbol{{\Omega^{(t)}_1} \Omega_2^{(t)}} \in \mathbb{R}^{r\times k}$ should align the prior output-input relationship matrix $\boldsymbol{P}$. The score is based on a basic reasoning that also underlies GRN inference techniques like [1], which argues that proxies of TF-TF interactions ($P$) or gene-gene interactions ($C$) are most readily available. The DASH pruning score aligns these proxies $P$ and $C$ with corresponding functions of the edges in the neural network. We have also provided more mathematical motivation below to Reviewer #4.

Finding new/undiscovered biological relations: Such a finding would be ideal and our model has been designed for this purpose. Inspired by your and other reviewer's question, we further analysed the obtained networks. Focusing on the Heme signaling pathway, which was uniquely identified by our suggested approaches. It turns out that this signaling pathway is highly relevant in breast cancer, we then proceeded to extract interaction partners from our model, i.e., the factors most highly affecting the Heme signaling molecules' dynamics, that could potentially be used as a new drug target. We provide further information on this finding in the general rebuttal response. We hope to discover further relations in the future, when we apply this approach in the field.

Pathway analysis: As described in Appendix B8, we first take a trained and sparsified model and compute gene influence scores, which we then use to compute pathway influence scores via permutation tests. The $z$-scores of the permutation tests are a measure of pathway enrichment.

Labelling consistency : Thank you for pointing this out. We will ensure consistency in a revised version.

[1] D Weighill et al. Gene Regulatory Network Inference as Relaxed Graph Matching, AAAI, 2021

Mathematical formulation: We apologize for the typos in dimensionality. Below we provide a brief note that includes all the corrected notation, definitions of key variables, and motivations behind the mathematical steps, especially pseudo-inverses.

DASH for $L=1$. For a single layer NN, with $k$ input and $r$ output neurons and corresponding weight matrix $\boldsymbol{W}$ should be $\in \mathbb{R}^{r\times k}$, and we compute pruning scores $\boldsymbol{\Omega} \in \mathbb{R}^{r\times k}$ using domain knowledge $\boldsymbol{P} \in \mathbb{R}^{r\times k}$. At any epoch $t$, we simply use $\boldsymbol{P}$ as the prior-based portion of $\boldsymbol{\Omega^{(t)}}$, the pruning scores at that epoch.

DASH for $L=2$. We are sorry for the confusion, during a revision of the paper we changed the notation inconsistently. The dimensionalities you provided are correct. The prior should be $C\in \mathbb{R}^{k\times k},$ where $k$ is the number of genes. Similarly $\boldsymbol{{\Omega^{(t)}_1}^{\intercal} \Omega_1^{(t)}} \in \mathbb{R}^{k\times k}$. Now, since $\boldsymbol{W_1^{(t)}}\in \mathbb{R}^{m\times k}$ learns how the $k$ inputs are encoded by $m$ hidden neurons, we surmise that the matrix product $\boldsymbol{{\Omega^{(t)}_1}^{\intercal} \Omega_1^{(t)}} \in \mathbb{R}^{k\times k}$ should approximately align with our prior knowledge of how the inputs co-vary, i.e. with $\boldsymbol{C}$. Since solving $\boldsymbol{{\Omega^{(t)}_1}^{\intercal} \Omega_1^{(t)}} = C$ is not feasible, we initialize $\boldsymbol{\Omega_1^{(0)}}$ as all 1s, and solve $\boldsymbol{{\Omega^{(t-1)}_1}^{\intercal} \Omega_1^{(t)}} = C$ instead. This is what motivates using the pseudoinverse.

Synthetic vs real world results: We agree that the synthetic data is only one part of the story, yet, it highlights distinct differences between the methods in a setting with available ground truth. If methods already fail to perform well on this "oversimplified" data, or show specific biases in terms of recovered structure, then such experiments are valuable. Moreover, we disagree that the real world results are "not as compelling", we get close to 90% accuracy in network recovery, which is much better than existing work (see Table 2). Naturally, there is a sparsity difference to the PathReg paper as the neural network architecture is different, thus different sparsity levels are required, as you correctly pointed out. We do recover a ranking of methods similar to the PathReg paper, which speaks for consistency. We are happy to dedicate more space to bone marrow results with the additional page in case of acceptance.

Balanced accuracy: Thank you for pointing that out, we will make this more concise and available in the beginning of the experiment section. In brief, the balanced accuracy measures whether an edge is correctly reconstructed weighted by the sparsity of the ground truth graph. The extraction is the same for all models, as each use the base PHOENIX architecture. The metric is computed with respect to a ground truth (in case of synthetic data) or gold standard (in case of real data). The gold standard ChIP-seq experiments are not included in and are distinct from the prior.

DASH+PHOENIX confounding: We believe that DASH should remain performant even when using an base model that is different from the PHOENIX. Hence we performed additional experiments (Table 1 in rebuttal PDF), with a simple two-layer MLP with ELU activation function as base model. The choice of ELU activation is motivated by the PathReg paper. We test a few sparsification strategies on this new base model applied to both synthetic data and the bone marrow data. We found DASH to still be performant.

Questions about PHOENIX: We apply all the sparsification strategies to the PHOENIX base model, that is the PHOENIX model without the biological regularization, a two-layer MLP with activation functions that resemble Hill kinetics. BioPrune is a pruning-based strategy that uses prior knowledge of the GRN to explicitly sparsify the neural network by setting weights to zero. Biological regularization takes an indirect route and implicitly encourages the neural network parameters towards zero through a penalty-term in the regularizer.

PINN+MP vs DASH: PINN+MP and BioPrune are novel models, as they have not been proposed in this context and we primarily discuss them as baselines compared to DASH. On synthetic data, PINN+MP is indeed a strong contender, with an MSE comparable to DASH (within error margin), yet deliver much less accurate representation of the underlying gene regulatory system (cf. balanced accuracy). Regarding Fig 3, the difference is barely visible in the small image. We provide a high-res version (Fig 1 in rebuttal PDF) where we also display the error between the inferred and true relationship. DASH outperforms PINN+MP which recovers lots of spurious structure.

Table 2: Yes, DASH and BioPrune do encourage the alignment of the model with the prior knowledge. But, the balanced accuracy is measured by comparing the sparsified model to ChipSeq validation data which is experimentally independent of the priors in question.

Training time: This is not a primary focus of our work, as runtime is generally negligible in comparison with applications in vision or language, e.g. DASH on PHOENIX takes approximately 40 minutes on the breast cancer data. The main motivation for pruning here is interpretability and alignment of sparsified models with known biology.

New knowledge: In principle, this approach can be used to generate new knowledge. We here focused on insights we can validate, which naturally means known knowledge. As a response to the question, we added a finding on Heme signaling based on new insights and its potential role in therapeutic design, which we will also add to the discussion. Due to the character limit, please refer to the general rebuttal for the details of this finding.

Noise levels: The noise was applied to both the expression data as well as the GRN prior, which we further describe in Appendix B1 and B3. We tested on three different noise levels (0%, 5%, 10%). We included results from the 5% setting in the main paper and the results for the remaining noise levels are in A1. Generally, we find that DASH is more robust to noise than its competitors and the performance gap increases with the noise level. This insight is also supported by the fact that DASH identifies more meaningful biology on real world data, which is indeed very noisy.

Evaluation with TF perturbation data: This is a great observation and we fully agree. However, this is prohibitively expensive and difficult to generate for time-course data and is hence not available up to our knowledge, especially on the genome-scale ($>10^4$ genes). We thus resorted to the next best thing, which is TF ChIP-seq data. If we could obtain perturbation data in the future, we would be eager to extend DASH to it.

DASH vs BioPrune: It is correct that BioPrune is simply using the prior GRN to calculate pruning scores, unlike DASH which takes both the prior GRN and the learned neural network weights into account when calculating pruning scores. However, the achieved sparsity is determined by the validation set. Thus, BioPrune can still contain edges that are not in the prior. The pruning process in concretely described in Appendix B.2.5, where we describe how each training epoch consisted of the entire training set being fed to the model, preceded by any pruning step that is prescribed by the pruning schedule. Training is terminated if the validation set performance fails to improve in 40 consecutive epochs. Upon training termination, we have obtained a model that has been iteratively sparsified to an extent that fails to improve the validation set performance. We discovered in our experiments that DASH's combination of data-based information (from the learned weights) and prior information (from the prior domain knowledge or GRN), leads to this achieved sparsity being much better than that of BioPrune. When we compare the GRN encoded by this sparse model structure back to the ground truth GRN, we see that the sparser model obtained by DASH encodes a much more accurate GRN than the denser model obtained by BioPrune.

Base model: This is the PHOENIX base architecture without any form of prior regularization. It is a fully connected two layer MLP with activation functions resembling Hill-like kinetics and gene specific multipliers for improved trainability. This base model is then subjected to the different sparsification strategies being benchmarked.

Why prior-informed regularization falls behind prior-informed pruning: This is because a pruning-based strategy (such as DASH) explicitly sparsifies the neural network by setting weights to zero, while regularization takes an indirect route and implicitly encourages the neural network parameters towards zero. The former leads to better sparsification in pruning based strategies, which is similar to the traditional pruning literature, where explicit pruning such as IMP, SNIP, SynFlow, etc usually outperform implicit $L_p$-based approaches.

Regarding Figure 3: The figure aids a qualitative comparison, where the overall structure of the matrix recovered by each method should be compared to the ground truth matrix (leftmost). It becomes evident that DASH as well as PINN+MP are much more aligned with the ground truth. A good quantitative measure of alignment here is the mean squared error (MSE) between the inferred and the ground truth relationships. We calculate this and include it in Additional Figure 1 (see rebuttal PDF), which clearly shows that DASH outperforms PINN+MP. We truncated Additional Figure 1 to only DASH and PINN+MP and increased the contrast to visualize the difference between these two, as they can appear to be very similar in Figure 3 of the paper. We will add this quantitative measure to all sparsification strategies in Figure 3 in the final version if the paper is accepted.

Thank you for staying engaged in the discussion!

As we tried to convey in our rebuttal, in principle, if large-scale perturbation data (e.g., by perturb-seq experiments) would exist for the tissues of interest such as breast cancer tissue, we would be happy to use it. But such experiments are usually (1) prohibitively expensive and (2) can not measure interventional perturbation effects within the tissue, but only on a cell-level - the tissue has to be extracted for gene-level perturbation experiments as they are in-vitro. Hence, we lose the ability of determining true effects in the whole tissue when using a perturb-seq-like experiment, while ChIP-seq gives a snapshot of the actual tissue state.

There also seems to be a misunderstanding, the provided reference [1] talks about expression difference after a global perturbation, for example an infection of cells. This gives expression difference of individual genes, but not how this was affected by other genes. In particular, this does not measure "causal" gene--gene effects such as the famous perturb-seq experiments based on CRISPR, which could be used as gold standard ground truth.

Considering perturb-seq and assuming it exists for our data, there is also an inherent limitation using (gene expression) perturbation data: we might see which genes $X$ have a causal effect (strength of up- or down-regulation) on a particular other gene $Y$, but these could be secondary or tertiary effects of the form $X \rightarrow Z_1 \rightarrow \ldots \rightarrow Y$, which should not correspond to links in the GRN. Only the immediate parent in the causal graph would ideally be represented there. With ChIP-seq, we do get the direct causal relationship (a TF $X$ is binding to the promoter of a gene $Y$), but lose the ability to tell what the exact causal effect is (up- or down-regulation, how strong, etc.).

That being said, in response to your review, we considered a comparison to a tf-inducing experiment (i.e., a bio-engineering to induce a particular TF's expression) on yeast, which was then manually curated to derive a "true" causal GRN [2], i.e., aiming to remove secondary and tertiary effects as discussed above. We note that while the underlying organism is the same yeast, the data was derived in different conditions, hence we expect differences in the GRNs and focus on relative performance differences between methods. To summarize the results, which we provide in the table below, the ranking of methods is similar to the comparison to the ChIP-seq gold standard, and there is still a large ($>10$ percentage points) improvement of DASH over existing state-of-the-art methods.

We would be happy to include further comparison in the future and appreciate any direct pointer to gene-level perturbation data relevant to us. Additional noise experiments will be provided in a separate comment.

[1] Kamal, A. et al. GRaNIE and GRaNPA: inference and evaluation of enhancer‐mediated gene regulatory networks. Mol Syst Biol e11627 (2023).

[2] Hacket, SR. et al. Learning causal networks using inducible transcription factors and transcriptome‐wide time series. Mol Sys Bio 16: e9174 (2020).

We understand your concern and added additional experiments with 20% respectively 40% noise on the prior, given in the table below. As expected, we do see a slight decrease in performance for prior-based methods correlated with the increase in prior noise. Yet, DASH still outperforms all existing work in terms of GRN accuracy even for 40% noise in the prior. We will add this additional analysis to the manuscript to provide a better discussion of robustness to prior noise. We thank the reviewer for their constructive feedback.

Large scale networks: For our particular domain, this problem was addressed by the design of the PHOENIX architecture, which scales up to large and complex networks commensurate to the whole human genome (on the order of $10^4$ genes/dimensions). Here we demonstrate that DASH works effectively on the PHOENIX architecture applied to such large scale datasets (e.g. the breast cancer dataset covers 11165 genes).

Regarding Figure 3: For the noise levels of 0%, 5%, and 10%, the percentage of true relationships captured by DASH are 98%, 95%, and 88% respectively, while the corresponding numbers for PINN+MP are 98%, 98% and 90%. However, this is not a good quantitative measure of the how well each method recovers the ground truth, since PINN+MP recovers a lot of spurious features. This is shown in additional Figure 1 (see rebuttal PDF) where we zoomed in to just these two methods and increased the contrast. A better quantitative measure of alignment is the mean squared error (MSE) between the inferred and the ground truth relationships. We calculate this and include it in additional Figure 1, which clearly shows that DASH outperforms PINN+MP.

Integration of domain knowledge for downstream disease related outcome tasks: We agree that incorporating domain knowledge leads to more meaningful learned dynamics, which is why we anticipate our approach DASH to find many relevant practical applications. In Genomics, also in the biomedical domain, the goal is often not prediction but inference, which is also the case here: We would like to learn more about the inherent structure of the disease in terms of the complex gene regulations. By learning about these structures, such as learning about specific genes that strongly affect other genes in a specific cancer, we can better understand the disease process and ultimately can design better treatments. More concretely, we can for example try to understand the effect of a drug on the gene regulatory network by investigating the change in gene regulatory dynamics between cells exposed to the drug versus those that were not exposed, which is what the authors of PHOENIX -- the model we use as basis for our experiments -- did in their original work. We also added an additional result, where we found an interesting biological insight, which we provide in the general rebuttal response due to character limit. In particular, we show how the insights derived from the GRN can be used to generate a potential new therapeutic approach.

"New" knowledge from DASH: For both breast cancer as well as yeast cell cycle, we use gold standard ChIP-seq data (a biological experiment to locate binding of TFs to the genome) to validate the biology of the inferred GRN. I.e., ChIP-seq experiments on the particular cells gives us a gold standard GRN for validation and we additionally note that this data has not been used for the construction of the prior. For higher-level knowledge, we have included a biological pathway analysis and found that the most relevant pathways identified in our model correspond to key paths in breast cancer progression, yeast cell cycle progression, or, respectively, bone marrow hematopoesis that align with the known biology of these processes. More concretely, we provide the example of the following pathways: "... TP53 activity and FOXO267 mediated cell death, both of which are highly relevant in cancer [ 37 , 24 ]". In response to your question, we further investigated the discovered biology focusing on the Heme signaling pathway, which was uniquely identified by our method. We provide more details on these interesting findings in the general author rebuttal due to space constraints.

A learnable $\lambda$: Here $\lambda$ is a hyperparameter of DASH. As mentioned in Appendix B4, the $\lambda$ values are determined using a $K$-fold cross validation approach. Specifically, we have automated a process that fits multiple models to the same dataset each with a different set of lambda values chosen from a grid. This automated grid-search is used to optimize the $\lambda$ values, based on predicted MSE on the validation set. This means that $\lambda$ adapts to the complexity of the data. Learning $\lambda$ in a differentiable manner instead would involve differentiating a complicated loss that is defined on cross-validation data, which is usually not efficient.

Test set size: For synthetic data, we know the ground truth generative model and evaluate on that. For breast cancer, we use 6% of data to evaluate the MSE, which we picked as data is really scarce and we need sufficient number of samples to train the model. We do, however, have data of an independent biological experiment (the ChipSeq data) to evaluate the inferred model in terms of reconstructed GRN, hence the performance is also generalizable. For the yeast cell cycle data, we do test on an entire biological replicate (i.e., have one replicate for training, one for testing). We do see that this information might got lost in the details and will clarify in the revised manuscript.

Complexity of methodology: The score is based on a basic reasoning that also underlies GRN inference techniques like [3], which argues that proxies of TF-TF interactions ($P$) or gene-gene interactions ($C$) are most readily available. The DASH pruning score aligns these proxies $P$ and $C$ with corresponding functions of the edges in the neural network. We have also provided more mathematical motivation below to Reviewer #4. We are happy to provide even more intuition on other aspects that the reviewer wants more clarity on.

Biological validation: In terms of expert validation, one of our co-authors has deep expertise in molecular biology, and they have helped in the biological analysis of the GRNs inferred by DASH. We further provide the example of the following pathways: "... TP53 activity and FOXO267 mediated cell death, both of which are highly relevant in cancer [ 37 , 24 ]". Unfortunately, due to the strong space constraint, we had to present further biological validation in the appendix (see Section A). Furthermore, note that yeast cell cycle as well as breast cancer results were validated using ChIP-seq, which is an independent experimental validation of biological plausibility. As an additional result, we found an interesting biological insight, which we provide in the general rebuttal response due to character limit.

Numerical examples of $P$: For synthetic data, we use a corrupted prior based on the data-generating model, which we elaborate on in App B.3 and B3.1. For real data, we use general information of transcription factor binding to gene promoter regions as prior information, which can be computed from binding motif matches with the corresponding genome (human respectively yeast). The result is a matching score that can be thresholded to get the {0,1} matrix encoding which (TF-encoding) gene has a relationship with which other gene. We follow the approach of Guebila et al. [1] to get matrix P. As prior C, we use the STRING database [2], which gives a general (i.e., not tissue-specific) graph of protein-protein interaction. Here, we use the interactions based on experimental evidence only and employ a cutoff of .6 to get a binary adjacency matrix. We will add the additional information to the discussion of priors in Appendix B.3.

Parameter tuning: As mentioned in Appendix B4, the $\lambda$ values are determined using a $K$-fold cross validation approach. Specifically, we have automated a process that fits multiple models to the same dataset each with a different set of lambda values chosen from a grid. This automated grid-search is used to optimize the $\lambda$ values, based on predicted MSE on the validation set.

Generalizability: We focused here on the considerably large field of gene regulatory networks, which has many applications in the biomedical domain. That said, we anticipate that DASH could be applicable to other core science domains, given the task's neural network structure follows a similar MLP layout and we have domain knowledge of similar structure (like matrices about input-output relationships $P$ and output-output relationships $C$). In principle, the general definition of the score does not depend on the activation functions. The alignment of prior information with the pruning can be seen as an alignment of the prior with the amount of information that flows through the specific edges. Regarding CNNs and transformers, the application of DASH might be less straight-forward and requires further thinking on which part of the network (connections through filters, parts of attention heads) should be aligned with the prior. This is not directly evident, as such components are not directly interpretable in this context and are not necessarily associated with genes. We would be happy to add a detailed discussion about this aspect.

Comparative analysis: We compared DASH to 8 other strategies for neural network sparsification, and additionally BioPrune and PINN+MP which are our own suggested strong baselines for comparison against DASH. Some of these eight strategies are classic benchmarks for comparison in the pruning literature, including IMP, which is considered as the standard baseline in the field, while others such as PathReg have been recently proposed for this exact problem setting. The PINN based method for inducing sparsity in the baseline PHOENIX model is as recent as 2024. We would appreciate specific references to relevant recent sparsification methods that include prior knowledge as we are not aware of any further ones.

Bias: This is a great point, and we would be happy to discuss this as part of ethical considerations.

[1] MB Guebila et al. GRAND: a database of gene regulatory network models across human conditions. Nucleic Acids Res. 2022 [2] D Szklarczyk et al. The STRING database in 2023: protein-protein association networks and functional enrichment analyses for any sequenced genome of interest. Nucleic Acids Res. 2023 [3] D Weighill et al. Gene Regulatory Network Inference as Relaxed Graph Matching, AAAI, 2021