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

Authors: intekhab  hossain
harvard t h chan
jonas  fischer
rebekka  burkholz
john  quackenbush
Venue: NeurIPS_cc_2024_Conference
Published: 2024-05-15 19:21:16

Abstract: The practical utility of machine learning models in the sciences often hinges on their interpretability. It is common to assess a model's merit for scientific discovery, and thus novel insights, by how well it aligns with already available domain knowledge - a dimension that is currently largely disregarded in the comparison of neural network models. While pruning can simplify deep neural network architectures and excels in identifying sparse models, as we show in the context of gene regulatory network inference, state-of-the-art techniques struggle with biologically meaningful structure learning. To address this issue, we propose DASH, a generalizable framework that guides network pruning by using domain-specific structural information in model fitting and leads to sparser, better interpretable models that are more robust to noise. Using both synthetic data with ground truth information, as well as real-world gene expression data, we show that DASH, using knowledge about gene interaction partners within the putative regulatory network, outperforms general pruning methods by a large margin and yields deeper insights into the biological systems being studied.

Main paper: 
 introduction: 
 With ever-growing neural network architectures encouraged by the success of overparametrization, with over a trillion parameters in a single model such as GPT4, there is a similarly growing demand for sparser, more parameter-efficient neural networks that are more resource-friendly and interpretable. The Lottery Ticket Hypothesis (LTH) provides an empirical existence proof of sparse, trainable network architectures [18] , that eventually achieve a similar performance as their dense counterpart. Subsequent work introduced structured pruning approaches, facilitating group-wise neuron- [32, 63] , or channel-sparsity [34, 25, 26] , which are, however, focused on the structure of the architecture design, aiming for better alignment with hardware implementations to eliminate operations, rather than structure that reflects relevant domain information. Especially for scientific discovery, an alignment of learned structure with such domain knowledge is, however, essential for interpretability, as only then the model represents meaningful domain-relevant relations. Such problem settings often occur for example in physics or biology where a learned model should give an explanation to be able to form a hypothesis. This poses the question: How can we select among multiple predictive models and promote the search for meaningful neural network structures?To guide the learning process, we argue that we need additional problem-specific structural information, and should leverage any available -and reliable -domain knowledge. One of the fundamental tasks of molecular biology is to understand the gene regulatory dynamics in health and disease. Gene regulatory dynamics describe the changes of the expression of a gene-the generation of small copies of a DNA segment that can serve among others as blueprint for proteins-dependent on other regulatory factors such as transcription factors, which are proteins that bind next to the gene (DNA segment) to modulate its expression. Yet, the exact dynamics are far from understood and changes in these dynamics can be drivers for diseases such as cancer. As such, improving an understanding of the mechanics behind these dynamics increases the understanding of the disease and can ultimately inform therapy design. This problem setting of estimating gene regulatory dynamics requires high interpretability, as an understanding of the true biological mechanics-the relationship between regulatory factors and a gene's expression-is needed, as well as sample-efficiency, as generating time-course data even for a few patients is extremely expensive. While models to estimate gene regulatory dynamics have been suggested [14, 1, 61, 27] , none of these is particularly sparse or interpretable. We propose a new approach of network sparsification that guides pruning by domain knowledge implemented for a neural model for estimating gene regulatory dynamics, which yields networks that are very sparse, align with underlying biology, while accurately predicting dynamics. This idea of prior-informed or domain-aware pruning is at the heart of this paper. In particular, we propose DASH (Domain-Aware Sparsity Heuristic, Fig. 1 ), an iterative pruning approach that scores parameters taking into account structural domain-knowledge. With DASH, it is possible to control the level of prior information taken into account for pruning and it automatically finds an optimal sparsity level aligned with both the prior and the data. Considering the task of estimating gene regulatory dynamics, we first show in synthetic experiments that DASH generally outperforms standard (taskand architecture agnostic) pruning as well as task-specific pruning approaches. On real data with a reference gene regulatory network (GRN) derived from gold standard biological experiments, we show that DASH better recovers the reference GRN and reflects more biologically plausible information. On recent single cell data on blood differentiation, we show that DASH, in contrast to existing work, identifies biologically relevant pathways that can be used to generate new insights and inform domain experts. We anticipate that our work serves both for future benchmarking on how well pruning approaches are in structure learning, as well as a blueprint for guided network pruning in fields where domain knowledge is readily available, such as in other hard sciences including physics or material science, where knowledge about variables, e.g., associations between atoms or molecules or equations relating quantities in a system, is available. 
