Using Fano factors to determine certain types of gene autoregulation

In general, genes are transcribed to mRNAs and then translated to proteins. We can use the abundance of mRNA or protein to represent the expression levels of genes. Both the synthesis and degradation of mRNAs/proteins can be affected (activated or inhibited) by the expression levels of other genes [42], which is called (mutual) gene regulation. Genes and their regulatory relations form a gene regulatory network (GRN) [16], generally represented as a directed graph: each vertex is a gene, and each directed edge is a regulatory relationship. See Fig. 1 for an example of a GRN.

The expression of one gene could promote/repress its own expression, which is called positive/negative autoregulation [11]. Autoregulation is very common in E. coli [63]. Positive autoregulation is also called autocatalysis or autoactivation, and negative autoregulation is also called autorepression [4, 21]. For instance, HOX proteins form and maintain spatially inhomogeneous expression of HOX genes [64]. For genes with position-specific expressions during development, it is common that the increase of one gene can further increase or decrease its level [81]. Autoregulation has the effect of stabilizing transposons in genomes [6], which can affect cell behavior [41, 79]. Autoregulation can also stabilize the cell phenotype [2], which is related to cancer development [96, 52, 15, 90].

While countless works infer the regulatory relationships between different genes (the GRN structure) [86], determining the existence of autoregulation is an equally important yet less-studied field. Due to technical limitations, it is difficult and sometimes impossible to directly detect autoregulation in experiments. Instead, we can measure gene expression profiles and infer the existence of autoregulation. In this paper, we consider a specific data type: measure the expression levels of certain genes without intervention for a single cell (which reaches stationarity) at a single time point, and repeat for many different cells to obtain a probability distribution for expression levels. Such single-cell non-interventional one-time gene expression data can be obtained with a relatively low cost [48].

With such single-cell level data for one gene $V$, we can calculate the ratio of variance and mean of the expression level (mRNA or protein count). This quantity is called the variance-to-mean ratio (VMR) or the Fano factor. Many papers that study gene expression systems with autoregulations have found that negative autoregulation can decrease noise (smaller VMR), and positive autoregulation can increase noise (larger VMR) [70, 68, 30, 51, 25, 18, 17]. This means VMR can be used to infer the existence of autoregulation.

We generalize the above observation and develop two mathematical results that use VMR to determine the existence of autoregulation. They apply to some genes that have autoregulation. For genes without autoregulation, these results cannot determine that autoregulation does not exist. We apply these results to four experimental gene expression data sets and detect some genes that might have autoregulation.

We start with some setup and introduce our main results (Section 2). Then we cite some previous works on this topic and compare them with our results (Section 3). For a single gene that is not regulated by other genes (Section 4) and multiple genes that regulate each other (Section 5), we develop mathematical results to identify the existence of autoregulation. These two mathematical sections can be skipped. We summarize the procedure of our method and apply it to experimental data (Section 6). We finish with some conclusions and discussions (Section 7).

One possible mechanism of “the increase of one gene’s expression level further increases its expression level” is a positive feedback loop between two genes [31]. Here $V_{1}$ and $V_{2}$ promote each other, so that the increase of $V_{1}$ increases $V_{2}$, which in return further increases $V_{1}$. We should not regard this feedback loop as autoregulation. When we define autoregulation for a gene $V$, we should fix environmental factors and other genes that regulate $V$, and observe whether the expression level of $V$ can affect itself. If $V$ is in a feedback loop that contains other genes, then those genes (which regulate $V$ and are regulated by $V$) cannot be fixed when we change $V$. Therefore, it is essentially difficult to determine whether $V$ has autoregulation in this scenario. In the following, we need to assume that $V$ is not contained in a feedback loop that involves other genes.

The actual gene expression mechanism might be complicated. Besides other genes/factors that can regulate a gene, for a gene $V$ itself, it might switch between inactivated (off) and activated (on) states [10]. These states correspond to different transcription rates to produce mRNAs. When mRNAs are translated into proteins, those proteins might affect the transition of gene activation states, which forms autoregulation [23]. See Fig. 2 for an illustration. Therefore, for a gene $V$, we should regard the gene activation state, mRNA count, and protein count as a triplet of random variables $(G,M,P)$, which depend on each other.

