In this paper, we conduct theoretical analyses on inferring the structure of gene regulatory networks. Depending on the experimental method and data type, the inference problem is classified into 20 different scenarios. For each scenario, we discuss the problem that with enough data, under what assumptions, what can be inferred about the structure. For scenarios that have been covered in the literature, we provide a brief review. For scenarios that have not been covered in literature, if the structure can be inferred, we propose new mathematical inference methods. Otherwise, we prove that the structure cannot be inferred.

Inference; Gene regulatory network; Independence; Differential equation.

Frequently used abbreviations:

GRN: gene regulatory network.

DAG: directed acyclic graph.

ODE: ordinary differential equation

§ AUTHOR SUMMARY

Many molecular regulators modulate gene expression, and these regulators form a gene regulatory network. In this paper, we study how to infer the gene regulatory network (or other networks) structure with experimental data. Specifically, we determine whether the expression of one gene can regulate (activate or inhibit) the expression of another gene. The network inference problem can be classified into different scenarios by the experimental data type. For each scenario (data type), we discuss whether the network structure can be inferred, and what assumptions are required. We introduce corresponding mathematical inference methods for each scenario, where some methods are known, and some are novel. This paper provides a unified framework to discuss various situations of the gene regulatory network structure inference problem.

§ INTRODUCTION

In living cells, most genes are transcribed to mRNAs and then translated to proteins. Some proteins serve as transcription factors that regulate the transcription of other gene(s). For example, in fruit flies, CLOCK (CLK) and CYCLE (CYC) are transcription factors that activate the transcription of period (per) and timeless (tim), the protein products of which in return inhibit CLK/CYC. Such transcription-translation negative feedback loop plays a central role in circadian rhythm [1]. Small molecules also participate in gene regulation. Retinoic acid inhibits Fgf8 in mice [2]. For E. coli, LacI protein inhibits the expression of lac operon, unless lactose level is high and glucose level is low [3]. The regulation of gene expression is a central topic in biology.

Broadly speaking, a gene regulatory network (GRN) consists of various molecular regulators that govern gene expression, namely, levels of mRNAs and corresponding proteins [4]. Such regulators interact with each other and form a complicated network. See Fig. <ref> for an example of GRN in E. coli. Here the expression of one gene can activate or inhibit the expression of another gene. In this paper, the major focus is to determine the regulatory relationship between genes. Regulations involving non-genetic factors, or even other regulation networks (not necessarily biological), can be determined with similar methods.

An example of GRN in E. coli [5]. Each vertex is a gene (some genes are omitted). A solid arrow means activation, and a dashed arrow means inhibition.

Most subjects of a GRN are large molecules confined within living cells. Therefore, it is difficult or even impossible to directly verify a regulatory network among a group of genes with biochemical methods. On the other hand, a large amount of data that are related to GRNs have been collected and documented, such as single-cell observational gene expression data [6] and bulk level interventional phenotype data [7]. The central question of this paper is to infer the GRN structure using experimental data (also called “reverse-engineering”). Besides gene regulation, the problem of inferring network structure with experimental data is extensively studied in other fields [8, 9, 10]. We consider two major situations:

(1) Consider genes $V_1,\ldots,V_n$. We can measure the expression levels of these genes (mRNA or protein count). The goal is to reconstruct the GRN that contains a specific gene, such as $V_1$. This means that we need to exclude irrelevant genes that do not belong to this GRN. Besides, for any two genes $V_i,V_j$ that belong to this GRN, we need to study whether the expression of $V_i$ regulates the expression of $V_j$. If yes, we determine whether the regulation is activation or inhibition.

(2) Consider genes $V_1,\ldots,V_n$ and a phenotype $V_0$ (growth rate etc.). Some genes can affect the level of this phenotype, and we assume this phenotype does not affect the expression of genes. We can measure the level of this phenotype, but not the expression levels of genes. Nevertheless, we can intervene with any genes, so as to change (decrease in general) their expression levels. By adding intervention on each gene $V_i$ and observe the phenotype $V_0$, we can exclude irrelevant genes. Therefore, we can stipulate that each $V_i$ can directly or indirectly affect $V_0$. Similarly, the goal is to reconstruct the GRN that contains $V_0$. We study whether the expression of $V_i$ activates or inhibits the expression of $V_j$.

Various data types can be used to infer GRN structures. Different data types require different mathematical inference tools. For certain data types, there have been numerous studies in determining GRN structures [5, 11, 12, 13, 14, 15, 16, 17]. We propose a unified framework to treat the GRN inference problem. Depending on possible data types that can be used to infer the GRN structure, we classify the GRN inference problem into 20 scenarios. For each scenario, we discuss what can be inferred about the GRN structure, and what assumptions are required. For scenarios that have been extensively studied, we only introduce a few representatives to show what we can infer in these scenarios. Readers may refer to some review papers for more details about these scenarios [11, 12, 18]. Other scenarios will be thoroughly discussed: whether and how the GRN can be inferred. We also develop novel mathematical methods if necessary. For identified regulatory relations, we determine whether they are activation or inhibition. Through this analysis, the existing mathematical methods, our novel methods, and potential methods that will be invented in the future can be treated in the same unified paradigm. This paradigm can be generalized to fit new data types.

The expression of one gene could promote its own expression, which is called autocatalysis [19]. We invent a mathematical method to identify the existence of autocatalysis in certain scenarios.

This paper is not a review of all important existing GRN inference methods. For example, machine learning methods and Bayesian type methods are not covered. Besides, this paper focuses on theoretical results, not practical algorithms.

Section <ref> briefly introduces possible biological data types regarding GRNs and related biological assumptions. Section <ref> discusses mathematical setup and related mathematical assumptions. Section <ref> presents the main results of this paper: with what data, under what assumptions, what we can infer about the GRN structures. Section <ref> discusses biological background about GRN and related measurements, and concerns about related assumptions. Section <ref> contains some classical mathematical results related to GRN structure inference. In Section <ref>, we develop some novel mathematical results that are used in the inference of GRN structures. Section <ref> provides a detailed explanation for the results in Section <ref>. Section <ref> discusses how to determine the existence of autocatalysis. We finish with some discussions in Section <ref>.

This paper is written to draw attention from readers with various backgrounds: experimental biologists who have related data and want to know what GRN structures can be inferred; mathematical biologists who need to implement inference methods on existing data; applied mathematicians who want to have a panoramic understanding of the GRN inference problem and develop novel mathematical methods. We encourage readers to skip certain sections if they are not interested in corresponding topics or do not have certain backgrounds. Readers who only want to understand the main results could stop at Section <ref>. Readers who are not interested in mathematical details could skip Subsections <ref>-<ref>, <ref>, and Sections <ref>-<ref>. Readers who are not interested in biological discussions could skip Section <ref>. Fig. <ref> describes the dependence relations of all sections.

\S2\ar[d]\ar[dr]&\S1\ar[l]\ar[r]\ar[r]&\S3\ar[dll]\ar[d]\ar[ld]\\
\S5\ar[rd]&\S4\ar[ld]\ar[rd]&\S6,\S7\ar[d]\\
\S9&\S10&\S8\ar[l]

Dependence relations of all sections in this paper.

§ BIOLOGICAL SETUP

§.§ Experimental data types

Consider genes $V_1,\ldots,V_n$. Various types of experimental data can be used to infer GRN structures. They can be classified in four major dimensions: (1) Modality. We can directly measure the expression levels of $V_1,\ldots,V_n$, or a phenotype $V_0$ that depends on the expression levels of $V_1,\ldots,V_n$. (2) Granularity. We can measure the levels of $V_1,\ldots,V_n$ or $V_0$ for a single cell and repeat many times, so as to obtain a group of random variables $X_1,\ldots,X_n$ or $X_0$ that represent the random values of $V_1,\ldots,V_n$ or $V_0$. We also obtain the corresponding probability distribution of these variables. Alternatively, we can measure these quantities over a large population of cells (bulk level), so that the randomness is averaged out. Then we obtain deterministic results of $X_1,\ldots,X_n$ or $X_0$, denoted by $x_1,\ldots,x_n$ or $x_0$. (3) Intervention. We can measure certain quantities after interfering with the expression of some genes, so that these genes and downstream genes are affected (not fully muted in general). Then we can compare measured quantities before intervention ($x_1,\ldots,x_n$) and quantities after intervention ($x_1',\ldots,x_n'$). We can also observe without any intervention. (4) Temporal resolution. We can measure at a single time point, $x_i(0)$, or measure at multiple time points as a time series, $x_i(0),x_i(1),x_i(2),$ and so on. (⁎) When we measure at single-cell level at multiple time points, there is an extra dimension: Can we measure the same cell multiple times, so that we can obtain the joint distribution for multiple time points, $\mathbb{P}[X_1(0)=c_0,X_1(1)=c_1,X_1(2)=c_2]$, or we need to measure different cells at different time points, so that we just have multiple marginal distributions for each time point, $\mathbb{P}[X_1(0)=c_0],\mathbb{P}[X_1(1)=c_1],\mathbb{P}[X_1(2)=c_2]$?

For these four major dimensions, we have $2^4=16$ different scenarios. Each scenario is named with a combination of four labels: (1) Gene expression or Phenotype; (2) Single-cell or Bulk; (3) Non-interventional or Interventional; (4) One-time or Time series. In four scenarios (Single-cell + Time series), there is an extra dimension of Joint distribution or Marginal distribution, meaning a total of 20 scenarios.

For each dimension, there is one type that is strictly more informative than the other: Gene expression $>$ Phenotype, Single-cell $>$ Bulk, Interventional $>$ Non-interventional, Time series $>$ One-time, Joint distribution $>$ Marginal distribution. Nevertheless, for more informative data types, generally the experiments are more difficult, more expensive, and less accurate.

There are intermediate states for each dimension. (1) When we can only measure the expression level of one gene (or a few genes), we can regard the measured gene as a phenotype, and apply methods developed for phenotype data. (2) When we increase the number of cells measured in one experiment, due to the central limit theorem [20], the data variance decreases fast. When the cell number is small enough, so that the stochasticity is significant, we can treat this case the same as the single-cell case. (3) When we can only add interventions on limited genes, some inference methods for interventional data might fail, and we need to treat such data as non-interventional. (4) For time series data, some methods only require measuring at two time points, while some methods require measuring at many more time points, depending on the number of genes involved. Therefore, when we have data at a few time points, we might need to use methods for one-time data.

§.§ Biological assumptions

Gene regulation in reality is complicated. To infer the GRN structure with limited experimental data, we need some biological assumptions.

(A) Besides genes, other non-genetic factors also regulate gene expression (retinoic acid for Fgf8). We assume that such non-genetic factors are kept at constants, so that they do not interfere with determining the regulation between genes.

(B) The regulation between genes might require specific conditions (lactose and glucose for Lacl and lac operon). We assume that all such conditions are satisfied, so that we can observe all possible regulations among considered genes.

(C) We assume all genes in the GRN that we want to reconstruct have been considered, and assume that the set of considered genes is not too large. Therefore, the set of genes considered consists of all genes in a small GRN and at most a few irrelevant genes.

(D) If $V_i$ regulates $V_j$, we assume the increase of $V_i$ always leads to the increase (activation) or decrease (inhibition) of $V_j$. We exclude the case that the regulation switches between activation and inhibition in different situations.

With these assumptions, the GRN structure inference problem can be transformed into a well defined mathematical problem. Nevertheless, these assumptions are not realistic in biology, and they will be further discussed in Subsection <ref>.

§ MATHEMATICAL BACKGROUND AND SETUP

§.§ Introduction of related mathematical assumptions

In the following, we introduce some necessary mathematics for understanding the problem of inferring GRN structures. We also introduce four mathematical assumptions that will be used in certain scenarios. Under these assumptions, the underlying GRN is simple enough, or the experimental data are regular enough, so that they follow certain mathematical models. In Section <ref>, we will introduce under different assumptions, what GRN structures can be inferred from experimental data. Nevertheless, all assumptions are more or less unrealistic in biology. We will discuss concerns about these assumptions in Subsection <ref>.

Four assumption are used in certain scenarios: (1) Path Blocking (PB): the intervention on one gene has no effect on another gene (or a phenotype), if and only if other intervened genes block all paths; (2) Directed Acyclic Graph (DAG): the GRN can be described by a directed graph without cycles; (3) Markov and Faithful (MF): the distribution of gene expression properly reflects the underlying DAG through conditional independence relations; (4) Linear System (LS): the gene expression (and possibly phenotype) time series data satisfy a linear ODE system. Most contents of this section, except the path blocking property, are from standard textbooks [21, 22, 23]. Readers can remember the names of these assumptions and skip the rest of this section, if they are not interested in mathematical details.