related work: 
 The Lottery Ticket Hypothesis (LTH) [18] provides an empirical existence proof of sparse, trainable neural network architectures. It conjectures that dense, randomly initialized neural networks contain subnetworks that can be trained in isolation with the same training algorithm that is successful for the dense networks. However, a strong version of this hypothesis [52, 66] , which has also been proven theoretically [43, 49, 47, 17, 5, 12, 15] , suggests that the identified initial parameters are not only specific to the sparse structure but also the learning task and benefit from information about the larger dense network that has been pruned [48] . Acknowledging this strong relation, other works have proposed to combine mask and parameter learning directly in continuous sparsification approaches [54] that employ regularization strategies that approximate L0 penalties [54, 41, 31, 56] . In the following, we recap the key ideas behind the methods most relevant to ours.Magnitude pruning (MP) In magnitude pruning (MP) a neural network is trained and then (post-hoc) pruned to a desired sparsity level by masking the corresponding proportion of smallest magnitude weights. This smaller, masked network is then further trained, reminiscent of fine-tuning [22] .Figure 1 : DASH. A NN, here a neural ODE for gene regulatory dynamics, is traditionally sparsified in a data-centric way (top). Pruning is done based on data alone, the pruning score Ω is a function of the learned weights W . Such sparsified models often do not learn plausible relationships in the data domain. We propose DASH (bottom), which additionally incorporates domain knowledge P into the pruning score Ω, yielding sparse networks giving meaningful and useful insights into the domain.More formally, we start with a fully connected neural network NN Θ with L layers parametrized by Θ := {(W l , b l )} L l=1 . MP is performed after training is complete (i.e. post hoc). For a suitable threshold of p%, MP "prunes" the trained Θ by setting the smallest (absolute value) p% of weights in Θ to 0. The choice of parameters to prune can either be made in an unstructured way, by choosing the the lowest p% across all {W l } L l=1 , semi-structured, by pruning the lowest weights per layer, or in a structured way such that, e.g., neurons with lowest p% average outgoing weight are pruned. Once pruning is complete the p%-sparse NN Θ ′ is fine-tuned on the data so that Θ ′ is learned appropriately. Iterative magnitude pruning (IMP) [18] suggest to alternate between training and magnitude pruning, iteratively sparsifying the network, to date still of the most successful pruning strategies. In practice, a pruning schedule is used to iteratively sparsify Θ, until a Θ ′ with desired target sparsity or predictive performance plateau is reached. Importantly, weights are reset to initial pre-training values, either after each round of pruning or once after the target sparsity has been reached. Pruning with rewinding [53] have demonstrated that rewinding weights to an earlier training point-a compromise between fine-tuning of MP and reset to initialization of IMP-provides good performance, which also has been suggested in the context of Neural ODEs as SparseFlow [37] .Implicit penalty-based approaches C-NODE [2] seek to reduce the overall number of input-output dependencies (i.e. paths of contribution from input neuron i to output neuron j) through {W l } L l=1 . The approach can result in both feature and weight sparsity in NeuralODEs. L 0 [40] incorporate a differentiable L 0 norm regularizer term in the objective. It implicitly prunes the network by encouraging weights to get exactly to zero. The L 0 regularizer is operationalized using non-negative stochastic gates which act as masks on the weights. PathReg [1] innovatively combines the strengths of both C-NODE and the L 0 approach to promote both weight and feature sparsity in NeuralODEs. It uses stochastic gates similar to [40] and add a C-NODE-inspired penalty terms that constrain the overall number of input-output paths by regularizing the probability of any path from input i to output j being non-zero.As application, we consider estimation of gene regulatory dynamics. Early work, such as COPASI [45] use a fully parametric modeling approach that are limited in their prediction capabilities. With recent advances in machine learning, tools such as Dynamo [51] , PROB [58] , and RNA-ODE [39] aim to learn regulatory ODEs using sparse kernel regression, Bayesian Lasso, and random forests, respectively. Leveraging high flexibility and performance of neural models, PRESCIENT [65] uses a simple NN to learn regulatory ODEs, whereas tools such Figure 2 : Results on simulated data. We visualize performance of pruning strategies in comparison to original PHOENIX (baseline) in terms of achieved sparsity (x-axis) and balanced accuracy (y-axis) of the recovered gene regulatory network against the ground truth on the SIM350 data with 5% noise. Error bars are omitted when error is smaller than depicted symbol. ✓ indicate methods that leverage prior information. Top left is best: recovering true, inherently sparse biological relationships.as DeepVelo [9] and sctour [35] have a variational autoencoder as backbone. The latest line of research [14, 1, 61, 27] uses neural ordinary differential equations [8] . However, a key limitation of these methods is the lack of interpretability arising from non-sparse dynamics that do not align with ground truth biology. Consequently, the induction of sparsity in gene regulatory ODEs has been an active area of research with C-NODE and PathReg as most recent advances [2, 1] , the achieved sparsity levels are, however, not yet sufficient to capture the relevant biology.3 Domain-aware pruning with DASH While the above sparsification strategies have shown to perform well in various settings, the resulting models are often either not particularly sparse or do not reflect meaningful domain knowledge. We hypothesize that this is due to two reasons: (1) the difficulty of identifying a good sparse network and (2) the current focus on hardware-centric rather than task-centric pruning, valuing structural pruning of a model in terms of groups of neurons (layers, channels) over structural pruning reflecting task-specific knowledge. To overcome these problems, we suggest to ground the model search (here, the network training) with existing domain-specific knowledge, which eases identifiability due to the introduced constraints and enables task-aware pruning to identify meaningful domain knowledge.In the following, we propose DASH (Domain-Aware Sparsity Heuristic), an iterative pruning-based approach that accounts for prior knowledge by scoring parameters in terms of their alignment with this prior, and show its usefulness for a neural model for the inference of gene regulatory dynamics. We assume the domain knowledge to be given as input-output relations (e.g., known protein-gene interactions in molecular biology), for which we want the network flow between any input and output to align with. Suppose our domain knowledge for a task from this domain is given as a real-valued relationship graph G = (V, E), where nodes are relevant entities from the domain, e.g. genes or proteins, and edges are a strength of association between these entities, e.g. evidence of association derived from literature or experiments. Examples for such relational information in case of protein-gene associations can be derived from protein binding profiles [57] . For the rest of the paper, we will assume that this domain knowledge is given as a matrix P ∈ R k×r , which encodes the strength of association between the k inputs and r outputs for our task of interest, such as known proxies of protein-gene interactions . Intuitively, for a one-layer neural network, we encourage pruning scores for a (neural) network edge to be proportional to the corresponding edge in the prior knowledge graph G while still taking into account the data-specific knowledge, thus enabling learning of new knowledge and robustness to wrong or missing information in the prior. We begin with this simple base case of task-aware pruning of a fully connected neural network with L = 1 layer and extend to more layers below.DASH for L = 1. For a single layer NN, with k input and r output neurons and corresponding weight matrix W ∈ R r×k , we compute non-negative pruning scores Ω ∈ R r×k by leveraging the domain knowledge P ∈ R r×k . In practice, we allow balancing between data-driven and prior-knowledge-driven pruning, implemented through a convex combination of the learned weights W and prior domain knowledge P controlled by the parameter λ ∈ [0, 1]. Alternating between training and pruning akin to Iterative Magnitude Pruning [18] , we set the following pruning score during a pruning phase:Ω (t) := (1 -λ) |W (t) | + λ|P | ,where |W (t) | represents the appropriately normalized matrix (details in Appendix B.2) of absolute weights as learned up to epoch t. We then prune the parameters in W (t) corresponding to the lowest absolute p t % of entries in Ω (t) , where p t is the desired sparsity level at time t given by a schedule.DASH for L = 2. For two-layered NNs with weights W 1 ∈ R m×k , W 2 ∈ R r×m , i.e k inputs, r outputs, and m hidden neurons, we consider knowledge about input-output relationships P ∈ R r×k as before. We can additionally use further knowledge about input-input relationships C ∈ R k×k . In molecular biology this could be information about binding or interaction partners available in databases such as STRINGDB [59] , which has also been employed to guide static gene regulatory network inference [62] , or co-regulators, derived from co-occurrence of proteins. Intuitively, we now project pruning scores for the first layer to the prior knowledge about input-input relations, encouraging closeness to this prior, while projecting the product of pruning scores of first and second layer to known input-output relations, reflecting the flow of information from input to output through these two layers. Given that W are the corresponding pruning scores, we surmise that the matrix productΩ (t) 1 ⊺ Ω (t)1 ∈ R k×k should approximately align with the prior knowledge C. Since solvingΩ (t) 1 ⊺ Ω (t) 1 = C is not directly feasible we initialize Ω (0)1 randomly, and resort to solving a recurrence relation version of problem, that isΩ (t-1) 1 ⊺ Ω (t) 1 = C.Using the left and right pseudo-inverse to obtain a solution to the above, defined as PInv L (X) = (X ⊺ X) -1 X ⊺ and PInv R (X) = X ⊺ (XX ⊺ ) -1 respectively: Ω (t) 1 := (1 -λ 1 ) |W (t) 1 | + λ 1 PInv L Ω (t-1) 1 ⊺ • C . PInv L Ω (t-1) 1 • C encourages Ω (t) 1 ⊺ Ω (tΩ (t) 2 := (1 -λ 2 ) |W (t) 2 | + λ 2 P • PInv R Ω (t) 1 .Similar to the case for DASH for L = 1 layer, we can now prune the parameters of W DASH for L > 2 For many interpretability-centric tasks, including our application to gene regulatory networks, small architectures of L = 2 are common, as domain experts are interested in understanding the exact flow of information through the network. Furthermore, we know that two-layer neural networks exhibit universal approximation [11] . We however hypothesize that the technique of computing pruning scores by fixing those of preceding layers can be extended to a larger number L of fully connected layers and elaborate in App. B.6.Flexibility While the λ l can be tuned using cross-validation (see App. B.4), we note that it allows for flexibly encoding different pruning philosophies. Specifically, when λ l = 0 ∀l, DASH corresponds to SparseFlow, and when λ l = 1 ∀l, DASH represents fully prior-based sparsification (which we term "BioPrune" and consider as experimental baseline). Here, we consider the task of predicting the regulatory dynamics of gene expression. To be able to understand the model and transfer it to clinical practice, interpretability is key. The most recent developments in modeling gene regulatory systems allow to model actual (temporal) regulatory dynamics, but require complex models, such as NeuralODEs, that hinder interpretability. While state-of-the-art results are now achieved with shallow architectures [27] that are more tractable than deep, heavily over-parameterized networks, these models still encode information across many thousands of weights and we show experimentally that such information does not reflect true biology well. In fact, true gene regulatory networks and hence their underlying dynamics are inherently sparse [6] . This sparsity should be properly reflected by neural dynamics models. The PHOENIX NeuralODE model will serve as our base model for applying sparsification strategies and we show that pruning aligned with prior domain knowledge improves interpretability as well as quality of inferred (new) knowledge.In a nutshell, given a time series gene expression sample for k genes, PHOENIX uses NeuralODEs to construct the predicted trajectory between gene expression g(t) ∈ R k (inputs) at time t = t i to any future expression g(t i+1 ) (outputs), by implicitly modeling the RNA velocity (dg/dt). PHOENIX uses biokinetics-inspired activation functions to separately model additive and multiplicative coregulatory effects. The trained model encodes the ODEs governing the dynamics of gene expression, which can be directly extracted for biological insights. We apply DASH to PHOENIX and give a brief review of the PHOENIX architecture in App. B.10 and a detailed account on how to apply DASH to this architecture in App. B.11. Next, we provide experiments on synthetic and real data showing the advantages of prior-informed pruning on the task of predicting gene-regulatory dynamics.Figure 3 : Reconstruction of ground truth relationships. Estimated effect of gene g j (x-axis) on the dynamics of gene g i (y-axis) in SIM350 for different levels of noise (rows). Ground truth is given on the left, our suggested approach and baselines (DASH, BioPrune, and PINN+MP) on the right with mean squared error between inferred regulatory relationships and ground truth in purple.TF perturbation experiments [21] . 
theory assumptions and proofs: 
 Question: For each theoretical result, does the paper provide the full set of assumptions and a complete (and correct) proof? Answer: [NA] Justification: There are no theorems nor proofs in this paper. Guidelines:• The answer NA means that the paper does not include theoretical results.• All the theorems, formulas, and proofs in the paper should be numbered and crossreferenced. • All assumptions should be clearly stated or referenced in the statement of any theorems.• The proofs can either appear in the main paper or the supplemental material, but if they appear in the supplemental material, the authors are encouraged to provide a short proof sketch to provide intuition. • Inversely, any informal proof provided in the core of the paper should be complemented by formal proofs provided in appendix or supplemental material. • Theorems and Lemmas that the proof relies upon should be properly referenced. 
experimental result reproducibility: 
 Question: Does the paper fully disclose all the information needed to reproduce the main experimental results of the paper to the extent that it affects the main claims and/or conclusions of the paper (regardless of whether the code and data are provided or not)? Answer: [Yes] Justification: We provide discussion of synthetic data generation (App. B.1), setup of training (App. B.2), and information about how prior knowledge was defined (App. B.3). We further provide information on which reference gold standard was used for evaluation of the inferred gene regulatory network (App. B.5), any pre-processing steps for the real world data (App. B.7) and how analyses were carried out (App. B.8). We further discuss all necessary details on how different pruning methods were implemented (App. B.9) and how regulatory networks were extracted (App. B.10.4). We provide an implementation of DASH as supplementary material.Guidelines:• The answer NA means that the paper does not include experiments.• If the paper includes experiments, a No answer to this question will not be perceived well by the reviewers: Making the paper reproducible is important, regardless of whether the code and data are provided or not. • If the contribution is a dataset and/or model, the authors should describe the steps taken to make their results reproducible or verifiable. • Depending on the contribution, reproducibility can be accomplished in various ways.For example, if the contribution is a novel architecture, describing the architecture fully might suffice, or if the contribution is a specific model and empirical evaluation, it may be necessary to either make it possible for others to replicate the model with the same dataset, or provide access to the model. In general. releasing code and data is often one good way to accomplish this, but reproducibility can also be provided via detailed instructions for how to replicate the results, access to a hosted model (e.g., in the case of a large language model), releasing of a model checkpoint, or other means that are appropriate to the research performed. • While NeurIPS does not require releasing code, the conference does require all submissions to provide some reasonable avenue for reproducibility, which may depend on the nature of the contribution. Justification: We do not see a direct harmful application of our suggested approach, it is a foundational line of research and our suggested application is in molecular biology. As such, it will rather bring a positive societal impact by improving understanding of health and disease and thereby improving therapy design.Guidelines:• The answer NA means that there is no societal impact of the work performed.• If the authors answer NA or No, they should explain why their work has no societal impact or why the paper does not address societal impact. • Examples of negative societal impacts include potential malicious or unintended uses (e.g., disinformation, generating fake profiles, surveillance), fairness considerations (e.g., deployment of technologies that could make decisions that unfairly impact specific groups), privacy considerations, and security considerations. • The conference expects that many papers will be foundational research and not tied to particular applications, let alone deployments. However, if there is a direct path to any negative applications, the authors should point it out. For example, it is legitimate to point out that an improvement in the quality of generative models could be used to generate deepfakes for disinformation. On the other hand, it is not needed to point out that a generic algorithm for optimizing neural networks could enable people to train models that generate Deepfakes faster. • The authors should consider possible harms that could arise when the technology is being used as intended and functioning correctly, harms that could arise when the technology is being used as intended but gives incorrect results, and harms following from (intentional or unintentional) misuse of the technology. • If there are negative societal impacts, the authors could also discuss possible mitigation strategies (e.g., gated release of models, providing defenses in addition to attacks, mechanisms for monitoring misuse, mechanisms to monitor how a system learns from feedback over time, improving the efficiency and accessibility of ML). 
experiments: 
 For evaluation we consider synthetic data from an established simulator tool [4] , as well as real world data of gene expression from breast cancer tissue [13] , from yeast with synchronized cell cycle [50] , and from human bone marrow [1] . In case of synthetic data, we use the ground truth regulatory system from the generating model for validation. For breast cancer and yest cell cycle data we use additional experimental data (ChIP-seq) from the corresponding studies, which are independent goldstandard biological experiment measuring sample-specific TF-gene interactions, to evaluate inferred regulatory relationships. We measure the correctness of a GRN learned by a model (see App. B.10.4) in terms of balanced accuracy, which measures whether an edge is correctly reconstructed weighted by the sparsity of the aforementioned ground truth graph. To evaluate predictive performance for real data, we use a 6% hold-out test set for breast cancer and one of the biological replicates hold out from training for the yeast data. As prior knowledge we leverage 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 a 0, 1-based matrix encoding which (TF-encoding) gene has a relationship with which other gene. We follow the approach of Guebila et al. [3] to obtain matrix P .As prior C, we use the STRING database [59] , 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. (for more details, see App. B.3).To compare pruning strategies, we consider the PHOENIX model as a basis, which is the state-of-theart NeuralODE for estimating gene regulatory dynamics [27] and provide an ablation on a standard MLP architecture (see App. Tab. 7, App. Tab. 10, and App B.12). We compare DASH against the PHOENIX model without additional pruning as performance reference, and suggest two simple yet powerful baselines, which is post-hoc magnitude pruning of weights followed by finetuning (PINN+MP), and BioPrune, a fully prior-based pruning (cf. Sec. 3). From the literature, we consider L 0 -regularized pruning [40] , C-NODE [2] , and PathReg [1] , which have been recently proposed for the inference of sparse gene-regulatory relationships, PHOENIX with biological regularization [27] , and dynamic sparse training (DST) [38] , all of which are implicit pruning approaches. We further consider explicit, iterative score-based pruning approaches including Iterative Magnitude Pruning (IMP) [18] , the flow-based model-agnostic pruning method SynFlow [60] , and the flow-based Neural ODE pruning SparseFlow [37] . We tune hyperparameters, including λ for DASH, on a validation set. Unlike other methods (e.g., L 0 ) DASH does not prolong the runtime of training much. A common pruning schedule where pruning scores are computed once every 10 epochs only increases runtime by <2% for the full model fitting process.We simulate gene expression time-series data from a fixed dynamical system, the ground truth was thus known (see App. B.1). In short, we adapt SimulatorGRN [4] to generate noisy time-series expression data from two synthetic gene regulatory systems (SIM350 and SIM690, consisting of 350 and 690 genes, respectively). We split trajectories into training (88%), validation (6% for tuning λ), and testing (6%).We evaluate all methods in terms of achieved sparsity and MSE of predicted gene expression values on the test set (App. B.2, B.4) and investigate biological plausibility by calculating balanced accuracy of regulatory relationships extracted from the PHOENIX model (for details, see App. B.5, B.10.4), We here report the results for the data of 350 genes and 5% noise, noting that results are consistent across different noise levels and with more number of genes (see App. Sec. A.1).A general trend across all experiments that aligns with our initial motivation is that dense models (sparsity < 50%) have a significantly worse reconstruction of the underlying biology -the ground truth GRN -than sparse models (sparsity > 80%) (see Fig. 2 ). Furthermore, we see that DASH retrieves not only among the sparsest networks, but also reflects the underlying GRN best across all methods, outperforming comparably sparse IMP by about 20 percentage points accuracy in different settings, even with decrease in quality of the prior (see sensitivity analysis in App. Tab. 6). Due to the prior-informed structured pruning, it is able to occupy the sweet spot of highly sparse at the same time biologically meaningful models.Consistent with the literature, PathReg outperforms L 0 as well as C-NODE in terms of sparsity [1] , we additionally find evidence that it also delivers more biologically meaningful results. Yet, IMP as well as prior-informed pruning approaches outperform PathReg by a large margin. The MSE of predicted gene expression of DASH is among the best, within one standard error of the best overall method.The only better approach is our suggested baseline, a combination of posthoc magnitude pruning of PHOENIX combined with additional finetuning (PINN+MP), which is, however, impractical as it requires to train and prune many PHOENIX models along a grid of sparsity levels (see App. B.9.3).Visualizing the estimated against ground truth regulatory effects (i.e., functional relationships between variables), we observe that DASH captures the effects much better than competitors (see Fig. 3 ). Virtually all existing approaches discover spurious regulatory effects, whereas prior-informed pruning identify the main regulatory effects correctly. Moreover, with increasing levels of noise in the simulation, we observe that both BioPrune as well as PINN+MP start finding spurious dependencies, while DASH still recovers the overall structure well. While not perfect, as seemingly there are more dependencies than in the sparse ground truth, potentially introduced by correlations between features, DASH provides a sparse estimation of regulatory effects that most closely resembles the ground truth relationships among existing work.To investigate the performance of DASH on real data, we consider gene expression measurements from a cross-sectional breast cancer study [13] . This data of 198 breast cancer patients with measurements for 22000 genes has been preprocessed and ordered in pseudotime [58] , which we use as basis for our experiments (cf App. B.7). Across methods, we observe that implicit sparsification methods generally perform poorly in terms of sparsity and accuracy of recovered relationships (see Tab. 2). While pruning-based sparsification approaches achieve greater sparsity and performance in predicted gene expression, with SparseFlow reaching the highest sparsity (95.7%) among all methods, the recovered biological relations are not better than random chance, which renders the underlying models useless for scientific discovery. DASH in contrast finds a comparably sparse network (92.7% sparsity) while having top of the line performance in terms of test MSE and high alignment with true biology (95.7% balanced accuracy). For this particular dataset, we observe that DASH primarily builds on the prior knowledge, not surprisingly performing similarly as BioPrune, which is our suggested baseline pruning approach taking only the prior into account. We will see for other real-world data that this weight of domain knowledge is highly task-specific and BioPrune yields sub-optimal results on different data.To better understand whether the inferred gene regulatory dynamics align with meaningful biology, we additionally perform a pathway analysis (see App. B.8). Such pathway analysis are a standard approach for domain experts to distill information for example for therapeutic design. The genes that show the highest impact on the dynamics within the derived model are tested whether they enrich in a specific higher level biological pathways. For the top-20 most significantly enriched pathway per model (App. Fig. 5 ), we observe that in contrast to prior-informed methods, the existing pruning approaches show only very few significant pathways, consistent with our quantitative results on [44, 28] , are only visible in models pruned with prior information. This provides evidence that pruning informed by a biological prior recovers biological signals that are relevant in the disease and which can not be picked up otherwise. Furthermore, we find Heme-signaling as a pathway uniquely identified as relevant in our approaches (cf. App. Fig. 5 ). Heme as a signaling molecule has key roles in the gene regulatory system [46] , and turns out to have an anti-tumor role in breast cancer specifically [19] . Subsequent approaches pharmaceutically targeting Heme signaling showed success [30] , with one of the key regulators affected being Bach1. To suggest further targets for e.g. combination treatment, we hence examined the top-5 regulatory factors in terms of weights in our estimated gene regulatory dynamics. These factors include PBX1 and FOXM1, for which a drug repurposing of existing compounds, such as [55] , could lead to a potential new treatment for this specific cancer.Yest cell-cycle data We next consider real data of synchronized yeast cell [50] (see App. B.2 for training setup). We observe an overall trend similar to the breast cancer study in terms of achieved sparsity and balanced accuracy (cf. Tab. 2), with implicit sparsification methods generally finding significantly less sparse models and all methods that do not incorporate prior knowledge having inferred relationships that are not better than random chance. For this data, however, DASH finds an optimal lambda value that incorporates more data-specific knowledge (λ = 0.75) compared to the breast cancer study above. This shows the advantage of DASH over our BioPrune baseline model (prior-only pruning), as here we gain about 10% points in balanced accuracy over BioPrune for ChIP-seq validation data [23] and 2% points over BioPrune for an independent TF perturbation network curated to derive a "true" causal GRN [21] (see App. Tab. 8). We also retrieve a 3% points sparser model. Comparing inferred biological knowledge between BioPrune and DASH through a pathway analysis, we see that DASH recovers more significantly enriched pathways related to cell cycle processes (cf. App. Fig. 6 ).Cell differentiation in human bone marrow Lastly, we investigate the performance of DASH in an exploratory setting with single cell data of human bone marrow ordered in pseudotime [1] . Here, we are interested in better understanding the gene regulatory dynamics of blood cell differentiation, the process of hematopoietic stem cells specializing into cells taking over roles such as immune response (e.g., B-and T-cells). This process is called hematopoiesis. We follow the steps of [1] to first split samples (i.e., cells) into the three different lineages (paths of differentiation), and train separate models for each (see App. B.7). We will here focus on the analysis of the Erythroid lineage.As before, DASH yields highly sparse (95%) networks, the most sparse among all competitors (App. Tab. 9). IMP shows similarly strong sparsification as DASH while PathReg achieves much less sparsity (14%). In terms of performance, all methods achieve similar MSE of predicted gene expression dynamics on the test set, meaning that even though much sparser, both DASH and IMP predict equally well as an order of magnitude more dense network. For this data, there are no gold-standard experiments for regulatory relations available, we hence focus on analysing the network topology. From the literature, we would expect sparser networks to be better align with biology [6] .DASH indeed shows the lowest out-degree in the inferred regulatory network, less than half of what IMP recovers. PathReg shows an order of magnitude larger average out-degree. To confirm the biological plausability, we again do a pathway analysis. DASH seems to find significant enrichment in biologically relevant pathways (App. Fig. 7 ) that can directly be linked to hematopoiesis, such as heme signaling or RUNX1 regulates differentiation of hematopoietic stem cells, which neither BioPrune nor SparseFlow-the only other method yielding a proper sparse model-could recover. 
discussion & conclusion: 
 We considered the problem of identifying sparse neural networks in the context of interpretability with a focus on the application to gene regulatory systems modeling. In domains such as biology and contexts when the true underlying systems are sparse, interpretability is key for experts, rendering the use of the common over-parametrized and complex neural network architectures difficult. Although NNs do not directly encode e.g. the regulatory relationship between genes its deep architecture is necessary to model complex functional relationships while ensureing stable learning. Yet, we can ensure that Recent advances in neural network pruning, such as those around the Lottery Ticket Hypothesis [18] , promise sparse and well-performing models, yet, hardness results prove finding optimally sparse models to be challenging [43] , which is also reflected by recent benchmarking results [16] . Our experiments confirmed that general pruning strategies provide sub-optimal sparsity, moreover, the underlying biological relationships are not properly reflected in the model. We proposed to guide pruning by domain knowledge, leveraging existing prior information to improve the interpretability and meaningfulness of pruned models. 
safeguards: 
 Question: Does the paper describe safeguards that have been put in place for responsible release of data or models that have a high risk for misuse (e.g., pretrained language models, image generators, or scraped datasets)?Answer: [NA] Justification: Not applicable Guidelines:• The answer NA means that the paper poses no such risks.• Released models that have a high risk for misuse or dual-use should be released with necessary safeguards to allow for controlled use of the model, for example by requiring that users adhere to usage guidelines or restrictions to access the model or implementing safety filters. • Datasets that have been scraped from the Internet could pose safety risks. The authors should describe how they avoided releasing unsafe images. • We recognize that providing effective safeguards is challenging, and many papers do not require this, but we encourage authors to take this into account and make a best faith effort.12. Licenses for existing assets Question: Are the creators or original owners of assets (e.g., code, data, models), used in the paper, properly credited and are the license and terms of use explicitly mentioned and properly respected?Answer: [No] Justification: We properly cite the work for the compared methods and all data used within the scope of this work. We also checked that each license is applicable. But we do not list each license and terms explicitly.Guidelines:• The answer NA means that the paper does not use existing assets.• The authors should cite the original paper that produced the code package or dataset.• The authors should state which version of the asset is used and, if possible, include a URL. • The name of the license (e.g., CC-BY 4.0) should be included for each asset.• For scraped data from a particular source (e.g., website), the copyright and terms of service of that source should be provided. • If assets are released, the license, copyright information, and terms of use in the package should be provided. For popular datasets, paperswithcode.com/datasets has curated licenses for some datasets. Their licensing guide can help determine the license of a dataset. • For existing datasets that are re-packaged, both the original license and the license of the derived asset (if it has changed) should be provided. • If this information is not available online, the authors are encouraged to reach out to the asset's creators. 13. New Assets Question: Are new assets introduced in the paper well documented and is the documentation provided alongside the assets? Answer: [NA] Justification: The primary contribution of this work is foundational, yet we do provide code and data-generating scripts as supplementary material.Guidelines:• The answer NA means that the paper does not release new assets.• Researchers should communicate the details of the dataset/code/model as part of their submissions via structured templates. This includes details about training, license, limitations, etc. • The paper should discuss whether and how consent was obtained from people whose asset is used.  
 
 Appendix: A.1 Synthetic data 8 : Balanced accuracies for experiments on yeast data. Balanced accuracy is based on 1) transcription factor binding ChIP-seq available for this data [23] and 2) A TF perturbation network created by Hackett et al. based on their TF perturbation experiments [21] . * marks our suggested baselines and method, ✓ marks methods that use prior information for sparsification.Table 9 : Results on Hematopoesis data for the Erythroid lineage. We give sparsity of pruned model and test MSE on predicted gene expression dynamics. No reference gene regulatory network is available to compute the accuracy of the recovered network we hence resort to reporting the average out-degree of nodes in the recovered network. DASH found an optimal λ-value of (0.80, 0.80). * marks our suggested baselines and method, ✓ marks methods that use prior information for sparsification. B Supplementary methodsThe purpose of simulation based data is so that the the underlying dynamical system that produced the this gene expression was known. Do this end, we closely follow the steps outlined by the simulation pipeline provided by [27] to generate reliable synthetic time-series gene expression data from two ground truth gene regulatory networks G 350 and G 690 consisting of 350 and 690 genes, respectively.The pipeline adapts SimulatorGRN [4] to generate from two synthetic S. cerevisiae gene regulatory systems (SIM350 and SIM690, consisting of 350 and 690 genes respectively). For every noise setting ∈ {0%, 5%, 10%}, the connectivity structure of each in silico system is used to synthesize 160 noisy expression trajectories for each gene across t ∈ T = {0, 2, 3, 7, 9}. We split up the trajectories into training (88%), validation (6% for tuning λ), and testing (6%). Since the average simulated expression value is ≈ 0.5, adding Gaussian noise of N (0, σ 2 ) using σ ∈ {0, 1 40 , 1 20 } corresponds roughly to average noise levels of {0%, 5%, 10%}.Since the number of genes k in each problem is different, the number of neurons m in PHOENIX's hidden layer is chosen to roughly scale with this k according to the original paper [27] : • SIM350: k = 350, m = 40For initialization values for each of NN sums , NN prods , NN Σcombine , and NN Πcombine , as well as that of υ i s we choose the default provided by the PHOENIX implementation [27] . The ODESolver (dopri5) and optimizer (Adam) are also chosen as the PHOENIX defaults across all experiments.We use iterative pruning schedules that are initially very aggressive and then become much more gradual. We found this approach to achieve high sparsity without adversely affecting the training dynamics (and subsequently the validation performance). Π | we approach again similarly, with the important modification that the gene-specific multipliers (from Section B.10.2) are row-wise multiplied into the weight matrices prior to normalization. This allows the effect of gene multipliers to appropriately be considered when performing pruning.The learning rate is used as the PHOENIX default of 10 -3 . We reduce the learning rate by 10% every 3 epochs, unless the validation set performance shows reasonable improvement. Importantly, we reset the learning rate back to 10 -3 immediately after a pruning step is completed, thereby allowing the newly sparsified model to start learning with a higher learning rate.We train for up to 500 epochs on an AWS c5.9xlarge instance, where each 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. Hence this training procedure automatically finds an optimal sparsity level using the validation set.As mentioned in Section 3, DASH can leverage prior matrices P and C to inform its pruning score. We use the following in our experiments:• SIM350 and SIM690:for synthetic experiments we choose P = A σ % to be noisy/corrupted versions (see B.3.1) of the adjacency matrices of ground truth networks G 350 and G 690 to reflect that transcription factor binding to target genes can itself be a noisy process in real life. A 1 A σ % represents prior knowledge of an interaction existing between two genes, and a 0 represented no interaction. -For C we use the outer product C = A σ % (A σ % ) ⊺ , to represent prior knowledge of coregulation. We again applied the corruption/missepecification procedure from B.3.1 so that C is also noisy.• Breast cancer data:-For the prior domain knowledge, we set P = W 0 , where W 0 was a motif map derived from the human reference genome, for the breast tissue specifically, which we obtained through GRAND [3] . W 0 is a binary matrix with W 0i,j ∈ {0, 1} where 1 indicates a likely occurence of a TF sequence motif in the promoter of the target gene, and hence indicating a putative interaction. More simply put, it indicates that whether there is (likely) a binding interface for the protein close to the target gene. -Based on information from the STRING database [59] , we obtained a protein-protein interaction matrix (PPI) which could use to operationalize our C matrix, since a PPI is a again a binary matrix that is suggestive of which transcription factors have combined (or coregulatory) effects. In a nutshell, the STRING database is a graph with proteins as vertices and knowledge about interactions between two proteins in the graph specifying edges. We set an entry P i j to 1 if the experimental evidence score on the edge between protein i and j is larger than 0.6, and set P i j to 0 otherwise.-For the prior domain knowledge, we followed a similar strategy as the breast cancer analysis, and set P and C based on the motif map and PPI matrix used in [62] . We appropriately subsetted P and C to only be limited to the k = 529 genes that were selected by the PathReg authors [1] in the analysis. • Yeast data:-For prior domain knowledge model we set P to reflect the regulatory network structure of a motif map. The map is based on predicted binding sites for 204 yeast transcription factors (TFs) [23] . These data include 4360 genes with tandem promoters. 3551 of these genes are also covered on the yeast cell cycle gene expression array. 105 total TFs in this data set target the promoter of one of these 3551 genes. The motif map between these 105 TFs and 3551 target genes provides the adjacency matrix A of 0s and 1s, representing whether or not a prior interaction is likely between TF and gene. -We set C to be the PPI matrix used for the same data in the PANDA paper [20] . We appropriately subsetted C to only be limited to the k = 3551 genes that were in the data. ) with corresponding adjacency matrix A σ % . We additionally performed sensitivity analyses using σ % ∈ {20%, 40%} to investigate the effect of even higher levels of prior corruption. We then used A σ % to obtain "corrupted" P and C as described in B.3.B.The choice of λ = (λ 1 , λ 2 ) is important for optimally combining prior information with model weights. Hence we implement a K-fold cross validation approach to choose λ. Test set performance is measured as the mean squared error between predictions and held-out expression values in the test set.To validate biological alignment of trained and sparsified models, we extracted GRNs from each models (as explained in B.10.4), and compared back to the validation networks. Specifically, once we extracted a GRN from the trained model, we looked at how well 0s vs non-zeros in that network aligned with 0s vs non-zeros in the validation network. Our comparison metric was balanced accuracy, which is the average of the true positive and true negative rates. The validation networks were as follows:• SIM350, SIM690: We used the ground truth networks G 350 and G 690 .• Breast cancer data: We used ChIP-seq data from the MCF7 cell line (breast cancer) in the ReMap2018 database [10] to create a validation network of TF-target interactions. • Bone marrow data: As also noted by [1] , validation network was not available, so we resorted to the pathway analysis. • Yeast data: We used two kinds of validation networks 1. ChIP-seq data [23] to create a network of TF-target interactions, and used this as a validation network to test explainability. The targets of transcription factors in this ChIP-chip data set were filtered using the criterion p < 0.001.Supposing we only have access to prior knowledge in the form of putative prior effect sizes between the n inputs and o outputs P ∈ R o×n . Then, for an NN with L -1 layers W 1 ∈ R m1×n , W 2 ∈ R m2×m1 , . . . , W L ∈ R o×m L ,we can adopt a strategy where we consider the pruning scores to be fixed for all but one layer.Since the productW (t) L . . . W (t) 2 • W (t) 1∈ R o×n represents the overall flow of information from inputs to outputs at epoch t, we surmise that Ω (t)L . . . Ω (t) 2 • Ω (t)1 ∈ R o×n should reflect P . We can thus prune as follows:1. Starting with the last layer, we fix the pruning scores of all other layers and compute as follows:Ω (t) L := (1 -λ L ) |W (t) L | + λ L P • PInv R Ω (t) L-1 . . . Ω (t)2 Ω (t) 1.2. For the middle layers l ∈ {2, 3, . . . , L -1}, we do:Ω (t) l := (1-λ l ) |W (t) l |+λ l PInv L Ω (t) L . . . Ω (t) l+2 Ω (t) l+1 •P •PInv R Ω (t) l-1 . . . Ω (t)2 Ω(t) 1.3. The first layer can be pruned using:Ω (t) 1 := (1 -λ 1 ) |W (t) 1 | + λ 1 PInv L Ω (t) L . . . Ω (t)3 Ω (t) 2• P .The original data set comes from a cross-sectional breast cancer study (GEO accession GSE7390 [13] ) consisting of microarray expression values for 22000 genes from 198 breast cancer patients, that is sorted along a pseudotime axis. We note that the same data set was also ordered in pseudotime by [58] in the PROB paper. For consistency in pseudotime inference, we obtained the same version of this data that was already preprocessed and sorted by PROB. We normalized the expression values to be between 0 and 1. We limited our analysis to the genes that had measurable expression and appeared in the aforementioned motif map and PPI matrices. This resulted in a pseudotrajectory of expression values for 11165 genes across 186 patients. We removed a contiguous interval of expression across 8 time points for testing (5%), and split up the remaining 178 time points into training (170, 90%) and validation for tuning λ prior (8, 5%).GPR files were downloaded from the Gene Expression Omnibus (accession GSE4987 [50] ), and consisted of two dye-swap technical replicates measured every five minutes for 120 minutes. Each of two replicates were separately ma-normalized using the maNorm() function in the marray library in R/Bioconductor [64] . The data were batch-corrected [29] using the ComBat() function in the sva library [33] and probe-sets mapping to the same gene were averaged, resulting in expression values for 5088 genes across fifty conditions. Two samples (corresponding to the 105 minute time point) were excluded for data-quality reasons, as noted in the original publication, and genes without motif information were then removed, resulting in an expression data set containing 48 samples (24 time points in each replicate) and 3551 genes.The data is originally from [42] (GEO accession code = GSE194122). The cleaning, preprocessing, and pseudotime analysis and was appropriately performed by [1] in the PathReg paper, and made publicly available, allowing us to access the processed version. Importantly, [1] split up the data into 3 different lineages (Erythroid, Monocyte, and B-Cell), and we fit a separate PHOENIX model on each lineage. The set contains 5 separate batches of data for each lineage, we used 1 for training (batch S1D2), 1 for validation (batch S1D1) and 3 for testing (batches S1D1, S2D4, and S3D6).We followed very closely the steps below from the Methods section of the PHOENIX paper [27] in order to compute pathway scores, with the only difference that we compute scores between different sparsification strategies Pru.Given M Pru a PHOENIX model trained on a dataset consisting of k genes, and sparsified using the pruning strategy Pru (for Pru ∈ {DASH, IMP, C-NODE, . . . , PathReg}), we performed perturbation analyses to compute gene influence scores IS Pru,j . We randomly generated 200 initial (t = 0) expression vectors via i.i.d standard uniform sampling {g(0) r ∈ R k } 200 r=1 . Next, for each gene j in M Pru , we created a perturbed version of these initial value vectors {g j (0) r } 200 r=1 , where only gene j was perturbed in each unperturbed vector of {g(0) r } 200 r=1 . We then fed both sets of initial values into M Pru to obtain two sets of predicted trajectories {{ g(t) r ∈ R k } t∈T } 200 r=1 and {{ g j (t) r ∈ R k } t∈T } 200 r=1 across a set of time points T . We calculated influence as the average absolute difference between the two sets of predictions, that represented how changes in initial (t = 0) expression of gene j affected subsequent (t > 0) predicted expression of all other genes in the Pru-dimensional systemIS Pru,j = 1 200 200 r=1 1 |T | t∈T t̸ =0 1 k k i=1 i̸ =j g i (t) r -g i j (t) r .Having computed gene influence scores IS Pru,j for each gene j in each dynamical system of dimension k genes sparsified with method Pru, we translated these gene influence scores into pathway influence scores. We used the Reactome pathway data set, GO biological process terms, and GO molecular function terms from MSigDB [36] , that map each biological pathway/process, to the genes that are involved in it. For each system sparsified by Pru, we obtained the pathway (p) influence scores (PS Pru,p ) as the sum of the influence scores of all genes involved in pathway pPS Pru,p = j∈p IS Pru,j .We statistically tested whether each pathway influence score is higher than expected by chance using empirical null distributions. We randomly permuted the gene influence scores across the genes to recompute "null" values PS 0 Pru,p . For each pathway, we performed Q = 1000 permutations to obtain a null distribution {PS 0 Pru,p,q } Q q=1 that can be compared to PS Pru,p . We could then compute an empirical p-value as p = 1Q Q q=1 I PS 0Pru,p,q >PSPru,p , where I is the indicator function. Finally, we used the mean (µ 0(Pru,p) ) and variance (σ 2 0(Pru,p) ) of the null distribution {PS 0 Pru,p,q } Q q=1 to obtain and visualize pathway z-scores that are now comparable across pathways (p) and sparsification strategies (Pru)EQUATION.B.9 Implementation details for other sparsification strategies on the PHOENIX architectureAs discussed in Section 3, IMP can be operationalized as a special case of DASH by setting λ 1 = λ 2 = 0.This is simply the PHOENIX model equipped with the prior-informed loss term. This loss-term in the original PHOENIX paper [27] is inspired by Physics-informed neural networks (PINNs).Once a PHOENIX model (trained including the prior-informed loss term, i.e. the PINN term) is fully trained (without any pruning), we inspect the trained model and pruned the lowest p% of parameters in each of W Σ , W Π , U Σ , U Π based on on the normalized weights (see B.2) to 0. We then fine-tune (i.e retrain without training the pruned parameters) this p% sparsified model and calculate its performance on the validation set. We repeat this process for a grid of values for p ∈ {0.50, 0.75, 0.83, 0.87, 0.90, 0.92, 0.95, 0.97, 0.99}. The validation set can then inform the best value of p. We repeated this entire procedure 3 times, so that we could apply the 1 standard error rule [24] and choose the optimal p as the sparsest fine-tuned model whose validation MSE is within 1 standard error of lowest obtained average validation MSE.C-NODE, PathReg, and L 0 implementations were obtained from the code associated with the PathReg paper [1] . We adapted the code so that the base NN architecture was exactly that of the PHOENIX model, including an implementation of the gene-specific multipliers (from Section B.10.2). Finally, we tuned the relevant parameters λ 0 and λ 1 in the objective function using the validation set.The code for DST was obtained from: https://github.com/junjieliu2910/DynamicSparseTrainingThe following are adapted from [27] and provided here for the reader's convenience.NeuralODEs [8] learn dynamical systems by parameterizing the underlying derivatives with neural networks: dg(t) dt = f (g(t), t) ≈ NN Θ (g(t), t). Given an initial condition g(t 0 ), the output g(t i ) at any given time-point t i can now be approximated using a numerical ODE solver of adaptive step size:g(t 1 ) = g(t 0 ) + t1 t0 NN Θ (g(t), t) dt.A loss function L g(t 1 ) ; g(t 0 )+ t1 t0 NN Θ (g(t), t) dt is then optimized for Θ via back propagation, using the adjoint sensitivity method [7] to carry the backpropagation through the integration steps of the ODESolver.PHOENIX models gene expression dynamics using NeuralODEs. Notably, for an expression vector of r genes, PHOENIX models both additive and multiplicative regulatory effects using two parallel linear layers with m neurons each: NN sums (with weights W Σ ∈ R m×r , and biases b Σ ∈ R m ) and NN prods (W Π ∈ R m×r , b Π ∈ R m ). Here, NN sums and NN prods are equipped with activation functions that model the Hill equationϕ Σ (x) = x -0.5 1+ | x -0.5 | and ϕ Π (x) = log ϕ Σ (x) + 1 .The Hill equation is a classical formula in biochemistry that models molecular binding in dependence of concentration. This results in outputs of the two parallel layersc Σ (g(t)) = W Σ ϕ Σ (g(t)) + b Σ and c Π (g(t)) = exp •(W Π ϕ Π (g(t)) + b Π ).As shown above, ϕ Π (x) yields the output of NN prods in the log space and is subsequently exponentiated in c Π to represent multiplicative effects in the linear space. The outputs c Σ (g(t)) and c Π (g(t)) are then separately fed into two more parallel linear layers NN Σcombine (with weights U Σ ∈ R r×m ) and NN Πcombine (U Π ∈ R r×m ), respectively. The outputs of NN Σcombine and NN Πcombine are summed to obtainc ∪ (g(t)) = U Σ c Σ (g(t)) + U Π c Π (g(t)). Finally, PHOENIX includes gene-specific multipliers υ ∈ R r for modeling steady states of genes that do not exhibit any temporal variation dgi(t) dt = 0, ∀t . Accordingly, the output for each gene i in c ∪ (g(t)) is multiplied with ReLU(υ i ) in the final estimate for the local derivativeNN Θ (g(t), t) = ReLU(υ) ⊙ c ∪ (g(t)) -g(t) .Although PHOENIX achieves some sparsity in its weight matrices (W Σ , W Π , U Σ , U Π ) without any external sparsification strategy, the achieved sparsity level is, however, at most 12% (see Figure 2 , and Tables 1, 3, 4).PHOENIX has the option to promote the NeuralODE to flexibly align with structural domain knowledge, while still explaining the observed gene expression data. This is operationalized via a modified loss function L mod g(t 1 ), g(t 1 ) =τ loss based on matching observed gene expression data L g(t 1 ) ; g(t 0 ) + t1 t0 NN Θ (g(t 1 ), t) dt+ (1 -τ )loss based on matching domain-knowledge L P * g(t 1 ) ; NN Θ (g(t 1 ), t) that incorporates the effect of any user-provided prior model P * , using a tuning parameter τ , and the original loss function L(x, x). PHOENIX implements P * as a simple linear model P * γ = A • γγ, where A is the adjacency matrix of likely connectivity structure based on prior domain knowledge (such as experimentally validated interactions) with A ij ∈ {+1, -1, 0} representing an activating, repressive, or no prior interaction, respectively.For synthetic experiments, we used the simple linear model: P * γ = A σ % • γγ, where we chose A σ % to be noisy/corrupted versions of the adjacency matrices of ground truth networks G 350 and G 690 (details in B.3.1). We set activating and repressive edges in A σ % to 1, while "no interaction" was represented using 0.We start with PHOENIX's prediction for the local derivative given a gene expression vector g(t) ∈ R r in an r-gene system:dg(t) dt = ReLU(υ) ⊙ W ∪ {c Σ (g(t)) ⊕ c Π (g(t))} -g(t), wherec Σ (g(t)) = W Σ ϕ Σ (g(t)) + b Σ and c Π (g(t)) = exp •(W Π ϕ Π (g(t)) + b Π )A trained PHOENIX model encodes interactions between genes primarily within the gene-specific multipliers υ ∈ R r , and the weight parameters from its neural network blocks W Π , W Σ ∈ R m×r and W ∪ ∈ R r×2m . This inspired an efficient means of projecting the estimated dynamical system down to a gene regulatory network (GRN) G n . In particular a matrix D ∈ R r×r is calculated, where D ij approximated the absolute contribution of gene j to the derivative of gene i's expressionD = W ∪ W Σ W Π .Gene-specific multipliers υ are applied, before adapting the marginal attribution approach described by Hackett et al. [21] . This resulted in the dynamics matrix D where D ij was scaled according to the relative contribution of gene j to the rate of change in gene i's expression:D ij = υ i D ij n j ′ =1 |υ i D ij ′ | .With its simple yet powerful architecture, PHOENIX provides an ideal base setting for applying and testing the merits of the discussed neural network sparsification strategies (Section 2), including DASH. We discuss the DASH implementation here and provide implementation details for other sparsification strategies in Appendix B.9.For DASH, we follow the steps of the L = 2 case in Section 3. Note that in PHOENIX we have k = r, since we are modeling how a set of k genes (inputs) affects each others derivatives (i.e. k outputs). We apply the steps from Section 3 on each parallel member of the first {W Σ , W Π } and second {U Σ , U Π } layers. For the first layer, we leverage prior knowledge of transcription factor coregulation (often available in the form of protein-protein interaction matrices) to formulate C ∈ R k×k . This is then used to calculate pruning scores Ω Σ , Ω Π for W Σ , W Π . Similarly, for U Σ , U Π in the second layer, we utilize P ∈ R r×k which are easily obtainable motif map matrices encoding prior knowledge of transcription factor binding sites around genes to calculate pruning scores Ψ Σ , Ψ Π .EQUATIONΠ , ΨΣ , ΨInputs: weightsEQUATIONΠ ,UΣ , UΠ ; priors C, P ; epoch t; previous scores Ω (t-1) Σ,Ω (t-1) Π ; tuning λ 1 , λ 2 • Normalize W (t) Σ to get | W (t) Σ | (see Appendix B.2) • Similarly obtain | W (t) Π |, | U (t) Σ |, and | U (t) Π | • If t = 0 -Initialize Ω (t)Σ and Ω (t)Π randomly with values from a standard Gaussian • Else-Ω (t) Σ := (1 -λ 1 ) |W (t) Σ | + λ 1 PInv L Ω (t-1) Σ ⊺ • C -Ω (t) Π := (1 -λ 1 ) |W (t) Π | + λ 1 PInv L Ω (t-1) Π ⊺ • C • Ψ (t) Σ := (1 -λ 2 ) |U (t) Σ | + λ 2 P • PInv R Ω (t) Σ • Ψ (t) Π := (1 -λ 2 ) |U (t) Π | + λ 2 P • PInv R Ω (t) ΠB.12 Ablation study using a plain MLP instead of PHOENIX as the base modelWe believe that DASH should remain performant even when using a base model that is different from PHOENIX. To be more explicit, the base model would give a new expression for NN Θ (g(t), t) in B.10.2.For PHOENIX as base model we have: NN Θ (g(t), t) = ReLU(υ) ⊙ c ∪ (g(t)) -g(t) .For alternative base model we have: NN Θ (g(t), t) = alternative base model(g(t))-g(t)In our ablation experiments, with used a simple two-layer MLP with ELU activation functions as the alternative base model . The choice of ELU activation is motivated by the PathReg paper [1] . We chose the number of hidden neurons such that the total number of parameters was comparable between the MLP and the PHOENIX base models. We tested a few sparsification strategies against DASH on this new base model, for both synthetic data and the bone marrow data. Guidelines:• The answer NA means that the paper has no limitation while the answer No means that the paper has limitations, but those are not discussed in the paper. • The authors are encouraged to create a separate "Limitations" section in their paper.• The paper should point out any strong assumptions and how robust the results are to violations of these assumptions (e.g., independence assumptions, noiseless settings, model well-specification, asymptotic approximations only holding locally). The authors should reflect on how these assumptions might be violated in practice and what the implications would be. • The authors should reflect on the scope of the claims made, e.g., if the approach was only tested on a few datasets or with a few runs. In general, empirical results often depend on implicit assumptions, which should be articulated. • The authors should reflect on the factors that influence the performance of the approach.For example, a facial recognition algorithm may perform poorly when image resolution is low or images are taken in low lighting. Or a speech-to-text system might not be used reliably to provide closed captions for online lectures because it fails to handle technical jargon. • The authors should discuss the computational efficiency of the proposed algorithms and how they scale with dataset size. • If applicable, the authors should discuss possible limitations of their approach to address problems of privacy and fairness. • While the authors might fear that complete honesty about limitations might be used by reviewers as grounds for rejection, a worse outcome might be that reviewers discover limitations that aren't acknowledged in the paper. The authors should use their best judgment and recognize that individual actions in favor of transparency play an important role in developing norms that preserve the integrity of the community. Reviewers will be specifically instructed to not penalize honesty concerning limitations.In case studies on gene regulatory dynamic inference, a key task in molecular biology with high relevance in cancer research, we showed based on simulated as well as real world data that our method, DASH, in contrast to a wide range of state-of-the-art methods, is able to recover neural networks that are both very sparse and at the same time biologically meaningful, allowing for direct extraction of a sparse gene regulatory network. On real data, DASH not only better aligns with gold-standard experimental evidence of regulatory interactions, but also uniquely reflects the dataspecific biological pathways, which can be used by domain experts to generate new insights.. It thus serves as a proof of concept that in critical domains, where interpretability is essential and domain knowledge exists, pruning can be heavily improved by alignment with prior knowledge. While our guided pruning approach is in principle agnostic to the type of neural network and task, we here focused on a specific case study that we deemed important. In the future, it would be interesting to apply DASH to different cases and domains, including other biomedical tasks, but also to physics or material sciences, where interpretability is also key and domain knowledge exists in the form of physical constraints and models. For any application, an important consideration to apply DASH is on the one hand the availability of prior knowledge, but on the other hand its quality; while we show here that even with incomplete and noisy prior knowledge we receive good results, factually wrong priors could steer the solution towards a wrong model. We hence assume that DASH will be of primary use in classical hard sciences mentioned above, with priors that stood the test of time over several decades. Another line of future work includes different architectural designs, such as convolution or attention mechanisms, where input-output relationships are less straightforward to project to across several layers.In summary, we make a case for pruning informed by domain knowledge and provide evidence that such approaches can massively improve sparsity along with domain specific interpretability. 


