Inference on autoregulation in gene expression with variance-to-mean ratio

Some genes can promote or repress their own expressions, which is called autoregulation. Although gene regulation is a central topic in biology, autoregulation is much less studied. In general, it is extremely difficult to determine the existence of autoregulation with direct biochemical approaches. Nevertheless, some papers have observed that certain types of autoregulations are linked to noise levels in gene expression. We generalize these results by two propositions on discrete-state continuous-time Markov chains. These two propositions form a simple but robust method to infer the existence of autoregulation from gene expression data. This method only needs to compare the mean and variance of the gene expression level. Compared to other methods for inferring autoregulation, our method only requires non-interventional one-time data, and does not need to estimate parameters. Besides, our method has few restrictions on the model. We apply this method to four groups of experimental data and find some genes that might have autoregulation. Some inferred autoregulations have been verified by experiments or other theoretical works.


PLCG
JNK P38 Figure 1: An example of GRN in human T cells [3].Each vertex is a gene.Each arrow is a regulatory relationship.Notice that it has no directed cycle.

Introduction
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).

Setup and main results
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 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 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 inactivated gene The mechanism of gene expression.A gene might switch between inactivated state and activated state.Gene is transcribed into mRNAs, which are translated into proteins.Proteins might (auto)regulate the state transition of the corresponding gene.
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 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 ′ < h n ′ −1 for some other n ′ , 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 P(X = n)/P(X = n − 1) = h/n, meaning that the distribution is Poisson with parameter h, P(X = n) = h n e −h /n!, and VMR = 1.If there exists positive autoregulation of certain forms, VMR > 1; if there exists negative autoregulation of certain forms, 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 VMR > 1 or 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 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 VMR ≥ 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 ′ in a known GRN.Assume V ′ is not contained in a feedback loop, and assume g n , the per molecule degradation rate of V ′ , 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 ′ has no autoregulation, then VMR ≥ 1.Therefore, VMR < 1 means autoregulation for V ′ .The conclusion "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 VMR ≥ 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 VMR ≥ 1 [38].In the non-autonomous scenario, we only focus on the less complicated case of 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 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 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.

Related works
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 VMR ≥ 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.
4 Scenario of a single isolated gene

Setup
We first consider the expression level (e.g., mRNA count) of one gene V in a single cell.At the single-cell level, gene expression is essentially stochastic, and we use a random variable X to represent the mRNA count of V .We assume V is not in a feedback loop.We also assume all environmental factors and other genes that can affect X are kept at constant levels, so that we can focus on V alone.This can be achieved if no other genes point to gene V in the GRN, such as PIP3 in Fig. 1.Then we assume that the expression of V is autonomous, thus X satisfies a time-homogeneous Markov chain defined on Z * .
Assume that the mRNA synthesis rate at X(t) = n − 1, namely the transition rate from X = n−1 to X = n, is f n > 0. Assume that with n mRNA molecules, the degradation rate for each mRNA molecule is g n > 0. Then the overall degradation rate at X(t) = n, namely the transition rate from X = n to X = n − 1, is g n n.The associated master equation is Define the relative growth rate h n = f n /g n .We assume that h n has a finite upper bound.This means that as time tends to infinity, the process reaches equilibrium.Thus at equilibrium, (1) the stationary probability distribution P n = lim t→∞ P[X(t) = n] exists, and (2) the mean E(X) and the variance σ 2 (X) are finite [48].If h n > h n−1 for some n, then there exists positive autoregulation.If h n < h n−1 for some n, then there exists negative autoregulation.If there is no autoregulation, then h n is a constant h, and the stationary distribution is Poisson with parameter h.In this setting, positive autoregulation and negative autoregulation might coexist, meaning that h n+1 < h n for some n and h n ′ +1 > h n ′ for some n ′ .