§.§ Graph theory and path blocking property

In graph theory, a directed graph consists of several vertices $\{V_1,\ldots,V_n\}$ and some directed edges between two vertices, such as $V_i\to V_j$.

A GRN can be represented as a directed graph. Each vertex $V_i$ represents a gene, and each directed edge $V_i\to V_j$ means that the expression of gene $V_i$ can regulate the expression of gene $V_j$. We require that the graph is connected.

If there are edges $V_i\to V_j,V_j\to V_k,\ldots,V_l\to V_m$, then $V_i\to V_j\to V_k\cdots\to V_l\to V_m$ is called a path from $V_i$ to $V_l$. Although the GRN might have cycles, we do not allow a path to pass a vertex more than once. The number of edges in a path is called the length of this path. If there exist paths from $V_i$ to $V_j$, the minimal length of such paths is called the distance from $V_i$ to $V_j$.

If there is an edge $V_i\to V_j$, then $V_i$ is a parent of $V_j$, and $V_j$ is a child of $V_i$. If there is a path from $V_i$ to $V_j$,  then $V_i$ is an ancestor of $V_j$, and $V_j$ is a descendant of $V_i$.

If we can only measure a phenotype $V_0$, then the concerned GRN consists of $V_0$ and genes $V_1,\ldots,V_n$ that have paths leading to $V_0$. We assume there is no edge starting from $V_0$.

Consider any two vertices $V_i,V_j$ and a set $\mathcal{S}\subset\{V_1,\ldots,V_n\}\backslash\{V_i,V_j\}$. If all paths from $V_i$ to $V_j$ need to pass through $\mathcal{S}$, we say that $\mathcal{S}$ blocks $V_i$ to $V_j$. Specifically, if there is no path from $V_i$ to $V_j$, then any set, including the empty set $\emptyset$, blocks $V_i$ to $V_j$.

Assume we can interfere with any genes and measure the expression levels of all genes. This means we can interfere with different sets of genes, and compare if they have the same effect on certain genes. The path blocking property means that for any $V_i,V_j$ and any set $\mathcal{S}\subset\{V_1,\ldots,V_n\}\backslash\{V_i,V_j\}$, after interfering with all genes in $\mathcal{S}$, intervention on $V_i$ cannot provide extra influence on the expression of $V_j$, if and only if $\mathcal{S}$ blocks $V_i$ to $V_j$. An equivalent expression is that the intervention on $\mathcal{S}$ and the intervention on $\mathcal{S}\cup\{V_i\}$ have the same effect on $V_j$. The path blocking property in this case means a set of genes that blocks $V_i$ to $V_j$ in the directed graph (GRN) also blocks the influence of $V_i$ to $V_j$ in gene expression, and vice versa. Therefore, by comparing the gene expression levels under different interventions, we can obtain partial information about the graph structure. For example, assume the GRN is $V_1\to V_2\to V_3$, and the path blocking property holds. After interfering with $V_2$, adding intervention on $V_1$ cannot further affect $V_3$, meaning that $\{V_2\}$ blocks $V_1$ to $V_3$. After interfering with $V_1$, adding intervention on $V_2$ can further affect $V_3$, meaning that $\{V_1\}$ does not block $V_2$ to $V_3$.

Assume we can interfere with any genes and measure one phenotype $V_0$. This means we can interfere with different sets of genes, and compare if they have the same effect on the measured phenotype. The path blocking property means that for any $V_i$ and any set $\mathcal{S}\subset\{V_1,\ldots,V_n\}\backslash\{V_i\}$, after interfering with all genes in $\mathcal{S}$, intervention on $V_i$ cannot provide extra influence on the expression of $V_0$, if and only if $\mathcal{S}$ blocks $V_i$ to $V_0$. Similarly, by comparing the phenotype under different interventions, we can obtain partial information about the graph structure.

§.§ Directed acyclic graph, Markov property, and faithful property

A directed acyclic graph (DAG) is a directed graph that has no directed cycles. For example, if we have edges $V_1\to V_2$ and $V_2\to V_3$, then we cannot have edge $V_3\to V_1$, but edge $V_1\to V_3$ is allowed. Specifically, an edge from a vertex $V_i$ to itself ($V_i\to V_i$) is regarded as a cycle.

In a causal DAG (also called Bayesian network) $\mathcal{G}$ with $n$ vertices, each vertex $V_i$ has an associated random variable $X_i$, which represents the stochastic expression level of gene $V_i$ in a single cell. The joint distribution of these variables is denoted as $\mathbb{P}$.

$\mathbb{P}$ is said to be Markov to $\mathcal{G}$, if $\mathbb{P}(X_1,\ldots,X_n)=\prod^{n}_{i=1}\mathbb{P}(X_i\mid \pi_{X_i})$, where $\pi_{X_i}$ means the expression levels of $V_i$'s parents. If all the conditional independence relations in $\mathbb{P}$ also appear in any other $\mathbb{P}'$ that is Markov to $\mathcal{G}$, then we say that $\mathbb{P}$ is faithful to $\mathcal{G}$. If a distribution $\mathbb{P}$ is Markov and faithful to $\mathcal{G}$, then the structure of $\mathcal{G}$ properly reflects the causal relations in $\mathbb{P}$.

Consider a DAG $\mathcal{G}$. A common model that fits this DAG is $X_i=f_i(\pi_{X_i})+\epsilon_i$. This means the value of $V_i$ depends on the values of its parents, plus an independent noise. In general, the joint distribution generated by this model is Markov and faithful to $\mathcal{G}$.

§.§ Linear system

If we measure the bulk level gene expression time series, the expression level of a gene $V_i$ can be regarded as a variable $x_i$ that changes continuously. If we measure the single-cell level gene expression $X_i$, then its expectation is $x_i$. The expression levels $\{x_i\}$ of all genes along time satisfy a linear system, or have the linearity property, if the following linear ordinary differential equations (ODEs) hold.

\begin{equation}
\label{eq1}
\begin{split}
\mathrm{d}x_1/\mathrm{d}t &=a_{11}x_1+a_{12}x_2+\cdots+a_{1n}x_n+b_1,\\
\mathrm{d}x_n/\mathrm{d}t &=a_{n1}x_1+a_{n2}x_2+\cdots+a_{nn}x_n+b_n.
\end{split}
\end{equation}

This system can be rewritten as $\mathrm{d}\vec{x}/\mathrm{d}t=A\vec{x}+\vec{b}$, where $\vec{x}=(x_1,\ldots,x_n)'$, $A=[a_{ij}]$, $\vec{b}=(b_1,\ldots,b_n)'$. Here $b_i$ describes the base synthesis rate of gene $V_i$, $a_{ii}$ describes the degradation (and possibly autocatalysis) of $V_i$, and $a_{ij}$ describes the effect of $V_j$ on the expression of $V_i$. If $a_{ij}>0$, $V_j$ activates $V_i$; if $a_{ij}<0$, $V_j$ inhibits $V_i$; if $a_{ij}=0$, $V_j$ does not affect $V_i$.

§ INFERRING GRN STRUCTURE WITH DIFFERENT DATA TYPES

We have discussed different data types and different assumptions in Section <ref> and Section <ref>. In this section, we present the main results of this paper: with what data type, under what assumptions, what can we infer about the GRN structure. See Table <ref> for a summary of all the results. For example, in Scenario 6 (Gene, Bulk, One-time, Interventional), we can fully infer the GRN structure under the path blocking assumption, and we can partially infer the GRN structure under the directed acyclic graph assumption. These results are explained in Section <ref>, based on the mathematical methods in Section <ref> and Section <ref>.

2|l|2*                                                                                2l|One-Time                                                                                2l|Time Series
2|l|                                                                                c]@l@Non-
Intervention                                    Intervention                                                                                c]@l@Non-
Intervention                                                                       Intervention
Expression     c]@l@Single-
Cell     c]@l@Scenario 1

c]@l@Scenario 2

PB: full
DAG: partial
c]@l@Scenario 3 a/b

3a Joint:
UC: full
3b Marginal:
c]@l@Scenario 4 a/b

4a Joint:
UC: full
4b Marginal:
LS: full
PB: full
DAG: partial
MF+DAG: full

Bulk                                                       c]@l@Scenario 5

No                                     c]@l@Scenario 6

PB: full
DAG: partial                        c]@l@Scenario 7

No                                                                        c]@l@Scenario 8

LS: full
PB: full
DAG: partial
2*Phenotype                                                     c]@l@Single-
Cell     c]@l@Scenario 9

No                                     c]@l@Scenario 10

PB: partial                                   c]@l@Scenario 11 a/b

No                                                                       c]@l@Scenario 12 a/b

PB: partial
Bulk                                                       c]@l@Scenario 13

No                                    c]@l@Scenario 14

PB: partial                                   c]@l@Scenario 15

No                                                                       c]@l@Scenario 16

PB: partial

GRN structure inference with different data types: under what assumptions, what structures can be inferred. There are 16 scenarios classified by the following dimensions of data types: Gene Expression vs. Phenotype; Single-Cell vs. Bulk; One-Time vs. Time Series; Non-Interventional vs. Interventional. In Scenarios 3/4/11/12, there is an extra dimension of Joint vs. Marginal. There are different assumptions: PB: path blocking; DAG: directed acyclic graph; MF: Markov and faithful; LS: linear system; UC: unconditional. Full/partial/no means all/some/no GRN structures can be inferred. For example, “MF+DAG: partial” means under MF assumption and DAG assumption, GRN structure can be partially inferred. The asterisk * in Scenarios 12/16 means for some identified regulatory relations, we cannot determine whether they are activation or inhibition.

The basic method of establishing causal relation is to vary one variable and examine if the other variable changes. When we cannot interfere with gene expression, the only chance is to profile gene expression on single-cell level. Within a single cell, the stochastic fluctuation of gene expression plays a similar role as intervention. Nevertheless, by comparing Scenario 1 vs. Scenario 2 and comparing Scenario 3 vs. Scenario 4, we can see that the stochastic fluctuation is less informative than intervention.

With interventional phenotype data (Scenarios 10, 12, 14, 16), the GRN is analogous to a black box. Although the gene expression levels are unknown, we can interfere by manipulating certain genes and use the resulting phenotype as output. Under different assumptions, we can partially reconstruct the GRN structure.

We can evaluate the performance of each scenario based on the level of GRN structures that can be inferred. Fully recovered GRN means 2 points, and partially recovered GRN means 1 point. For example, Scenario 1 (partial GRN) has 1 point. The summation over eight scenarios (1-8) with gene expression data is 11 points, and the summation over eight scenarios with phenotype data (9-16) is 4 points. We compare the overall score for each data dimension. Gene expression vs. Phenotype: 11 vs. 4; Single-cell vs. Bulk: 9 vs. 6; Time series vs. One-time: 8 vs. 7; Interventional vs. Non-interventional: 12 vs. 3. We can see that the gene expression data and interventional data are more informative.

Scenario 4 (Gene, Single-cell, Time series, Interventional) is the most informative case. We can see that in Scenario 4, the GRN structure can be fully inferred under various assumptions.

§ BIOLOGICAL DETAILS AND CONCERNS

§.§ Measurements related to GRN

We briefly introduce biological experiments related to inferring GRN structures and their restrictions, as an extension for Subsection <ref>. We do not aim at covering all important related papers. Readers may refer to some review papers for detailed summaries on related biological techniques [24, 13, 25].

The expression of a protein-coding gene consists of transcription (DNA to mRNA) and translation (mRNA to protein). Therefore, to measure the expression level of a given gene, we can either measure the amount of the corresponding mRNA or measure the amount of the corresponding protein. There are many methods with various reliabilities and accuracies to measure the amount of mRNA and/or protein on the bulk level and the single-cell level [26, 27, 28]. By now (2021), these methods are not 100% accurate, especially on single-cell level [29]. A more important problem is that most mRNAs and proteins are confined within living cells (except for secreted proteins). This means cells have to be killed before these molecules can be harvested, and then some analytical methods can be applied to quantify the abundance of these molecules. Therefore, with these methods, a cell or a cell population can be measured only once [26].