Decision: The scores from reviewers are 5, 5, 4, 3, 5. The reviews make many valuable points, but the authors have responded very well, so as the area chair, I believe that the numerical scores deserve to be higher, and this submission should be accepted.

The authors' responses include additional experiments with positive results, and a novel biological finding that is plausible and may have clinical utility.

Here is a high level question for the authors to discuss in the final version. They write "While pruning can simplify deep neural network architectures ... in the context of gene regulatory network inference, state-of-the-art techniques struggle with biologically meaningful structure learning." There are two different meanings for the word "network" here: the DNN versus the actual biological mechanisms. An ideal algorithm would output the actual regulatory structure, not just a pruned DNN that can answer questions about this. A DNN is an indirect representation of the actual mechanisms. The question is, to what extent can the literal regulatory structure be represented and discovered directly?

Another high-level suggestion: Mention in the abstract, and discuss more in the text, applicability of the DASH method outside molecular biology, expanding lines 61, 62 and 326, 327. In which domains does or does not prior knowledge tend to be of the type required by DASH?

Reviews: review: paper_summary: paper_summary: The paper presents DASH (Domain-Aware Sparsity Heuristic), a new framework for pruning neural network models by incorporating domain-specific knowledge. The primary goal is to improve the interpretability and biological relevance of models used for gene regulatory network (GRN) inference. Traditional pruning methods often fail to reflect biologically meaningful structures, leading to less interpretable models.