Theoretical results
With single-cell non-interventional one-time gene expression data for one gene, we have the stationary distribution of the Markov chain X.We can infer the existence of autoregulation with the VMR of X, defined as VMR(X) = σ 2 (X)/E(X).The idea is that if we let f n increase/decrease with n, and control g n to make E(X) invariant, then the variance σ 2 (X) increases/decreases [49].We shall prove that VMR > 1 means positive autoregulation, and VMR < 1 means negative autoregulation.Notice that VMR > 1 does not exclude the possibility that negative autoregulation exists for some expression level.This also applies to VMR < 1 and positive autoregulation.
We can illustrate this result with a linear model: set Lemma 1.Consider the Markov chain model for one gene with general transition coefficients f n , g n , described by Eq. 1. Calculate VMR(X) at stationary.(1) Assume h n+1 ≥ h n for all n.We have VMR(X) ≥ 1; moreover, VMR(X) = 1 if and only if h n+1 = h n for all n.(2) Assume h n+1 ≤ h n for all n.We have VMR(X) ≤ 1; moreover, VMR(X) = 1 if and only if h n+1 = h n for all n.
From Lemma 1, we can directly obtain the following proposition.Proposition 1.In the setting of Lemma 1, (1) If VMR(X) > 1, then there exist values of n for which h n+1 > h n ; thus this gene has positive autoregulation.(2) If VMR(X) < 1, then there exist values of n for which h n+1 < h n ; thus this gene has negative autoregulation.
(3) If VMR(X) = 1, then either (A) h n+1 = h n for all n, meaning that this gene has no autoregulation; or (B) h n+1 < h n for some n and h n ′ +1 > h n ′ for some n ′ , meaning that this gene has both positive and negative autoregulation (at different expression levels).
Remark 1. Proposition 1 requires that the gene expression is autonomous.In reality, many genes are non-autonomous, and transcriptional/translational bursting can make the VMR to be larger than 100 [22].
Remark 2. Results similar to Proposition 1 have been proven in a non-autonomous model of gene expression [50].
Proof of Lemma 1. Define λ = − log P 0 , so that P 0 = exp(−λ).Define d n = n i=1 h i > 0 and stipulate that d 0 = 1.We can see that Also, Besides, ( Here the first inequality is from the Cauchy inequality, and the second inequality is from for all n (the second inequality of Eq. 2).The equality condition is equivalent to h n+1 = h n for all n.
(2) Assume h n+1 ≤ h n for all n.Then d n+2 /d n+1 ≤ d n+1 /d n , and Since 0 < d n ≤ h n 1 , this series converges for all t ∈ C, so that H(t) is a well-defined analytical function on C, and In the following, we only consider H(t), H ′ (t), H ′′ (t) as real functions for t ∈ R.
Consider any t 1 , t 2 ∈ I with t 1 ≤ t 2 and any p, q ∈ N with p ≥ q.Since d p+1 /d p ≤ d q+1 /d q , and t p−q 1 ≤ t p−q 2 , we have About the condition for the equality to hold, assume h n ′ +1 < h n ′ for a given n ′ .Then for all t 1 , t 2 ∈ I with t 1 ≤ t 2 and a constant C that does not depend on t 1 , t 2 .Therefore, Since H(t) has a finite positive upper bound A and a positive lower bound B on I, we have and thus ∀t ∈ I, ] 2 } < 0, and VMR(X) < 1.
We have proved in (1) that if h n+1 = h n for all n, then VMR(X) = 1.Thus when h n+1 ≤ h n for all n, VMR(X) = 1 if and only if h n+1 = h n for all n.
In sum, for the Markov chain model of one gene (by assuming the expression to be autonomous), when we have the stationary distribution from single-cell non-interventional one-time gene expression data, we can calculate the VMR of X. VMR(X) > 1 means the existence of positive autoregulation, and VMR(X) < 1 means the existence of negative autoregulation.VMR(X) = 1 means either (1) no autoregulation exists; or (2) both positive autoregulation and negative autoregulation exist (at different expression levels).

Scenario of multiple entangled genes 5.1 Setup
We consider m genes V 1 , . . ., V m for a single cell.Denote their expression levels by random variables X 1 , . . ., X m .The change of X i can depend on X j (mutual regulation) and X i itself (autoregulation).Since these genes regulate each other, and their expression levels are not fixed, we cannot consider them separately.If the expression of gene V k is non-autonomous, we also need to add its interior factors (gene activation state, etc.) into X 1 , . . ., X m .
We can use a continuous-time Markov chain on (Z * ) m to describe the dynamics.Each state of this Markov chain, (X 1 = n 1 , . . ., X i = n i , . . ., X m = n m ), can be abbreviated as n = (n 1 , . . ., n i , . . ., n m ).For gene V i , the transition rate of n i − 1 → n i is f i (n), and the the transition rate of n (3) Define nī = (n 1 , . . ., n i−1 , n i+1 , . . ., n m ).Define h i (n) = f i (n)/g i (n) to be the relative growth rate of gene V i .Autoregulation means for some fixed nī, h i (n) is (locally) increasing/decreasing with n i , thus f i (n) increases/decreases and/or g i (n) decreases/increases with n i .For the non-autonomous scenario, another possibility for autoregulation is that V i can affect its interior factors.