Many mRNAs have less than 20 copies in one cell [30]. Thus the gene expression of a single cell is too stochastic to be described by a deterministic model. We can repeat the gene expression measurement experiment over multiple single cells and obtain a probability distribution of single-cell gene expression [31]. When the observation is based on a large number of cells (bulk level), stochasticity is averaged out, so that the dynamics should be deterministic. Theoretically, repeating the gene expression measurement experiment over a large number of cells should produce the same result. In reality, various technical issues weaken the repeatability of measuring gene expression [32]. For example, external factors that affect gene expression, such as nutrition concentrations, are hardly the same for two cell populations. Therefore, when we sample two populations from the same cell line (especially if we sample at different time points), their gene expression profiles might not be exactly the same.

Traditional approaches in the studies of gene expression are observational, meaning that we directly measure cell(s) at stationary without interfering with gene expression. Recent techniques like gene knockdown or gene knockout allow us to temporarily or permanently alter the expression of genes of interest [33, 34, 35]. Most interventions decrease the expression of corresponding genes, but increasing is also possible [36]. After the intervention, we can measure the gene expression data at a certain time point. By now (2021), on bulk level, commonly used knockdown interventions cannot robustly decrease the expression of targeted genes to zero [37]. For example, Hurley et al. applied siRNA to disrupt certain genes. For only around 70% genes, the bulk-level expression decreases to less than 40% of the original level. For around 10% genes, siRNAs even increase their bulk expression levels [14]. On the single-cell level, we can measure and select cells that are successfully intervened. However, the measurement is not 100% accurate, so that the selection is not always correct. Besides, selection might be biased towards cells with low baseline expression levels of targeted genes. Interventions cannot be used to maintain gene expression to fixed non-zero levels.

After the intervention, we hope to measure how the gene expression evolves over time. Under current technologies, we can only measure a cell or a cell population once. Therefore, at the bulk level, we can measure different populations at different time points after the same intervention. The gene expression level at the bulk level is a deterministic number. Thus theoretically, measuring different populations at different time points is equivalent to measuring the same population at different time points. At the single-cell level, we can repeat the measurement for different cells at the same time point to obtain the marginal distribution of gene expression at each time point. However, since we can only measure a cell once, we cannot obtain the joint distribution of gene expression at different time points. The joint distribution can provide extra information (correlation coefficient, etc.) than marginal distributions. Regardless of technology restrictions, we discuss both cases, depending on whether the joint distribution can be measured.

Some phenotypes, such as growth rate, number of molecules released, and drug resistance, could be used to reflect gene expression levels [7]. These phenotypes can be measured on bulk level or single-cell level. Besides, some phenotypes can be measured without perturbing or even killing the cells. This means such phenotypes can be measured at different time points for the same cell(s). Nevertheless, phenotypic transitions involve complicated phenomena and mechanisms [38]. Therefore, we need to be cautious when utilizing phenotypes in determining GRN structures.

§.§ Concerns about biological GRN assumptions

In Subsection <ref>, we propose four biological assumptions to simplify the GRN structure inference problem. Here we discuss biological concerns about these assumptions.

Various non-genetic factors, such as vitamins and minerals, directly regulate gene expression. Since this paper focuses on regulations between genes, we assume such non-genetic factors are kept at constants. Nevertheless, we do not fully understand what factors can affect gene expression, and it is difficult to control the abundance of so many kinds of molecules. Besides, when we manipulate the expression of genes, the concentrations of some non-genetic factors might also be changed [39].

The regulation between genes might occur only under specific conditions. We assume that such conditions are satisfied, so that all regulations can be observed. Since we need to infer whether a regulation between genes exists, we cannot guarantee that the conditions for this regulation are satisfied. We can also remove this assumption, and just aim at reconstructing the GRN that operates under naive conditions.

To correctly infer the GRN structure, we need to assume that all genes in the target GRN are considered. For example, if the true GRN is $V_1\to V_2\to V_3$, but $V_2$ cannot be measured or intervened with, then we can observe that varying $V_1$ in any case would lead to the change of $V_3$. Therefore, we would obtain a false direct edge $V_1\to V_3$. We also need to assume that the set of considered genes is not too large. The amount of data needed and the computational time cost for GRN structure inference methods grow fast with the number of genes considered. If we cannot exclude most genes that are irrelevant to the target GRN beforehand, then either the inference methods behave poorly, or the cost (time and money) is unbearable. Nevertheless, before inferring the GRN structure, it is almost impossible to accurately distinguish which genes are involved or not involved in the target GRN.

To simplify the classification of activation and inhibition, we assume the regulation is monotonic, meaning that if $V_i$ regulates $V_j$, then the increase of $V_i$ always leads to the increase or decrease of $V_j$. In reality, Crn1 at different concentrations can either activate or inhibit the Arp2/3 complex in yeast [40]. This illustrates that gene regulation might be non-monotonic.

§.§ Concerns about mathematical GRN assumptions

In Section <ref>, we propose four mathematical assumptions that will be used in certain scenarios. Here we discuss biological concerns about these assumptions.

The Path Blocking (PB) assumption states that the intervention on one gene has no effect on another gene (or a phenotype), if and only if other intervened genes block all paths. Consider a GRN $V_1\to V_2\to V_3$. Some interventions (e.g., siRNA) are not strong enough, so that after interfering with $V_2$, $V_1$ can still affect $V_2$, then $V_3$. On the other hand, if we assume that $V_1$ is necessary for the expression of $V_2$, then as we knock out $V_1$, $V_2$ cannot express. In this case, interventions on $V_2$ have no extra effect on $V_3$, although the path $V_2\to V_3$ is not blocked. Even though the path blocking assumption is powerful in the GRN structure inference, it is not very realistic.

The directed acyclic graph (DAG) assumption states that the GRN has no cycle, which is crucial in many scenarios. Nevertheless, in reality, it is common for a GRN to have directed cycles. In fact, there are some pairs of genes that could mutually regulate each other [5].

The Markov and faithful (MF) assumption states that the conditional independence relations in the distribution of gene expression are consistent with the causal DAG. Under this assumption, conditional independence should be correctly identified. However, due to the heterogeneity of cells [41], there might be multiple cell types, each of which has its own gene expression distribution, and the measured gene expression distribution might be the combination of multiple distributions. Such mixing might make independent genes show pseudo dependence. Mathematically, assume $X_1,X_2$ (expression of two genes for cell type 1) are independent, and $X_1',X_2'$ (expression of two genes for cell type 2) are independent. Define $\bar{X}_1=X_1$ with probability $0.5$, and $\bar{X}_1=X_1'$ with probability $0.5$, and $\bar{X}_2$ is defined similarly (overall expression of two genes for mixed cell population). Then $\bar{X}_1$ and $\bar{X}_2$ might not be independent.

The linear system (LS) assumption states that the net change rate of one gene (corresponding mRNA or protein) linearly depends on the quantity of another gene, regardless of other genes. In reality, the regulation of gene expression might be highly nonlinear. Specifically, if gene $V_1$ inhibits gene $V_2$, then when the level of $V_2$ is $0$, the net change rate of $V_2$ might still be negative, a contradiction.

§ EXISTING RELATED MATHEMATICAL RESULTS

In this section, each subsection contains a known mathematical method that will be used to infer GRN structures in certain scenarios. At the end of each subsection, there is an “A/I Remark” that explains how to determine whether an identified regulatory relation is activation or inhibition.

§.§ Use conditional independence to infer DAG structure

This subsection is a standard topic in causal inference. Notice that the causal DAG is also called the Bayesian network. Readers may refer to related monographs for background and details [22].

For a distribution $\mathbb{P}$ that is Markov and faithful to a causal DAG, there is an edge between $V_i$ and $V_j$ if and only if $V_i$ and $V_j$ are independent conditioned on some other variables. However, we still need to determine the direction of this edge.

For two causal DAGs $\mathcal{G}_1$, $\mathcal{G}_2$ with a bijection between vertices, if any distribution $\mathbb{P}$ that is Markov and faithful to $\mathcal{G}_1$ is also Markov and faithful to $\mathcal{G}_2$, and vice versa, then we say that $\mathcal{G}_1$ and $\mathcal{G}_2$ are Markov equivalent. With this relation, all causal DAGs can be classified into different Markov equivalent classes. For example, the following three DAGs form an equivalent class: $V_1\to V_2\to V_3$, $V_1\leftarrow V_2 \leftarrow V_3$, $V_1\leftarrow V_2 \to V_3$. If $\mathbb{P}$ is Markov and faithful to an unknown causal DAG $\mathcal{G}$, then we can only find the Markov equivalent class that contains $\mathcal{G}$. DAGs in this equivalent class are indistinguishable under $\mathbb{P}$. Markov equivalent causal DAGs have the same edges, except that some edges have opposite directions. Therefore, with the information of $\mathbb{P}$, we can determine all edges of the unknown $\mathcal{G}$, but some edge directions are unknown.

Notice that we cannot directly use correlation coefficient to determine the existence of an edge. Consider a DAG $V_2\leftarrow V_1 \to V_3$. Although there is no edge between $V_2$ and $V_3$, the correlation coefficient of $X_2,X_3$ is not $0$ in general. This means “correlation does not imply causation”.

Consider a GRN that $V_i$ strongly activates $V_k$, $V_k$ strongly activates $V_j$, and $V_i$ weakly inhibits $V_j$. If $V_i$ increases, the weak inhibition of $V_i\to V_j$ is covered by the strong activation of $V_i\to V_k\to V_j$, meaning that $V_j$ also increases. If we naively calculate the correlation coefficient of $V_i$ and $V_j$, we would obtain the false conclusion that $V_i$ activates $V_j$. Since there are multiple paths from $V_i$ to $V_j$, we need to control other variables to determine the direct regulation from $V_i$ to $V_j$. Therefore, we propose to calculate the conditional covariance (or conditional correlation, partial correlation, etc.) between $V_i$ and $V_j$, conditioned on other possible parents of $V_j$. Positive conditional covariance means activation, and negative conditional covariance means inhibition. If there is only one path from $V_i$ to $V_j$, namely the edge $V_i\to V_j$, then we can directly check the effect of $V_i$ on $V_j$, such as calculating the correlation coefficient. Using these methods, we can also quantify the strength of a regulatory relation.

§.§ Inferring GRN structure with ODE model

For the linear ODE system Eq. <ref>, if we know the value of $x_i(t)$ at different time points, we can calculate $\mathrm{d}x_i/\mathrm{d}t$, and then calculate the parameters $\{a_{ij}\}$ and $\{b_i\}$ by solving a linear algebraic equation system [42]. Here we need to assume that the interventions added on different genes only change the initial values $x_i(0)$, but not system parameters $\{a_{ij}\}$ and $\{b_i\}$. If the genes in a GRN satisfy such a linear system, then $a_{ij}\ne 0$ if and only if there is an edge $V_j\to V_i$.

In the ODE model, if $a_{ij}>0$, $V_j$ activates $V_i$; if $a_{ij}<0$, $V_j$ inhibits $V_i$. The value of $a_{ij}$ quantifies the regulation strength.

§.§ Determining GRN in stochastic process model

On single-cell level, gene expression is stochastic. If we can only measure each cell once, and repeat for many different cells, we obtain some random variables, where each variable corresponds to the expression level of a gene. If we can measure the same cell at different time points, then we obtain the joint distribution of a stochastic process. We consider a discrete-time stochastic process that describes the evolution of multiple variables. This model is also called the dynamic Bayesian network. For example, we consider three variables at different time points, $[X_1(0),X_2(0),X_3(0)]$, $[X_1(1),X_2(1),X_3(1)]$, $[X_1(2),X_2(2),X_3(2)]\ldots$. Here $X_1(t),X_2(t),X_3(t)$ are the expression levels of genes $V_1,V_2,V_3$ at time $t$. We assume there are no other genes that affect the expression of $V_1,V_2,V_3$. One interpretation is to define $[X_1(t+1),X_2(t+1),X_3(t+1)]$ as a stochastic function of $[X_1(t),X_2(t),X_3(t)]$:

\begin{equation}
\label{eq2}
\begin{split}
\end{split}
\end{equation}

Here different $\epsilon_{i,t}$ are independent noise terms. This interpretation implies that the process is Markov, meaning that that if we condition on the variables at one time point ($[X_1(2),X_2(2),X_3(2)]$), then the past ($[X_1(1),X_2(1),X_3(1)]$) and the future ($[X_1(3),X_2(3),X_3(3)]$) are independent. For a continuous-time Markov process, if we can only observe at discrete time points, then the observed variables form a discrete-time Markov process. Thus we can focus on discrete-time processes.