### Key Contributions:
- Introduction of DASH: A pruning method that uses both learned weights and prior domain knowledge to iteratively score and prune network parameters, ensuring the resulting models are both sparse and biologically meaningful.
- Improved Model Interpretability: By guiding pruning with structural information about gene interactions, DASH produces models that align better with biological insights compared to traditional methods.
- Experimental Validation:
On synthetic data, DASH outperformed general pruning methods in accurately recovering the underlying GRN.
On real-world gene expression data, DASH identified biologically relevant pathways that other methods missed.
- Robustness to Noise: The framework maintained model robustness even in noisy data environments, showcasing its practical utility.
strength_weakness: strength_weakness: ### Testing Set Size
One notable limitation of the paper is the relatively small size of the testing set, which comprises only 6% of the total data. This raises concerns about the generalizability and robustness of the reported results. A small testing set can lead to overfitting and may not adequately capture the model's performance on unseen data, especially in cases where the performance differences between methods are quite small.

### Complexity of the Methodology
While the methodology is detailed and thorough, it may appear overly complex for some readers. Simplifying certain aspects or providing more intuitive explanations could improve accessibility and comprehension.

### Biological Validation
The paper could further substantiate the claim that DASH improves biological relevance by involving domain experts and providing concrete examples.

Suggestions for Improvement:

- **Expert Validation**: Include validation from domain experts to verify the biological significance of the results.
- **Concrete Examples**: Provide specific examples where DASH's results align with known biological mechanisms or lead to new hypotheses.
questions: questions: ### 1. Format and Numerical Examples of $P$:
It appears that using prior knowledge alone (i.e., the "BioPrune" baseline) can perform reasonably well for both synthetic and real datasets. However, the description of the domain knowledge used is not entirely clear. It seems that different sets of domain knowledge are employed for different datasets. While we understand the generic form of $P$, can you provide numerical examples for both synthetic and real datasets to illustrate the exact format of $P$?

### 2. Parameter Tuning
How are the lambda values $\lambda$ determined for each dataset? Are they manually tuned, or is there an automated process for selecting these values? Can you provide the results and details on the validation process for tuning these parameters?

### 3. Generalizability (to other neural network, to other domain)
The paper demonstrates the effectiveness of DASH on gene regulatory networks. How generalizable is this approach to other domains or types of neural networks (CNN, transformer)? Also, what if there are activation functions in the current used models (MLPs)? Have there been any preliminary tests or considerations for applying DASH to fields such as physics, material science, or other areas of computational biology?

### 4. Comparative Analysis
While DASH is compared with several pruning methods, are there any other recent state-of-the-art techniques or methods that should be included in the comparison? How does DASH stand against the very latest advancements in neural network sparsification?
limitations: limitations: The paper does discuss certain limitations, but more details could be added to enhance transparency and improve the assessment of the work:

**Scope of Application**: The generalizability of DASH to other domains outside gene regulatory networks isn't thoroughly explored. Adding more discussion on the potential limitations when applying DASH to different types of neural networks or scientific fields could be beneficial.

**Bias and Fairness**: The datasets used in biological research can sometimes contain biases, such as underrepresentation of certain populations. Discussing how the method handles such biases and ensuring that the models do not reinforce existing disparities is important. Suggesting ways to incorporate fairness checks and balance the representation in datasets would strengthen the ethical considerations.
paper_summary: paper_summary: The proposed DASH method underscores the importance of interpretability in network pruning for biological discoveries, emphasizing the need for alignment with domain knowledge. Using both synthetic and real data, DASH demonstrates superior performance beyond baselines and offers insights into biological systems.
strength_weakness: strength_weakness: Depending on the size and complexity of the domain knowledge, DASH might be hard to apply to large-scale or highly complex networks.
questions: questions: 1. In Figure 3, what is the percentage of true relationships comparing DASH and PINN + MP?
2. It’s not surprising that incorporating the known domain knowledge enhance the interpretability and learned more meaningful dynamics. Moreover, how does the learned gene regulatory dynamics itself can benefit the downstream disease related outcome predictions?
3. The paper mentioned that DASH has the better quality of inferred (new) knowledge. Is there any discussion or literatures to support the inferred “new” knowledge?
4. How about making lambda a learnable parameter to mimic the complex and dynamic biological system?
limitations: limitations: The authors clearly addressed the limitations.
paper_summary: paper_summary: Gene regulatory network inference is an important, but difficult problem.

