Inference on autoregulation in gene expression

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 are affected (activated or inhibited) by the expression levels of other genes [1], which is called (mutual) gene regulation. Genes and their regulatory relations form a gene regulatory network (GRN) [2], 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 GRN.

The expression of one gene could promote/repress its own expression, which is called positive/negative autoregulation [4]. Autoregulation is very common in E. coli [5]. Positive autoregulation is also called autocatalysis or autoactivation, and negative autoregulation is also called autorepression [6, 7]. For instance, HOX proteins form and maintain spatially inhomogeneous expression of HOX genes [8]. For genes with position-specific expressions during development, it is common that the increase of one gene can further increase or decrease its level [9].

While countless works infer the regulatory relationships between different genes (GRN structure) [10], 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 stationary) 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 lower cost [11].

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) [12, 13, 14, 15, 16, 17, 18]. 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 [19]. 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 [20]. 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 [21]. 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. The state space is on/off (for $G$) or the mRNA/protein count on $\mathbb{Z}$ (for $M,P$). When we consider the expression level $M$ or $P$ (but do not control $G$), sometimes itself still follows a Markov chain, and we call this scenario “autonomous”. In other cases, $M$ or $P$ itself is no longer Markovian, and we call this scenario “non-autonomous”. We need to consider the triplet $G,M,P$ in the non-autonomous scenario.

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.

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$, 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 Poisson 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 [22, 23], quadratic functions [24], or Hill functions [25]. There are other papers that consider Markov chain models in gene expression/regulation [26, 27, 28, 29, 30, 31], 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 [32, 33, 34, 26]. 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, which cannot be fixed. Due to this non-controllable factor, there might be transcriptional bursting [35, 36] or translational bursting [37], 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. Assume $V^{\prime}$ is not contained in a feedback loop, and assume $g_{n}$, the per molecule degradation rate of $V^{\prime}$, is not regulated by other genes or its interior factors. 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}\geq 1$. Therefore, $\text{VMR}<1$ means autoregulation for $V^{\prime}$. The conclusion “$\text{VMR}<1$ means autoregulation” has been observed by Munsky et al. [15] for a single gene that is non-autonomous.

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$ [38]. 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.

There are other mathematical approaches to infer the existence of autoregulation in gene expression [39, 40, 41, 42, 43, 44, 45, 46]. We introduce some works and compare them with our method. (A) Sanchez-Castillo et al. [39] 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. [40] 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. [41] applied a Bayesian method for model selection. This method (1) needs time series data; (2) estimates a group of parameters. (D) Veerman et al. [42] 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. [43] 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 methods, our method has some advantages: (1) Our method uses 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 our method are much easier and cheaper to obtain. (2) Our method does 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. (3) Our method has few restrictions on the model, making them applicable to various scenarios with different dynamics. In sum, our method is simple and universal, and have lower requirements on data quality.

Compared to other 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 stationary distribution, sometimes it is impossible to build the causal relationship (including autoregulation) [47]. Thus with this data type, some disadvantages are inevitable.