When we fix environmental factors and other genes that affect $V$, the triplet $(G,M,P)$ should follow a continuous-time Markov chain. A possible state is the gene activation state on/off (for $G$), the mRNA count on $\mathbb{Z}$ (for $M$), and the protein count on $\mathbb{Z}$ (for $P$). Thus the total space is $\{0,1\}\times\mathbb{Z}\times\mathbb{Z}$. When we consider the expression level $M$ or $P$ (but have no access to the value of $G$), sometimes itself is Markovian (its dynamics can be fully determined by itself, without the knowledge of $G$), and we call this scenario “autonomous”. In other cases, $M$ or $P$ itself is no longer Markovian (its dynamics explicitly depends on $G$), and we call this scenario “non-autonomous”. We need to consider the triplet $(G,M,P)$ in the non-autonomous scenario. This is similar to a hidden Markov model, where a two-dimensional Markov chain is no longer Markovian if projected to one dimension (since this dimension depends on the other dimension).

For the autonomous scenario, we can fully classify autoregulation for a gene $V$. Assume environmental factors and other genes that affect the expression of $V$ are kept at constants. Define the expression level (mRNA count for example) of one cell to be $X=n$, the mRNA synthesis rate at $X=n-1$ to be $f_{n}$, and the degradation rate for each mRNA molecule at $X=n$ to be $g_{n}$. This is a standard continuous-time Markov chain on $\mathbb{Z}$ with transition rates

Define the relative growth rate $h_{n}=f_{n}/g_{n}$. If there is no autoregulation, then $h_{n}$ is a constant. Positive autoregulation means $h_{n}>h_{n-1}$ for some $n$, so that $f_{n}>f_{n-1}$ and/or $g_{n}<g_{n-1}$; negative autoregulation means $h_{n}<h_{n-1}$ for some $n$, so that $f_{n}<f_{n-1}$ and/or $g_{n}>g_{n-1}$. Notice that we can have $h_{n}>h_{n-1}$ for some $n$ and $h_{n^{\prime}}<h_{n^{\prime}-1}$ for some other $n^{\prime}$, meaning that positive autoregulation and negative autoregulation can both exist for the same gene, but occur at different expression levels. Such non-monotonicity in regulating gene expression often appear in reality [1].

For the non-autonomous scenario, we can still define autoregulation. Consider the expression level $X$ of $V$ (mRNA count or protein count) and its interior factor $I$. If $X$ is the mRNA count, then $I$ is the gene state; if $X$ is the protein count, then $I$ is the gene state and the mRNA count. If there is no autoregulation, then $X$ cannot affect $I$, and for each value of $I$, the relative growth rate $h_{n}$ of $X$ is a constant. If $X$ can affect $I$, or $h_{n}$ is not a constant, then there is autoregulation. When $X$ can affect $I$, there is a directed cycle ($X\to I\to X$), and the change of $X$ can affect itself through $I$. In this case, it is not always easy to distinguish between positive autoregulation and negative autoregulation.

Quantitatively, for the autonomous scenario, when we fix other factors that might regulate this gene $V$, if $V$ has no autoregulation, then $h_{n}=f_{n}/g_{n}$ is a constant $h$ for all $n$. In this case, the stationary distribution of $V$ satisfies $\mathbb{P}(X=n)/\mathbb{P}(X=n-1)=h/n$, meaning that the distribution is Poissonian with parameter $h$, $\mathbb{P}(X=n)=h^{n}e^{-h}/n!$, and $\text{VMR}=1$. If there exists positive autoregulation of certain forms, $\text{VMR}>1$; if there exists negative autoregulation of certain forms, $\text{VMR}<1$. However, such results are derived by assuming that $f_{n},g_{n}$ take certain functional forms, such as linear functions [54, 57], quadratic functions [24], or Hill functions [67]. There are other papers that consider Markov chain models in gene expression/regulation [32, 61, 65, 14, 62, 45], but the role of VMR is not thoroughly studied.

In this paper, we generalize the above result of inferring autoregulation with VMR by dropping the restrictions on parameters. Consider a gene $V$ in a known GRN, and assume it is not regulated by other genes, or assume other factors that regulate $V$ are fixed. Assume we have the autonomous scenario, meaning that its expression level $X=n$ satisfies a general Markov chain with synthesis rate $f_{n}$ and per molecule degradation rate $g_{n}$. We do not add any restrictions on $f_{n}$ and $g_{n}$. Use the single-cell non-interventional one-time gene expression data to calculate the VMR of $V$. Proposition 1 states that $\text{VMR}>1$ or $\text{VMR}<1$ means the existence of positive/negative autoregulation.

