DynGFN: Towards Bayesian Inference of Gene Regulatory Networks with GFlowNets
Authors: Lazar Atanackovic
Alexander Tong
BO WANG
Leo J Lee
Yoshua Bengio
Jason Hartford
Venue: NeurIPS 2023 poster (2023)
Keywords: Bayesian Structure Learning
Generative Flow Networks
Single-cell
Dynamical Systems
TL;DR: A GFlowNet approach for Bayesian dynamic structure learning of gene regulatory networks.

Abstract: One of the grand challenges of cell biology is inferring the gene regulatory network (GRN) which describes interactions between genes and their products that control gene expression and cellular function. We can treat this as a causal discovery problem but with two non-standard challenges: (1) regulatory networks are inherently cyclic so we should not model a GRN as a directed acyclic graph (DAG), and (2) observations have significant measurement noise so for typical sample sizes, there will always be a large equivalence class of graphs that are likely given the data, and we want methods that capture this uncertainty. Existing methods either focus on challenge (1), identifying cyclic structure from dynamics, or on challenge (2) learning complex Bayesian posteriors over directed acyclic graphs, but not both. In this paper we leverage the fact that it is possible to estimate the ``velocity'' of the expression of a gene with RNA velocity techniques to develop an approach that addresses both challenges. Because we have access to velocity information, we can treat the Bayesian structure learning problem as a problem of sparse identification of a dynamical system, capturing cyclic feedback loops through time. We leverage Generative Flow Networks (GFlowNets) to estimate the posterior distribution over the combinatorial space of possible sparse dependencies. Our results indicate that our method learns posteriors that better encapsulate the distributions of cyclic structures compared to counterpart state-of-the-art Bayesian structure learning approaches.

DynGFN: Towards Bayesian Inference of Gene Regulatory Networks with GFlowNetts
Lazar Atanackovic
Alexander Tong2.4+
Bo Wang
4.2. 12,5
Leo J. Lee
Yoshua Bengio 3,4,6
Jason Hartford 3.4.7
University of Toronto, Vector Institute Mila Quebec Al Institute, *Université de Montréal University Health Network, "CIAAR Fellow, VValence Labs
Abstract
One of the grand challenges of cell biology is inferring the gene regulatory network (GRN) which describes interactions between genes and their products that control gene expression and cellular function. We treat this causal discovery problem but with two non-standard challenges: (1) regulatory networks are inher- should not model GRN directed acyclic graph (DAG), and (2) observations have significant measurement noise, so for typical sample sizes there will always be large equivalence class of graphs that are likely given the data, and we want methods that capture this uncertainty. Existing methods either focus on challenge (1), identifying cyclic structure from dynamics, or on challenge (2) learning complex Bayesian posteriors over DAGs, but not both. In this paper we leverage the fact that it is possible estimate the "velocity" of gene expression with RNA velocity techniques to develop an approach that addresses both challenges. Because we have access velocity information, we can treat the Bayesian structure learning problem as problem of sparse identification of a dynamical system, capturing cyclic feedback loops through time. Since our objec- tive is to model uncertainty over discrete structures, we leverage Generative Flow Networks (GFlowNets) estimate the posterior distribution over the combinatorial space of possible sparse dependencies. Our results indicate that our method learns posteriors that better encapsulate the distributions cyclic structures compared counterpart state-of-the-a Bayesian structure learning approaches.
ently cyclic
Introduction
Inferring gene regulatory networks (GRNs) a long standing problem in cell biology [25, 44]. If we knew the GRN, it would dramatically simplify the design of biological experiments and the development of drugs because it would serve as map of which genes to perturb manipulate protein and gene expression. GRNs concisely represent the complex system of directed interactions between genes and their regulatory products that govern cellular function through control RNA causal discovery
(gene) expression and protein concentration. We can treat GRN inference problem by treating the regulatory structure between genes (variables) as causal dependencies (edges) that we infer / rule out by using gene expression data. Structure learning methods aim to automate this task by inferring set of directed acyelic graphs (DAGs) that are consistent with the conditional independencies that
measure among the variables [14, 41, 42]. While there may be multiple
"Equal Contribution Correspondence to: a att.aacckoo.@@mmll...cooo
37th Conference on Neural Information Processing Systems (NcurIPS 2023).

---