Assume the underlying GRN is $V_1\to V_2\to V_3$, then this means current $X_1$ can directly affect future $X_1$, current $X_1,X_2$ can directly affect future $X_2$, and current $X_2,X_3$ can directly affect future $X_3$. The equations of $[X_1(t+1),X_2(t+1),X_3(t+1)]$ becomes

\[X_1(t+1)=f[X_1(t)]+\epsilon_{1,t},\]

\[X_2(t+1)=g[X_1(t),X_2(t)]+\epsilon_{2,t},\]

\[X_3(t+1)=h[X_2(t),X_3(t)]+\epsilon_{3,t}.\]

The causal relations between different variables in this process are shown in Fig. <ref>. An arrow can only point from a variable at time $t$ to a variable at time $t+1$. Therefore, this illustration is a DAG. We can project this process DAG along the time axis to obtain the underlying GRN. In Fig. <ref>, the underlying GRN has a cycle. We can still use the process DAG to illustrate the causal relations.

Stochastic process model. $X_1(0)$ can affect $X_1(1)$ and $X_2(1)$; $X_2(0)$ can affect $X_2(1)$ and $X_3(1)$; $X_3(0)$ can affect $X_3(1)$. The right most $V_1\to V_2\to V_3$ is the corresponding GRN.

Stochastic process model. $X_1(0)$ can affect $X_1(1)$ and $X_2(1)$; $X_2(0)$ can affect $X_2(1)$ and $X_3(1)$; $X_3(0)$ can affect $X_3(1)$ and $X_1(1)$. The right most $V_1\to V_2\to V_3$ with an arrow $V_3\to V_1$ is the corresponding GRN.

If we have the joint distribution of two neighboring time points, then we can use conditional independence to infer the structure of the process DAG, and thus the underlying GRN. Here we need to introduce the notion of Markov boundary. For a set of random variables $\mathcal{S}$ and another random variable $X$, a Markov boundary $\mathcal{S}_1$ of $X$ in $\mathcal{S}$ is a subset of $\mathcal{S}$, such that (1) $X\ind \mathcal{S}\mid \mathcal{S}_1$; (2) for any $\mathcal{S}_2\subsetneqq \mathcal{S}_1$, $X\nind \mathcal{S}\mid \mathcal{S}_2$. Here $X_1\ind X_2\mid X_3$ or $X_1\nind X_2\mid X_3$ means $X_1$ and $X_2$ are (not) independent conditioned on $X_3$. A Markov boundary is a minimal subset that fully preserves the information of $V$ in $\mathcal{S}$. It is possible to have multiple Markov boundaries in reality [43], and such multiplicity causes severe problems in identifying causal relations [44]. For simplicity, we shall assume the Markov boundary is unique. In Fig. <ref>, if we set $\mathcal{S}=\{X_1(0),X_2(0),X_3(0)\}$, then the Markov boundary of $X=X_3(1)$ in $\mathcal{S}$ is $\{X_2(0),X_3(0)\}$. This means that we should have arrows from the Markov boundary to $X$, namely, $X_2(0)\to X_3(1)$ and $X_3(0)\to X_3(1)$. Similarly, the Markov boundary of $X_2(1)$ is $\{X_1(0),X_2(0)\}$, meaning arrows $X_1(0)\to X_2(1)$ and $X_2(0)\to X_2(1)$. The Markov boundary of $X_1(1)$ is $\{X_1(0)\}$, meaning an arrow $X_1(0)\to X_1(1)$. When we project this process DAG along the time axis, we obtain the underlying GRN $V_1\to V_2 \to V_3$. Since $X_i(t+1)$ naturally depends on $X_i(t)$, we do not project this edge as $V_i\to V_i$.

We have shown that for a causal DAG, we can use conditional independence to discover all the edges, except that the directions of some edges are unknown. For the stochastic process model, each arrow in the process DAG has a predetermined direction (from an earlier time to a later time). Therefore, with the help of the time arrow, we can fully determine the structure of the process DAG, and thus the underlying GRN.

In order to correctly recover the GRN, we need to require that the observation time step is small enough, so that the conditional independence correctly reflects the process (DAG). For example, in Fig. <ref>, if the observation time step is 2, so that $[X_1(1),X_2(1),X_3(1)]$, $[X_1(3),X_2(3),X_3(3)],\ldots$ are invisible, then the structure of the process DAG is shown in Fig. <ref>. Here an indirect relation from $X_1$ to $X_3$, $X_1(0)\to X_2(1) \to X_3(2)$, is misinterpreted as a direct relation $X_1(0)\to X_3(2)$, so that the projected GRN has an extra arrow $V_1\to V_3$. Similarly, if the process in Fig. <ref> has observation time step 2, then as shown in Fig. <ref>, the inferred GRN has several extra arrows. In this case, we can determine whether one gene affects another gene, but we cannot determine whether this influence is direct or indirect.

We denote $\mathbb{E}X_1(t),\mathbb{E}X_2(t),\mathbb{E}X_3(t)$ by $x_1(t),x_2(t),x_3(t)$. In Eq. <ref>, if the functions $f,g,h$ are all linear, and $\mathbb{E}\epsilon_{i,t}=0$, then we have

\begin{equation*}
\begin{split}
\end{split}
\end{equation*}

This is a discretized version of Eq. <ref>. Therefore, when the dynamics is linear, we can regard the ODE model as the expectation of the stochastic process model. When the dynamics is nonlinear, we might not have $\mathbb{E}f(X_1)=f(\mathbb{E}X_1)$, and the evolution of $x_1(t+1)$ is essentially more complicated, such that it does not only depend on $x_1(t),x_2(t),x_3(t)$, and the ODE model fails.

For GRN $V_1\to V_2\to V_3$, from Fig. <ref>, we can see that conditioned on $[X_2(0),X_2(1),\ldots]$, $[X_1(0),X_1(1),\ldots]$ and $[X_3(0),X_3(1),\ldots]$ are independent. However, for any $t$, even if the process has reached stationary, we might not have $X_1(t)\ind X_3(t)\mid X_2(t)$. Therefore, in the stochastic process model, if we only have the stationary distribution, we cannot use conditional independence to infer the GRN structure. The causal DAG model discussed in Subsection <ref> might not be regarded as the stationary situation of the stochastic process model.

Stochastic process model. It is the same as Fig. <ref>, except that the observation time step is 2, so that $X_1(1),X_2(1),X_3(1),X_1(3),X_2(3),X_3(3),\ldots$ are invisible. Therefore, there is a direct arrow $X_1(0)\to X_3(2)$. The right most $V_1\to V_2\to V_3$ with an arrow $V_1\to V_3$ is the corresponding (fake) GRN.

Stochastic process model. It is the same as Fig. <ref>, except that the observation time step is 2, so that $X_1(1),X_2(1),X_3(1),X_1(3),X_2(3),X_3(3),\ldots$ are invisible. Therefore, there are direct arrows $X_1(0)\to X_3(2)$, $X_2(0)\to X_1(2)$, and $X_3(0)\to X_2(2)$. The right most $V_1\leftrightarrow V_2\leftrightarrow V_3$ with an arrow $V_1\leftrightarrow V_3$ is the corresponding (fake) GRN.

In the stochastic process model, similar to the discussion in Subsection <ref>, we can use conditional covariance to determine whether an arrow represents activation or inhibition.

§ NOVEL RELATED MATHEMATICAL RESULTS

In this section, each subsection contains a novel mathematical method that will be used to infer GRN structures in certain scenarios. At the end of each subsection, there is an “A/I Remark” that explains how to determine whether an identified regulatory relation is activation or inhibition.

§.§ Use path blocking relation to infer GRN structure with gene expression data

Consider an unknown GRN with genes $V_1,\ldots,V_n$. We can intervene with any genes and measure the expression profile of all genes. Assume the path blocking property holds. Then for each $V_i,V_j$ and each subset $\mathcal{S}$ of $\{V_1,\ldots,V_n\}\backslash\{V_i,V_j\}$, we know whether $\mathcal{S}$ blocks $V_i$ to $V_j$, meaning that any directed path from $V_i$ to $V_j$ passes at least one vertex in $\mathcal{S}$. In this case, we can infer the whole GRN: there is an edge $V_i\to V_j$ if and only if other vertices cannot block $V_i$ to $V_j$.

For an identified edge $V_i\to V_j$, add interventions on certain genes, so that other paths from $V_i$ to $V_j$ are blocked. Then $V_i$ activates $V_j$ if the decrease of $V_i$ leads to the decrease of $V_j$, and $V_i$ inhibits $V_j$ if the decrease of $V_i$ leads to the increase of $V_j$.

§.§ Use path blocking relation to infer GRN structure with phenotype data

Consider an unknown GRN with genes $V_1,\ldots,V_n$ and a phenotype $V_0$. We can intervene with any genes, but we can only measure the level of the phenotype, not the genes. We assume there is a directed path from each $V_i$ to $V_0$, and there is no edge starting from $V_0$. Assume the path blocking property holds. Then for each $V_i$ and each subset $\mathcal{S}$ of $\{V_1,\ldots,V_n\}\backslash\{V_i\}$, we know whether $\mathcal{S}$ blocks $V_i$ to $V_0$, meaning that any directed path from $V_i$ to $V_0$ passes at least one vertex in $\mathcal{S}$. In this case, the GRN structure can be partially inferred (at least $n$ edges).

If a subset $\mathcal{S}$ of $\{V_1,\ldots,V_n\}\backslash\{V_i\}$ blocks $V_i$ to $V_0$, but any proper subset of $\mathcal{S}$ cannot block $V_i$ to $V_0$, then $\mathcal{S}$ is called a minimal blocking set. If a blocking set $\mathcal{S}$ is not minimal, then $\mathcal{S}$ contains a blocking subset that is minimal. Define $\beta(V_i)$ to be all minimal subsets that block $V_i$ to $V_0$. If two GRNs have the same path blocking relations, then they have the same $\beta(V_i)$ for each $V_i$, and vice versa.

For the two GRNs in Fig. <ref>, $\beta(V_1)=\{\{V_3\}\}$, $\beta(V_2)=\{\{V_3\}\}$, $\beta(V_1)=\emptyset$. Thus they are equivalent in the sense of path blocking. If a directed edge appears in all equivalent GRNs, we can determine that this edge exists; if a directed edge appears in none of equivalent GRNs, we can determine that this edge does not exist; if a directed edge appears in some but not all equivalent GRNs, we cannot determine whether this edge exists.

Two equivalent GRNs in the sense of path blocking.

The following procedure describes how to determine certain edges. (1) There is an edge $V_i\to V_0$ if and only if $\beta(V_i)=\emptyset$. (2) If there exists $\mathcal{S}\in\beta(V_i)$, so that $V_j\notin \mathcal{S}$, and $\mathcal{S}$ cannot block $V_j$ to $V_0$, then there is no edge $V_i\to V_j$. (3) If (2) is not satisfied, but there exists $\mathcal{S}\in\beta(V_i)$, so that $V_j\in \mathcal{S}$, then there is an edge $V_i\to V_j$. (4) If for any $\mathcal{S}\in\beta(V_i)$, we have $V_j\notin \mathcal{S}$, and $\mathcal{S}$ blocks $V_j$ to $V_0$, then we cannot determine whether $V_i\to V_j$ exists.

(1) If there is an edge $V_i\to V_0$, then other vertices cannot block this path. If there is no edge $V_i\to V_0$, then $\{V_1,\ldots,V_n\}\backslash\{V_i\}$ can block $V_i$ to $V_0$, and it contains a minimal blocking subset.

(2) Choose a path $V_j\to\cdots\to V_0$ that cannot be blocked by $\mathcal{S}$. If there exists an edge $V_i\to V_j$, then the path $V_i\to V_j\to\cdots\to V_0$ cannot be blocked by $\mathcal{S}$, a contradiction.

(3) Define $\mathcal{S}'=\mathcal{S}\backslash\{V_j\}$. Assume there is no edge $V_i\to V_j$. Since $\mathcal{S}'$ cannot block $V_i$ to $V_0$, there is a path $V_i\to \cdots \to V_k\to V_j$ that is not blocked by $\mathcal{S}'$, and there is a path $V_j\to \cdots \to V_0$ that is not blocked by $\mathcal{S}'$.

$V_k$ cannot appear in the path $V_j\to \cdots \to V_0$. Otherwise, we can combine these two paths and shrink out $V_j$, to obtain a path $V_i\to \cdots\to V_k\to \cdots V_0$. This path from $V_i$ to $V_0$ does not pass through $V_j$, and it cannot be blocked by $\mathcal{S}'$, which means it cannot be blocked by $\mathcal{S}$, a contradiction.