Nevertheless, the autonomous condition requires some assumptions, and often does not hold in reality [5, 39, 33, 32]. Consider a gene $V$ that is not regulated by other genes, and has no autoregulation. The mRNA count or the protein count is regulated by the gene activation state (an interior factor), which cannot be fixed. Due to this non-controllable factor, there might be transcriptional bursting [60, 19] or translational bursting [8], where transcription or translation can occur in bursts, and we have $\text{VMR}>1$. This does not mean that Proposition 1 is wrong. Instead, it means that the expression level itself is not Markovian, and the scenario is non-autonomous. In this scenario, we should apply Proposition 2, described below, which states that no autoregulation means $\text{VMR}\geq 1$.

We extend the idea of inferring autoregulation with VMR to a gene that is regulated by other genes, or with non-autonomous expression. Consider a gene $V^{\prime}$ in a known GRN. Denote other genes that regulate $V^{\prime}$ and the interior factors (gene state and/or mRNA count) of $V^{\prime}$ by $\boldsymbol{F}$. Denote the values of $V^{\prime},\boldsymbol{F}$ as $X,\boldsymbol{Y}$. Assume $V^{\prime}$ is not contained in a feedback loop, and assume $g_{n}$, the per molecule degradation rate of $V^{\prime}$ at $X=n$, is not regulated by other genes or its interior factors (gene state and/or mRNA count). We do not add any restrictions on the synthesis rate $f_{n}$. Proposition 2 states that if $V^{\prime}$ has no autoregulation, then $\text{VMR}(X)\geq 1$. Therefore, $\text{VMR}(X)<1$ means autoregulation for $V^{\prime}$.

Proposition 2 is derived in a “one-step” Markov chain model, where at one time point, only transitions to the nearest neighbors are allowed: $(X=n,\boldsymbol{Y}=\boldsymbol{a})\to(X=n+1,\boldsymbol{Y}=\boldsymbol{a})$, $(X=n,\boldsymbol{Y}=\boldsymbol{a})\to(X=n-1,\boldsymbol{Y}=\boldsymbol{a})$, and $(X=n,\boldsymbol{Y}=\boldsymbol{a})\to(X=n,\boldsymbol{Y}=\boldsymbol{a}^{\prime})$. This one-step Markov chain model is the most common approach in stochastic representations of gene regulation [70, 30, 54, 51, 17]. Recently, there are some studies that consider “multi-step” Markov chain models, where at one time point, the change of mRNA/protein count can be accompanied with the change of other factors, such as the gene state [7, 43, 74]. For example, the following transition is allowed: $(G^{*},M=n)\to(G,M=n+1)$. In this multi-step model, Proposition 2 is no longer valid: even without autoregulation, it is possible that $\text{VMR}(X)<1$. Consider an example that the production of one mRNA molecule needs many steps of gene state transition, and the gene returns to the initial step after producing one mRNA molecule: $G_{1}\to G_{2}\to\cdots\to G_{k}\to G_{1}+M$, $M\to\emptyset$. Since there are many steps, the total time for one cycle of $G_{1}\to\cdots\to G_{k}\to G_{1}+M$ can be highly deterministic, such as $1$ second. Assume the degradation probability for each mRNA molecule in $1$ second is $0.01$. Then the mRNA count is highly concentrated near $100$, and $\text{VMR}(X)<1$ (close to $0.5$ in numerical simulations).

Since multi-step models allow more transitions, they are more general than one-step models. However, it is still a question that whether such generalizations are necessary, since one-step models have good fitting with experimental data [38, 18, 9]. Proposition 2 provides a method to verify this problem: If a gene has $\text{VMR}(X)<1$, but we use other methods to determine that it has no autoregulation, then Proposition 2 states that one-step models deviate from reality, and multi-step models should be adopted. Therefore, when one-step models hold, Proposition 2 is a valid method to determine the existence of autoregulation; when one-step models do not hold, combined with other methods to determine autoregulation, Proposition 2 can detect the failure of one-step models.

