Networks with positive and negative edges (signed networks) are ubiquitous across various domains. They work well for situations when the objects modeled as network vertices have positive and negative interactions. Accounting for edge types helps substantially with many important network-related problems, such as missing link prediction (Li, Fang, and Zhang 2017;Li et al. 2020), node ranking (Li, Fang, and Zhang 2019), network synchronization (Monteiro et al. 2022). Signed networks were successfully applied to model social interactions (trust/distrust in Epinions (Xu et al. 2019), epidemic spreading (Li et al. 2021), political interactions between US Congressmen (Thomas, Pang, and Lee 2006), and gene regulation (Mason et al. 2009).
One of the most interesting things about signed networks is the distribution of signs on the edges -i.e., the signage. Systems modeled by signed networks rarely exhibit random interactions between their elements -instead, they follow certain interaction patterns, that could provide insights into the underlying formation process. Generative models considering signage formation are essential for testing hypotheses about the emergence of interactions, as well as creating synthetic datasets for algorithm benchmarking (especially in areas where obtaining real-world datasets is difficult).
One of the fields affected by the scarcity of realistic signed network datasets is gene regulation. The diversity of cells, organs, and eventually organisms arises from regulatory interactions between genes through their biochemical products (proteins or RNAs). If a gene product increases the synthesis of a target gene product, the corresponding directed edge is modeled as positive; if it decreases, the edge is modeled as negative. The set of all regulatory interactions along with a set of genes is called a gene regulatory network (GRN). GRNs play a key role in the process of organism development, and when broken can lead to serious organism abnormalities and diseases (Alon 2019). Understanding gene regulation mechanisms as well as identifying drug interventions often rely on reconstructing GRNs from the dynamics of the corresponding gene transcripts (Zhang et al. 2017;Lopes-Ramos et al. 2020).
In this work, we designed a new class of signage models for gene regulatory networks along with a maximum-likelihoodbased framework assessing their goodness-of-fit. Our signage models work for scenarios in which the graph topology is formed first and later refined with the edge signs. Motivated by underlying biological processes, our models involve latent non-overlapping groups of nodes and the edge signs are generated based on the endpoints' groups. A real-life example is forming social interactions in a closed community (e.g. dormitories) in which one gets acquainted with the people one lives and interacts with and later decides on their attitude to them. Some evolutionary hypotheses suggest a similar origin of gene regulation, with the genome structure defining the gene interactions in the course of species divergence (Bylino, Ibragimov, and Shidlovskii 2020;Wittkopp and Kalay 2012). Signage models facilitate the comparison of node interactions regardless of the underlying topology -this might be useful for graphs generated with different topologies but fol-lowing the same signage model, allowing topology-agnostic comparisons and testing the plausibility of node interactions hypotheses.
We evaluated the fit of our models using maximum likelihood as an objective criterion, see, for example (Bezáková, Kalai, and Santhanam 2006). Most existing methods evaluating a model's goodness-of-fit rely on high-level network characteristics (degree distribution, diameter, frequency of specific subgraphs, etc.), which might be prone to overfitting and subjectivity in criteria selection especially when real-world data is scarce (as in the case of GRNs). In contrast, likelihood-based goodness-of-fit is a universal approach which does not rely on an arbitrarily chosen network or signage properties. The main challenge of this approach is the computation of the likelihood.
Contributions 1) We design a class of graph signage models with a latent node group assignment, which gets generated first, followed by the edge signage, the distribution of which is conditional on the node group assignment. 2) We pose a new optimization problem for modeling graph signage for gene regulatory networks, GRASMOS 1 , which aims to select the model and its parameters that best fit the observed data using Maximum Likelihood. Our algorithmic framework can be used as a guide to design likelihood estimators for other models, facilitating direct comparison of their goodness-of-fit. 3) We analyze the GRASMOS problem space and identify four cases of parameter combinations that have different intuitions and induce different computational complexity of the likelihood computation. We give efficient algorithms for estimating the likelihood corresponding to the given parameters. Depending on the parameter class they belong to, our algorithms range from closed-form expression to MCMC sampling. 4) We obtain a 16-fold reduction in the MCMC algorithm runtime, by identifying parts of the datasets where MCMC can be replaced by exact computation and by accounting for symmetry in the parameter space. 5) We evaluate our framework on two real-world bacteria GRN datasets -E.coli and B.subtilis (both with thousands of nodes and edges), as well as on synthetic datasets of comparable size. (Derr, Aggarwal, and Tang 2018;Jung, Park, and Kang 2020) designed generators of signed networks based on the structural balance theory, that capture the distribution of signed triangles. Such generators model a signage jointly with a topology, while in our work, the signage is formed on the basis of the existing topology. These generators do not appear to generate self-loops, while our work admits topologies with self-loops, which play an essential part in GRNs (Burda et al. 2011;Alon 2019). Other works designed algorithms generating a node embedding in signed networks and tested its performance for missing link prediction and node classification (Li, Fang, and Zhang 2017;Li et al. 2020). Such algorithms aim to have high predictive power while in our work we aim for a model with high explanatory power.
1 https://github.com/Restel/grasmos Preliminaries A signed graph G ± = (V, E, A) consists of a directed graph G = (V, E) with node set V and edge set E ⊆ V × V , and a signage function A : E → {+, -} determining a positive or a negative sign for each edge. Let n = |V | be the number of vertices and m = |E| = m + + m -, where m + and m - is the number of positive and negative edges, respectively.
For any v ∈ V let out s (v) be the number of outgoing edges from v of sign s ∈ {+, -}, and define in s (v) analogously for incoming edges to v. For a sign s, we define s as its complementary sign. The total number of outgoing edges from v is out(v) = out s (v) + out s(v). Let self s (v) represent the number of s-self-loops of node v.
We propose a graph signage model, where the signs of the edges of a given directed graph G = (V, E) are driven by a random latent node partition C and a parameter matrix ξ. In particular, let S be a set of symbolic node group labels, q be a distribution over S (i.e., q : S → [0, 1] where s∈S q(s) = 1), and ξ be an |S| × |S| matrix of probabilities ξ x,y , where x, y ∈ S. The signage model first randomly creates a latent node partition C : V → S by assigning each node independently to one of the groups in S according to the distribution q. Then, for each edge (u, v) ∈ E, the model assigns the sign + to this edge with probability ξ C(u),C(v) , and the signotherwise, obtaining a signed graph G ± . The parameters of the model are combined in the tuple Θ = (ξ, q ) . We define |S| × |S| probability matrices P + and P -denoting the probabilities that, based on the node assignment of the end-points of an edge, the edge gets sign + or -: For s 1 , s 2 ∈ S, let P + s1,s2 := ξ s1,s2 and P - s1,s2 := 1 -ξ s1,s2 .
In this section, we formalize our optimization problem, GRASMOS, aimed at finding the parameters of the generative graph signage model with the best explanatory power of the observed data, based on maximum likelihood. Recall that the signage model only generates the signs of the edges, not the topology; in a sense, it "augments" the edges with signs based on the existing graph topology. Formally, given a signed graph G ± = (V, A), we are looking for such combination of parameters Θ that has the highest probability of generating edge signs A, denoted P (A|Θ). For a fixed node group assignment, C, we can find the corresponding probability of the signs as:
Since the node group assignment is generated randomly, the probability of interest is a sum of conditional probabilities over all possible node assignments:
(2)
Our model and the above expression are reminiscent of the well-known Ising model from statistical physics and its partition function (Ising 1925;Jerrum and Sinclair 1993). However, there are crucial differences between the two models:
the Ising model is defined for undirected graphs and its Hamiltonian does not capture our setting. GRASMOS aims to find model parameters Θ MLE with the highest probability (or likelihood) of realizing the edge signs A. As is standard for the maximum likelihood approach, for numerical purposes we frame the optimization problem as finding the parameters Θ with the lowest negative loglikelihood:
(3)
GRASMOS is an optimization problem over an (|S| 2 + |S| -1)-dimensional parameter space, since ξ is of dimensions |S| × |S| and q is determined by |S| -1 probabilities.
While our signage model is general, in this work we focus on modeling GRNs where S = {A, R}. Therefore, from now on we assume S = {A, R} and, therefore, the optimization is over a 5-dimensional parameter space. For each parameter setting Θ, the likelihood computation is potentially an exponential summation over 2 n node group assignments. However, on closer scrutiny, it turns out that points in some areas of the problem space are easier to estimate than the rest and it is the relationship between ξ probabilities that determines the computational difficulty of the point estimation. In this section, we subdivide our main optimization problem (3) into several cases (models) according to the relationship between probabilities in ξ and provide intuition for each model. We present the models from the simplest to the most general.
For the problem space points with identical ξ entries, the node group partition does not affect the likelihood of Θ, so we can estimate the likelihood of such points analytically. We call this case the Node-Oblivious model because the signs of edges are independent of the nodes' groups. Define ξ = ξ A,A = ξ A,R = ξ R,A = ξ R,R . Then, P + s1,s2 = ξ and P - s1,s2 = 1 -ξ for every s 1 , s 2 ∈ S. This model can be viewed as a signage analog of the well-studied Erdös-Rényi random graph model.
The edge sign probabilities depend only on the source node's group. In the context of gene regulation, that means that a gene mostly performs the regulation of the same sign: an activator gene tends to activate, and a repressor gene represses. We call the parameter instances belonging to this case the Source-Consistent model. This observation allows us to reformulate (1) as a product of probabilities of signed outgoing edges over the set of vertices:
Target-Consistent model (TC) Similarly to the previous model, in the Target-Consistent model, the following edge probabilities are pairwise equal, i.e. ξ A,A = ξ R,A =: ξ * A and ξ A,R = ξ R,R =: ξ * R and the calculation of P (A|Θ, C) is analogous to the SC model, using incoming edges.
Bi-Node-Consistent model (BNC) In the Bi-Node-Consistent model, the edge signs depend on the node group assignment of both the source and the target nodes and the probabilities in ξ can be arbitrary. In this parameter subspace, we can see ξ instances inducing assortativity: e.g., whether nodes from the same group tend to have positive edges between each other, and negative edges to the nodes from the other group (Mussmann et al. 2015), (Newman 2002). We reformulate (1) for the likelihood estimation of the BNC parameters as a product of the likelihood contributions of non-self-loop (outgoing) signed edges and self-loop signed edges for each node given node assignment C. Let l(v, C) out and l(v, C) self be the likelihood contribution of the outgoing and self-loop edges of node v to the Θ likelihood given C in (1). Thus,
We rewrite (1) with respect to the contribution of each node to the overall likelihood:
In this section, we show how to estimate the likelihood of the GRASMOS model parameters, depending on which of the above-stated models they belong to. The methods we design range in their complexity depending on the computational difficulty of the problem: we provide a closed-form expression for the Node-Oblivious model, a polynomial-time exact algorithm for the source-and Target-Consistent models, and an MCMC sampling algorithm for the Bi-Node-Consistent model. We prove that the models (or, more precisely, their natural generalizations) are self-reducible (Jerrum, Valiant, and Vazirani 1986), which allows us to estimate the likelihood for a given Θ through a product of likelihood ratios of progressively smaller instances. For the BNC model, we use sampling to get an approximation of each of the likelihood ratios.
For the NO model, from ( 1) and ( 2) we get
Therefore, the optimization problem (3) can be solved analytically: it is maximized for ξ max := m+ m++m-. Since q does not contribute to P (A|Θ), Θ MLE = (ξ MLE , q MLE ), where ξ MLE has all entries equal to ξ max and q MLE is arbitrary.
The total likelihood of Θ consists of the sum of 2 n likelihoods of Θ conditional on the node assignment (see ( 2)). We express L(Θ) as the product of ratios of pairs of likelihoods on decreasingly smaller spaces of node group assignments: the reduction in space is achieved by fixing the groups for some nodes.
Let V = {v 1 , v 2 , . . . , v n }. For s 1 , . . . , s j ∈ S, we define C [s1,...,sj ] j as the set of node group assignments C : V → {A, R} such that C(v i ) = s i for every i ≤ j (vertices v 1 , . . . , v j have their group determined by s 1 , . . . , s j ). Notice that C 0 is the set of all node group assignments (with no restrictions) and that |C [s1,...,sj ] j | = 2 n-j for any s 1 , . . . , s j . For our self-reducibility approach, let us fix a "master" node assignment C : V → {A, R} that gradually more and more vertices will adhere to. Let
where P (C j |j) is the probability that vertices v j+1 , . . . , v n get the node assignment given by C j . Notice that P (C j |j) is a probability distribution over the node group assignment subspace C
[ C(v1),..., C(vj )] j
, and, therefore, Z j is the likelihood of Θ restricted to this subspace.
Our overall goal is to estimate L(Θ) = Z 0 . We do this via the following product of likelihood ratios:
Notice that Z n = P (A|Θ, C) can be easily computed via (1). Therefore, if we estimate each of the ratios σ j := Zj Zj-1 , we can compute Z 0 as follows:
For numerical stability, we will choose C so that each ratio σ j is reasonably far from 0.
Theorem 1. There exists C such that σ j ≥ 1/2 for every j ∈ {1, . . . , n}.
Proof. We will define C inductively. First, notice that the ratio σ j does not rely on the entire C but only on its restriction to {v 1 , . . . , v j }. For j ∈ {1, . . . , n}, suppose C(v 1 ), . . . , C(v j-1 ) have been chosen so that σ i ≥ 1/2 for every i < j. We will show how to choose C(v j ) so that σ j ≥ 1/2 by considering the two possible cases: C(v j ) can be either an activator or a repressor. For s ∈ {A, R}, define Z s j using ( 8) with
The proof builds C constructively but the algorithmic efficiency is unclear (it involves summations over exponentially many terms). In the following sections, we will show how to construct C and compute all the σ * j := max{σ A j , σ R j } (and thus L(Θ)) efficiently.
We show how to estimate the ratio of likelihoods σ * j for the SC models in linear time. Our algorithms rely on the fact that, in the SC model, the group assignment of any node, C(v j ), only affects the probabilities of edges connected to v j . Therefore, the calculation of σ * j depends only on node v j and its incident edges. The proof and the analogous computation for the TC model are included in the full version (Brilliantova, Miller, and Bezáková 2022). Theorem 2. When Θ belongs to the Source-Consistent model, for every j ∈ {1, . . . , n} and every C : V → {A, R}:
where
Corollary 3. For the Source-Consistent model, the likelihood of Θ, L(Θ), can be calculated in O(n + m) time.
In contrast to the SC and TC models, in which an edge sign depends only on one endpoint and the ratio of likelihoods σ * j can then be computed analytically, we are not aware of any analytical approach for σ * j in the BNC model where both edge end-points influence the sign of the edge. The main difficulty arises from the fact that the calculation of σ * j under BNC might depend on the colors of the remaining n -j -1 nodes, which, at that point, are still unassigned in the course of the algorithm. Instead of exact computation, we estimate each σ * j via a Markov Chain Monte-Carlo (MCMC) sampling of node group assignments on the subspace corresponding to Z j , with first j vertices assigned to groups.
Suppose C(v 1 ), . . . , C(v j ) has been already defined. Let us define w j (C j ) := P (A|Θ, C j ) • P (C j |j) as the weight of the node group assignment
we will randomly generate C j , with probability proportional to its weight. Therefore, the stationary distribution of the node group assignments should be µ j (C j ) := w j (C j )/Z j . To obtain this stationary distribution, we use the Metropolis-Hasting technique (Metropolis et al. 1953;Hastings 1970), with state space Ω j := C [ C(v1),..., C(vj )] j and the Markov chain transitions τ : Ω j × Ω j → [0, 1] defined as follows. Let C j be the current state.
• Choose a uniformly random z ∈ {j + 1, . . . , n}. Let C j be identical C j , except C j (v z ) = Cj (v z ) (the assignment of v z is opposite in C j and C j ).
• With probability min{1,
(2) (the self-loops τ (C j , C j ) correspond to the remaining probability, so that τ is a stochastic matrix).
The underlying transition graph is analogous to the hypercube with n -j dimensions, and, therefore, the state space is connected (i.e., we can get from every state to every other state using transitions of the Markov chain). This Markov chain is also aperiodic due to the presence of selfloop transitions and, therefore, it has a unique stationary distribution. The Metropolis-Hastings technique ensures that this stationary distribution is exactly the distribution µ j , i.e., proportional to the weights w j (C j ).
Following the earlier outline (see ( 9) and ( 10)), we will be estimating σ j = Zj Zj-1 . Recall that, for numerical stability, we wanted σ j ≥ 1/2. We showed that C(v j ) = arg max s∈{A,R} σ s j yields σ j = max{σ A j , σ R j } ≥ 1/2. We will estimate both σ A j and σ R j simultaneously by generating samples from
Let C j-1 be a random sample from Ω j-1 , drawn according to µ j-1 . Define f A : Ω j-1 → {0, 1} as the indicator function that the corresponding node assignment assigned v j as an activator:
Therefore, we can use the expectation of f A to estimate σ A j (and, similarly, we can define f R and estimate σ R j since
The accuracy of the estimate increases with the number of samples: we will draw k independent samples C
(1)
). We use the same set of samples for the average f A and the average f R computation (each sample contributes to either f A or f R ). Since the expectation E µj-1 [f A (C j-1 )] corresponds to the probability of drawing a sample with v j assigned as an activator, we draw such sample with probability qZ A j Zj-1 , and we draw a sample with v j as a repressor with probability
(1-q)Z R j Zj-1 . Since these two types of samples cover the state space, at least one of the types will be drawn with probability ≥ 1/2. This means that we will very likely have a sufficient number of samples of at least one of the types, which will allow us to estimate the corresponding expected value closely. Since we do not know which of the types is more likely, we generate k samples and estimate both the average f A and the average f R . Then, we choose the larger average and define the corresponding A or R as C(v j ). Then, divide the larger average by q or (1 -q) as appropriate, obtaining an estimate on σ j . As a side remark, this σ j might potentially be different than max{σ A j , σ R j } but we still have σ j ≥ 1/2 since qσ A j ≥ 1/2 or (1 -q)σ R j ≥ 1/2. Next, we discuss the accuracy of this estimate. Suppose we aim to be within a (1 + ε)-factor of the true L(Θ), for some small ε. Since we have n self-reducibility steps (i.e., n ratios σ j to estimate, see ( 9) and ( 10)) and the estimated quantities are larger than 1/2, we can follow the derivation in (Jerrum 2003)[page 26] almost verbatim to obtain the desired number of samples per estimate. However, this number of samples is rather high for the datasets of our size, so we (successfully) employ convergence heuristics for the estimates.
We combine the MCMC approach with exact calculation to achieve significant speed-up for real-world datasets. GRNs, like many other real-world networks, are not dense: they usually have a few highly connected nodes (called hubs in network science) and many low-degree nodes. For such lowdegree nodes, the use of MCMC is unnecessary and computationally wasteful: as soon as all of their adjacent vertices are assigned, the ratio of likelihoods for such nodes can be calculated analytically. For such v j , this theorem shows how to compute σ * j exactly in time proportional to v j 's degree. The proof is in the full version (Brilliantova, Miller, and Bezáková 2022). Theorem 4. Let C : V → {A, R, undefined} be a partial node assignment. Reorder the vertices so that v 1 , . . . , v j-1 are assigned by C to either A or R, and all other vertices are undefined. Suppose that all of v j 's neighbors are already assigned by C. Then σ * j can be calculated as:
Synthetic datasets To validate that our algorithms can correctly identify the GRASMOS parameters, we created synthetic datasets with known Θ, resembling real-world GRNs by their network parameters and size. The characteristics of all synthetic and real-world datasets we used are shown in Table 1. Our evaluation pipeline had two steps: generation and reconstruction. During the generation part, we generated network topologies from the Directed Scale-Free (DSF) model, known to capture the characteristics of GRNs well (Van den Bulcke et al. 2006), then using several different Θ parameters we sampled the node group partition and assigned + andsigns to the edges.
For illustrative purposes, in the Results section below we present our validation for the Source-Consistent model: We chose several Θ (we refer to them as Θ-generators) from this model and generated corresponding signages for our graph. During the reconstruction part, we tested multiple candidate Θ (Θ-candidates) belonging to the SC, TC, and NO models and checked the proportion of samples in which the candidate with the best likelihood Θ MLE matched the Θ-generator. We also calculated the L 1 -norm between the Θ-generator and Θ MLE . The parameters within the Θ-candidates varied between 0.1 and 0.9, with an increment step of 0.1. We tested 4 different Θ-generators in the generation part and for each of 1458 of Θ-candidates in the reconstruction part. To account for stochasticity we repeated the pipeline 10 times. Tests for additional Θ-generators can be found in the full version (Brilliantova, Miller, and Bezáková 2022).
Real-world GRNs We used two public databases: Regu-lonDB (Santos-Zavaleta et al. 2019) and SubtiWiki (Pedreira, Elfmann, and Stülke 2022;Flórez et al. 2009). Both databases contain experimentally validated information about gene regulatory interactions and their type (activation/repression) for bacteria species: Regulon for Escherichia coli and SubtiWiki for Bacillus We only left the entries associated with transcriptional gene regulation. From both datasets, we filtered out regulation edges of unknown types and duplicated edges. Regulatory interactions can be context-dependent: under different conditions, the same gene can either activate or repress the target gene (Ong and Corces 2011). Our model does not account for such scenarios, so we had to filter out duplicated edges.
For fitting the GRASMOS parameters of real-world GRNs, we: used a fine-grained exploration of parameters varying in [0.1, 0.9] with an increment of 0.05 belonging to the SC, TC, and NO models; and, for complexity reasons, we evaluated the BNC model subspace using a coarse-grid with each parameter varying in {0.25, 0.5, 0.75}. For Θs belonging to multiple models, we compared the results from these approaches, made sure that they were consistent, and identified a good candidate to run a "refined" BNC search with 432 additional Θ-candidates in the vicinity (searching through Θs roughly within ±0.125 from the identified candidate).
We estimated L(Θ) of the BNC parameters using a parallel implementation on the RIT research computing cluster. We used up to 432 nodes with each core estimating a likelihood of a single Θ-candidate via MCMC sampling. Each core is equipped with Intel ® Xeon ® Gold 6150 CPU @ 2.70GHz. The RAM upper limit for our computation was 2048 MB.
For synthetic datasets, we tested how well our framework based on the total likelihood reconstructs the Θ values that were used to generate an edge signage instance. We found that, overall, it was able to reconstruct the Θ-generator well.
Figure 1(left) shows the percentage of test cases when the Θgenerator equaled the reconstructed Θ MLE exactly, and when the ξ-portion of Θ and Θ MLE were identical, respectively. Reconstruction of the ξ-parameters was particularly successful. Moreover, even when Θ and Θ MLE differed, they were very close in their L 1 -norm, see Figure 1(right).
Table 2 presents 5 best Θ candidates per model for the Sub-tiWiki and Regulon datasets. For both datasets, the BNC parameters had the best performance among all models. However, some of the top BNC parameters belonged to the SC model for SubtiWiki. BNC's winning score is not surprising as it generalizes the other models. However, when another model overtakes BNC (as in the case of the SC model and Subtiwiki), that might lead to biological hypotheses about the underlying processes, which is our main reason for considering the simpler models. The source-consistency hypothesis corresponds well to the existence of operons -groups/clusters of bacterial genes that are co-located on the DNA strand and controlled by the same gene-regulator (Salgado et al. 2001). Interestingly, in the case of Regulon, none of the top-5 Θ belonged to the SC model. Moreover, all of the top parameters implied that genes tend to activate the genes from groups other than they belong to. This observation might corroborate the existence of feed-forward transcriptional control with interchanging edge signs suggested for certain GRNs (Sasse and Gerber 2015) and raises the question of whether the bacteria species corresponding to these two datasets indeed have different gene regulation mechanisms, or the discrepancies should be explained by other factors.
We assessed the accuracy of the MCMC sampling for the likelihood estimation of those Θ instances for which we have the exact solution (i.e., Θ in the SC, TC, and NO models). We aimed to have a (1 + 1/n)-multiplicative accuracy, i.e., at most 0.2% likelihood error for our datasets. We used this target accuracy for our computation of the number of needed samples. In most cases, our MCMC computations were within the target accuracy. In those cases when we were off by more than (1 + 1/n), the corresponding Θ-candidates had L(Θ) far from L(Θ MLE ), and in the coarse search, the Θ-candidate with the smallest L(Θ) was close to L(Θ MLE ). This meant that despite having more inaccurate L(Θ) estimates for these (few) Θs than we hoped for, the coarse search eliminated them from consideration and thus eliminated the inaccuracies.
As for the running time, the analytical calculation of some vertices (Section Algorithmic speed-up of MCMC) and the reduction of the parameter space by half (discussed in the full version of this paper) resulted in 16x actual speed-up for BNC Θ parameters compared to the naive approach. Due to the parallelization of the grid search, the total runtime of fitting the BNC parameters was around 2 days for Regulon and 4 days for the SubtiWiki.
In this work, we present a novel Maximum-Likelihood-based graph signage model selection problem for gene regulation, GRASMOS, developed a fitting framework for the problem and showcased its usage for two gene regulatory networks of B.subtilis and E.coli. Our graph signage models and the model selection framework opened up a plethora of directions for future research. Among them: 1) Does the high explanatory power of the Θ M LE translate into high predictive power? For our models to have high predictive power, we need to be able to reconstruct the hidden node assignment. We used the master node assignment C (visually shown for SubtiWiki dataset in Figure 2) to estimate the likelihood, but how is it related to the "optimal" node assignment C M LE = arg max C∈C L(Θ|C)? 2) Is our framework robust when there are missing edges? What is the percentage of the missing edges it can handle?
Additionally, we acknowledge that the signage can depend on the topology. A natural next step for this work is to compare the likelihood of models that generate the signage jointly with the topology, and an approach that uses existing topology generators followed by our signage model.