Theoretical results
With expression data for multiple genes, there are various methods to infer the regulatory relationships between different genes, so that the GRN can be reconstructed [10].In the GRN, if there is a directed path from gene V i to gene V j , meaning that V i can directly or indirectly regulate V j , then V i is an ancestor of V j , and V j is a descendant of V i .
Fix a gene V k in a GRN.We first consider a simple case that V k is not contained in any directed cycle (feedback loop), which means no gene is both an ancestor and a descendant of V k , such as PIP2 in Fig. 1.This means V k itself is a strongly connected component of the GRN.This condition is automatically satisfied if the GRN has no directed cycle.If the expression of V k is non-autonomous, we need to add the interior factors of V k into V 1 , . . ., V m , and it is acceptable that V k regulates its interior factors.In this case, we can prove that if V k does not regulate itself, meaning that h k (n) is a constant for fixed nk and different n k , and X k does not affect its interior factors (if non-autonomous), then VMR(X k ) ≥ 1.The reason is that VMR < 1 requires either a feedback loop or autoregulation.We need to assume that the per molecule degradation rate g k (•) for V k is not affected by V 1 , . . ., V m , which is not always true in reality [1].With this result, when VMR < 1, there might be autoregulation.Proposition 2. Consider the Markov chain model for multiple genes, described by Eq. 3.
Assume the GRN has no directed cycle, or at least there is no directed cycle that contains gene V k .Assume g k (•) is a constant for all n.If V k has no autoregulation, meaning that h k (•) and f k (•) do not depend on n k , and V k does not regulate its interior factors, then This Proposition is in an unpublished work by Paulsson et al., who study a similar problem [51,52].It also appears in a preprint by Mahajan et al. [53], but the proof is based on a linear approximation.We propose a rigorous proof independently.

Proof. Denote the expression level of
For simplicity, denote the expression levels of V 1 , . . ., V l by a (high-dimensional) random variable Y .Assume V k has no autoregulation.Since V k does not regulate V 1 , . . ., V l , W does not affect Y .Denote the transition rate from Y = i to Y = j by q ij ≥ 0. Stipulate that q ii = − j =i q ij .When Y = i, the transition rate from W = n to W = n + 1 is F i (does not depend on n), and the transition rate from W = n to W = n − 1 is G.

The master equation of this process is dP
Assume there is a unique stationary probability distribution P n,i = lim t→∞ P[W (t) = n, Y (t) = i].Then we have Define P i = n P n,i .Sum over n for Eq. 4 to obtain meaning that P i is the stationary probability distribution of Y .Define W i to be W conditioned on Y = i at stationary.Then P(W i = n) = P(W = n | Y = i) = P n,i /P i , and E(W i ) = n nP n,i /P i .Multiply Eq. 4 by n and sum over n to obtain G + j Sum over i for Eq. 6 to obtain Multiply Eq. 4 by n 2 and sum over n to obtain Sum over i for Eq. 8 to obtain Multiply Eq. 6 by E(W i ) and sum over i to obtain Then we have i Here the first equality is from Eq. 10, the third equality is from Eq. 5, and other equalities are equivalent transformations.Now we have where the first equality is by definition, the second equality is from Eqs. 7,9, the first inequality is from Eq. 11, the third equality is from i P i = 1, and the second inequality is the Cauchy inequality.
We hypothesize that the requirement for g k (•) in Proposition 2 can be dropped.
Conjecture 1. Assume V k is not contained in a directed cycle in the GRN, and V k does not regulate its interior factors.If V k has no autoregulation, meaning that h k (•) does not depend on n k (but might depend on nk), then V k has VMR ≥ 1.
If the GRN has directed cycles, there is a conjecture by Paulsson et al. [51,52], which has been numerically verified but not proved yet.Conjecture 2. Assume for each V i , g i (•) does not depend on n, and f i (•) does not depend on n i (no autoregulation).Then for at least one gene V j , we have VMR ≥ 1.
Notice that Conjecture 2 does not hold if g i depends on nī.One counterexample is Then VMR = 2e/(4e − 1) ≈ 0.55 for both genes.
Assume Conjecture 2 is correct.For m genes, if we find that VMR for each gene is less than 1, then we can infer that autoregulation exists, although we do not know which gene has autoregulation.
6 Applying theoretical results to experimental data We summarize our theoretical results into Algorithm 1. Proposition 1 applies to a gene that has no ancestor in the GRN.However, it requires the corresponding gene has autonomous expression, which is difficult to validate and often does not hold in reality.Thus the inference result by Proposition 1 for VMR > 1 (positive autoregulation) is not very reliable.When VMR < 1 and Proposition 1 could apply, we should instead apply Proposition 2 to determine the existence of autoregulation, since Proposition 2 does not require the expression to be autonomous, thus being much more reliable.Proposition 2 applies when the gene is not in a feedback loop and has VMR < 1.Notice that our result cannot determine that a gene has no autoregulation.
For a given gene without autoregulation, its expression level satisfies a Poisson distribution, and VMR is 1.If we have n samples of its expression level, then the sample VMR (sample variance divided by sample mean) asymptotically satisfies a Gamma distribution Γ[(n−1)/2, 2/(n−1)], and we can determine the confidence interval of sample VMR [54].If the sample VMR is out of this confidence interval, then we know that VMR is significantly different from 1, and Propositions 1,2 might apply.
We apply our method to four groups of single-cell non-interventional one-time gene expression data from experiments, where the corresponding GRNs are known.Notice that we need to convert indirect measurements into protein/mRNA count.See Table 1 for our inference results and theoretical/experimental evidence that partially validates our results.See Appendix A for details.There are 186 genes in these four data sets, and we can only determine that 12 genes have autoregulation (7 genes determined by Proposition 1, and 5 genes determined by Proposition 2).Not every VMR is less than 1, so that Conjecture 2 does not apply.For the other 174 genes, Proposition 1 and Proposition 2 do not apply, and we do not know whether they have autoregulation.
In some cases, we have experimental evidence that some genes have autoregulation, so that we can partially validate our inference results.Nevertheless, as discussed in the Introduction, there is no gold standard to evaluate our inference results.
In the data set by Guo et al. [55], Sanchez-Castillo et al. [39] inferred that 17 of 39 genes have autoregulation, and 22 genes do not have autoregulation.We infer that 5 genes have autoregulation, and 34 genes cannot be determined.Here 3 genes are shared by both inference results to have autoregulation.Consider a random classifier that randomly picks 5 genes and claims they have autoregulation.Using Sanchez-Castillo et al. as the standard, this random classifier has probability 62.55% to be worse than our result, and 10.17% to be better than our result.Thus our inference result is better than a random classifier, but the advantage is not significant.