In the scenario that Proposition 2 may apply, if $\text{VMR}\geq 1$, Proposition 2 cannot determine whether autoregulation exists. In fact, with VMR, or even the full probability distribution, we might not distinguish a non-autonomous system with autoregulation from a non-autonomous system without autoregulation, which both have $\text{VMR}\geq 1$ [9]. In the non-autonomous scenario, we only focus on the less complicated case of $\text{VMR}<1$, and derive Proposition 2 that firmly links VMR and autoregulation.

In reality, Proposition 1 and Proposition 2 can only apply to a few genes (which are not regulated by other genes or have $\text{VMR}<1$), and they cannot determine negative results. Thus the inference results about autoregulation are a few “yes” and many “we do not know”. Besides, for the results inferred by Proposition 1, especially those with $\text{VMR}>1$ (positive autoregulation), we cannot verify whether their expression is autonomous, and the inference results are less reliable.

Current experimental methods can hardly determine the existence of autoregulation, and to determine that a gene does not have autoregulation is even more difficult. Therefore, about whether genes in a GRN have autoregulation, experimentally, we do not have “yes” or “no”, but a few “yes” and many “we do not know”. Thus there is no gold standard to thoroughly evaluate the performance of our inference results. We can only report that some genes inferred by our method to have autoregulation are also verified by experiments or other inference methods to have autoregulation. Besides, if the result by Proposition 2 does not match with other methods, it is possible that the one-step model fails. Instead, in Section A, we test our methods with numerical simulations, and the performances of both Propositions are satisfactory.

There have been some results of inferring autoregulation with VMR [70, 68, 30, 51, 25, 18, 17]. However, these VMR-based methods have various restrictions on the model, and some of them are derived by applying linear noise approximations, which are not always reliable in gene regulatory networks [72].

Besides VMR-based methods, there are other mathematical approaches to infer the existence of autoregulation in gene expression [59, 92, 22, 73, 37, 97, 35, 34]. We introduce some works and compare them with our method. (A) Sanchez-Castillo et al. [59] considered an autoregressive model for multiple genes. This method (1) needs time series data; (2) requires the dynamics to be linear; (3) estimates a group of parameters. (B) Xing et al. [92] applied causal inference to a complicated gene expression model. This method (1) needs promoter sequences and information on transcription factor binding sites; (2) requires linearity for certain steps; (3) estimates a group of parameters. (C) Feigelman et al. [22] applied a Bayesian method for model selection. This method (1) needs time series data; (2) estimates a group of parameters. (D) Veerman et al. [73] considered the probability-generating function of a propagator model. This method (1) needs time series data; (2) estimates a group of parameters; (3) needs to approximate a Cauchy integral. (E) Jia et al. [37] compared the relaxation rate with degradation rate. This method (1) needs interventional data; (2) only works for a single gene that is not regulated by other genes; (3) requires that the per molecule degradation rate is a constant.

Compared to other more complicated methods, VMR-based methods (including ours) have two advantages: (1) VMR-based methods use non-interventional one-time data. Time series data require measuring the same cell multiple times without killing it, and interventional data require some techniques to interfere with gene expression, such as gene knockdown. Therefore, non-interventional one-time data used in VMR-based methods are much easier and cheaper to obtain. (2) VMR-based methods do not estimate parameters, and only calculate the mean and variance of the expression level. Some other methods need to estimate many parameters or approximate some complicated quantities, meaning that they need large data size and high data accuracy. Therefore, our method is easy to calculate, and need lower data accuracy and smaller data size.

Compared to other VMR-based methods, our method has few restrictions on the model, making them applicable to various scenarios with different dynamics. Besides, our derivations do not use any approximations.

In sum, compared to other VMR-based methods, our method is universal. Compared to other non-VMR-based methods, our method is simple, and has lower requirements on data quality.

Compared to other non-VMR-based methods, our method has some disadvantages: (1) The GRN structure needs to be known. (2) Our method does not work for certain genes, depending on regulatory relationships. Proposition 1 only works for a gene that is not regulated by other genes, and we require its expression to be autonomous; Proposition 2 only works for a gene that is not in a feedback loop. (3) Proposition 2 requires the per molecule degradation rate to be a constant, and it cannot provide information about autoregulation if $\text{VMR}\geq 1$. (4) Our method only works for cells at equilibrium. Thus time series data that contain time-specific information cannot be utilized other than treated as one-time data. With just the stationary distribution, sometimes it is impossible to build the causal relationship (including autoregulation) [85]. Thus with this data type, some disadvantages are inevitable. Such impossibility results might be generalized to other data types or even other fields [77].