The manuscript explores a novel approach to build domain knowledge into a general NODE model for this problem via pruning. The approach could work in other areas.
strength_weakness: strength_weakness: It’s very hard to match real data biases for the prior GRN when studying simulated data. I’m struggling to understand whether the noise % refers to this GRN prior or the expression data. 5% noise for the prior seems completely unrealistic and would render these experiments irrelevant to real data analysis.

Analyzing real gene expression data, the authors use priors based on TF motifs mapped to promoter sequences, which enable DASH to identify genes with true TF ChIP binding with apparently high accuracy. Regulatory network inference method evaluations generally find that genes with nearby TF binding events do not correspond well to genes whose expression changes upon TF perturbation. Evaluation with TF perturbation data is the gold standard here and would make for a more compelling DASH evaluation.
questions: questions: In Fig 2, what does it mean for DASH to have 2x fewer parameters than BioPrune but achieve greater GRN accuracy? Isn’t BioPrune simply using the GRN?

In Fig 2, what is ”Base model”? 

It would be helpful to understand why PHOENIX with the same GRN prior information, but using regularization instead of pruning, falls behind the DASH pruning strategy.

I’m not really sure what I’m supposed to be able to understand from Figure 3. I can’t see any of the specific matrix entries well enough to compare them.
limitations: limitations: Yes
paper_summary: paper_summary: The authors propose a network pruning approach which is guided by prior biological domain knowledge, which they call Domain-Aware Sparsity Heuristic (DASH). They aim to obtain highly sparse networks, which align with known biology and gene regulatory dynamics. To do so they propose computing pruning scores which combine learned weights with prior domain knowledge, balanced by a tunable parameter $\lambda$. The authors apply DASH to PHOENIX, a neuralODE designed to model gene regulatory dynamics and conduct experiments on four datasets, one synthetic and three real world. They also compare DASH to other existing pruning techniques, such as $L_0$, PathReg and SparseFlow. The experiments compare DASH to these techniques in terms of sparsity, accuracy in recovering known biological relationships and MSE of predicted gene expression values, as well as pathway analysis.
strength_weakness: strength_weakness: I think this is relevant and important research, but in my opinion there are four main issues to address in this paper: the mathematical formulation for DASH is very confusing, there is an over emphasis of results obtained on synthetic data, as well as a reliance on a metric (balanced accuracy) which might potentially lead to inflating results, and a lack of clarity about potential confounding effects between DASH/PHOENIX. I will detail each of these concerns below. 

1 - Mathematical formulation of DASH, in the case of L=2 (at the top of p.5 - I have no numbering on the PDF here for some reason): 

  - In line 139 you define your weight matrix as  $W \in \mathbb{R}^{k \times r}$, where $k$ is your input corresponding to your genes - yet a paragraph later you transpose your notation and start using $W_1\in \mathbb{R}^{m\times k}$, $W_2\in \mathbb{R}^{r\times m}$, which is confusing.
   - You then define an input-input matrix $C \in \mathbb{R}^{n\times n}$ which should be $\mathbb{R}^{k\times k}$ as you defined your inputs a few lines above, you follow by saying that $W_1$ has $n$ inputs and that the matrix product $\Omega_1^{(t)^T} \Omega_1^{(t)}  \in \mathbb{R}^{n\times n}$ -  this doesn't make sense dimensionally. 
   - You posit $\Omega_1^{(t)^T} \Omega_1^{(t)}  \in \mathbb{R}^{n\times n}$  approximates to $C$ without first clearly defining what $\Omega_1^{(t)}$ is and without clearly enunciating why this is most likely a good approximation. 
   - You then define a recurrent relation between $\Omega_1^{(t)}$ and $\Omega_1^{(t-1)}$, which needs an initial condition to be well defined. Furthermore, this temporal dependency is not present in the L=1 case. 
   - Why are you using the left and right pseudoinverses on $\Omega_1^{(t-1)}$ and $\Omega_1^t$ respectively? It would be good to take time to explain this approach in more detail, as this is the central premise of this paper. You don't need to define what the left and right pseudoinverse of a matrix is. You do need to explain why you are using them here and to what end. 
   - The product of $W_2^{(t)} \cdot W_1^{(t)}$ should be $\in \mathbb{R}^{k\times r}$, not $\mathbb{R}^{o\times n}$
       

The definition of DASH is the central premise of the paper, yet there are many notation errors, variables which are not properly defined, and most importantly only vaguely motivates why the proposed substitution makes sense. I suggest you introduce all the notation you will use, even if it seems trivial, and be more careful with errors in notation. It would be great if you expanded on this paragraph by giving more intuitions behind the physical reasoning, as well as explaining the mathematical steps described. 

2 - Over emphasis of results obtained on synthetic data:

A significant portion of results (Table 1, Figure 2 and 3) and discussion focus on synthetic data. Although useful for initial validation and exploration, an in depth examination of the three real-world datasets studied would be more relevant. The data generation process likely oversimplifies the complexities of real regulatory networks, which could lead to an unrealistic assessment of the model's performance.  The results obtained on the real-world datasets are not as compelling as those obtained on the synthetic dataset suggesting DASH is not generalising as well as it could and might be fitting to characteristics of the synthetic dataset. Expanding on the results obtained for the real-world datasets would help address these concerns. For example, the results presented in A.3 Table 6 for the bone marrow data differ significantly from the Sparisity and MSE results obtained for $L_0$, C-NODE and PathReg in the PathReg paper. I would be interested to know more about these differences. Could it be coming from using PHOENIX as a baseline, rather than the NN model used in PathReg?  In my opinion the results on the bone marrow dataset should be in the main body of the work, as it reproduces an experiment from the PathReg paper, permitting a more direct comparison between the methods. 

3/4 - Balanced accuracy to validate biological alignment \& confounding effects between DASH/PHOENIX

The balanced accuracy for biological alignment metric is not clearly defined in the main body of the article, yet is one of the main results in Table 1, 2 and 6, where the authors argue that because DASH is obtaining higher balanced accuracy scores it is better able to recover true biology. They point towards B.5, where they say "To validate biological alignment of trained and sparsified models, we extracted GRNs from each models (as explained in B.10.4), and compared back to the validation networks." In B.10.4 they describe retrieving GRN from the PHOENIX model. This lack of clarity regarding how the biological alignment is measured makes it difficult to evaluate these results. Is it possible that the extraction process itself could be biased towards producing results that align with the prior knowledge used in DASH? This could artificially inflate the balanced accuracy scores for DASH and BioPrune. Could it be there is a confounding effect in the ablation study coming from applying the sparsification methods only to the PHOENIX model? To check this the results would need to be expanded to include a different SOTA neuralODE model to verify DASH works across different models and is not a result of the DASH/PHOENIX combination.
questions: questions: In addition to the questions and concerns I raised above, I have the following questions/comments:

- For clarity, you apply all the pruning  methods mentioned - $L_0$, C-NODE, PathReg, DST, IMP, SynFlow, SparseFlow - to PHOENIX?
    
- What's the difference between PHOENIX with biological regularisation and BioPrune?
    
- In Results, you say PINN+MP and BioPrune are your models, yet they both correspond to post-hoc pruning if I understand correctly. Therefore, neither of them have undergone pruning with DASH. Yet, both perform strongly, with PINN+MP outperforming DASH on MSE in Table 1 and seemingly doing a better job at reconstructing ground truth relationships in Figure 3. Why do you think that is?
    
- The main difference between models in Table 2 seems to be the Balanced Accuracy, but is it not expected for DASH to recover similar levels of sparsity to BioPrune given you are "forcing" convergence towards $C$ in your model and hence achieving a higher Balanced Accuracy?
    
- If one of the takeaways of DASH is the reduced training time and/or memory usage it would be good to see training/inference time or GFLOPs.
    
- Can this be used to generate new biological knowledge? Some discussion on this topic would enrich the conclusions of this paper in my opinion.
limitations: limitations: I do not think the potential limitations of DASH are adequately discussed, see my comments above.
paper_summary: paper_summary: The authors proposed an interesting model called DASH. The model can recover the gene interact relationships in a sparsity way and make a high accuracy.
strength_weakness: strength_weakness: The paper is not well organized. E.g., The explanation of the equations and experiments. (see questions)

More biological background should be included for better understanding.
questions: questions: Can the knowledge be updated by the model? 
Can you model and simulate the different activities of the gene regulatory? E.g. activate the genes or deactivate the genes. 
How to align the parameter n with \Omiga in DASH for H=2? What is its biological meaning and align with the prior knowledge? 
How to measure the accuracy? Do you consider whether the reconstructed relations exist in biology but have not been discovered? E.g. regulatory network rewiring in the cells
How to do the pathway analysis? How to get the enriched pathways?

You use Fig. 1 and Fig. 3 but Figure 2. Please ensure they are in the same form.
limitations: limitations: The authors don't set a limitation section in the paper. The authors proposed a method that simulates the bio system. However, the biology system is more complex. More things can be discussion in the future.

Comments: title: Adding biological background
Response to Authors Rebuttal
Additional experiments with more noise in prior
Clarification ChIP-seq and perturbation experiments with additional experiments
comment: Thank you for the useful suggestion. It is important to note that a great deal of the motivation for this work is the need in the biological sciences to have methods that balance predictive power with explainability. A highly predictive model has some utility but we want to understand why the system evolves from State A to State B so that we can search for meaningful ways to perturb the system in its evolution and affect the outcome. For example, if we have a tumor that is responsive to chemotherapy, we would want to block the tumor from evolving to a chemo-resistant state. A model with high predictability and high explainability would allow us to develop hypotheses about interventions. We would be happy to further clarify the motivation and to provide more insight into the biology of the systems we are modeling.
We would be happy to include further explanations and we do believe this would make for a simple revision -- after all the motivation, method, and experiments are almost the same, only additional background of the particular application is added, which is usually considered a minor change.
We are happy that you appreciate the additional results. Note that the results are not directly comparable as the number of layers is different. Our main goal is to show that our method is also applicable to other architectures (MLP + ELU activations as in the pathreg paper) rather than replication of previous results. Regarding the requested changes, we do believe that the suggested changes about a restructuring of text is usually considered a minor change as motivation, method, as well as experiments are still similar.
Thanks for your response. For the biological background, I think it is more important to describe what is the meaning of some biological names. For example, what is a regulatory factor? You also give the transcription factors as one example of regular factors, which is harder to understand. I believe the biological background can provide a more easy way for reviewers to understand the importance of your solutions for the problems. The paper needs further polishing for publication. I will keep my score.
I thank the authors for the time and effort put into their rebuttal, especially adding the ablation results comparing PHOENIX to another baseline model in Table 1 (R) for bone marrow data and simulated yeast time-series data. I also appreciate the authors addressing the mathematical notation errors in the text. 

However, I still do not understand the discrepancy between obtained results and those shown in PathReg (cf. Appendix C where they show extended results, where the sparsity levels and MSE obtained differ markedly from those reported here). Moreover, I strongly feel the overall structure of the paper needs substantial reorganising and editing to make it more readable and accessible, with an emphasis on properly motivating the DASH model both mathematically and biologically, introducing PHOENIX early on and clearly enunciating the difference between PINN+MP/BioPrune and DASH. Given this would require substantial modification of the work as it stands, I am maintaining my original rating.
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.


| Strategy | Prior corruption | Sparsity(\%) | Bal. Acc.(\%) | MSE ($10^{-3}$)  |
| :----- | :---: | :---: | :---: | :---: |
| None/Baseline | - | 11.5 | 54.8 | 3.6 |
| $L_0$ | -  | 34.7 | 61.3 | 6.1 | 
| C-NODE | -  | 10.7 | 60.5 | 1.9  |
| PathReg | -  | 59.7 | 64.2| 6.1 |
| DST |- | 94.3 | 72.3 | 4.2 |
| IMP | - | 86.1 | 63.2 | 4.1  |
| Iter. SynFlow | - | 79.1 | 60.0 | 2.3 |
| SparseFlow | - | 95.8 | 72.8 | 2.9 |
| PINN |  0\%  | 11.3 | 60.3 | 2.3  | 
| PINN |   20\%  | 12.4 | 60.8 | 3.1 |
| PINN |   40\%  | 11.2 | 60.6 | 2.7 |
| BioPrune |  0\%  | 83.5 | 88.0 | 3.6 |
| BioPrune|   20\%  | 80.9 | 81.5 | 7.6 |
| BioPrune |   40\% | 86.8 | 80.1 | 11.1 |
| DASH |  0\%  | 92.6 | 91.1 | 1.9  |
| DASH |   20\%  | 92.4 | 86.2 | 6.7 |
| DASH |   40\%  | 85.9 | 79.5 | 6.1 |
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.

| Strategy | Sparsity | Bal. Acc. (ChIP-seq) | Bal. Acc. (TF Perturb.) |
| :------- | :------: | :-----------------: | :----------------------: |
| None/Baseline   | 0.10\% | 49.87\% | 49.92\% |
| $L_0$   | 34.43\% | 48.43\% | 49.28\% |
| C-NODE   | 10.89\% | 50.04\% | 50.17\% |
| PathReg   | 12.09\% | 50.11\% | 49.92\% |
| PINN   | 0.17\% | 49.93\% | 50.01\% |
| DST  | 77.80\% | 49.92\% | 50.33\% |
| IMP| 83.22\% | 49.99\%  | 48.45\% |
|Iter. SynFlow  | 85.65\% | 49.57\% | 49.77\% |
| SparseFlow  | 95.22\% | 49.89\% | 51.58\% |
| BioPrune   | 94.69\% | 79.23\% | 64.50\% |
| DASH | **97.18\%** | **88.43\%** | **66.79\%**|
| PINN + MP  | 95.01\% | 55.39\%  | 52.95\% |

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).
I have read the authors’ rebuttal. I hope you’ll include these clarifications and additional information in your revision.

First, I apologize that I missed copy pasted my summary of your paper. Here was the original version I wrote in my own notes:
*This paper proposes a new method called DASH (Domain-Aware Sparsity Heuristic) for pruning neural ordinary differential equation (NODE) models to infer gene regulatory dynamics from time series gene expression data. The key innovation is incorporating prior biological knowledge as soft constraints during pruning to achieve more interpretable and biologically plausible sparse models. The authors evaluate DASH against existing pruning methods on both synthetic and real gene expression datasets, showing it can achieve high sparsity while maintaining predictive accuracy and recovering “known” regulatory interactions. Overall, the paper demonstrates value from biology-informed pruning for inferring interpretable gene regulatory networks from dynamic data.*

I don't understand why you can't use TF perturbation data to evaluate your method. There's a large literature on this. E.g. see Kamal, A. et al. GRaNIE and GRaNPA: inference and evaluation of enhancer‐mediated gene regulatory networks. Mol Syst Biol e11627 (2023).

I remain concerned that the simulations use impractically small noise levels on the prior.

Altogether, I’ll maintain my current scores.