Find all possible $V_l$ with the property that there exists a path $V_i\to \cdots \to V_l\to V_j$ that is not blocked by $\mathcal{S}'$, then $V_l$ also cannot appear in the path $V_j\to \cdots \to V_0$. Define the set of all such $V_l$ (including $V_k$) to be $\mathcal{T}$. $\mathcal{T}\cup\mathcal{S}'$ blocks any path from $V_i$ to $V_0$, and it cannot block the path $V_j\to \cdots \to V_0$. We can find a minimal blocking subset in $\mathcal{T}\cup\mathcal{S}'$, which belongs to $\beta(V_i)$, and it cannot block $V_j$ to $V_0$. Thus the condition of (2) is satisfied, a contradiction.

(4) If the condition of (4) is satisfied, we show that when the edge $V_i\to V_j$ does not exist, adding this edge into the GRN does not change path blocking relations, and when the edge $V_i\to V_j$ exists, deleting this edge from the GRN does not change path blocking relations. Therefore, we cannot determine whether the edge $V_i\to V_j$ exists.

Assume the edge $V_i\to V_j$ is not in the GRN. Consider $\mathcal{S}$ that blocks $V_i$ to $V_0$, and $\mathcal{R}$ that cannot block $V_i$ to $V_0$. Since $\mathcal{S}$ contains a minimal subset, and this subset blocks $V_j$ to $V_0$, after adding the edge $V_i\to V_j$, it can still block $V_i$ to $V_0$. After adding the edge $V_i\to V_j$ in the GRN, $\mathcal{R}$ still cannot block $V_i$ to $V_0$. For $V_k$ that has a path to $V_i$, we can use the same argument to show that blocking subsets of $V_k$ to $V_0$ are not changed.

Assume the edge $V_i\to V_j$ is in the GRN. We prove that after deleting this edge, the condition of (4) still holds. Then using the above argument, we can see that after deleting this edge, adding it back does not change the path blocking relations. After deleting this edge, if there is a minimal subset $\mathcal{T}$ that contains $V_j$ and blocks $V_i$ to $V_0$, then before deleting the edge $V_i\to V_j$, $\mathcal{T}$ is still a minimal subset that blocks $V_i$ to $V_0$, contradicting to the condition of (4). After deleting this edge, assume there is a minimal subset $\mathcal{R}$ that blocks $V_i$ to $V_0$ but does not contain $V_j$ and does not block $V_j$ to $V_0$. Then before deleting the edge $V_i\to V_j$, $\mathcal{R}\cup\{V_j\}$ blocks $V_i$ to $V_0$, and it contains a minimal blocking subset $\mathcal{R}'$ that contains $V_j$. The reason for $V_j\in\mathcal{R}'$ is that $\mathcal{R}$ cannot block $V_j$ to $V_0$ after deleting the edge $V_i\to V_j$. Such $\mathcal{R}'$ violates the condition of (4). In sum, after deleting the edge $V_i\to V_j$, $\beta(V_i)$ is not changed. Similarly, for another $V_k$, $\beta(V_k)$ is not changed.

Consider the shortest path from $V_i$ to $V_0$, $V_i\to V_k \to \cdots\to V_0$. Define $\mathcal{S}$ to be the set of all children vertices of $V_i$, then $\mathcal{S}$ blocks $V_i$ to $V_0$. Since vertices in $\mathcal{S}\backslash \{V_k\}$ cannot be closer to $V_0$ than $V_k$, $\mathcal{S}\backslash \{V_k\}$ cannot block the path $V_i\to V_k \to \cdots\to V_0$. Therefore, a minimal subset of $\mathcal{S}$ that blocks $V_i$ to $V_0$ must contain $V_k$. This means that the edge $V_i\to V_k$ can be identified. Therefore, we can identify at least one edge starting from each $V_i$. For an unknown GRN, we can identify at least $n$ edges. Consider the GRN $V_n\to\cdots\to V_1\to V_0$. We can identify exactly $n$ edges in this GRN. Thus $n$ is the minimal number of edges that could be identified.

Consider a GRN that there is an edge from each vertex in $\mathcal{S}_2=\{V_{n/2+1},\ldots,V_{n}\}$ to each vertex in $\mathcal{S}_1=\{V_1,\ldots,V_{n/2}\}$, and there is an edge from each vertex in $\mathcal{S}_1$ to $V_0$. All $n^2/4+n/2$ edges in this GRN can be identified. We guess $n^2/4+n/2$ is the maximal number of edges that could be identified.

Assume we have identified an edge $V_i\to V_j$. From the above argument, all edges in the shortest path from $V_j$ to $V_0$ can be identified. We add interventions on all genes except $V_i$ and genes on this shortest path. If the intervention on $V_i$ and the intervention on $V_j$ have the same sign on the change of $V_0$, then $V_i$ activates $V_j$. If they have different signs, then $V_i$ inhibits $V_j$. Here we assume that the effect of intervention is known to be decreasing gene expression.

§.§ Use ancestor-descendant relation to infer GRN structure with gene expression data

Consider an unknown GRN with genes $V_1,\ldots,V_n$, and assume the GRN is a DAG. We can add interventions on any genes and measure their expression levels. If the intervention on $V_i$ leads to the change of $V_j$, then $V_i$ is an ancestor of $V_j$, and there is a directed path from $V_i$ to $V_j$ on the graph. With such ancestor-descendant relations, we can partially infer the GRN structure: at least $n-1$ edges can be identified.

For the two GRNs in Fig. <ref>, $V_1$ has descendants $V_2,V_3,V_4$; $V_2$ has descendants $V_3,V_4$; $V_3$ has descendant $V_4$; $V_4$ has no descendant. Thus they are equivalent in the sense of ancestor-descendant relations. If a directed edge appears in all equivalent GRNs, we can determine that this edge exists; if a directed edge appears in none of equivalent GRNs, we can determine that this edge does not exist; if a directed edge appears in some but not all equivalent GRNs, we cannot determine whether this edge exists.

Two equivalent GRNs in the sense of ancestor-descendant relations.

The following procedure describes how to determine certain edges. (1) If $V_j$ is not a descendant of $V_i$, then we can determine that the edge $V_i\to V_j$ does not exist. (2) If $V_j$ is a descendant of $V_i$, and $V_i$ has another descendant $V_k$, which is an ancestor of $V_j$, then we cannot determine the existence of the edge $V_i\to V_j$. (3) If $V_j$ is a descendant of $V_i$, and $V_i$ does not have another descendant $V_k$, which is an ancestor of $V_j$, then we can determine that the edge $V_i\to V_j$ exists.

(1) If there is an edge $V_i\to V_j$, then $V_j$ is a descendant of $V_i$.

(2) Since $V_i$ is an ancestor of $V_k$, we can add an edge $V_i\to V_k$ if it does not exist, and the ancestor-descendant relations are not affected. Since $V_k$ is an ancestor of $V_j$, we can add an edge $V_k\to V_j$ if it does not exist, and the ancestor-descendant relations are not affected. Since we have a path $V_i\to V_k \to V_j$, if there is an edge $V_i\to V_j$, we can delete it without affecting the ancestor-descendant relations. Since $V_i$ is an ancestor of $V_k$, if the edge $V_i\to V_k$ does not exist, we can add it, and the ancestor-descendant relations are not affected.

(3) If there is no edge $V_i\to V_j$, then there is a path $V_i\to V_k\to \cdots\to V_j$. Thus $V_k$ is a descendant of $V_i$, and $V_j$ is a descendant of $V_k$, a contradiction.

When we can confirm the existence of edge $V_i\to V_j$, there is no other path from $V_i$ to $V_j$. Therefore, when we interfere with $V_i$, the change of $V_j$ is directly from $V_i$. Assume we can measure gene expression profiles. After the intervention on $V_i$, if the change of $V_i$ and the change of $V_j$ have the same sign, $V_i$ activates $V_j$; otherwise, $V_i$ inhibits $V_j$.

Consider $\mathcal{S}_1=\{V_1,\ldots,V_{n/2}\}$, $\mathcal{S}_2=\{V_{n/2+1},\ldots,V_n\}$. There is an edge from each vertex in $\mathcal{S}_1$ to each vertex in $\mathcal{S}_2$, and no other edge exists. Then all these $n^2/4$ edges can be identified. We guess this is the maximal number of edges that could be identified. When the GRN is a DAG, we can see that identified edges cannot form a triangle: if $V_i\to V_j$ and $V_j\to V_k$ can be identified, then $V_i\to V_k$ cannot be identified. Due to Turán's Theorem [21], a triangle-free graph with $n$ vertices can have at most $n^2/4$ edges. Therefore, in the DAG case, $n^2/4$ is the maximal number of edges that could be identified.

If the GRN has cycles, we might recover no edge. For example, consider a GRN $V_1\to V_2\to V_3\to V_4 \to V_5$ with an edge $V_5\to V_1$, and a GRN $V_1\to V_3\to V_5\to V_2 \to V_4$ with an edge $V_4\to V_1$. In both GRNs, any vertex is a descendant of any vertex. Therefore, these two GRNs have the same ancestor-descendant relations, but share no edge.

If the GRN is a DAG with $n$ vertices, then we can use ancestor-descendant relations to identify at least $n-1$ edges.

The DAG is associated with a natural partial order: if there is an edge $V_i\to V_j$, then $V_i<V_j$. This partial order can be extended to a total order. If $V_i<V_j$ in the total order, then $V_i$ cannot be a descendant of $V_j$.

When $n=2$, this proposition is trivial. Assume this proposition is true for any $n<K$. When $n=K$, find a vertex $V_i$ which has no child. Assume that after deleting $V_i$ and edges linked to $V_i$, the DAG is divided into $m$ connected components $\mathcal{S}_1,\ldots,\mathcal{S}_m$. For the parents of $V_i$ in each $\mathcal{S}_j$, choose the largest one $V_j'$ under the total order. There is no other path from $V_j'$ to $V_i$. Thus the edge $V_j'\to V_i$ can be identified. With $m$ connected components, we can identify $m$ edges that lead to $V_i$. An edge that can be identified in each connected component should be able to be identified in the original DAG. Thus deleting $V_i$ does not affect identifying edges in each connected component. Apply this proposition to each connected component with size $n_j$, then we can identify at least $n_j-1$ edges. In sum, we can identify at least $\sum_j n_j-m+m=K-1$ edges.

Assume we have identified an edge $V_i\to V_j$. This means no $V_k$ is a descendant of $V_i$ and an ancestor of $V_j$. Thus $V_i\to V_j$ is the only directed path from $V_i$ to $V_j$. When we intervene with $V_i$, so that the expression level of $V_i$ decreases, we can measure the expression level of $V_j$. If the expression level of $V_j$ decreases, then $V_i$ activates $V_j$. Otherwise, $V_i$ inhibits $V_j$.

§.§ Use conditional independence and ancestor-descendant relation to infer GRN structure

Assume the GRN is a DAG. We can measure the joint distribution of gene expression levels, and this distribution is Markov and faithful to the DAG. We can also intervene with genes to obtain the ancestor-descendant relations. With the conditional independence relations, we can find all edges of this DAG, but some edge directions cannot be determined. With the ancestor-descendant relations, we can determine the direction of each edge. Therefore, combining these two data types, we can reconstruct the full DAG.

When the full DAG is inferred, we can use the same methods in Subsection <ref> to determine whether an edge $V_i\to V_j$ is activation or inhibition.

§.§ Use partially observed ODE to infer GRN structure with phenotype data

Consider an unknown GRN with phenotype $V_0$ and genes $V_1,\ldots,V_n$. Assume the GRN is a DAG. We can intervene with any genes, but we can only measure the phenotype on bulk level, not the genes. We assume there is a directed path from each $V_i$ to $V_0$, and there is no edge starting from $V_0$. Denote the levels of $V_0,V_1,\ldots,V_n$ by $x_0(t),x_1(t),\ldots,x_n(t)$, and assume they satisfy a linear ODE system. There is a partial order on $V_0,V_1,\ldots,V_n$, associated with edges in this DAG. We can extend it into a total order and reorder the genes, so that in the linear ODE system $\mathrm{d}\vec{x}/\mathrm{d}t=A\vec{x}+\vec{b}$, $A$ is upper-triangular. This means there is no edge from $V_i$ to $V_j$ if $i<j$. The diagonal elements of $A$ are the degradation rates, which are negative. Therefore, the linear system has a unique fixed point $-A^{-1}\vec{b}$, which is stable.