Conclusions
For a single gene that is not affected by other genes, or a group of genes that form a connected GRN, we develop theoretical results to determine the existence of autoregulation.These results generalize known relationships between autoregulation and VMR by dropping restrictions on parameters.Our results only depend on VMR, which is easy to compute and more robust than other complicated statistics.We also apply our method to experimental data and detect some genes that might have autoregulation.We prove two propositions for Markov chains, which might have theoretical values.We introduce two conjectures that have been numerically verified but not yet proved.They are of theoretical interest and worth further consideration.The Markov chain models in this paper can be studied via lifting into a higher-dimensional space [63], treating as a random dynamical system [64], or as a branching process [65].With the expression profiles of different genes, we can construct a similarity graph [66].If we know the existence of autoregulation for some genes, we can use the similarity graph to infer other genes.
Our method requires independent and identically distributed samples from the exact stationary distribution of a fully observed Markov chain, plus a known GRN.Proposition 1 requires the expression is autonomous.Proposition 2 requires that the GRN has no directed cycle, and degradation is not regulated.If our inference fails, then some requirements are not met: (1) cells might affect each other, making the samples dependent; (2) cells are heterogeneous; (3) the measurements have extra errors; (4) the cells are not at stationary; (5) there exist unobserved variables that affect gene expression; (6) the GRN is inferred by a theoretical method, which can be interfered by the existence of autoregulation; (7) the expression is non-autonomous; (8) the GRN has unknown directed cycles; (9) the degradation rate is regulated by other genes.Such situations, especially the unobserved variables, are unavoidable.Therefore, current data might not satisfy these requirements, and our inference results should be interpreted as informative findings, not ground truths.In fact, other theoretical works that determine gene autoregulation, or general gene regulation, also need similar assumptions and might fail.Nevertheless, with the development of experimental technologies, there will be more data with higher quality that fit the requirements of our method.Thus we believe that our method will be more applicable in the future.
Cells keep growing and dividing, and the gene expression fluctuates along the cell cycle.Discussions on such non-stationary situations can be found in other papers [20,67,68,69].
About cell heterogeneity, we prove a result in Appendix B that if several cell types have VMR ≥ 1, then for a mixed population of such cell types, we still have VMR ≥ 1.Therefore, cell heterogeneity does not fail Proposition 2, since VMR < 1 for the mixture of several cell types means VMR < 1 for at least one cell type.