p(G.0. = P.1...000111
GFN Graph Sampler
Structural Model
Parameter HyperNetwork
dxp
QytG)
(G) oc PG))
p [0|
p(D (DDGG,0 -
hg(G) -
Go(x. - Tx
Figure I; Architecture for Bayesian structure learning of dynamical systems. DynGFN consists of three main components: GFlowNet modeling posterior distribution over graphs Qu (GID), HyperNetwork modeling posterior over parameters given graph (0 D), and the structural equation model scoring G and 0 according how well they fit the data. Although the figure shows the case where (G) is modelled with GFlowNet, this can be any arbitrary graph sampler that can sample discrete structures G (G)D)
DAGs in this set-the "Markov equivalence class"-when we are able to perturb the variables with enough experimental interventions, it is possible to uniquely identify causal graph [18].
However, structure learning for inferring GRNs comes with two non-standard challenges: (1) gene regulation contains inherent cyclic feedback mechanisms, hence we should not model a GRN as DAG, and (2) observations are limited and have significant measurement noise, hence there exists a large equivalence class of graphs that are likely given datasets with typical sample sizes. Existing methods either focus on (1) 100-idene identifying graphs with cyclic structure by leveraging dynamics [15, 13] or assuming the system is in equilibrium [36], or (2) learning complex Bayesian posteriors over explanatory DAGs [12], but not both. In this work, we address both challenges concurrently in a fully differentiable end-to-end pipeline (see Figure 1).
To accomplish this, we treat structure learning system. From a dynamical systems perspective, one can model both causal structure between variables as well as their tome-dependen system response with the drift function [37, 43]. We leverage the fact that we can estimate the rate of change of gene's ' expression (velocity) with RNA velocity methods [9]. This data takes the form of dynamic tuple pairs (x, dx). which we can use to pose the dynamical system learning problem
problem of sparse identification dynamical
regression task (see Figure 1). This significantly simplifies the learning objective as we can model system dynamics while also learning structure without the need for numerically intensive differrntial equation solvers. We view this step towards Bayesian structure learning from continuous dynamics- we term this Bayesian dynamic structure learning.
Our approach estimates the posterior over the sparse dependencies and parameters of the dynamical system. This is important scientific applications because it is usually prohibitively expensive to acquire a enough data to uniquely identify the true graph underlying data generating process. Capturing the complex distribution over candidate structure is critical for downstream scientific applications and is an essential step in active causal discovery [39, 52, 19]. This is especially important settings where experiments are expensive, e.g. conducting genetic perturbations for inference of GRNs. Bayesian structure learning is a class of methods that try to model this distribution over structure from observed data. These methods model posteriors over admissible structures P(G that explain the observations [32, 11, 4, 12, 31], but focus on modelling distributions over DAGs.
Our approach leverages Generative Flow Networks (GFlowNets) to model complex distributions over cyelic structures. GFlowNets [7, 8] parameterize the distribution over any discrete object (e.g. graphs) through a sequential policy, and as a result avoid needing to make restrictive parametric assumptions on the distribution. This makes them a useful tool in structure learning, particularly in cases where D) discrete and complex [12). In this work, we use GFlowNets to learn posteriors over the sparse structure in a dynamical system, and separately learn the posteriors over the parameters of the drift function via HyperNetwork [16] that conditions on inferred structures. Our main contributions summarized as follows:

---

We develop a novel framework for Bayesian structure learning under the lens of dynamical system identification for modelling complex posteriors over eyclic graphs. We consider flexible parameterizations for the structural model such that we can capture both linear and non-linear dynamic relationships.
We design novel Get architecture, Dynamic Gel (DynGFN), tailored for modelling posteriors over cyclic structures. We propose per-node factorization within DynGFN that enables efficient search over the discrete space of cyclic graphs.
We empirically evaluated DynGFN on synthetic dynamic data designed to induce highly multi- modal posteriors over graphs.
We showease the use of DynGFN on a real biological system using single-cell RNA- ANA-velocity data for learning posteriors of GRNs.
2 Related Work
There are many works on the problem of identifying causal structure G from either observational [e.g.
49, 54, 35] or interventional [e.g. 26, 29, 38] data, but the majority of existing methods return only the most likely DAG under the observed data. By returning only the most likely graph, these methods are often overconfident in their predictions. Bayesian approaches attempt to explicitly model posterior distribution over DAGs given the data and model specification.
Bayesian Structure Learning: Recently, there has been significant interest in fully differentiable
Bayesian methods for structure learning in the static case. DiBS [32], BCD-Nets [11], VCN [4], and DAG-GFlowNet [12] all attempt to learn distribution over structural models from fully observed system. The key difference is how these methods parameterize the graph. DiBS is a particle variational inference method that uses two matrices U and V where G = sigmoid(UU) ((UTV) where the sigmoid is applied elementwise which is similar to graph autoencoders. BCD-Nets and DP-DAG use the Humbell-Snkhom distribution to parameterize permutation and direct parameterization of lower triangular matrix. VCN uses an autoregressive LSTT to generate the graph as this gets rid of the standard uni-modal constraint of Gaussian distributed parameters. DAG-GFN has shown success for modelling P(G 'D) [12]. However, it remains restrictive to assume the underlying structure of the observed system is DAG as natural dynamical systems typically contain regulating feedback mechanisms. This can be particularly challenging for GelowNets since including cycles in the underlying structure exponentially increases the discrete search space. We show that under certain assumptions we can in part alleviate this shortcoming for learning Bayesian posteriors over cyclic structures for dynamical systems. In small graphs, these methods can model the uncertainty over possible models (including over Markov equivalence classes).
Dynamic and Cyclic Structure Learning: There has been comparatively little work towards
Bayesian structure learning from dynamics. Recent works in this direction based on NeuralODEs [10]
propose a single explanatory structure [50, 6, 1, 2]. CD-NOD leverages heterogeneous non-stationary
data for causal discovery when the underlying generative process changes over time [53, 21). A similar approach uses non-stationary time-serics data for causal discovery and forecasting [20]. DYNOTEARS is score-based approach that uses time-series to learn structure [40]. However, these methods do not attempt to explicitly model a distribution over the explanatory structure. Other methods aim to learn cyclic dependencies in the underlying graph [24, 36, 28, 3]. For instance, [24]
propose an iterative method that leverages interventional data to learn directed cyclic graphs. is suggested that CD-NOD is also extendable to learn cyclic structure [21]. But these methods do not model a posterior over structure. In general, there remains a gap for the problem of Bayesian structure learning over cyclic graphs.
We include further discussion on related work for GRN inference from single-cell transcriptomic data and cell dynamics in Appendix C.
3

---

3 Preliminaries
3.1 Bayesian Dynamic Structure Learning
We consider a finite dataset, D, of dynamic pairs de) E I2 x R where x respresents the state of the system sampled from an underlying time-invariant stochastic dynamical system governed by latent drift de f(z, e) where is noise term that parameterizes the SDE; and are mutually independent. The latent drift has some fixed sparsity pattern i.e. 0f Bxz + for small set of variables, which can be parameterized by a graph G such that 9-j 1 ar, * 0). where Gi G. i l,.... .....dd l.....,d. d. The variables 2, for which 0s 0z, o can be interpreted as the causal parents of z,, denoted Pa(æ;). This lets us define equivalent dynamic structural model [37, 43] of the form,
Problem Setup:
drç(1) = ()) dt
for i =1=,,, 1.....,,,, For the graph G to be identifiable, we assume that all relevant variables are observed, such that causal sufficiency is satisfied.
Our goal to model our posterior over explanatory graphs Q(G[D) given the data. We aim to jointly learn distribution over parameters 0 that parameterize the latent drift f(x); these parameters will typically depend on the sparsity pattern such that p(0]G) # p(0). We can factorize this generative model follows,
p 0, D) =
This factorization forms the basis of our inference procedure. We learn parameterized function fo(x) : TRd + Rd that approximates the structural model defined in (1). To model this joint the combinatorial space of
distribution, need a way of representing P(G), a distribution possible sparsity patterns, and P()[G), the posterior over the parameters of fa. We use GFlowNets [7] represent P(G), and HyperNetwork to parameterize P(O
3.2 Generative Flow Networks
GFlowNets are an approach for learning generative models over spaces of discrete objects [7, 8]. GelowNets learn stochastic policy (P(T) to sequentially sample object from a discrete space , Sn)) represents a full Markovian trajectory over plausible discrete states, where s, is the terminating state (i.e. end of trajectory) [34]. The GFlowNet is trained such that at convergence, sequential samples from the stochastic policy over trajectory, x ~ PF(T), x Sn, R(x) equal in distribution to samples from the normalized reward distribution P(x) Ex'ex R(XT)
X. Here 7 = ($g. $1....,
The GFlowNet policies are typically trained by optimizing either the Trajectory Balance (TB) loss [34], Subtrajectory Balance (Sub-TB) loss [33], or the Detailed Balance (DB) loss [12]. In this work, we exploit the DB loss to learn a stochastic policy for directed graph structure.
Detailed Balance Loss: The DB loss [12] leverages the fact that the reward function can be evaluated for any partially constructed graph (i.e. any prefix of 7), and hence get intermediate reward signals for training the GFlowNet policy. The DB loss is defined as:
log R(s:)) (n v Pr (8m18111:1
LDn (Dn(St,
R( Pr
where Pr (s ()) and PB(s;-1[s;; v)) represent the forward transition probability and backward transition probability, and a trainable normalizing constant, respectively. Under this formulation, during GFlowNet training the reward is evaluated every state. For this reason, the DB formulation is in general advantageous for the structure learning problem where any sampled graph can be viewed as a complete state, hence more robustly inform gradients when training the stochastic policy than counterpart losses. Previous work has shown GFlowNets are useful in settings with multi-modal posteriors. This is of particular interest us where many admissible structures can explain the observed data equally well. We model (G) using PF (p 8 v) and learn the parameters ş.

---

Algorithm Batch update training of DynGFN
1: Input: Data batch (ra, drs), initial NN weights v. b, I sparsity prior Ao and learning rate
2: 80 + Oxxdxd
b Training is paralleled over B graph trajectories Sample initial actions vector
3: (s( so: so:1/).
4: while a not & do
5:
Compute PF(s: (5(5s(5--1i1)) Pa(Si-1(s5;))
6:
0 + (s(ss)
7:
dæp + sz(
(s:):) + e e-lda --Iddy-dull 1t, -
8:
9:
# v -V CV 8;-1)
LDDB(,,, 8;-1) computed in Equation b Take action step to go to next state
10:
an Pr(s;(a(1;;;)) 1; + 8i+1
@ ++ + log R
return Updated GFN weights vp and updated HyperNetwork weights Ø
DynGFN for Bayesian Dynamic Structure Learning
We present a general framework for Bayesian dynamic structure learning and propose a GFlowNet architecture, DyrGFN, tailored for modelling posterior over discrete cyclic graphical structures. We summarize our framework Figure and Algorithm 1. DynGFN consists of key modules:
1. A graph sampler that samples graphical structures that encode the structural dependencies among the observed variables. This is parameterized with GFlowNet that iteratively adds edges to graph.
2. A model that approximates the structural equations defined in (1) model the functional relationships between the observed variables, indexed by parameters 6. This is a class of functions that respect the conditional independencies implied by the graph sampled in step
1. We enforce this by masking inputs according to the graph.
3. Because the functional relationships between variables may be different depending on which graph is sampled, we HyperNetwork architecture that outputs the parameters of the structural equations as a function of the graph. We also show that under linear assumptions of the structural modules, we solve for optimal 0 analytically (i.e. do not need the HyperNetwork).
For training, we assume 7 sparsity of graphs G to constrain the large discrete search space over possible structures. We use reward R for a graph G and I penalty of the form: R(G) (illde-di1 We motivate this set-up so we can estimate dr close to da in an end-to-end learning pipeline. Since estimates for da are dependent on G and 0, this reward informs gradients to learn policy that can approximate Q(G) given dynamic data.
The main advantage of DynGFN comes when modelling complex posteriors with many modes. Prior work has shown GFlowNets are able to efficiently model distributions where
share information between different modes [34]. The challenge we tackle is how to do this with a changing objective function, the GelowNet objective is function of the current parameter HyperNetwork and the structural equations. We use multilayer perceptrons (MLPs) to parameterize the stochastic GFlowNet policy, HyperNetwork architecture, and the dynamic structural model!.
4.1 Graph Sampler
DynGFN models a posterior distribution over graphs Q(G[D) given a finite set of observations. To learn (G)))) DynGFN needs explore over large discrete state space. Since we aim learn a bipartite graph between r and dr, DynGFN needs to search over 2d possible structures, where d denotes the dimensionality of the system and 2 the number of possible edges in G. For even moderate d, this discrete space is very large (e.g. for 20 we have 2400 possible graphs).
'hen we assume linear dynamic structural relationships, can solve for the parameters analytically, thus do not need MLPs for the HyperNetwork and dynamic structural model. This is further discussed in section 4.2
5

---

However, under the assumption of causal sufficiency, we can significantly reduce this search space, by taking advantage of the fact that QGG(D) factorizes as follows,
Q(G)D) - II Q((G!,,jID) iE[1.....]]
By using this model, we reduce the search space from 2d + d2d For d = 20 this reduces the search space from 2400 to ze 2213, While still intractable to search over, it is still vast improvement over the unfactorized case. We call this model a per-node posterior, and we use per-node GFlowNet going forward. We discuss details regarding encouraging forward policy exploration during training in Appendix B.
4.2 HyperNetwork and Structural Model
We aim to jointly learn the structural encoding G and parameters 0 that together model the structural relationships dz = fo (z. G) of the dynamical system variables. To accomplish this, we propose learning an individual set of parameters I for each graph G. independent of the input data z. This approach encapsulates G() in (2). We use HyperNetwork architecture that takes G as input and outputs the structural equation model parameters 0, i.e. 0 = h,) hence P(O 8(0 allowing us to learn a separate for each G. This HyperNetwork model does not capture uncertainty in the parameters, however the formulation may be extended to the Bayesian setting by placing prior on the HyperNetwork parameters Φ. Although h& allows for expressive parameterizations for 0, may not be easy to learn". HyperNetworks have shown success in learning parameters for more complex models (e.g. LSTMs and CNNs) [16], hence motivates their fit for our application.
Linear Assumption on Dynamic Structural Model: In some cases may suffice to assume a linear differential form dr Ar to approximate dynamics. In this setting, given sampled graph G ~ Q(G)) and n i.i.d. observations of (Ir dx) we can solve for 0 A analytically. To induce dependence the graph structure, we use the sampled G' as mask and construct ã; = G." Then we can solve for 0 on a per-node basis as
0: = (G 3. + XT) 57 dxi,
where i 1........ d, A > 0 the precision of an independent Gaussian prior over the parameters, and is the identity matrix. We use À = 0.01 throughout this work.
5 A Useful Model of Indeterminacy
In order to evaluate the ability of DynGFN to model complex posteriors over graphs, we need structure learning problem with large equivalence class of admissible graphs. We present simple way to augment set of identifiable dynamics under some model to create a combinatorial number of equally likely dynamics under the same model. More specifically, this creates a ground truth posterior Q ( D) x TT())) where T(.) 9 + G an analytically computable transformation over graphs and G" is the identified graph under the original dynamics. We use this system to test how well we can learn a posterior over structures that matches what we see in single-cell data.
Specifically, given a dataset of (Ir, dx) E Rd x Rd pairs, we create a new dataset with d ++ variables where the 'new" variable v' is perfectly correlated with an existing variable v. In causal terms, this new variable inherits the same parents as v, that is Pa(v')= := Pa(v) and the same structural equations as v. that du' dv. This is depicted in Figure 2. This creates number of new possible explanatory graphs, which we generalize with the following proposition.
Proposition 1. Given any d dimensional ODE system with G identifiable under f E F. the D = d + dimensional system dæ di As, denote the vector of multiplicities m E Nd with ma as the number of repetitions of each variable. Then this construction creates an admissible family of graphs g' where |g'| = (2 - 1)8.1. Furthermore, under an T penalty on G. this reduces to Im,),
See Appendix A for full proof. The intuition behind this proposition can be seen from the case of adding a single copied variable. This corresponds to A = [o. where 6, is a vector with a on
We discuss training dynamics when using ho in Appendix B.7.

---

Dynamic Graph
Static Graph
x
dx
dr
P2
G
Figure 2: Visualization of modeling cyclic dependencies over time (left). To create a family of admissible graphs (right), we add new variable v/g which has the same values as 03 and creates three possible explanations for the data (green arrows). If we apply sparsity penalty, then we can eliminate the last possibility (which has two additional edges) for only two possible graphs. G" denotes the ground truth graph while G' denotes the admissible family of graphs induced by 2'g.
node e and zeros elsewhere, and Ia is the d-dimensional identity matrix. Let v have e children, such that Pa(c) in the identifiable system, then any of those c child nodes could depend either on v or on the new node o' or both. This creates 3 possible explanatory graphs. If we restrict ourselves to the set of graphs with minimal 77 norm, then we eliminate the possibility of child node depending on both v and e', this gives 2 possible graphs, choosing either or v' as parent.
6 Experimental Results
In this section we evaluate the performance of DynGFN against counterpart Bayesian structure learning methods (see Appendix B.2 for details). Since our primary objective is to learn Bayesian posteriors over discrete structure G, we compare to Bayesian methods that can also accomplish this task, i.e. versions of BCD-Nets [11] and DiBS [32]. We show in certain cases, DynGFN is able to better capture the true posterior when there are a large number of modes. We evaluate methods according to four metries: Bayes-SHD, area under the receiver operator characteristic curve (AUC), Kullback-Leibler (KL) divergence between learned posteriors (G) and the distribution over true graphs P(G"), and the negative log-likelihood (NLL.) P(D[G, 0) (in our setting this reduces to the mean squared етrоr between dr and dr. given 0 and sampled G)s). Since the analytic linear solver requires data run-ting run-time to compute optimal parameters for the structural model, we include the NLL metric only for models using the HyperNetwork solver. Bayes-SHD measures the average distance to the closest structure in the admissible set of graphs according to the structural hamming distance, which in this case is simply the hamming distance of the adjacency matrix representation to the closest admissible graph. We assume P(G*) is uniform over G and include further details about evaluating the quality of learned posteriors in Appendix B.8.
6.1 A Toy Example with Synthetic Data
To validate the ability of DynGEN to learn cyclic dependencies we consider identifiable acyclic and cyclic 3 variable toy systems and provide a comparison with a DAG structure learning method (NOTEARS [54]). We show results for this toy example in Figure 3. NOTEARS does not model cyclic dependencies and therefore struggles yield accurate predictions of A in the cyclic setting. We can also verify this result by considering conditional independence test over cyclic dependencies. is easy to see that the conditional independence test fails the cyclic setting: in the acyclic case, we can identify the v-structure by observing that æ1 1. 23 and L .ra Z2, which implies that 22 a collider (i.e. æn and v3 are marginally independent and conditionally dependent); while in the cyclic example, we introduce time dependencies such that there are cycles in the summary graph that render these variables marginally dependent. We show that DynGFN able to identify the true cyclic dependencies in this toy example. We note that in this toy example DynGFN exhibits degree
7

---

Grache
DynoFN
Acydic
AUC-ROC
-3007
NOTEARS
DynSFN
Aoyalic
D.57
0.21
Cyclic
6M
Cyclic
0.58
1.0
6:11
Figure 3: Toy example for learning DAG structure (top) and cyclic graph structure (bottom) from observational data. In this example compare NOTEARS and DynGFN for learning acyclic and cyclic dependencies from observational data. Data is generated from linear system dr = Ar with a corresponding acyclic and cyclic A as defined in the figure. We use 500 observational samples for each system. NOTEARS implemented using the CausalNex [5] library.
Bayesian dynamic structure learning of linear and non-linear systems with d = 20 variables. The graphs representing the structural dynamic relationships of the linear and non-linear systems have 50 edges out of possible 400. The ground truth discrete distribution P(G*) contains 1024 admissible graphs for each respective system. The l and pre-fix denote usage of the analytic linear solver and HyperNetwork solver for structural model parameters, respectively. Results are reported on held out test data over 5 model seeds.
Table I:
Linear System AUC
Model
Bayes-SHD
KL 1
NLL
B-DynBCD
32.0 = 0.27
0.71 +: 0.0
1707.45 + 9.66
4 4-DynDiBS
29.2 = 0.78
0.71
0.0
6622.43 + 171.67
B-DynGFN
22.8 + 1.4
0.75 + 0.01
1091.60 + 35.72
R-DynBBD
5.5 1.1
1 0.04
701.19 46.99
(9.83 1 0.59)
h-DynDiBS
28.5 4.2
0.51 + 0.07
7934.90
381.80
(8.17 1 1.30
h-DynGFN
6.7 0.0
0.94
350.92 30.15
(8.35 + 0.02)
Non-linear System
Model
Bayes-SHD
AUC
KL |
NLL 4
{-DynBCD
77.5 + 8.3
0.42 + 0.03
3814.86 + 354.56
D-ynDiBS
75.7 = 7.7
0,59 + 0.01
5893.65 + 59.66
DDDnGFF
45.7 0.6
0.55 0.55 0.0
226.25 + 6.58
h-DynBCD
192.9 -. 0.7
0.50 1.501:00 0.0
9108.69 1 51.34
(3.83 1 0.32) 4
h-DynDiBS
48.1 1 9.0
0.53 + 0.10
8716.64 1: 265.29
(4.06 1 0.10
h-DynGFN
32.6 + 0.9
0.67 0.01
193.28 + 8.53
(1.47 + 0.11)E
of convergence sensitivity on a per-run basis. the following sections we provide comprehensive results over 5 random seeds and consider larger systems.
6.2 Experiments with Synthetic Data
We generated synthetic data from two systems using our indeterminacy model presented in section 5: (1) linear dynamical system dr Ax, and (2) non-linear dynamical system dr sigmoid(A) (Arr). We consider E-DynGFN and h- h-DynGFN, i.e. DynGFN with the linear analytic parameter solver as shown in (5), and DynGFN with the HyperNetwork parameter solver h3- Likewise, we compare DynGFN and h- h-DynGFN to counterpart Bayesian baselines which we call Ê-DynBCD, é-DynDiBS, h- DynBCD, and h-DynDDBS. To constrain the discrete search procedure, we assume a sparse prior on structure (i.e. the graphs G), using the L prior. Due to challenging iterative optimization dynamics present when using (G) for DynGFN, to train initialize the forward policy Pr(s: 88-1: v) using the & Iearned in E-DynGFN to provide more admissible starting point for learning hi (we

---

discuss further details in Appendix B.7). We do not need to do this for h-DynBCD and h-DynDDBS as we are able to train both models end-to-end without iterative optimization. In Table we show results of our synthetic experiments for learning posteriors over multi-modal distributions of cyclic graphs. We observe the DynGFN is most competitive on both synthetic systems for modelling the true posterior over structure. Details about DynGFN, baselines, and accompanying hyper-paramenn parameters can be found in Appendix B.
6.3 Ablations Over Sparsity and Linearity of Dynamic Systems
We conduct two ablations: (I) ablation over sparsity of the dynamic system structure, and (2) ablation over At, the time difference between data points of dynamic simulation. For a sparsity level of 0.9, the ground truth graphs have 50 edges out of d2 possible edges. In these experiments, P(G*) for d = 20 and sparsity 0.9 has 1024 modes. We conduct the ablations over 5 random seeds for each set of experiments.
Sparsity: DynGFN uses the 7 prior on G throughout training. Under this setting, system sparsity carries significant weight on the ability to learn posteriors over the structured dynamics of a system. We show
Table Ablation for 4-DynGFN on d 20 systems with varying levels of sparsity and fixed At = 0.03.
Sparsity
Bayes-SHD
AUC
KL |
0.95
16.4 + 1.71
0.79 0.0
889.57 = 31.24
this trend in Table 2. We note that comput-
0.90
22.8 +
1.41 1.41
0.75
1091.60 + 35.72
ing the KL-divergence for DynGFN, specif-
0.85
32.8 1.88
0.71
ically computing the probability of gener-
0.80
39.2 4: 0.69
0.71
ating a true G, becomes computationally intractable as G is less sparse'. For sys- tems of 0.9 and 0.95 sparsity, we observe a decreasing trend in KL and Bayes-SHD and an increasing trend AUC. This result is expected as DynGFN can better traverse sparse graphs as the combinatorial space over possible trajectories is smaller relative
0.75
60.2 = 1.17
0.66
Table 3: Ablation for (-DynGFN on d = 20 non-linear systems with varying At and fixed sparsity 0.9.
At
Bayes-SHD +
AUC
KL 4
0.001
38.7 0.80
0.61 0.0
202.41 + 9.95
0.005
39.0 0.81
0.60
206.83 + 11.55
to denser systems.
0.01
1.13
0.59 0.0
202.71 1 7.74
0.05
1.662
1.00
226.25 1 6.58
Linearity: Training DynGFN via the lin- ear solver for the structural model parame- ters is an easier objective due to simplified training dynamics. Because of this, we explore the performance of E-DynGFN assuming fo for modelling equation (1) to be linear in the non-linear system. We do this by conducting an ablation over At and find that the performance of é-DynGFN on the non-linear system improves At - 0. We show a portion of this trend in Table 3.
0.1
0.50
2641.86 1 2.17
6.4 Experiments on Single-Cell RNA-velocio Data
To show how DynGFN can be applied to single cell data we use a cell cycle dataset of human Fibroblasts [46]. As a motivating example we show the correlation structure of single-cell RNA-seq data from human Fibroblast cells [46] Figure 4. We show both the raw correlation and the correlation over cell cycle time, which is significantly higher. With such pure cell population whose primary axis of variation state in the cell cycle by aggregating over cell cycle time we expect observation noise to be averaged out, leading to "truer" view of the correlation between latent variables. Further details for this experimental set-up are provided in Appendix C.I.3. Since there are many genes which affected by the cell cycle phase, there many correlated variables that are downstream of the true cell cycle regulators. This provides a natural way of using cell cycle data to evaluate a model's ability to capture the Bayesian posterior. In Table 4 we show results for learning posteriors over an undetermined GRN using RNA velocity data. We find that E-DynGFN and h-DynGFN yield low KL
$ FFor example, since DynGFN constructs one object G sequentially over state space distribution, we must compute probabilities of all combinatorial state trajectories for constructing G = (S.,..., s.).. The space of combinatorial state trajectories is n! in nature, hence this computation is only possible for small graphs and/or sparse graphs.

---

(a)
(b)
(c)
(d)
Figure 4: (a) Correlation structure in the raw single cell data 5000 cells and 2000 genes selected
by selo [8] pre- pre-processing. (b) Correlation structure among genes over (inferred) cell cycle times. This stronger correlation structure is more reflective of the correlation in the underlying system. (c) Histogram of pairwise Pearson correlations between all genes passing pre-processing, comparing the absolute values of the elements in (a) and (b). (d) Shows the ground truth GRN extracted subset of the KEGG cell cycle pathway. Cdc25A is known to inhibit Cdk1 which is known to inhibit Cdc25C, while the Mem complex is highly correlated with Cdc25A, they do directly interact with Cdk1 [24].
Table 4: Bayesian dynamic structure learning 5-D cellular system using scRNA velocity data. The dynamics of this system are unknowns, however identify 81 admissible graphs between variables (genes) that describe the data. We train models over 5 seeds. The graphs of this system contain of 7
true edges.
Cellular System RNA Velocity
Model
Bayes-SHD
AUC 1
KL
NLL +
E-DynBCD
2.6 : 0.1
0.56 + 0.01
321.95
3.34
E-DynDiBS
6.5 0.4
0.47 0.01
550.17
16.63
i-DyaGFN
3.3 0.4
0.59 1 0.03
44.98 1 18.60
h- h-DynBCD
10.1 0.8
0.53 1 0.03
587.11 + 24.00
0.094 + 0.003
DynDiBS
9.6 4.2
0.51 : 0.13
560.85 = 83.83
0.084 = 0.0
h-DyGGFN
+ 1.2
0.58 + 0.05
39.82 + 28.05
0.109 + 0.001
and moderate Mayes-SHD. While (-DynBCD performs well in terms of identify a small distribution of true G's, it falls short in modelling the true posterior (this can be seen from low Bayes-SHD. high KL).
7 Conclusion
We presented DynGFN, a method for Bayesian dynamic structure learning. In low dimensions we found that DynGFN is able better model the distribution over possible explanatory structures than counterpart Bayesian structure learning baseline methods. As proof of concept, we presented an example of learning the distribution over likely explanatory graphs for linear and non-linear synthetic systems where complex uncertainty over explanatory structure is prevalent. We demonstrate the use of DynGFN for learning gene regulatory structure from single-cell transcriptomic data where there are many possible graphs, showing DynGFN can better model the uncertainty over possible explanations of this data rather than capturing a single explanation.
Limitations and Future Work:
We have demonstrated degree of efficacy when using DynGFN for Bayesian structure learning with dynamic observational data. A key limitation of DynGFN is scaling to larger systems. To effectively model P(G, 0, D), DynGFN needs to search over environment state space of possible graphs. This state space grows exponentially with the number of
possible edges, i.e. 2d or d2d for nod-nnddGGNN where d is the number of variables in the system. Therefore, DynGFN is currently limited smaller systems. Nevertheless, there are many applications where Bayesian structure learning, even over 5-20 dimensional examples that we explore here, could be extraordinarily useful. We include further discussion of scaling DynGFN Appendix C.3 with some ideas on how to approach this challenge. We found that training DynGFN requires some selection of hyper-parane perre-pprseeters and in particular parameters that shape the reward function. Selecting paramearssssess for the baseline methods prove more difficult for this task.
10

---

8 Acknowledgments and Disclosure of Funding
This research was enabled in part by the computational researches provided by Mila (mila. quebec) and Compute Canada (ccdb, computecandd. ca). In addition, resources used in preparing this research were in part provided by the Province of Ontario and companies sponsoring the Vector Institute (vectorinsti tute ai/partners/), All authors are funded by their primary academic institution. We also acknowledged funding from the Natural Sciences and Engineering Research Council of Canada, Recursion Pharmaceuticals, CIFAR, Samsung, and IBM. We are grateful to Cristian Dragos Manta for catching some typos the manuscript.
References
[1] Aliee, H., Theis, F. J., and Kilbertus, N. Beyond predictions in neural odes: Identification and interventions. arXiv preprint 1206.1243, 2021.
[2] Alice, H., Richter, T., Solonin, M., Ibarra, I., Theis, F., and Kilbertus, N. Sparsity in continuous- depth neural networks. Advances in Neural Information Processing Systems (NeurlPS), 2022.
[3] Améndola, C., Dettling, P., Drton, M., Onori, F., and Wu, J. Structure learning for cyclic linear causal models. Uncertainty in Artificial Intelligence (UAI), 2020.
[4] Annadani, Y., Rothfuss, J., Lacoste, A., Scherrer, N., Goyal, A., Bengio, Y., and Bauer, Variational causal networks: Approximate bayesian inference over causal structures. arXiv preprint, 2021.
[5] Beaumont, P., Horsburgh, B., Pilgerstorfer, P., Droth, A., Oentaryo, R., Ler, S., Nguyen, H., Ferreira, G. A., Patel, Z., and Leong, W. CausalNex, 2021. URL https://githu/ com/
quantumblacklabs/caaaa causalexx
[6] Bellot, A. and Branson, K. Neural Graphical Modelling in Continuous Time: Consistency Guarantees and Algorithms. International Conference on Learning Representations (ICLR), 2022.
[7] Bengio, E., Jain, M., Korablyov, M., Precup. D.,. and Bengio, Y. Flow Network based Generative Models for Non-Iterative Diverse Candidate Generation. Advances in Neural Information Processing Systems (NeurIPS), 2021.
[8] Bengio, Y., Deleu, T., Hu, J., Lahlou, S., Tiwari, M., and Bengio, E. GFlowNet Foundations. arXiv preprint 2111.09266, 2022.
[9] Bergen, V., Lange, M., Peidli, S., Wolf, F. A., and Theis, J. Generalizing RNA velocity to transient cell states through dynamical modeling. BioRxiv preprint 820936, 2019.
[10] Chen, R. T. Q., Rubanova, Y., Bettencourt, and Duvenaud, D. Neural Ordinary Differential Equations. Advances in Neural Information Processing Systems (NeurIPS), 2018.
[11] Cundy, C., Grover, A., and Ermon, S. BCD Nets: Scalable Variational Approaches for Bayesian Causal Discovery. Advances in Neural Information Processing Systems (NeurIPS). 2021.
[12] Deleu, T., Góis, A., Emezue, C., Rankawat, M., Lacoste-Julien, S., Bauer, S., and Bengio, Y. Bayesian Structure Learning with Generative Flow Networks. Uncertainty in Artificial Intelligence (UAT), 2022.
[13] Friedman, N., Murphy, K., and Russell, S. Learning the structure of dynamic probabilistic networks. Uncertainty in Artificial Intelligence (UAI), 1998.
[14] Glymour, C., Zhang, K., and Spirtes, P. Review of Causal Discovery Methods Based on Graphical Models. Frontiers in Genetics, 2019.
[15] Granger, C. W. Investigating causal relations by econometric models and cross-spectral methods. Econometrica: journal of the Econometric Society, 1969.
[16] Ha, D., Dai, M., and Le, Q. Hypernetworks. International Conference on Learning Representations (ICLR), 2017.
11

---

[17] Hashimoto, T. B.. Gifford, D. K., and Jaakkola, T. S. Learning Population-Level Diffusions with Generative Recurrent Networks. International Conference on Machine Learning (ICML), 2016.
[18] Hauser, A. and Bühlmann, P. Characterization and greedy learning of interventional markov equivalence classes of directed acyclic graphs. Journal of Machine Learning Research (JMLR), 13(1), 2012.
[19] He, Y.-B. and Geng, Z. Active learning of causal networks with intervention experiments and optimal designs, Journal of Machine Learning Research (JMLR), 9, 2008.
[20] Huang, B., Zhang, K., Gong, M., and Glymour, C. Causal discovery and forecasting in nonstationary environments with state-space models. International Conference on Machine Learning (ICML), 2019.
[21] Huang, B., Zhang, K., Zhang. J., Ramsey, J., Manchez-Romero R., Glymour, C., and Schölkopf, B. Causal discovery from hecerogeccoous/noostaatoorr data. The Journal of Machine Learning Research (JMLR), 21, 2020.
[22] Huguet, G., Magruder, D. S., Tong, A., Fasina, O., Kuchroo, M., Wolf, G., and Krishnaswamy, S. Manifold interpolating optimal-transpae flows for trajectory inference. Advances Neural Information Processing Systems (NeurIPS), 2022.
[23] Huguet, G., Tong, A., Zapatero, M. R., Wolf, G., and Krishnaswamy, S. Geodesic Sinkhorm: Optimal transport for high-dimensional datasets. arXiv preprint 2211.00805, 2022.
[24] Itani, S., Ohannessian, M., Sachs, K., Nolan, G. P., and Dahleh, M. A. Structure learning in causal cyclic networks. Proceedings of Workshop on Causality: Objectives and Assessment at NIPS 2008, 2010.
[25] Karlebach, G. and Shamir, R. Modelling and analysis of gene regulatory networks. Nature reviews Molecular cell biology, 9, 2008.
[26] Ke, N. R., Bilaniuk, O., Goyal, A., Bauer, S., Larochelle, H., Pal, C., and Bengio, Y. Learning neural causal models from unknown interventions. arXiv preprint 1910.01075, 2019.
[27] La Manno, G., Soldatov, R., Zeisel, A., Braun, F., Hochgerner, H., Petukhov, V., Lidschreiber, K., Kastriti, M. E., Lönnerberg, P., Furlan, A., Fan, J., Borm, L.. E., Liu, Z., van Bruggen, D., Guo, J., He, X., Barker, R., Sundström, E., Eestelo-Brance G., Cramer, P., Adameyko, I., Linnarsson, S., and Kharchenko, P. V. RNA velocity of single cells. Nature, 560, 2018.
[28] Lacerda, G., Spirtes, P. L., Ramsey, J., and Hoyer, P. O. Discovering eyclic causal models by independent components analysis. arXiv preprint 1206.3273, 2012.
[29] Lachapelle, S., Brouillard, P., Deleu, T., and Lacoste-Julien, S. Gradient-Basee Neural DAG Learning. International Conference on Learning Representations (ICLR), 2020.
[30] Liu, Q. and Wang, D. Stein variational gradient descent: general purpose bayesian inference algorithm. Advances in Neural Information Processing Systems (NeurIPS), 2016.
[31] Lopez, R., Hütter, J.-C., Pritchard, J. K., and Regev, A. Large-scale differentiable causal discovery of factor graphs. Advances in Neural Information Processing Systems (NeurIPS), 2022.
[32] Lorch, L., Rothfuss, J., Schölkopf, B., and Krause, A. DIBS: Differentiable Bayesian Structure Learning. Advances in Neural Information Processing Systems (NeurIPS), 2021.
[33] Madan, K., Secttr-Brooks. J., Korablyov, M., Bengio, E., Jain, M., Nica, A., Bosc, T., Bengio, Y.. and Malkin, N. Learning gflownets from partial episodes for improved convergence and stability. arXiv preprint 2209.12782 2022.
[34] Malkin, N., Jain, M., Bengio, E.,, Sun, C., and Bengio, Y. Trajectory Balance: Improved Credit Assignment in GFlowNets. Advances in Neural Information Processing Systems (NeurIPS), 2022.
12

---

[35] Monti, R. P. Zhang, K., and Hyvärinen. A. Causal discovery with general non-linear relation- ships using non-linear ica. Uncertainty in Artificial Intelligence (UAI). 2020.
[36] Mooij, J. M., Janzing, D., Heskes, T., and Schölkopf, B. On causal discovery with cyclic additive noise models. Advances Neural Information Processing Systems (NeurIPS), 2011.
[37] Mooij, J. M., Janzing, D., and Schölkopf, B. From Ordinary Differential Equations to Structural Causal Models: The deterministic case. Uncertainty in Artificial Intelligence (UAI), 2013.
[38] Mooij. M., Magliacane, S., and Claassen. T. Joint causal inference from multiple contexts. The Journal of Machine Learning Research (JMLR), 21, 2020.
[39] Murphy, K. P. Active learning of causal bayes net structure. Technical report, UC Berkeley, 2001.
[40] Pamfil, R., Sriwattanaworachhii, N., Desai, S., Pilgerstorfer, P., Georgatzis, K., Beaumont, P., and Aragam, B. Dynotears: Structure learning from time- sen-ssries data. International Conference on Artificial Intelligence and Statistics (AISTATS), 2020.
[41] Pearl, Causality, Caimbridge University Press, second edition, 2009.
[42] Peters, J., Janzing, D., and Schölkopf, B. Elements of causal inference: foundations and learning algorithms. The MIT Press, 2017.
[43] Peters. J., Bauer. S., and Pfister, N. Causal models for dynamical systems. Probabilistic and Causal Inference: The Works of Judea Pearl, 2022.
[44] Pratapa, A., Jalihal, A. P., Law, J. N., Bharadwaj. A., and Murali, T. Benchmarking algorithms for gene regulatory network inference from single-cell transcriptomic data. Nature methods, 17, 2020.
[45] Qiu, X., Zhang, Y., Martin-Rufing D., Weng. C., Hosseinzadeh, S., Yang, D., Pogson, N., Hein, M. Y., Hoi (Joseph) Min, K., Wang, L., Grody, E. I., Shurtleff, M. J., Yuan, R., Xu, S., Ma. Y., Replogle, J. M., Lander, E. S., Darmanis, S., Bahar, I., Sankaran, V. G., Xing, J., and
Weissman, J. S. Mapping transcriptomic vector fields of single cells. Cell, 185, 2022.
[46] Riba, A., Oravecz, A., Durik, M., Jiménez, S., Alunni, V., Cerciat, M., Jung, M., Keime, C., Keyes, W. M., and Molina, N. Cell cycle gene regulation dynamics revealed by RNA velocity and deep-learning. Nature Communications, 13, 2022.
[47] Saelens, W., Cannoodt, R., Todorov, H., and Saeys, Y. A comparison of single-cell trajectory inference methods. Nature Biotechnology, 37, 2019.
[48] Schiebinger, G., Shu, J., Tabaka, M., Cleary. B., Subramanian, V., Solomon, A., Gould, J., Liu, S., Lin, S., Berube, P., Lee, L.,, Chen, J., Brumbaugh, J., Rigollet, P., Hochedlinger, K., Jaenisch, R., Regev, A., and Lander, E. S. Transpor Analysis of Single-Cell Gene Expression Identifies Developmental Trajectories in Reprogramming. Cell, 176, 2019.
[49] Spirtes, P. and Glymour, C. An algorithm for fast recovery of sparse causal graphs. Social Science Computer Review, 9, 1991.
[50] Tank, A., Covert, I., Foti, N., Shojaie, A., and Fox, E. Neural Granger Causality. IEEE Transactions on Pattern Analysis and Machine Intelligence, 2021.
[51] Tong, A., Huang, J., Wolf, G., van Dijk, D., and Krishnaswamy, S. TrajectoryNet: A Dynamic Optimal Transport Network for Modeling Cellular Dynamics. International Conference on Machine Learning (ICML), 2020.
[52] Tong, S. and Koller, Active learning for structure bayesian networks. International Joint Conference Artificial Intelligence (IJCAI), 17, 2001.
[53] Zhang, K., Huang, B., Zhang, J., Glymour, C., and Schölkopf, B. Causal discovery from nonsta- tionary/heterogeeeous data: Skeleton estimation and orientation determination. International Joint Conference on Artificial Intelligence (IJCAI), 2017.
[54] Zheng, X., Aragam, B., Ravikumar, P., and Xing, E. P. DAGs with NO TEARS: Continuous Optimization for Structure Learning. Advances Neural Information Processing Systems (NeurIPS), 2018.
13

---

Supplementary Material
A Proof of Proposition 1
Proposition calculates the number of admissible strueture graphs for a linear ODE system with correlated variables. We will first show the general case this is Ii 1)).1, then analyze the case of an L penalty on the edges of G, which reduces the size of the of admissible graphs to (m,,
Proof. Consider an identifiable linear system de = Az where we directly observe (x, de with G admissible graph by definition. For each
dt identifiable. Then the system with 1 has exactly node, we analyze its set of child nodes in G. i.c. c(u) = {u E 1/ s.t. a - v E g). For an identifiable system, each child must have an incoming edge from its parent.
at
Next, we consider the process of adding correlated variable, i.e. consider the situation of w.l.o.g. consider m = (s, 1...... 1) for some > 1. Then for each child of e, (21), there are now possible parents. This has multiplied the number of possible graphs by 2 1. Since each element m independent, this leads to the first statement, i.e. | = Ileea
1)0.1.
Under an L penalty, then we constrain the possible graphs to different graphs, where each child node picks exactly one of the possible parents. This leads the second statement, [g'| (m,)).
B Experimental Details
B.1 Single Cell Dataset Preprocessing
We start with the processed data from [46]. We first filter it applying steps from the Sclelo tutorial. We then sub-select the genes of interest and use the "Ms" and "velocity" layers, which we normalize to mean zero standard deviation one for the states and scale the dz with the same parameters.
B.2 Baselines for Bayesian Dynamic Structure Learning
Existing Bayesian structure learning methods are typically constrained to learning DAGs. Tem- poral information about the the dynamic relationships amongst variables in a system can help alleviate this constraint. DiBS and BCD-Nets ing approaches for static systems. We apply versions of DiBS and BCD-Nets such that they applicable in our Bayesian dynamic structure learning framework for learning cyclic graph structure from dynamic data. We use the approach taken in DiBS and parameterize the distribu- tion over graphs Pa. Z) II. II, Po, (Gig Zig), where Z =UVV U, B&xd Here
state-of-the-ar Bayesian structure learn-
Pa, (G., = 1[ = (0x ((x) 1/(1 +e aaarr and a = ac(t) (t denotes the training iteration. We use c(l) = Vl. As 1 1 xc,
xc, Pa. (GIZ) - 8(G] In DiBS, Stein variational gradient decent (SVGD) [30] used iteratively transport particles Z learn the target distribu- tion. Following from the above parameterization, we implement a version of BCD-Nets by treating U I~ NM(p, R.), V. (Ni,, oe), and learning Yuu /ly; au and ov. Since our framework uses dynamic data, we incorporate DiBS and BCD-Nets within the framework (labeled DynDiBS and DynBCD, respectively) leverage dynamic information for Bayesian structure learning of cyclic graphs.
B.3 Hyper- prrper-prrameers for Baselines
For both DynBCD and DynDiBS use A = d across datasets. Since DynDiBS is ensemble based method, we use 1024 samples for the linear and non non-linear synthetic systems and 1000 samples for the cellular system (both training and evaluation). Since DynBCD is variational approach and doesn't 't require parallelized model ensembles, we use a large quantity of samples for training and evaluation. In the of the cellular system, since there is significantly smaller quantity of admissible graphs, we use less samples for DynBCD. We use graph sparsity regularization denoted by
An and temperature parameter I' that scales the magnitude of the likelihood (e.g. f f (dz,, dr)).
14

---

In Table 5 and Table 6 we outline the hyper-paramen parameters we found to yield the most competitive results. We use grid search tune DynBCD and DynDiBS. All baselines are trained for 1000 epochs.
Table 5: Hyper-parameters for DynBCD. We define learning rate as E.
Linear System
Model
Ao
T
Qt
samples
B-DynBCD
0.0001
0.001
0.01
0.1
5000
h-DynBCD
0.0001
0.0025
0.01
2
2000
Non-linear System
Model
Ao
a
samples
&-DynBCD
0.00005
0.01
0.01
2
5000
h-DynBCD
0.0001
0.001
0.01
1
2000
Cellular System
Model
X
T
Ox
samples
B-DynBCD
0.0001
0.001
0.05
0.05
1000
h-DynBCD
0.00001
0.0005
0.1
2
1000
Table 6:
for DynDiBS. We define learning rate Linear System
Model
A
T
Y
e-DynDiBS
0.0025
500
0.01
0.0001
3000
h-DynDiBS
0.0001
3
0.01
0.0001
10000
Non-linear System
Model
Ao
T
Y
f-DynDiBS
0.001
10
0.01
0.0001
3000
N-DynDiBS
0.0001
0.1
0.01
0.0001
10000
Cellular System
Model
A
T
Y
(-DynDiBS
0.0025
1
0.05
0.0001
3000
h-DynDiBS
0.00001
0.1
0.01
0.01
3000
We note that when evaluation on validation and test data for Bayes-SHD and AUC metrics, we hard threshold Pa, ZG) We find that through training this the final 04 is typically small enough in magnitude such that Pa., (G) Z) does not yield full threshold of Z. To this end, we select large when computing the KL metric mimic hard threshold behaviour as experienced during training. We use 0₁ = 1 x 10 In DynDiBS the parameter y helps control separation of particles Z during training. In general, we found DynBCD and DynDiBS baselines are challenging to train and to find hyper-parameter settings with good performance. In part, we believe this is due to the numerous hyper- papem-taraeeters required to tune as well as the general difficulty of the objective.
B.4 Neural Network Architectures and Hyper-parameters
We parameterize (s:sss::-;;;)) and he with MLP architectures. (s takes the current state as input and first computes common representations using a layer MLP. Then a 2 layer MLP with a softmax output activation takes the representations as input and outputs a distribution over possible actions. The latter MLP used to parameterize one head for each distribution Pp(si[8-11 v). We use hidden unit dimension of 128 and leaky rectified linear unit (Leaky ReI.U) activation functions for the 1=(8:11-1;;8 v) MI. architecture. We use a uniform backward policy for Ps(s.-1 1ss v) To parameterize hgg we use a layer MLP with hidden layer dimensions of
15

---

Table Hlper-parameterr for DynGFN. We define learning rate as c, Mirrin as number of training
samples, and Meval the number of evaluation samples.
Linear System
Model
A
T
Maraón
Meval
B-DynGFN
0.0001
100
0.01
1024
5000
h-DynGFN
0.00001
100
0.005
256
3000
Non-linear System
Model
Ao
T
Marain
Mexal
2-DynGFN
0.0001
150
0.01
1024
5000
h-DynGFN
0.00001
150
0.005
256
3000
Cellular System
Model
Ao
T
Marain
Meval
B-DynGFN
0.00005
45
0.01
1024
1000
h-DynGFN 0.0001
10
0.1
1024
1000
(64, 64, 64) and exponential linear unit activations (ELU). We consider two parametrizations for fo: single linear parameters, i.e dæ: = Or, and a single hidden layer neural network da = =f)()). We use these parametrizations model linear and non-linear node-wise parent-child structural models, where x E Idd are the node-wise input observations.
B.5 Hyper-parameee for DynGFN
DynGFN requires setting variety of hyper-parameters that lead different trade off's in model performance. In particular, Ao (sparsity encouragement for identified graphs), a temperature parameter
T' that scales the magnitude of the reward likelihood (e.g. + MSE(dz, dr) )) learning rate softmax tempering e (we always use cosine schedule for c, with discrete period of 5 epochs), and number of training epochs. In our experiments we select hyper-parameter values that we find lead competitive performance (this pertains to 2-DynGFN and it- li-DynGFN models) by observing performance over a few hyper-parameter values. We outline the selected syper-parameters for each respective model in Table For GFlowNet sampling, due to computational limits, we use less training samples Mtrain than evaluation samples meval for DynGEN. We train DynGFN for 1000 epochs.
B.6 GFlowNet Exploration vs. Exploitation
The general procedure for training GelowNets is inspired from reinforcement learning where the primary objective is to learn a stochastic policy Tuass) to sample actions from an action space given current state. In our setting, the action space represents possible locations where an additional edge can be placed to an existing graph and each state is represented by a current graph. Since under this training procedure attributing a reward associated to the sampled state/graph, the policy is susceptible to exploitation: Pr(s:|s;-1:ve v)) samples graph(s) with a high reward, it becomes easy for the policy to focus sampling said graphs since they yield high reward. To alleviate this we encourage exploration using softmax tempering
are sampling from the GFlowNet policy (f(s)|s--i
at every iteration then
stochastic policy, by multiplying the logits of our forward policy by 1/c before applying the softmax function. A larger flattens the stochastic policy such that exploration within the action space encouraged. However, setting the parameters is challenging and there exists a trade-off between exploring and exploiting the stochastic policy during optimization. We address this by using a cosine schedule for such that 1 1000. We treat the period of the cosine schedule as a hyper-parameter
B.7 Discussion of Training Dynamics
GFlowNets are a relatively recent class of models that can be challenging to optimize. We discuss some of the challenges with training them especially in the context of a learned energy function.
16

---

0.08
MSE
Bayes-SHD
175
0.07
0.06
150
0.05
125
0.04
100
0.03
75
0.02
50
0.01
200
400
600
800
1000
200
400
600
800
1000
epoch
epoch
Figure 5: Validation curves throughout training of (-DynGFN on linear system. Mean squared error (left) and Bayes-SHD (right).
We observed that in settings where the energy reward is fixed and we could proportionally penalize missing edges as well as the addition of incorrect edges (c.g. Ł-DynGFN), we were able to better learn posteriors over admissible graphs over models that require sparse priors and/or trainable energies. This suggests that DynGFN may be limited by an inadequate energy reward. However, we found training DynGFN with a trainable energy function challenging since the GFlowNet stochastic policy depends on the rewards, and vice versa. Further investigation and experimentation into this alternating optimization procedure is required. To in part alleviate the constraints of this challenge, in this work we presented results for h h-DynGFF which used initialized GFlowNet parameters from trained effective
(-DynGFN model the same task. We found that this initialization technique enabled iterative optimization of the forward policy and energy reward. In Figure we show validation curves throughout training of E-DynGFN linear system with 20 and sparsity value of 0.9.
B.8 Evaluating Quality Learned Posteriors
Using our indeterminacy model defined in section 5, can determine P(G*) for given set of correlated variables, where G" denotes the set of true equally admissible structures. Here we assume P(G*) uniform distribution over G and determine the KL between Q(G) and P(G"). We compute the KL considering only the probabilities of trained models to generate all structures in G+. ie. Q(G*). This due the computational constraints for calculating the KL for DynGFN since even for sparse graphs of moderate size this is a combinatorial computation. Nonetheless, this approach allows to directly compare the learned posteriors QG)) to ground truth discrete distribution over structure G to evaluate the effectiveness of Bayesian structure learning approaches.
B.9 Implementation Details
Our model is implemented in Pytorch and Pytorch Lightning and is available at https:: h/go::///it thub. aalaaaaaaaaaddnnaaaaa Models were trained on heterogeneous mix of HPC clusters for total of - 1,000 GPU hours primarily on NVIDIA RTX8000 GPUs.
17

---

C
Additional Details
C.1
Single-cell Biology and Gene Regulatory Network Inference
C.1.1 Gene Regulatory Networks and Cell Dynamies
One dynamical system of interest is that of cells. Cellular response to environmental stimuli or genetic perturbations can be modelled as a complex time-varying dynamical system [17. 1]. In general, dynamical system models are a useful tool for downstream scientific reasoning. In this work we are primarily interested in identifying the underlying cell dynamics from data. A reasonable model for cell dynamies is as a stochastic dynamical system with many, possibly unobserved, components. There are many data collection models for gaining insight into this system from single-cell RNA-sequencing data. We will primarily focus on RNA velocity type methods, where both x and an estimate of da are available in each cell, but note that there are other assumptions to infer dynamics and regulation such as pseudotime-
methods [47, 1J, and optimal transport methods [17, 48, 51, 22, 23]. After learning possible explanatory regulation, this is used in downstream tasks, but the resulting conclusions drawn from these models are nccessarily conditional on the inferred regulation. Motivated by gene regulatory networks, we explicitly model uncertainty over graphs which allows us to propagate the resulting uncertainty to downstream conclusions.
C.1.2 Learning Gene Regulatory Networks From Single-cell Data
Single-cell transcriptomics has an interesting property in that from a single measurement we can estimate both the current state and the current velocity dr. Because mechanistically RNA undergoes a splicing process, we can measure the quantities of both the unspliced (carly) and spliced (late) RNA in the cell. From these two quantities we can estimate the current RNA content for each gene and the current transcription rate. There exist many models for denoising and interpreting this data [27, 9, 45]. Furthermore, there exist more elaborate measurement techniques extract more accurate velocity estimates [45]. The fact that we have an estimate of the current velocity is exceptionally useful for continuous time structure discovery because it allows us to avoid explicitly unrolling the dynamical system.
Learning the underlying causal structure from data is one of the open problems in biology. There many works that attempt to learn the effect of a change in a gene, or the addition of a drug. These works often build models that directly predict the outcome of an intervention. This may be useful for certain applications, but often does not generalize well out of distribution. We would like to learn a model of the underlying instantaneous dynamics that give rise to effects longer time scales. This approach has a number of advantages. (1) it is closer to the mechanistic model; it may be easier to learn model of the instantaneous dynamics rather than the dynamics over long time scales (details in Appendix C.2). (2) One model can be trained and applied to data from many sources including RNA-velocity, Pseudotime, Single-cell time series, and steady state perturbational data. (3) The instantaneous graph may be significantly sparser (and therefore easier to learn) than the summary graph or the equilibrium graph.
C.1.3 Further Details on Experiments with Single-cell RNA-velocity Data
The process of Eukaryotic cell division can be divided into four well regulated stages based on the phenotype, Gap 1 (G1), Synthesis (S), Gap 2 (G2), and Mitosis (M). This process is a good starting point for GRN discovery as it is (I) relatively well understood, (2) deterministic, and (3) well studied with plentiful data. While there is an underlying control loop controlling the progression of the cell cycle, there are many other genes that also change during this cycle. To rediscover the true control process from data we must disentangle the true causal genes from the downstream correlated genes. This may become very difficult when we only observe dynamics at longer time scales. We hide a cell cycle regulator among two downstream genes that are highly correlated (Spearman p > 0.75) and test whether we can model the Bayesian posterior namely that we are uncertain about which of the three genes (Cdc25A, Mem2, Mcm5) is the true causal parent of Cdkl.
18

---

C.2 Instantaneous Graph and Long-horizon Graph
The graph recovered depends on the time scale considered. We make distinction between the conditional structure the graph based on the time scale. Consider a system governed by dt
da(
f(x(t)). We define the instantaneous graph as:
(2))))))
gi.j
U,1
i= ,.....,,, d.j....jj..... ;.....; d.
(6)
dx;
Here, 1(2) =1 if z * 0, otherwise ((2) 0. We define the temporal summary of the system drift / as:
I. f(x(t))dt.
F(T, x) To +
Then, we can define a long-horizon graph over long time-scales as:
:= (OF(T.z),
i= i= L= d,j I....,dd dd
(8)
Ox:
are integrating the drift f over long time-scales, the long-horizon graph may be less sparse than the instantaneous graph. In cellular systems, this equates to observing cell dynamics over long time-scales, in turn observing increase quantities of correlations between variables of the system. Thus, trying to delineate the instantaneous dynamics from long time-scales may be difficult depending on the underlying system dynamics.
C.3 Further Discussion on Future Work
Although DynGFN is currently limited to smaller systems, we foresce approaches that would enable some degree of scaling DynGFN larger systems. One approach is to leverage biological information of known gene-gene connections as more informative prior for DynGFN. Currently, DynGFN learns a forward stochastic policy for (G) starting from an initialized state so of all zeros. Instead, we can define šu using prior of high confidence biological connections and sequentially add edges starting from this new initial state. This would reduce the number of possible structures DynGFN would need to search over, thus improving the potentially scalability of DynGFN. Another approach to learn structure between sets of genes (variables) rather than single genes. Since GRNs are generally very sparse, it makes sense group genes in sets. Consequently, we can then learn structure between these grouped genes, rather than just individual genes. In turn, DynGFN can explore/learn the structure over smaller space while effectively capturing structure between a significantly larger set of genes. To group genes, we would use prior biological information, either form existing literature expert domain knowledge.
In this work we exploit the use of a minimal prior, i.e. L sparsity prior, for learning Bayesian dynamic structure between variables. In general, the aforementioned approaches for scaling DynGFN to larger systems involve use of more informative priors on G. Although we mention two ways we foresee approaching this in the biological context of GRN inference, the general approach of using more informative priors can help scale DynGFN larger systems across applications.
C.4 Broader Impacts
While it is important to acknowledge the potential risks of drawing incorrect scientific conclusions due incorrect assumptions, our work embraces Bayesian perspective for structure learning. A key component of our work is to account for uncertainty within our method, aiming to minimize the chances of incorrect conclusions. It is important to note that the accuracy of conclusions relies on applying the method in settings that align with the underlying assumptions, such as causal sufficiency and the use of dynamic observational data. By adhering these guidelines, our approach holds promise for producing robust and reliable scientific outcomes.
19


Reviews:
[{"review_id": "Wjz87nGyMI", "reviewer": "Reviewer_PVQy", "summary": "The authors propose to learn dynamical systems with cyclic dependencies. They factorise the generative model using a variant of the GFlowNet model, a HyperNetwork and MLPs.", "questions": "To illustrate the advantages of the proposed approach:\n1.   use simple systems with few nodes, with and without cycles, and show how well the true dependencies are recovered. \n2. show a clear example where failing to model the cyclic dependencies results in a poor fit.  \n3. show how to use the multimodal representation capacity when explaining the cause dependencies in a gene regulatory network.", "limitations": "Yes.", "rating": 7, "confidence": 4, "soundness": 2, "presentation": 2, "contribution": 2, "strengths": "Positives include:\n1. extending the formulation from DAGs to cyclic graphs\n2. introducing a per-node posterior formulation  to improve the computational complexity\n3. using GFlowNets to address the multimodality issues", "weaknesses": "Negatives include:\n1. due to the high density the paper is hard to understand: a lot is explained and detailed in the Appendix\n2. it is hard to gauge the quality of the results: even for the artificial case it is hard to assess intuitively if the approach solves the problem in a satisfactory way. Perhaps the authors should use simple systems with few nodes, with and without cycles, and show how well the true dependencies are recovered. \n3. it is unclear how important it is to being able to model cyclic structures: where are the cyclic dependencies in the examples shown in Fig. 2 and Fig. 3? The authors should show a clear example where failing to model the cyclic dependencies results in a poor fit.  \n4. it is unclear how to translate the declared advantage of being able to model multi-modal posteriors  in practical terms: how would we make use of this representation capacity when explaining the cause dependencies in a gene regulatory network?", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We thank the reviewer for their comments and insightful feedback. Please refer to our feedback summary and the attached document regarding questions about cyclic dependencies in Fig. 2 and Fig. 3. We also include further discussion regarding the toy example in the feedback summary. We now address each comment individually.  \n\n>Perhaps the authors should use simple systems with few nodes, with and without cycles, and show how well the true dependencies are recovered.\n\nThe reviewer suggests a valuable addition. To convey this, we introduced a toy example with a simple 3 variable system (with and without cycles), see Fig. 4 in the attached document. We consider a method for learning acyclic graphs (NOTEARS [1]) which does not model cyclic dependencies, unlike DynGFN. It is clear from this example that NOTEARS struggles to recover cyclic dependencies. We will add Fig. 4 to the manuscript. We can also verify this result by considering a conditional independence test over cyclic dependencies. With reference the attached Fig. 4 it is easy to see that the conditional independence test fails in the cyclic setting: In the acyclic case, we can identify the v-structure by observing that $x_1 \\perp x_3$ and $x_1 \\not\\perp x_3 | x_2$, which implies that $x_2$ is a collider (i.e. $x_1$ and $x_3$ are marginally independent and conditionally dependent); in the cyclic example, we introduce time dependencies such that there are cycles in the summary graph that render these variables marginally dependent.\n\n>it is unclear how to translate the declared advantage of being able to model multi-modal posteriors in practical terms\n\nLearning multi-modal posteriors over structure allows us to quantify uncertainty over a particular causal dependency. In turn, it makes it easier to answer the active learning question, “how should we select interventions such that we minimize uncertainty?” [2]. We use this as motivation for our work and cite relevant papers (see lines 54-55). We will add the aforementioned citation to this list. As a first step, we need to effectively learn the multi-modal posterior over structure, hence the motivation and focus of this work. \n\nWe hope our response fully addresses all the important and salient points raised by the reviewer. We believe the new example improves the clarity of our work, and we thank the reviewer for their excellent suggestion. Due to the increased clarity of these new examples, we ask that the reviewer consider raising their score. If you have any additional comments and suggestions we would be happy to discuss further. \n\n[1] Zheng, Xun, et al. \"Dags with no tears: Continuous optimization for structure learning.\" Advances in neural information processing systems 31 (2018)\n\n[2] Toth, Christian, et al. \"Active bayesian causal inference.\" Advances in Neural Information Processing Systems 35 (2022)"}, {"review_id": "95sghdCVZl", "reviewer": "Reviewer_b9zX", "summary": "Authors introduced a principled methodology called DynGFN, which effectively identified cyclic structures and concurrently modeled the Bayesian posteriors over directed acyclic graphs (DAGs). Leveraging RNA velocity, the authors formulated a dynamic system that unveiled the underlying gene regulatory networks. DynGFN was meticulously designed with three modules, and its superior performance was demonstrated through synthetic experiments and real-world analysis.", "questions": "Major concerns\n1. The synthetic experiments and real analysis did not explicitly present any results related to feedback loop identification, which serves as one of the motivations behind DynGFN. Could you elaborate on this point?\n2. How many epochs did it take to stabilize the training of l-DynGFN?\n3. Given the inference of RNA velocity involves multiple free parameters, how do these parameters impact the inference results obtained from DynGFN?\n4. In the synthetic experiments simulating dynamic systems, how were the ground truth regulatory networks constructed based on the simulated dynamic systems?\n5. Appendix B.3 mentions the construction of the validation and test sets to fine-tune the hyperparameters. Could you explain how these sets were constructed for the experiment?\n6. Regarding Table 4, could you clarify what the ground truth graph was used for calculating Bayes-SHD and AUC? Was it formulated based on external knowledge?\n7. In Figure 3, it is important to note that gene regulation is not solely determined by the direct interaction between corresponding proteins. In fact, it has been reported that \"CDK1 targets MCM2-7 complex...\" [1].\n\n[1] Enserink, Jorrit M., and Richard D. Kolodner. \"An overview of Cdk1-controlled targets and processes.\" Cell division 5.1 (2010): 1-41.\n\n\n\nMinor concerns\n1. Figure 1 lacks definitions for numerous parameters, which may cause confusion in understanding the illustration.", "limitations": "The authors mentioned two limitations: 1) scalability issues concerning larger systems and 2) hyperparameter tuning. To address these challenges, the authors suggested employing more informative priors or grouping genes together as nodes.", "rating": 6, "confidence": 4, "soundness": 4, "presentation": 2, "contribution": 3, "strengths": "The real biological systems can be rarely formulated as DAGs and there are always feedback loops to make the system work. DynGFN was appropriately motivated with a real biological thinking.", "weaknesses": "The overall manuscript was hard to follow as there was too much content included. There were still a lot of technical details that needed to be clarified. It would improve substantially if the manuscript can be carefully restructured.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We thank the reviewer for their comments and insightful feedback. Please refer to our summary to all reviewers and the attached document regarding questions about cyclic dependencies in our synthetic and real experiments, as well as for clarifications regarding Fig. 1. We also introduced a new toy example to exemplify our model’s capability for learning cyclic dependencies (See attached Fig. 4 and response to Reviewer [PVQy]). We now address each salient point individually.\n\n>How many epochs for l-DynGFN to stabilize?\n\nTo show when l-DynGFN stabilizes during training, we have introduced validation curves for mean squared error and Bayes-SHD over the course of training on the linear system with $d=20$ (See Fig. 5 in the attached). In this experiment I-DynGFN stabilizes in around 500 epochs.\n\n>How do RNA velocity parameters impact the inference results obtained by DynGFN?\n\nThe reviewer points at a great question. We agree that these free parameters may impact the results of DynGFN. DynGFN currently takes RNA velocity as input, but an ideal solution would incorporate this source of uncertainty (uncertainty due to RNA velocity inference) into the posterior over the graph and the parameters. Due to the difficulty of the cyclic structure learning problem and given that DynGFN already requires multiple posteriors, we leave investigation into a fully Bayesian RNA velocity and gene network recovery method to future work. We will add discussion of this to the limitations section. \n\n>How were the ground truth regulatory networks constructed for the synthetic dynamical systems?\n\nWe first construct the ground truth regulatory network (GRN) and then use this ground truth GRN to simulate the system dynamics. Consider the linear system $dx = Ax$. To construct the system, we randomly sample $A$ such that $A \\in \\mathbb{R}^{d \\times d}$ with a specified sparsity. We constrain $A$ to follow some properties to ensure the system is stable. Specifically, we subtract the maximum eigenvalues from the diagonal of $A$. This $A$ is then used to simulate $dx = Ax$. The procedure follows for the non-linear system with the addition of the non-linearity. We will include this description in the Appendix. For more specific implementation details see line 58 of `src/datamodules/simulated_datamodule.py` in the attached code.\n\n>How were the validation and test sets constructed for the experiments? How was the ground truth graph in the real data setting formulated?\n\nIn the synthetic experiments, we simulate a system and generate train, validation, and test observations $(x, dx)$ for a given system with true connections defined by $A$. In the real data setting, we split the observed data $(x, dx)$ into train, validation, and test datasets. $A$ in the real data setting is constructed using prior biological knowledge. Specifically, we extracted a subset of the gene network from the KEGG cell cycle pathway entry hsa04110. Fig. 3 (see attached document for updated version) shows the constructed ground truth graph. \n\n>It is important to note that gene regulation is not solely determined by the direct interaction between corresponding proteins.\n\nWe agree, we will add a mention of this in Section 6.3, and to the future work section. We find this quite interesting and gets back to what gene regulatory network we are trying to discover. Our model is searching for a posterior over sparse gene regulatory networks which explain the data dynamics. It could be argued that as the time goes to zero, we will find only direct interactions, however we leave further investigation of this to future work. \n\nWe thank the reviewer for their time and effort in reviewing our work and we hope the reviewer would kindly consider a fresh evaluation of our work given the main clarifying points outlined above. If you have any additional comments, we would be happy to discuss further. "}, {"review_id": "AW2RmD45D0", "reviewer": "Reviewer_uDCE", "summary": "The authors tackle the well-known problem of Gene Regulatory Network (GRN) inference from time series data. Starting from an estimate of RNA velocity, GRN inference is framed as a causal discovery task, with two specificities: the inferred graphs may be cyclic, and a probabilistic distribution over candidate graphs is sought rather than a point estimate, so as to account for the high amounts of noise inherently present in biological data. The method is articulated around two assumptions:\n\n- a dynamic structural model $\\frac{dx_i(t)}{dt} = f_i(\\text{Pa}(x_i), \\epsilon_i)$ for each gene. $\\text{Pa}(x_i)$ is directly obtained from the regulatory graph $G$, and the structural causal model (SCM) $f_i$ is parametrized by $\\theta$\n- a factorisation of the joint distribution of the graph $G$, the structural causal model parameters $\\theta$, and the observed data $\\mathcal{D}$ as $p(G,\\theta,\\mathcal{D}) = p(\\mathcal{D}|G,\\theta)p(\\theta|G)p(G)$. The graph sampler $p(G)$ is modeled as a GFlowNets. Next, since the parameters of the SCM depend on the graph structure, they are computed by a HyperNetwork, i.e. $p(\\theta|G)$ is substituted to a neural network which takes as input a graph $G$ and outputs the structural equation model parameters $\\theta$: $p(\\theta|G) = \\delta(\\theta|G)$.\n \nNext, the authors evaluate the proposed method, DynGFN, on several synthetic examples as well as a real-world one involving single-cell RNA velocity data. The authors introduce other baselines by using different approaches than a GFlowNets for the graph sampler $p(G)$. This being said, DynGFN offers substantial improvements when it comes to jointly recovering the ground truth structure while characterizing the uncertainty around it.", "questions": "- Any insights on why $h$-DynGFN yields significantly better results than $\\ell$-DynGFN on the linear system presented in table 1?\n\n-  In the single cell experiment, where does the ground truth network (used to compute AUC / Bayes-SHD) come from? Could you explain a little more what does \"correlation over cell cycle time\" means in Section 6.3?\n\nI have read the rebuttal done by the authors.", "limitations": "The authors mention as main limitation the difficulties to scale to large gene regulatory networks due to the combinatorial explosion. Causal sufficiency is also assumed here, which means that all relevant variables are observed. This assumption hardly ever holds in real-world biological problems. This issue is also mentioned by the authors in the Appendix.", "rating": 6, "confidence": 3, "soundness": 2, "presentation": 2, "contribution": 3, "strengths": "The big picture of the method is well explained. The work described here is likely to have a significant impact for the GRN inference community, as I believe that not many probabilistic methods are available in this field. I also believe that DynGFN could be used as a starting point for the more general problem of physical/biological dynamical system discovery, where one is not only concerned with the graph adjacency matrix but with providing a mechanistic description of the structural causal model $f_i$ [1,2]\n\n[1] Discovering governing equations from data by sparse identification of nonlinear dynamical systems, PNAS, 2016\n[2] Identification of dynamic mass-action biochemical reaction networks using sparse Bayesian methods, PLoS Comp. Biol., 2022", "weaknesses": "- The paper builds on GFlowNets for the graph sampler $P(G)$ and Hypernetworks to generate the structural causal model parameters $\\theta$ given a sampled graph $G$. I found it slightly difficult to understand everything $in$ $detail$ given that these tools are each introduced only in a small paragraph.\n\n- Even though the method focuses on Bayesian Inference, it would have been interesting to assess DynGFN's performances compared to SOTA methods for GRN inference from time series data such as DynGENIE3 [1], BINGO [2] or others. Lastly, these methods use time-resolved gene expression data, from which one can get an estimate of the RNA velocity using for instance finite differences, which may be a quite rough estimate for sparsely sampled measurements. Do I understand right that the RNA velocity data employed in this submission somehow constitutes a ``more principled'' way to estimate the RNA velocity?\n\n[1] dynGENIE3: dynamical GENIE3 for the inference of gene networks from time series expression data, Scientific Reports, 2018\n[2] Gene regulatory network inference from sparsely sampled noisy data, Nature Communications, 2020", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We thank the reviewer for their comments and insightful feedback. We address your comments in what follows.\n\n>The paper builds on GFlowNets for the graph sampler and Hypernetworks to generate the structural causal model parameters. I found it slightly difficult to understand everything given that these tools are each introduced only in a small paragraph.\n\nThank you for bringing up this issue, we realize our work is quite dense and the page limit was an issue. We will of course expand the exposition in the main text with the addition of an extra page. We are happy to incorporate any specific suggestions you have on how to improve the exposition.\n\n> Even though the method focuses on Bayesian Inference, it would have been interesting to assess DynGFNs performance to SOTA methods from time-series data.\n\nWe are also very interested in how DynGFN compares to and can be combined with SOTA methods on inferring GRNs from time-series data! As noted by the reviewer, methods such as DynGENIE3 and BINGO consider a different setting where we have access to multiple samples from the time series. They are generally not well suited to velocity measurements, which effectively amounts to having two extremely close time points, as these methods often using changes over time to infer the gene regulatory network dynamics. We believe adapting these state of the art methods to infer regulation from velocity data would be an interesting direction given the potential quantity of RNA-velocity analysis, but we believe this is out of scope for this project. Future work on how to adapt DynGFN to the time series setting, potentially using ideas developed in [3], is also extremely promising. \n\n>Does the RNA velocity used in this submission constitute a more ‘principled’ way of estimating RNA velocity?\n\nWe do not explicitly estimate RNA velocity in this work. Rather, we use estimated RNA velocity (using scVelo [1]) to help formulate the problem of Bayesian structure learning over cyclic graphs using observational data consisting of dynamic tuple pairs $(x, dx)$. We leave inference of RNA velocity using DynGFN for future work. See response to Reviewer [b9zX] for further discussion.\n\n>Insights on why h-DynGFN yields better results compared to l-DynGFN on the linear system?\n\nThe reviewer brings up an insightful question. We suspect that l-DynGFN may be more sensitive to the stochasticity between batches of data compared to h-DynGFN since the analytic linear solver used in l-DynGFN directly uses the minibatch to solve for the parameters (the factorization P(\\theta | G) is poorly enforced). This means that there is greater stochasticity in the reward depending on the batch selected. We believe this is why we see that l-DynGFN yielding worse results on the linear system compared to h-DynGFN.\n\n>Where does the ground truth network in the single-cell experiment come from?\n\nThis ground truth network is constructed using external prior biological knowledge. Specifically, we extracted a subset of the gene network from the KEGG cell cycle pathway entry hsa04110.\n\n>What is meant by ‘correlation over cell cycle time’?\n\nThis data is from Riba et al. 2022 [2]. Here they supply public data with a cell cycle pseudo-time label called `cell_cycle_theta` in the public data. We correlate bin this label into 10 bins and approximate the correlation of gene expression over time.\n\nWe would like to thank the reviewer for their review of our paper. We believe we have answered all the great points raised by the reviewer in our author response, and we kindly ask the reviewer to consider upgrading their score if the reviewer is satisfied with our responses. Please let us know if you have any additional feedback or comments. We would be happy to discuss. \n\n[1] Bergen, Volker, et al. \"Generalizing RNA velocity to transient cell states through dynamical modeling.\" Nature biotechnology 38.12 (2020)\n\n[2] Riba, Andrea, et al. \"Cell cycle gene regulation dynamics revealed by RNA velocity and deep-learning\" Nature Communications (2022)\n\n[3] Tong et al. \"Simulation-Free Schrödinger Bridges via Score and Flow Matching\" ArXiv (2023)"}, {"review_id": "XE43JdQd9p", "reviewer": "Reviewer_2wco", "summary": "In this work, the authors propose an extension of GFlowNets to enable posterior sampling of cyclic graph structures and apply their approach to infer gene regulatory networks (GRN). While existing methods are able to infer cyclic graph structures or sample from Bayesian posteriors over DAGS, they cannot do both simultaneously. The authors formulate the problem as a dynamical systems identification problem which can be factored into the priors and model likelihood. The priors are learned using either linear or non-linear Hyper Networks while the model likelihood is estimated using the GFlowNets framework with detailed balance loss. The authors show their model is robust to several parameters including edge sparsity and time intervals between data points. Furthermore, they introduce several optimization ideas enabling a tractable search over admissible graph structures. Finally, they benchmark their approach against several baseline methods and apply their approach to an RNA velocity dataset to show their method is able to recover known and putative gene-gene interactions. \n", "questions": "see above.", "limitations": "see above.", "rating": 6, "confidence": 3, "soundness": 3, "presentation": 2, "contribution": 3, "strengths": "The paper addresses an interesting/important question which has been the center of much research for decades: posterior sampling over directed graphs, this time with cycles and not DAGs. The authors appear knowledgable in this field and they take advantage of GFlowNets, building a method to assess posterior structures with cycles as they unfold them across time.", "weaknesses": "Our biggest grievance is that the paper is just not reader friendly for anyone who is not an expert in that specific line of work. \nAs noted the work seems technically very solid, one that offers specific new solution for the stated problem. Still, we as researchers not in this specific field, remain unsure after reading it of the overall innovation/significance of creating posterior sampling of structures over GFlowNets. We add more details below.\n\n(1)\nThe termination criteria of Algorithm 1 is not clear. Under what conditions does the transition probability become null? Although one of the main innovations of the paper is sampling from the posterior of GFlowNet structures, the authors don’t reference any standard Bayesian analysis. How are posterior samples summarized (for example, how is a point estimate of p obtained from posterior samples in Figure 3)? How and when does Algorithm 1 converge to the posterior? Is the approach sensitive to variable initializations? \n\n(2)\nAlthough the authors compare their approach to sensible baselines using sensible metrics, it is not clear the approach has a clear advantage over the original GFlowNet in this applied setting. For example, how does the MAP parameter estimate (or equivalent)  of this method compare to the original GFlowNet under the AUC metric (these are presumably comparable). Does GFlowNet produce different insight on gene-gene interactions in the RNA velocity dataset?\n\n(3)\nThe authors assume causal sufficiency as a criteria for model identifiability. Is this a reasonable assumption given the sparsity of single cell data and limitations in the scalability of the model? In other words, how can we be reasonably certain that all casual gene measurements are observed given the model effectively only functions on a gene space of size 20 or lower?\n\n(4)\nSection 5 on graph augmentation is very abstruse. Efforts could be made to improve readability/accessibility for an audience that doesn’t work directly in this field. The related Figure 2 is also very cryptic. \n\n(5)\nIt is not clear or obvious how the factorization into the “per-node” posterior in section 4.1 leads to the reduction in search space from 2^d^2 to d2^d. A brief explanation would be useful.\n\n(6)\nLine 154: “Previous work has shown GFlowNets are useful in settings with multi-modal posteriors.” This needs a citation.\n\n(7)\nThere is inconsistency in the notation in Figure 1 which leads to some confusion. Specifically, it appears the Q variable is overloaded. How is the graph prior (represented as an adjacency matrix) sampled? Is there a hyperprior that is used or are graph structures uniform? \n\n(8)\nThe patterns that the authors suggest that are evident in the Figure 3 heatmaps are not convincing. Perhaps it would be better to sort both heatmaps the same way instead of independent clustering.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We thank the reviewer for their comments, insightful feedback, and helpful suggestions to improve the manuscript. Please refer to our feedback summary and the attached document regarding questions/clarifications on Fig. 1, Section 5 and Fig. 2, and Fig. 3. In the following, we address individual comments/questions brought up by the reviewer to improve the overall exposition quality of the manuscript. We are happy to incorporate any additional suggestions and feedback. \n\n>The termination criteria of Algorithm 1 is not clear. How and when does Algorithm 1 converge to the posterior?\n\nWhile the stopping criteria is not explicit, it is part of the GFlowNet policy $P_F(s_i | s_{i-1}; \\psi)$. See [1] for further details and for sufficient conditions under which the optimal GFlowNet policy (which includes the stopping criterion) correctly converges to the posterior. \n\n>How are posterior samples summarized? How is a point estimate of p obtained from posterior samples in Fig. 3?\n\nWe are not certain which $p$ is being referred to, but any statistics of the posterior are computed by either Monte Carlo (sampling from the GFlowNet policy) or using log-likelihood evaluation of the model. In Fig 3 there is a $\\rho$ which quantifies correlation between genes used to determine the ground truth graph. Is this what is being referred to?\n\n>Is the approach sensitive to variable initialization?\n\nWe do not think DynGFN is very sensitive to variable initialization. Examining Tables 1-4 in Section 6 we see that the standard deviation over 5 seeds is much smaller than the differences between methods. This suggests that DynGFN is relatively stable in that it outputs similar performing posteriors over different initializations. \n\n>Section 5 on graph augmentation is very abstruse.\n\nTo clarify this section, we have amended Fig. 2 (see attached document) and included more informative labels and captions to correspond to the language used in Section 5. \n\n>It is not clear that the approach has a clear advantage over the original GFlowNet in the applied setting? \n\nA key contribution of this paper is on how to leverage time and velocity information to address Bayesian structure learning over **cyclic** graphs representing dynamics, while DAG-GFN [2] (we assume this is what is meant by original GFlowNet) only works in the acyclic case. Therefore, it is not directly applicable to graphs with cycles, as are known to frequently occur in gene regulatory networks. We refer the reviewer to Fig. 4 in the attached PDF where we show an explicit example where using the DAG assumption for cyclic structure learning fails. \n\n\n>Is the causal sufficiency assumption a reasonable assumption..?\n\nThe reviewer identifies a great point. We agree that the assumption of causal sufficiency likely fails for single-cell data with only 20 genes. However, as we grow the number of genes that we are able to infer over, this assumption becomes more plausible. Ideally, we want an approach that yields (1) good posteriors, (2) good identifiability, and (3) good scalability, but due to the significant challenges associated with each of them, we have devoted our effort to address the first two in this work. If GFlowNets (or some other inference procedure) improved scalability, DynGFN could immediately leverage this advancement to do inference over a larger number of genes, which would make the assumption more plausible. \n\n>It is not clear or obvious how the factorization into the “per-node” posterior in section 4.1 leads to the reduction in search space.\n\nThe per-node factorization in effect uses $d$ GFlowNets to learn $d$ graphs each of size $(d \\times 1)$ for $dx_i$ and its parents $x$ (ref equation 4). This is possible because the posterior is factorizable into independent components. This is an advantage of DAG-GFN [2] as this is not possible in the DAG setting. Each of the d graphs are then aggregated to form the full $(d \\times d)$ graph $G$. In this case, each of the $d$ GFlowNets trains over a search space of $2^d$ combinations. Therefore, the state space of the per-node factorized DynGFN is $d \\cdot 2^d$. In contrast, if a single GFlowNet is used, then there are $2^{(d^2)}$ possible states. We will include this explanation in section 4.1. Additionally, we thank the reviewer for noting this and bringing it to our attention. We correct line 191: $d \\cdot 2^d$ = $20 \\cdot 2^{20} ~= 2^{104}$ for $d=20$ should be $d \\cdot 2^d ~= 2^{(4.3 + 20)}$ for $d=20$.\n\n>Previous work has shown GFlowNets are useful in settings with multi-modal posteriors. This needs a citation.\n\nWe thank the reviewer for noting this. We agree, we will add citations [1-4] to this statement.\n\n>Figure 3\n\nSee the attached PDF for an updated Figure 3 with the same ordering. To make our point more obvious, we additionally include a histogram of correlation values for the full correlation and the correlation over cell cycle time. We can see that the (absolute value) of the correlation over cell cycle time is substantially higher on average in distribution.\n\nWe thank the reviewer for their insightful and helpful feedback on our paper. We hope that our rebuttal fully addresses all the important points raised by the reviewer. If our responses, alongside the accompanying additions, improve the overall quality of our work, we hope the reviewer would kindly raise their score. We are more than happy to answer any further questions.\n\n[1] Bengio, Yoshua, et al. \"Gflownet foundations.\" arXiv preprint arXiv:2111.09266 (2021).\n\n[2] Deleu, Tristan, et al. \"Bayesian structure learning with generative flow networks.\" Uncertainty in Artificial Intelligence. PMLR, 2022)\n\n[3] Malkin, Nikolay, et al. \"Trajectory balance: Improved credit assignment in gflownets.\" Advances in Neural Information Processing Systems 35 (2022)\n\n[4] Madan, Kanika, et al. \"Learning GFlowNets from partial episodes for improved convergence and stability.\" International Conference on Machine Learning. PMLR (2023)"}]

Author rebuttal:
We would like to thank all the reviewers for their time and insightful feedback when reviewing our paper. We appreciate the constructive criticisms and suggestions which will serve to improve the overall quality of our paper. 

Our method focuses on Bayesian structure learning over cyclic graphs for application in gene regulatory network (GRN) inference. We consider a dynamical systems perspective to model cycles through time while leveraging advancements in generative flow networks (GFlowNets) to learn complex posteriors over explanatory structure given the observational data. In general, reviewers found our work to address an important and long-standing problem in biology while employing a probabilistic/Bayesian perspective ([uDCE], [2wco]), to be well motivated for tackling the problem of cyclic Bayesian structure learning for GRNs ([PVQy, b9zX, uDCE]), and to provide a good quality contribution to the area ([b9zX, uDCE, 2wco]). Reviewers’ primary concerns consist of clarifying questions. Here we address some general concerns raised across reviewers. 

1. **Cyclic dependencies in synthetic and real experiments:** Reviewers asked about cyclic dependencies in the synthetic and real experiments ([PVQy], [b9zX]). We highlight that our experiments do all contain cyclic dependencies in the ground truth graphs and amend the accompanying figures (Fig. 2 and Fig. 3) to appropriately demonstrate this. Please see the attached document for reference. In Fig. 2 we have added a 3D realization showing the difference between the dynamic graph and static graph. We show how we model cyclic dependencies present in the static graph through directed edges in the dynamic graph. This corresponds to the motivation of using dynamic data of the form (x, dx). In Fig. 3 we have adjusted the presented GRN to include the cyclic dependencies of the ground truth GRN that were considered in the experiment.    
2. **Comparison to DAG learning method(s):** Reviewers raised points about demonstrating the effectiveness of DynGFN to learn cyclic dependencies relative to counterpart directed acyclic graph (DAG) learning methods ([PVQy], [2wco]). To show this, we have added an additional experiment (toy example) that considers a DAG-based system and accompanying system with cyclic dependencies. We demonstrate that the DAG structure learning method (NOTEARS [1]) cannot learn the cyclic dependencies in the cyclic system, compared to DynGFN which performs very well in the cyclic setting. We show this result in Fig. 4 (included in the attached document) which we will include in the revised manuscript.  
3. **Figure clarifications:** Some reviewers asked clarifying questions regarding presentation of Fig.1 ([b9zX], [2wco]) and Fig. 3 ([PVQy], [b9zX], [uDCE], [2wco]). To clarify Fig. 1, we will add subscripts to the $Q$ variable to separate usage for modeling the posterior over graphs and posterior over parameters to outline the difference between the 2 $Q$’s. Specifically, for posterior over graphs we will state $Q_\Psi(G | D)$ (this is consistent with the notation in the figure), while for the posterior over parameters we will state $Q_\phi(\theta | G, D)$. Since we model $Q_\phi(\theta | G, D)$ as a Dirac, this collapses to $\theta = h_\phi(G)$, as shown in the figure. For Fig. 3, we have amended the ground truth GRN (stated in item 1 above) and the accompanying heatmaps. Please see the attached document for the amended Fig.3. 

[1] Zheng, Xun, et al. "Dags with no tears: Continuous optimization for structure learning." Advances in neural information processing systems 31 (2018)

Meta-review: Reviewers were excited by the paper's motivation of developing methodology for structure learning in cyclic systems, in this case gene-regulatory networks (GNRs). They found the proposed methodology, based on GFlowNets, to be natural and innovative and were impressed by the method's empirical performance. While they were hesitant about some of the method's assumptions (most notably, the causal sufficiency assumption) as they pertained to GNRs, the reviewers nevertheless felt like ideas in the paper would have significant impact in the GNR literature. All of the reviewers commented on the paper's readability, which detracted from both their understanding and ultimate evaluation. While the reviewers were able to understand enough to advocate for the paper, I encourage the authors to revise the camera-ready for an audience that is less familiar with GFlowNets.
Decision: Accept (poster)