Assume the system starts from the fixed point, and each time, we intervene with the value of one gene $V_i$. This means that $\vec{x}(0)=-A^{-1}\vec{b}+\vec{\delta}_i$, where $\vec{\delta}_i$ is a zero vector except its $i$th component. After the intervention, we can observe how $x_0(t)$ changes. Here we need to assume that the interventions added on different genes only change the initial values $x_i(0)$, but not system parameters $\{a_{ij}\}$ and $\{b_i\}$. We consider a question: based on only $x_0(t)$ under different interventions, can we partially infer the DAG? Surprisingly, the answer is yes. We can determine whether a gene $V_i$ is an ancestor of another gene $V_j$.

Assume the diagonal elements of $A$ are different. Then the solution of $x_0$ is

\[x_0(t)=c_{0}e^{a_{00}t}+c_{1}e^{a_{11}t}+\cdots+c_{n}e^{a_{nn}t}+d,\]

where the coefficients $c_i$ and $d$ depend on $A,\vec{b},\vec{x}(0)$. Assume in $V_0,V_1,\ldots,V_n$, the descendants of $V_i$ are $V_0,V_1,\ldots,V_{i-1}$. Since we start from the fixed point, except that $V_i$ is perturbed, those $V_j$ with $j>i$ shall be fixed, and the solution of $x_0$ is

\[x_0(t)=c_{0}e^{a_{00}t}+c_{1}e^{a_{11}t}+\cdots+c_{i}e^{a_{ii}t}+d.\]

The coefficients $c_0,c_1,\ldots,c_i$ are all non-zero, unless $A$ satisfies an algebraic equation. However, such $A$ is a zero-measured set in $\mathbb{R}^{(n+1)\times (n+1)}$ [45]. Therefore, when we perturb $x_i(0)$, we can observe $x_0(t)$ as a linear combination of $i+1$ exponential functions. We can numerically determine the values of $a_{00},a_{11},\ldots,a_{ii}$. This means that when we perturb $x_i$, we can find $a_{00},a_{11},\ldots,a_{ii}$, which should correspond to $V_i$ and its descendants. The question is to find the correspondence between $a_{ii}$ and $V_i$. If a component $a_{nn}$ only appears after perturbing $x_n$, then we can make sure that $a_{n}$ corresponds to $V_n$. Then if $a_{jj}$ only appears after perturbing $x_n$ and $x_j$, then we know that $a_{jj}$ corresponds to $V_j$. Since there is a total order on $V_0,V_1,\ldots,V_n$, determined by the underlying DAG, we can determine the correspondence one by one. With the established correspondence, for each $V_i$, we can determine which genes are the descendants of $V_i$, by observing which components appear in $x_0(t)$.

For example, consider one phenotype $V_0$ and three genes $V_1,V_2,V_3$, and we can measure $x_0(t)$ after perturbing any one gene. If we perturb $x_1$, $x_0(t)$ has two exponential components $\lambda_0,\lambda_1$. If we perturb $x_2$, $x_0(t)$ has two exponential components $\lambda_0,\lambda_2$. If we perturb $x_3$, $x_0(t)$ has four exponential components $\lambda_0,\lambda_1,\lambda_2,\lambda_3$. Since component $\lambda_3$ only appears after perturbing $x_3$, $\lambda_3$ corresponds to $V_3$. Since component $\lambda_2$ only appears after perturbing $x_2$ or $x_3$, and $V_3$ already corresponds to $\lambda_3$, we know that $\lambda_2$ corresponds to $V_2$. Since component $\lambda_1$ only appears after perturbing $x_1$ or $x_3$, and $V_3$ already corresponds to $\lambda_3$, we know that $\lambda_1$ corresponds to $V_1$. Finally, $\lambda_0$ corresponds to $V_0$. Therefore, we know that $V_1$ has descendant $V_0$; $V_2$ has descendant $V_0$; $V_3$ has descendants $V_0,V_1,V_2$. We can determine that the DAG has edges $V_3\to V_2$, $V_3\to V_1$, $V_2\to V_0$, $V_1\to V_0$, with a possible edge $V_3\to V_0$.

In sum, assume that we can observe the phenotype variable $x_0(t)$ under the intervention on each ancestor gene $V_i$ of $V_0$. Then we can determine whether one gene $V_i$ is an ancestor of another gene $V_j$. Using the same method in Subsection <ref>, this ancestor-descendant relation can partially determine the DAG. Certainly, the calculation of $a_{00},a_{11},\ldots,a_{ii}$ from $x_0(t)$ is numerically infeasible, and requires unreasonable accuracy for $x_0(t)$. Therefore, this method is not applicable, unless the DAG is very small.

We can only observe the phenotype, not the genes, although we can intervene with genes. Consider a GRN $V_2\to V_1\to V_0$ with an edge $V_2\to V_0$. From the ancestor-descendant relations, we can only identify the edges $V_2\to V_1$ and $V_1\to V_0$. Assume $V_1$ activates $V_0$. If the level of $V_0$ decreases as the level of $V_2$ decreases (by an intervention), then we can explain that $V_2$ strongly activates $V_0$, and $V_2$ weakly activates or inhibits $V_1$. If the level of $V_0$ increases as the level of $V_2$ decreases (by an intervention), then we can explain that $V_2$ strongly inhibits $V_0$, and $V_2$ weakly activates or inhibits $V_1$. Therefore, we cannot determine whether $V_2\to V_1$ is activation or inhibition.

§.§ Use path blocking relation and partially observed ODE to infer GRN structure with phenotype data

Consider an unknown GRN with genes $V_1,\ldots,V_n$ and a phenotype $V_0$. Assume the GRN is a DAG. We can intervene with any genes, but we can only measure the phenotype on bulk level, not the genes. We assume there is a directed path from each $V_i$ to $V_0$, and there is no edge starting from $V_0$. Assume the path blocking property holds, and the levels of $V_0,V_1,\ldots,V_n$ satisfy a linear ODE system. Using the method in Subsection <ref>, we can infer some edges. Using the method in Subsection <ref>, we can also infer some edges. However, combining these two results might not fully recover the full GRN. For example, the two GRNs in Fig. <ref> have the same ancestor-descendant relations and path blocking relations. Thus we cannot determine whether $V_1\to V_3$ exists in this setting.

Two equivalent GRNs in the sense of ancestor-descendant relations and path blocking relations.

For the two GRNs in Fig. <ref>, we can identify edges $V_1\to V_2$ and $V_2\to V_3$, but not the edge $V_1$. Since $V_2\to V_3\to V_4\to V_0$ is the only directed path from $V_2$ to $V_0$, if we need to observe the influence of $V_2$ on $V_0$, or the influence of $V_1$ on $V_2$, we cannot intervene with $V_3$ or $V_4$ (otherwise, the path is blocked). Therefore, to determine the effect of $V_1$ on $V_2$, we can only interfere with $V_1$ or $V_2$. Similar to the discussions in Subsection <ref>, we cannot determine $V_1\to V_2$ is activation or inhibition.

§ DIFFERENT SCENARIOS OF INFERRING GRN STRUCTURE

In this section, we discuss the question that with each possible data type, whether the GRN structure can be inferred. As introduced in Section <ref>, each scenario (data type) is described in the following four dimensions: (1) Do we measure Gene expression or Phenotype? (2) Do we measure on Single-cell level or Bulk level? (3) Do we measure at a single time point (One-time) or at multiple time points (Time series)? (4) Is the measurement Interventional or Non-interventional? When the measurement is on the single-cell level and at multiple time points, we need to discuss whether we can measure the Joint distribution or Marginal distributions for multiple time points.

In order to place each discussion to the most proper scenario, this section might not be read fully in order. For example, Scenario 4 partially depends on the discussions in Scenario 6 and Scenario 8.

We shall mainly focus on methods that could use accurate data to produce a determined GRN, not various Bayesian type methods that produce a posterior distribution on possible GRNs. Nevertheless, Bayesian type methods have the advantage that they could easily utilize prior information [11].

For each scenario, there are different approaches to infer the GRN structure, corresponding to different models and different assumptions. We do not aim at exhausting all existing and potential approaches.

§.§ Scenario 1: gene, single-cell, one-time, non-interventional

In this scenario, we have the joint distribution of steady-state expression levels for all genes. As discussed in Subsection <ref>, if the GRN is a DAG, and the joint distribution is Markov and faithful to this DAG, then we can partially determine the GRN structure: all edges can be identified, except that the directions of some edges are unknown.

This method has been applied in determining GRN structures [12]. Various algorithms utilize experimental data to determine conditional independence relations [46, 47], and then infer DAG structures [48].

There are also other approaches to model this scenario. In a network, we can use neighboring relations to predict the characters associated with each node [49]. Reversely, we can use the node characters (gene expression profiles) to infer the network structures. This approach can only provide a rough guess among different possibilities. Structural equation model is also used to infer GRN structures in this scenario [15]. This method does not determine the true GRN, but evaluates all possible GRNs. Linearity is required for this method.

If the underlying dynamics satisfies the stochastic process model in Subsection <ref>, then at stationary, the joint distribution might not satisfy the Markov and faithful condition. For example, if the true GRN is $V_1\to V_2\to V_3$, then we only have $X_1(t-1)\ind X_3(t+1)\mid [X_2(0),\ldots,X_2(t)]$, or at most $X_1(t-1)\ind X_3(t+1)\mid X_2(t)$, but not $X_1(t)\ind X_3(t)\mid X_2(t)$. In this case, we will obtain fewer conditional independence relations.

§.§ Scenario 2: gene, single-cell, one-time, interventional

Scenario 2 is a combination of Scenario 1 and Scenario 6. Assume the GRN is a DAG, and the stationary joint distribution is Markov and faithful to this DAG. From the stationary distribution, we can obtain the conditional independence relations. From the interventions, we can obtain the ancestor-descendant relations. As discussed in Subsection <ref>, we can combine these two kinds of relations to fully reconstruct the DAG.

If the DAG condition holds, as discussed in Subsection <ref>, we can use ancestor-descendant relations to partially determine the GRN.

If the path blocking property holds under the interventions, similar to Scenario 6, we can fully determine the GRN.

§.§ Scenario 3: gene, single-cell, time series, non-interventional

In this scenario, we measure the single-cell gene expression profiles at multiple time points without intervention.

Scenario 3a: Joint. If we can measure the joint distribution for different time points, then as discussed in Subsection <ref>, we can fully reconstruct the GRN, even if there are cycles. If the measurement time step is too large, we can only obtain the ancestor-descendant relations, and we can only partially recover the GRN if it is a DAG.

This method has been applied in determining GRN structures [16]. There are some other approaches, such as Granger causality [50] and Boolean networks [17]. Readers may refer to a review paper [11] for different methods applicable for this scenario.

Scenario 3b: Marginal. If we can only measure the marginal distribution for each time point, then this is essentially the same as Scenario 1. If the GRN is a DAG, and the joint distribution is Markov and faithful to this DAG, then we can partially determine the GRN structure.

§.§ Scenario 4: gene, single-cell, time series, interventional

Scenario 4a: Joint. If we can measure the joint distribution for different time points, then similar to Scenario 3, we can fully reconstruct the GRN.

Scenario 4b: Marginal. If we can only measure the marginal distribution for each time point, there are different approaches. If we take expectations for expression profiles, this scenario moves to Scenario 8, and we can use the ODE model to fully recover the GRN under the linearity assumption. If the path blocking property holds under the interventions, similar to Scenario 6, we can use the blocking relations to fully recover the GRN. If the GRN is a DAG, and the stationary joint distribution is Markov and faithful to this DAG, similar to Scenario 2, we can use the conditional independence relations and the ancestor-descendant relations to fully recover the GRN. If the GRN is a DAG, similar to Scenario 6, we can use the ancestor-descendant relations to partially recover the GRN.

§.§ Scenario 5: gene, bulk, one-time, non-interventional

In this scenario, we can only observe the bulk level expression profile. Since there is no intervention, and the expression level of each gene is just a deterministic number, we cannot infer the structure of GRN. One might propose an approach that we can repeat the experiments on different cell populations, whose results might differ due to random noise. Then we can obtain a joint distribution of expression. This approach looks like Scenario 1, but it does not work. In Scenario 1, if there is an edge $V_1\to V_2$, the expression levels $X_1$ and $X_2$ might be $X_1=C+\epsilon_1$, $X_2=f(X_1)+\epsilon_2=f(C+\epsilon_1)+\epsilon_2$. Since the noise $\epsilon_1$ can transmit through the edge $V_1\to V_2$, we have $X_1\nind X_2$. In this scenario, the expression levels of $V_1$ and $V_2$ might be $X_1=C+\epsilon_1$, $X_2=f(C)+\epsilon_2$. The noise $\epsilon_1$ is external, so that it cannot transmit through the edge $V_1\to V_2$, and we have $X_1\ind X_2$.