A Details of applications on experimental data
In experiments, the expression levels of genes are not directly measured as mRNA or protein counts.Rather, they are measured as cycle threshold (Ct) values or fluorescence intensity values.Such indirected measurements need to be converted.Related details can be found in other papers [50].
Guo et al. [55] measured the expression (mRNA) levels of 48 genes for mouse embryo cells at different developmental stages.We consider three groups (16-cell stage, 32-cell stage, 64-cell stage) that have more than 50 samples.Sanchez-Castillo et al. [39] used such data to infer the GRN structure, including autoregulation, but the GRN only contains 39 genes.Thus we ignore the other 9 genes.In the inferred GRN, genes BMP4, CREB312, and TCFAP2C are not contained in directed cycles.In the 16 to infer that FN1 and HNF4A (VMR > 1 for all three cell groups) might have positive autoregulation.Nevertheless, it is more likely that the expressions of FN1 and HNF4A are non-autonomous, and there is no autoregulation.Sanchez-Castillo et al. [39] inferred that BMP4, HNF4A, TCFAP2C have autoregulation.Besides, there is experimental evidence that BMP4 [56], HNF4A [57], TCFAP2C [58] have autoregulation.Therefore, our inference results are partially validated.Psaila et al. [59] measured the expression (mRNA) levels of 90 genes for human megakaryocyte-erythroid progenitor cells.Chan et al. [70] inferred the GRN structure (autoregulation not included).In the inferred GRN, genes BIM, CCND1, ECT2, PFKP have no ancestors.BIM has 214 effective samples, and VMR is 187.7,outside of the 95% confidence interval [0.8191, 1.1987].CCND1 has 68 effective samples, and VMR is 111.3, outside of the 95% confidence interval [0.6905, 1.3660].ECT2 has 56 effective samples, and VMR is 8.2, outside of the 95% confidence interval [0.6618, 1.4069].PFKP has 134 effective samples, and VMR is 82.1, outside of the 95% confidence interval [0.7742, 1.2543].Thus we can apply Proposition 1 to infer that BIM, CCND1, ECT2, PFKP might have positive autoregulation.Nevertheless, it is more likely that the expressions of these four genes are non-autonomous, and there is no autoregulation.There is experimental evidence that ECT2 has autoregulation [60], which partially validates our inference results.No other gene fits the requirement of Proposition 2.

Fig. 1 ) 1 :
and VMR < 1 Output V k has autoregulation // Assume the degradation of V k is not regulated by V 1 , . . ., V m Else If V k has no ancestor in the GRN (like PIP3 in Fig. 1) and VMR > 1 Output V k has autoregulation //Assume the expression of V k is autonomous Else Output We cannot determine whether V k has autoregulation End of if End of if Algorithm Detailed workflow of inferring autoregulation with gene expression data.
) ≥ 1. Equality holds if and only if d n /d n+2 = d n+1 /d n+3 for all n (the first inequality of Eq. 2) and d n d n+2 = d 2

Table 1 :
The autoregulation inference results by our method on four data sets.Source column is the paper that contains this data set.Proposition 1 column is the genes that can be only inferred by Proposition 1 to have autoregulation.Proposition 2 column is the genes that can be inferred by Proposition 2 to have autoregulation.Theory column is the genes inferred by both our method and other theoretical works to have autoregulation.Experiment column is the genes inferred by both our method and other experimental works to have autoregulation.Bold font means the inferred gene with autoregulation is validated by other results.Details can be found in Appendix A.
-cell stage group with 75 samples, if there is no autoregulation, then the 95% confidence interval of VMR is [0.7041, 1.3470].BMP4 (VMR = 0.2139), CREB312 (VMR = 0.1971), and TCFAP2C (VMR = 0.3468) have significantly small VMR, and we can apply Proposition 2 to infer that BMP4, CREB312, and TCFAP2C might have autoregulation.In the other two groups, these genes do not have VMR < 1, and the results are relatively weak.Besides, in the inferred GRN, genes FN1 and HNF4A have no ancestors.For the 16-cell stage with 75 samples, the VMR of FN1 and HNF4A are 3.4522 and 1.3599, outside of the 95% confidence interval [0.7041, 1.3470]; for the 32-cell stage with 113 samples, the VMR of FN1 and HNF4A are 93.1070 and 46.7688, outside of the 95% confidence interval [0.7554, 1.2784]; for the 64-cell stage with 159 samples, the VMR of FN1 and HNF4A are 117.3059and 93.9589, outside of the 95% confidence interval [0.7917, 1.2322].Thus we can apply Proposition 1