§.§ Scenario 6: gene, bulk, one-time, interventional

In this scenario, we can interfere with any genes, then measure the bulk level expression profile. We can compare whether two different interventions produce the same result.

If the path blocking property holds under the interventions, as discussed in Subsection <ref>, there is an edge $V_i\to V_j$ if and only if they cannot be blocked by all other vertices. This means that we can fully determine the GRN even if it has cycles.

If the path blocking property does not hold, then through intervention, we can only obtain the ancestor-descendant relation: the intervention on $V_i$ leads to the change of $V_j$, if and only if $V_i$ is an ancestor of $V_j$. From Subsection <ref>, we can partially infer the GRN structure (at least $n-1$ edges for a GRN with $n$ vertices), if the GRN is a DAG. If the GRN has cycles, we might not infer any edge.

§.§ Scenario 7: gene, bulk, time series, non-interventional

Without intervention, the cell population is in equilibrium. Thus multiple observations at different time points should provide the same information. Scenario 7 is essentially the same as Scenario 5. We cannot infer the GRN structure.

§.§ Scenario 8: gene, bulk, time series, interventional

In this scenario, we can measure how the bulk expression level of each gene changes over time. Assume the dynamics is linear, meaning that the expression levels satisfy a linear ODE system. As discussed in Subsection <ref>, we can use the observed data to calculate the parameters, then determine the GRN structure. This method is extensively used [5, 36]. We can add some nonlinear terms in the ODE system to allow for more complicated dynamics. Certainly, this increases the quantity of data needed. When there are not enough data, we might need to assume that the coefficient matrix is sparse.

If we do not assume the dynamics is linear (or satisfying any other specified forms), then we meet the same situation as Scenario 6. If the path blocking property holds, we can fully determine the GRN. If the DAG assumption holds, we can partially infer the GRN structure through ancestor-descendant relations.

§.§ Scenario 9: phenotype, single-cell, one-time, non-interventional

In this scenario, we can neither perturb a gene nor measure the expression level of a gene. Since there is no information about genes in the data, we cannot infer the structure of GRN.

§.§ Scenario 10: phenotype, single-cell, one-time, interventional

We assume the path blocking property holds from any gene to the phenotype: the intervention on one gene $V_i$ does not change the phenotype $V_0$, if and only if all paths from $V_i$ to $V_0$ are blocked by intervened genes. For example, consider the GRN $V_1\to V_2\to V_3\to V_0$. If $V_3$ is already intervened, then the intervention on $V_2$ does not affect $V_0$. If just $V_1$ is intervened, then the intervention on $V_2$ could affect $V_0$.

Through intervening with different genes, we can determine the path blocking relation. As discussed in Subsection <ref>, we can use the path blocking relation to partially determine the GRN. We can determine at least $n$ edges if the GRN has $n$ genes.

§.§ Scenario 11: phenotype, single-cell, time series, non-interventional

Scenario 11a/b: Joint/Marginal. Similar to Scenario 9, there is no information about genes in the data, and we cannot infer the structure of GRN. In this scenario, whether we can measure the joint distribution for multiple time points does not change the conclusion.

§.§ Scenario 12: phenotype, single-cell, time series, interventional

Scenario 12a/b: Joint/Marginal. This scenario is almost the same as Scenario 16. The only difference is that we measure single-cell level phenotype data, which have higher accuracy in detecting differences. For a GRN with $n$ genes, we can identify at least $n$ edges if (1) the path blocking property holds; or (2) the dynamics is linear, and the GRN is a DAG; or (3) the path blocking property holds, the dynamics is linear, and the GRN is a DAG. Notice that in (2) and (3), for some identified edges, we cannot determine whether the regulatory relation is activation or inhibition. In this scenario, whether we can measure the joint distribution for multiple time points does not change the conclusion.

§.§ Scenario 13: phenotype, bulk, one-time, non-interventional

Similar to Scenario 9, there is no information about genes in the data, and we cannot infer the structure of GRN.

§.§ Scenario 14: phenotype, bulk, one-time, interventional

This scenario is almost the same as Scenario 10. The only difference is that we measure bulk level phenotype data, which have lower accuracy in detecting differences. We can partially determine the GRN (at least $n$ edges if the GRN has $n$ genes).

§.§ Scenario 15: phenotype, bulk, time series, non-interventional

Similar to Scenario 9, there is no information about genes in the data, and we cannot infer the structure of GRN.

§.§ Scenario 16: phenotype, bulk, time series, interventional

If the path blocking property holds, then as discussed in Subsection <ref>, we can partially determine the GRN structure (at least $n$ edges for a GRN with $n$ genes).

If the dynamics is linear, and the GRN is a DAG, then as discussed in Subsection <ref>, we can determine the ancestor-descendant relations, and partially determine the GRN structure (at least $n$ edges for a GRN with $n$ genes). Notice that for some identified edges, we cannot determine whether the regulatory relation is activation or inhibition.

If the path blocking property holds, the dynamics is linear, and the GRN is a DAG, then as discussed in Subsection <ref>, we can combine the above two methods and partially infer the GRN structure. Notice that for some identified edges, we cannot determine whether the regulatory relation is activation or inhibition.

§ HOW TO EXAMINE THE EXISTENCE OF AUTOCATALYSIS

Consider a certain gene $V$ and its expression level (e.g., mRNA count) $x$. In general, when the expression levels of other genes that affect $V$ are kept at constants, the synthesis rate of $V$ is a constant $k$. Besides, the total degradation rate $cx$ of $V$ is proportional to $x$, where $c$ is the coefficient of degradation. Therefore, the overall change rate of $x$ is $k-cx$.

In special cases, there exists autocatalysis, meaning that the synthesis rate of $V$ increases with the current expression level $x$. Thus the synthesis rate can be expressed as $k+f(x)$, where $f(x)$ is an increasing function of $x$. Autocatalysis can be observed if the overall change rate of $x$ increases with $x$ itself. Since $x$ cannot increase indefinitely, autocatalysis can happen only for $x$ in a certain range. This means that $f(x)$ has to be nonlinear. In nonlinear systems, if $V_1,V_2$ form a positive feedback loop, it is possible that the increase of $V_1$ leads to a higher stable point [51]. This can be regarded as generalized autocatalysis.

If we assume autocatalysis is linear, then the synthesis rate is $k+bx$, where $b$ is the coefficient of autocatalysis. Now the overall change rate of $x$ is $k+bx-cx$. Since $x$ cannot grow to infinity, we need to have $b<c$. Here we assume that all other genes have fixed expression levels, so that we can focus on one gene.

In ODE models, the change rate of $x$ is $\mathrm{d}x/\mathrm{d}t=k-(c-b)x$. With the bulk level gene expression data, we can use $\mathrm{d}x/\mathrm{d}t$ and $x$ to calculate $k$ and $(c-b)$. However, we cannot determine whether $b=0$ or $b>0$, since $c=10,b=7$ and $c=3,b=0$ are equivalent in this model. The existence of autocatalysis is masked by the existence of degradation, whose rate also increases with the expression level.

The identification of autocatalysis is possible only with single-cell gene expression data (Scenarios 1-4). In stochastic process models, we can examine the existence of autocatalysis just with the stationary distribution. The idea is that if we increase both $b$ and $c$, the variance would increase [41]. With the assumptions above, when the expression level $X=n$, the synthesis rate, namely the transition rate from $X=n$ to $X=n+1$, is $bn+k$; the degradation rate, namely the transition rate from $X=n$ to $X=n-1$, is $cn$. At stationary, define $P_n$ to be the probability that the expression level $X=n$. Since $\{P_n\}$ is stationary, we have

\begin{equation}
\label{eq3}
\end{equation}

In Eq. (<ref>), take the summation from $n=1$ to $\infty$ to obtain $b\mathbb{E}(X)+k=c\mathbb{E}(X)$. Thus $\mathbb{E}(X)=k/(c-b)$. Now multiply Eq. (<ref>) by $n-1$, and take the summation from $n=1$ to $\infty$. Then we have $b\mathbb{E}(X^2)+k\mathbb{E}(X)=c\mathbb{E}(X^2)-c\mathbb{E}(X)$. Thus $\mathbb{E}(X^2)=(k^2+ck)/(c-b)^2$, and $\sigma^2(X)=ck/(c-b)^2$. We can see that the mean and variance of $X$ satisfy $\sigma^2(X)/\mathbb{E}(X)=c/(c-b)$. If $\sigma^2(X)=\mathbb{E}(X)$, then $b=0$; if $\sigma^2(X)>\mathbb{E}(X)$, then $b>0$. We have reached a criterion for autocatalysis: at stationary, if the variance of $X$ is larger than the mean of $X$, then $X$ has autocatalysis. To apply this method, we need the single-cell level distribution of gene expression profiles.

This criterion is valid only if other genes are fixed, or do not interfere with the considered gene. In reality, all genes in the same GRN reach stationary as a network. Thus the marginal distribution of a gene at stationary might be different from the distribution discussed above. In Scenarios 1-3, if we can confirm that a gene is not affected by other genes (which is not true in general), then we can use the criterion to determine whether autocatalysis exists.

If we have single-cell level interventional gene expression time series data (Scenario 4), we can assume linearity, and calculate the values of $c_i-b_i$ and $k_i$ for each gene in ODE models, as discussed in Subsection <ref>. Then we can use the variance and covariance of gene expression to solve the value of each $b_i$ and $c_i$, similar to the one-dimensional case. Since the system is much more complicated, the condition of $b_i=0$ might explicitly depend on the values of $c_i-b_i$ and $k_i$. Therefore, we have to concretely calculate all $b_i$ to see if they equal $0$, not just comparing mean and variance.

In sum, when the gene expression profile changes linearly, we can determine the existence of autocatalysis for some genes with single-cell gene expression data (Scenarios 1-4). If the data are interventional time series (Scenario 4), we can determine autocatalysis for all genes. Here the extra dimension of Joint vs. Marginal in Scenario 4 does not matter, since we only need the interventional time series to calculate $b_i,c_i,k_i$ with ODE methods.

§ DISCUSSION

In this paper, we introduce the problem of GRN structure inference and provide complete answers to all situations. To handle some less-discussed scenarios, we develop some novel mathematical methods. We introduce a method to determine the existence of autocatalysis in certain scenarios.

This paper does not cover all mathematical approaches to the GRN structure inference problem, and there will be new mathematical methods for scenarios discussed in this paper. With the development of biological technologies, there will be new data types that are not included in our 20 scenarios. This work just provides a general paradigm of inferring GRN (or general networks) structures: given what type of data, under what assumptions, what we can infer about the GRN structure.

We only focus on theoretical inference methods, but not practical algorithms that work on real data. The designing of such algorithms is related to some mathematical fields not covered in this paper: statistical inference [52], numerical computing [53], etc.

For well-studied Scenarios 1, 3, 8, we mention multiple algorithms that work on finite experimental data. Although the principles of different algorithms are more or less the same, the efficiencies are different. The reason is that we only have measured data, not accurate values/distributions, so that different approaches might introduce different levels of errors and have different computational time costs [44]. (A classical result in numerical linear algebra: In theory, Gaussian elimination and Cramer's rule can both solve a linear equation system. In practice, Gaussian elimination is significantly faster and more accurate [54].) Besides, some algorithms have certain assumptions, under which the computation is much faster, and the data required is much less. For example, if we need to determine the independence of two random variables, and it is assumed that they are joint Gaussian, then we only need to check whether their correlation coefficient is $0$.

ODE-based models (including dynamical system methods) have an advantage: once the parameters are determined, we can directly predict the behavior under new perturbations. To derive ODE models, we assume there are infinitely many identical cells living in a stationary environment, and the underlying dynamics is time-homogeneous. These assumptions are not always true: we only have finitely many cells; cells have heterogeneity; the environment keeps changing; cells are mutating over time. Besides, this derivation might fail in nonlinear cases. For a nonlinear function $f$ and a random variable $X$, the expectation generally cannot cross $f$: $\mathbb{E}[f(X)]\ne f(\mathbb{E}X)$. Therefore, the ODE-based models cannot match the reality with extremely high accuracy, which conflicts with the fact that some ODE-based models (such as the method in Subsection <ref>) are numerically sensitive to perturbations.

When there are not enough data, one solution is to produce more data through interpolation [5]. With this method, the dynamics is already stipulated by the data interpolation method. The new data can only be used to confirm that their dynamics follow the interpolation method. Therefore, this approach is essentially equivalent to adding more assumptions (e.g., linearity) into the model. We need to be cautious not to regard assumptions as conclusions.

In Scenario 6 and Scenario 8, we have the bulk level gene expression profiles after different interventions. One approach of determining possible gene regulations is to directly calculate the correlation coefficient or mutual information of the expression levels of two genes [14]. Here data after different interventions are regarded as samples from the same distribution. We provide a simple example to show that this approach might build a false relation between two independent genes. Assume there are three independent genes $V_1,V_2,V_3$. After intervention on $V_1$, the expression levels of $V_1,V_2$ are $0.5$ and $1$; after intervention on $V_2$, the expression levels of $V_1,V_2$ are $1$ and $0.5$; after intervention on $V_3$, the expression levels of $V_1,V_2$ are $1$ and $1$. Then the correlation coefficient of $V_1,V_2$ is $-0.5$, although $V_1,V_2$ are independent.

When we measure bulk level gene expression and repeat several times, due to various external noise, the measured quantities might not be exactly equal. When we measure at multiple time points (time series), and each cell population can be measured only once, we cannot match the quantity at one time point to the quantity at another time point. For example, assume for the same population, one quantity fluctuates periodically like a sine wave. Different populations have the same frequency but different phases. Since we cannot measure the same population more than once, we can only observe that the marginal distribution at each time point is the same, and the fluctuation cannot be identified. This problem is even related to the thermodynamic side of gene expression dynamics: if we do not know the joint distribution of multiple time points, then we cannot distinguish equilibrium steady state and non-equilibrium steady state (which differ by entropy productions) [55]. This problem also appears for the single-cell level time series data.

For the methods in Subsection <ref> (causal DAG), Subsection <ref> (stochastic process), when the number of genes in the GRN increases, the conditional independence tests need to concern more variables, and the quantity of data needed grows exponentially. Besides, the computational time cost also grows exponentially. For the method in Subsections <ref>, <ref> (path blocking), the number of experiments and computational time cost also grow exponentially with the number of genes. To avoid such problems, we need to assume the GRN is sparse, or just consider a few genes. Nevertheless, such solutions might lead to unreliable results, since GRNs in reality might be large and dense [5].

§ ACKNOWLEDGEMENTS

Y.W. would like to thank Dr. Nadya Morozova and Dr. Andrey Minarsky for helpful discussions.

[1]
Patke A, Young MW, Axelrod S.
Molecular mechanisms and physiological importance of circadian
Nat Rev Mol Cell Biol. 2020;21(2):67-84.

[2]
Wang Y, Kropp J, Morozova N.
Biological notion of positional information/value in morphogenesis
Int J Dev Biol. 2020;64(10-11-12):453-63.

[3]
Müller-Hill B.
The lac Operon.
de Gruyter; 2011.

[4]
Cunningham TJ, Duester G.
Mechanisms of retinoic acid signalling and its roles in organ and
limb development.
Nat Rev Mol Cell Biol. 2015;16(2):110-23.

[5]
Bansal M, Gatta GD, Di Bernardo D.
Inference of gene regulatory networks and compound mode of action
from time course gene expression profiles.
Bioinformatics. 2006;22(7):815-22.

[6]
Cao J, O'Day DR, Pliner HA, Kingsley PD, Deng M, Daza RM, et al.
A human cell atlas of fetal gene expression.
Science. 2020;370(6518).

[7]
Galluzzi L, Vitale I, Senovilla L, Olaussen KA, Pinna G, Eisenberg T, et al.
Prognostic impact of vitamin B6 metabolism in lung cancer.
Cell Rep. 2012;2(2):257-69.

[8]
Brugere I, Gallagher B, Berger-Wolf TY.
Network structure inference, a survey: Motivations, methods, and
ACM Comput Surv. 2018;51(2):1-39.

[9]
Breza E, Chandrasekhar AG, McCormick TH, Pan M.
Using aggregated relational data to feasibly identify network
structure without network data.
Am Econ Rev. 2020;110(8):2454-84.

[10]
Mateos G, Segarra S, Marques AG, Ribeiro A.
Connecting the dots: Identifying network structure via graph signal
IEEE Signal Process Mag. 2019;36(3):16-43.

[11]
Penfold CA, Wild DL.
How to infer gene networks from expression profiles, revisited.
Interface Focus. 2011;1(6):857-70.

[12]
Emmert-Streib F, Glazko G, Altay G, De Matos Simoes R.
Statistical inference and reverse engineering of gene regulatory
networks from observational expression data.
Front Genet. 2012;3:8.

[13]
Wang T, Li B, Nelson CE, Nabavi S.
Comparative analysis of differential gene expression analysis tools
for single-cell RNA sequencing data.
BMC Bioinformatics. 2019;20(1):1-16.

[14]
Hurley D, Araki H, Tamada Y, Dunmore B, Sanders D, Humphreys S, et al.
Gene network inference and visualization tools for biologists:
application to new human transcriptome datasets.
Nucleic Acids Res. 2012;40(6):2377-98.

[15]
Xiong M, Li J, Fang X.
Identification of genetic networks.
Genetics. 2004;166(2):1037-52.

[16]
Yu J, Smith VA, Wang PP, Hartemink AJ, Jarvis ED.
Advances to Bayesian network inference for generating causal
networks from observational biological data.
Bioinformatics. 2004;20(18):3594-603.

[17]
Shmulevich I, Dougherty ER, Zhang W.
Gene perturbation and intervention in probabilistic Boolean
Bioinformatics. 2002;18(10):1319-31.

[18]
Huynh-Thu VA, Sanguinetti G.
Gene regulatory network inference: an introductory survey.
In: Gene Regulatory Networks. Springer; 2019. p. 1-23.

[19]
Carrier TA, Keasling JD.
Investigating autocatalytic gene expression systems through
mechanistic modeling.
J Theor Biol. 1999;201(1):25-36.

[20]
Durrett R.
Probability: theory and examples. vol. 49.
Cambridge university press; 2019.

[21]
Diestel R.
Graph Theory.
2nd ed. New York: Springer-Verlag; 2000.

[22]
Pearl J, Glymour M, Jewell NP.
Causal inference in statistics: A primer.
John Wiley & Sons; 2016.

[23]
Murray JD.
Mathematical Biology I: An Introduction.
Berlin, Heidelberg: Springer-Verlag; 2001.

[24]
Ding J, Adiconis X, Simmons SK, Kowalczyk MS, Hession CC, Marjanovic ND, et al.
Systematic comparison of single-cell and single-nucleus
RNA-sequencing methods.
Nat Biotechnol. 2020;38(6):737-46.

[25]
Perrimon N, Ni JQ, Perkins L.
In vivo RNAi: today and tomorrow.
Cold Spring Harb Perspect Biol. 2010;2(8):a003640.

[26]
Andrews TS, Hemberg M.
Identifying cell populations with scRNASeq.
Mol Aspects Med. 2018;59:114-22.

[27]
Wang Z, Gerstein M, Snyder M.
RNA-Seq: a revolutionary tool for transcriptomics.
Nat Rev Genet. 2009;10(1):57-63.

[28]
Sinkala E, Sollier-Christen E, Renier C, Rosas-Canyelles E, Che J, Heirich K,
et al.
Profiling protein expression in circulating tumour cells using
microfluidic western blotting.
Nat Commun. 2017;8(1):1-12.

[29]
Svensson V.
Droplet scRNA-seq is not zero-inflated.
Nat Biotechnol. 2020;38(2):147-50.

[30]
So LH, Ghosh A, Zong C, Sepúlveda LA, Segev R, Golding I.
General properties of transcriptional time series in Escherichia
Nat Genet. 2011;43(6):554-60.

[31]
Kaern M, Elston TC, Blake WJ, Collins JJ.
Stochasticity in gene expression: from theories to phenotypes.
Nat Rev Genet. 2005;6(6):451-64.

[32]
Rondina MT, Voora D, Simon LM, Schwertz H, Harper JF, Lee O, et al.
Longitudinal RNA-Seq analysis of the repeatability of gene
expression and splicing in human platelets identifies a platelet SELP
splice QTL.
Circ Res. 2020;126(4):501-16.

[33]
Nikam RR, Gore KR.
Journey of siRNA: clinical developments and targeted delivery.
Nucleic Acid Ther. 2018;28(4):209-24.

[34]
Pickar-Oliver A, Gersbach CA.
The next generation of CRISPR–Cas technologies and applications.
Nat Rev Mol Cell Biol. 2019;20(8):490-507.

[35]
Gujral TS, Peshkin L, Kirschner MW.
Exploiting polypharmacology for drug target deconvolution.
Proc Natl Acad Sci USA. 2014;111(13):5048-53.

[36]
Wang YK, Hurley DG, Schnell S, Print CG, Crampin EJ.
Integration of steady-state and temporal gene expression data for the
inference of gene regulatory networks.
PLOS ONE. 2013;8(8):e72103.

[37]
Munkácsy G, Sztupinszki Z, Herman P, Bán B, Pénzváltó Z,
Szarvas N, et al.
Validation of RNAi silencing efficiency using gene array data shows
18.5% failure rate across 429 independent experiments.
Mol Ther Nucleic Acids. 2016;5:e366.

[38]
Jiang DQ, Wang Y, Zhou D.
Phenotypic equilibrium as probabilistic convergence in
multi-phenotype cell population dynamics.
PLOS ONE. 2017;12(2):e0170916.

[39]
Reidling JC, Subramanian VS, Dahhan T, Sadat M, Said HM.
Mechanisms and regulation of vitamin C uptake: studies of the
hSVCT systems in human liver epithelial cells.
Am J Physiol Gastrointest Liver Physiol. 2008;295(6):G1217-27.

[40]
Liu SL, Needham KM, May JR, Nolen BJ.
Mechanism of a concentration-dependent switch between activation and
inhibition of Arp2/3 complex by coronin.
J Biol Chem. 2011;286(19):17039-46.

[41]
Wang Y.
Some Problems in Stochastic Dynamics and Statistical Analysis of
Single-Cell Biology of Cancer [Ph.D. thesis].
University of Washington; 2018.

[42]
Pollicott M, Wang H, Weiss H.
Extracting the time-dependent transmission rate from infection data
via solution of an inverse ODE problem.
J Biol Dyn. 2012;6(2):509-23.

[43]
Ein-Dor L, Kela I, Getz G, Givol D, Domany E.
Outcome signature genes in breast cancer: is there a unique set?
Bioinformatics. 2004;21(2):171-8.

[44]
Wang Y, Wang L.
Causal inference in degenerate systems: An impossibility result.
In: Proceedings of the Twenty Third International Conference on
Artificial Intelligence and Statistics. PMLR; 2020. p. 3383-92.

[45]
Zhou D, Wang Y, Wu B.
A multi-phenotypic cancer model with cell plasticity.
J Theor Biol. 2014;357:35-45.

[46]
Huang TM.
Testing conditional independence using maximal nonlinear conditional
Ann Stat. 2010;38(4):2047-91.

[47]
Zhang K, Peters J, Janzing D, Schölkopf B.
Kernel-based conditional independence test and application in causal
In: Proceedings of the Twenty-Seventh Conference on Uncertainty in
Artificial Intelligence; 2011. p. 804-13.

[48]
Pearl J.
Probabilistic Inference in Intelligent Systems.
San Mateo: Morgan Kaufmann; 1988.

[49]
Wang Y, Zhang B, Kropp J, Morozova N.
Inference on tissue transplantation experiments.
J Theor Biol. 2021;520:110645.

[50]
Barnett L, Barrett AB, Seth AK.
Granger causality and transfer entropy are equivalent for Gaussian
Phys Rev Lett. 2009;103(23):238701.

[51]
Santillán M, Mackey MC, Zeron ES.
Origin of bistability in the lac operon.
Biophys J. 2007;92(11):3830-42.

[52]
Casella G, Berger RL.
Statistical inference.
Cengage Learning; 2021.

[53]
Pozrikidis C.
Numerical computation in science and engineering.
Oxford university press New York; 1998.

[54]
Trefethen LN, Bau III D.
Numerical Linear Algebra.
SIAM; 1997.

[55]
Wang Y, Qian H.
Mathematical representation of Clausius' and Kelvin's statements
of the second law and irreversibility.
J Stat Phys. 2020;179:808-37.