Pseudospectral method for fourth-order fractional Sturm-Liouville problems

1.   Introduction

Enhancers, located upstream or downstream from the transcription start sites (TSSs), are important cis-acting elements, which can activate or promote the transcription levels of their target genes, independent of their position and orientation to the target genes ^[1,2,3,4,5]. It has been reported that there are more than 400,000 enhancers in human genome ^[6,7], and this number is far more than that of protein-encoding genes. Active enhancers can recruit a number of transcription factors (TFs) and RNAPII to control precise gene expression in tissue-specific patterns required for the directional development of cell types and lineages ^{[8,9,10,11,12]}.

Biological experiments have reported fundamental characteristics of enhancers: (1) enhancers with the promoters of protein-encoding genes can form chromatin loops, thus producing a large amount of eRNAs ^[13,14]; (2) Functionally active enhancers can be transcribed into different types of eRNAs in an unidirectional or a bidirectional manner ^{[9,15,16,17,18,19]}; (3) eRNAs produced by enhancers are extremely unstable because they are quickly degraded by the RNA exosome, a ribonucleic acid exonuclease complex involved in the degradation and processing of intracellular RNAs ^{[20,21,22,23]}; (4) The eRNA abundance is so low that eRNAs are often viewed as the byproducts of enhancer transcription ^[24]. These experimental facts provide guidelines for modeling and analysis of gene expression regulated by eRNAs, but biological functions of the eRNAs remain elusive.

Two theories have been proposed to interpret the roles of eRNAs. The first theory is the so-called "transcription noise theory", which considers that eRNAs are merely the byproducts of enhancer transcription without functional significance. Some biologists supported this idea, e.g., Albin Sandelin thought that enhancer and the regulatory promoter may form a stable DNA loop when high concentrations of RNAP and TFs accumulate ^[25]. The excess RNAPII machinery collides stochastically with enhancers to initiate transcription ^{[18,19,20,21,22,23,24,25,26]}. The second theory is the so-called "functional theory", which considers that eRNAs themselves may play a role of activating gene expression since transcription does not seem a random phenomenon but would be controlled by various factors ^{[27,28,29,30,31,32,33]}. For example, stationary enhancers have no transcriptional activities, meaning that chromatin accessibility is decided by DNase I hypersensitivity ^[34]. With regard to the "functional theory", several hypotheses have been made to describe the role of eRNAs at distinct stages: (1) Enhancer-promoter (E-P) looping. For example, in human breast cancer cells, the eRNAs expressed from estrogen receptors that bind enhancers promote specific E-P interactions by recruiting adhesion proteins to the enhancers ^[27]. Also for example, in prostate cancer cells, the eRNAs expressed from the Kallikrein-related peptidase 3 enhancer accelerate the interaction between the KLK3 enhancer and the KLK2 promoter, by forming a complex with the androgen receptor and subunit of Mediator complex, MED1 ^[29]; (2) Chromatin remodeling. During the myogenic differentiation of C2C12 skeleton muscle cells, eRNAs expressed from MYOD1 core enhancer facilitate the accessibility of chromatin and enhance the occupancy of RNAPII to the promoter region and the subsequent gene expression ^[31]; (3) Early transcription elongating. In the induction of immediate early genes in neuronal cells, eRNAs facilitate the efficient release of negative elongation factors (NELF) from their target gene promoters through direct association with the RNA recognition motif in the NELF-E subunit, thus promoting the synchronous expressions of the target genes ^[33]. A schematic diagram for the interaction between enhancer and promoter is shown in Figure 1. Which of the above three hypotheses is more reasonable is debating. In particular, the mechanism of how eRNAs affect gene expression is unclear.

So far, many molecular mechanisms have been proposed to explain the impact of eRNAs on gene expression, but to our knowledge, there seem no theoretical models to investigate eRNAs' potential functions. In this paper, we build a mechanistic model of stochastic gene expression, which considers both an E-P loop characterizing the E-P interaction and different feedbacks representing different ways that eRNAs impact gene expression. Although simple, this model still captures fundamental events taking place in gene expression, such as random synthesis and degradation of eRNA and mRNA, and stochastic switching between promoter activity states as well as eRNA regulation of the mRNA expression. It also considers three possible regulation ways: Type-I regulation by which eRNAs regulates transcriptional activity by facilitating the formation of the E-P loop's formation, Type-II regulation by which eRNAs directly increase the production rate of mRNA, and the combination of these two ways (called mixed regulation). Note that the "transcription noise theory" does not consider regulation.

We analytically derive approximate steady-state mRNA and eRNA probability distributions as well as their statistical quantities such as mean and expression noise (defined as the ratio of the variance over the squared mean). Interestingly, we find that with the increase of the E-P loop length, the mRNA probability distribution can switch from unimodality to bimodality and then be back to unimodality in all the three cases of eRNA regulation. We also find that in contrast to the other two regulations, the Type-II regulation can result in the highest mean mRNA level and the lowest mRNA noise, independent of the E-P loop length. Since these results are qualitative, we conclude that the E-P loop can induce phenotypic switching, which may be utilized by organisms or live cells. Our analysis validates the eRNA's "functional theory", whereas our results would have potential biological implications.

2.   Model description

2.1.   Mathematical model

Although the interaction between the enhancer and the promoter of a gene would be diverse (in particular, in eukaryotic cells), here we assume that they form a loop (a representative chromatin structure) so that RNAPII can access a specific binding site to initiate transcription. This assumption is reasonable. In fact, experimental evidences indicated that the formation of a DNA loop is beneficial to structural stability and thus to recruitment of TFs and RNA polymerases. For example, E. Savitskaya et al. verified by experiment that long-distance functional interactions between Su (Hw) insulators that can regulate the E-P communication in Drosophila melanogaster ^[35]. Looping is critical for process ranging from gene regulation to recombination and repair, and the boundaries of loop domains (insulators) are the determinants of chromosome form and function in multicellular eukaryotes ^[36]. M. Julien et al. suggested that the genome segment delimited by topological insulators and comprising the gene and its enhancer regulatory sequences is the plausible adaptor, which mediating the correspondence between the concentrations of specific transcription factors in the cell nucleus and gene expression level ^[37]. In addition, P. Charles et al. proposed that the Activity-by contact model of enhancer-promoter regulation from thousands of CRISPR perturbations to map and predict which enhancers regulate which genes ^[38]. We also assume that eRNA and mRNA are expressed only when the E-P loop is formed. This assumption is also reasonable since a promoter must recruit TFs and RNA polymerases to initiate transcription. Using the idea of the mapping method proposed in ref. ^[39], we build a mechanistic model of stochastic gene expression to model the effect of eRNA on gene expression, referring to Figure 2. This model is, in essence, similar to the extensively studied ON-OFF model of stochastic transcription ^{[40,41,42,43,44]}. That is, the gene promoter has one active (ON) state that corresponds to the case that the E-P loop is formed so that RNAPII can access a specific binding site to initiate transcription, where DNAs are transcribed into mRNAs and enhancers are transcribed into eRNAs, and one inactive (OFF) state that corresponds to the case that the E-P loop is not formed, where transcription is not existent, namely, neither the gene nor the enhancer is expressed. Furthermore, stochastic transitions exist between ON and OFF states.

In addition, the model also considers feedback. Usually, mRNA and eRNA need to use the same resource, polymerases, which, in principle, may couple two populations of the molecules. If a polymerase pool is sufficiently large, then we can neglect the mRNA-eRNA feedback bridged by polymerase. The produced eRNAs then serve as transcriptional factors to regulate the target gene expressions. Based on experimental evidences, we introduce three possible regulation ways (or mechanisms) of eRNA: Type-I regulation, Type-II regulation and mixed regulation, as described in the introduction. Note that the first regulation mechanism corresponds to the first hypothesis in the "functional theory", whereas the second regulation mechanism integrates the latter two hypotheses in the "functional theory". Nevertheless, it is unclear which of the three regulation mechanisms is more possibly utilized by live organisms or cells.

For clarity, we list all biochemical reactions corresponding to Figure 2A as follows

where ${k_0}$ and ${k_1}$ are respectively transition rates from inactive ($I$) to active ($A$) states and vice versa in the absence of eRNA regulation, ${\mu _1}$ and ${\mu _2}$ are respectively transcription rates, ${\delta _1}$and ${\delta _2}$ are respectively degradation rates of mRNA and eRNA. In Eq 1, $\left[{eRNA} \right]$ represents the number of eRNA molecules or the eRNA concentration, $\alpha$and $\beta$ represent two feedback strengths and are assumed as nonnegative constants. If $\alpha > 0$ but $\beta \, {\text{ = }}\, 0$, then the model reduces to the case of DNA looping. If $\alpha = 0$ but $\beta > 0$, then it reduces to the case of Type-II regulation. If $\alpha > 0$ and $\beta > 0$, then it corresponds to the case of mixed regulation. If $\alpha = 0$ and $\beta = 0$, then it corresponds to the case of no regulation. Thus, we can study the effect of eRNA on the mRNA expression in a united framework.

Let $m$and $n$ represent the numbers of mRNA and eRNA molecules, respectively. To better trace the time evolution of probability distributions of mRNA and eRNA molecule numbers, we introduce two distribution functions ${P_I}\left({m, n; t} \right)$ and ${P_A}\left({m, n; t} \right)$, which represent the probabilities that mRNA have $m$ molecules and eRNA have $n$molecules at time $t$when the gene is respectively at $I$ and $A$ states. Thus, the total probability, denoted by$P\left({m, n; t} \right)$, is given by $P = {P_A} + {P_I}$. Then, the chemical master equation corresponding to Eq 1 takes the following form

where $m, n = 0, 1, 2, \cdots$, $S$ is the unit operator, and ${E_l}$ is the common step operator and $E_l^{-1}$ is its inverse (where $l = m$ or $n$). The initial condition is set as ${P_I}\left({0, 0;0} \right) = 1$ and ${P_A}\left({0, 0;0} \right) = 0$. We are interested in stationary eRNA and mRNA distributions. Note that at steady state, Eq 2 becomes

Introduce two factorial marginal distributions for eRNA and mRNA, which are respectively defined as ${Q_i}\left(n \right) = \sum\limits_m {{P_i}\left({m, n} \right)}$ and ${R_i}\left(m \right) = \sum\limits_n {{P_i}\left({m, n} \right)}$, where $i = I, A.$ Thus, the total stationary distributions of eRNA and mRNA are $Q\left(n \right) = {Q_I}\left(n \right) + {Q_A}\left(n \right)$ and $R\left(m \right) = {R_I}\left(m \right) + {R_A}\left(m \right)$, respectively. One main aim of this paper is to derive the analytical expressions of $Q\left(n \right)$ and $R\left(m \right)$.

2.2.   Dependence of looping rate on E-P loop length

Apart from the eRNA's regulation mechanism, the E-P looping rate is also an important factor affecting gene expression. Since an enhancer's position relative to its promoter is not fixed, the E-P looping rate is assumed to change in a finite range. Note that the transition rate from OFF to ON states in Figure 2, ${k_0}$, actually represents the rate that the enhancer and the promoter form a DNA loop (denoted by ${k_{loop}}$), i.e., ${k_0} = {k_{loop}}$. In order to quantify how the E-P loop length (denoted by $d$) impacts gene expression, it is needed to know the dependence of ${k_0}$ on $d.$ For a single DNA loop, previous works studied the impact of the E-P loop length on the looping rate, and even gave empirical formulae ^[47,48]. In our case, the corresponding formula takes the form of ${k_{loop}} = {k_{on}}\exp \left({ - \frac{u}{d} - v\log d + w{\kern 1pt} d + z} \right)$, which is motivated by the worm-like chain model of DNA bending ^[49], where ${k_{loop}}$ represents the rate of the DNA looping, and $d$ represents the E-P loop length along the DNA line, i.e., the distance between the promoter and the enhancer. Accordingly, some parameter values can be set as ${k_{on}} = 1, u = 140.6, v = 2.52, w = 0.0014, z = 19.9$, which are obtained by fitting experimental data ^[47,48]. As mentioned on the Introduction, enhancers located upstream or downstream of the target gene. For target gene, the different position of enhancer will directly affect the E-P looping rate ${k_{loop}}.$In other words, once the target gene and the enhancer are fixed, then the ${k_{loop}}$is determined. However, for different target gene, the position of enhancer is different, so ${k_{loop}}$ is also different. Under this setting, given the value of ${k_0}$, we can determine the DNA loop length (also called distance) $d$from the above empirical formula. However, the case of ${k_1} = {k_{loop}}$ is also possible in theory. Therefore, based on the above empirical formula, we can also determine $d$ if the value of $k{\, _1}$ is given.

As is well known, enhancers' regulation has a long-distance effect. Although enhancers are cis-acting elements, the position relative to their target genes may be highly changeable. Generally, an enhancer is located ~200$bp$ upstream from the target promoter, but can also facilitate the transcription activity of a distant promoter even with the distance being > 10$kb$ ^[3,50]. In the following sections, we will set the E-P loop length between 200$bp$ and 8000$bp$ in our numerical simulations, to obtain qualitative results.

Finally, in this section, we point out: (1) Most eRNA's have been proposed to have very short half-lives (~7.5 mins) and only a few have long half-lives ^[33], implying that the degradation rate of eRNA is in general larger than that of mRNA. In our simulation, we will set ${\delta _1} = 1$ (one unit) and ${\delta _2} = 5$. (2) The amount of eRNAs is usually much less than that of mRNAs ^[51], implying that the production rate of the former is much smaller than that of the latter. In our simulation, we will set ${\mu _1} = 20$ and ${\mu _2} = 2$. (3) The choice of other parameter values is based on the setting of parameter values used often in the study of the common ON-OFF model at the transcription level ^[40].

3.   Results

First, there are two possible ways to determine the DNA loop length ($d$) if the other parameters are fixed: increasing ${k_0}$, and reducing ${k_1}$. Second, we numerically find that in whichever way (i.e., given the value of ${k_0}$ or ${k_1}$, which determines $d$according to the above empirical formula), the qualitative impact of increasing ${k_0}$ or reducing ${k_1}$ on the mRNA expression is the same (referring to an example shown in the Appendix of this paper). Therefore, we will only show results in the case that $d$ is determined by the previous setting of ${k_0}$ that is allowed to change in a finite range (implying that $d$ changes in another finite range).

3.1.   Analytical distributions and statistics

Analytical results can often provide the comprehensive understanding of a system's behavior. Here we present analytical results for eRNA and mRNA distributions, and statistical quantities such as expectations and variances as well as noise intensities. The detailed mathematical derivations are put in Appendix A of this paper.

First, consider the eRNA case. We can show that the exact, analytical eRNA distribution is given by

where all the parameters are normalized by ${\delta _2}$, i.e., ${\tilde \mu _2} = {{{\mu _2}} \mathord{\left/ {\vphantom {{{\mu _2}} {{\delta _2}}}} \right. } {{\delta _2}}}, {\tilde k_0} = {{{k_0}} \mathord{\left/ {\vphantom {{{k_0}} {{\delta _2}}}} \right. } {{\delta _2}}}$, ${\tilde k_1} = {{{k_1}} \mathord{\left/ {\vphantom {{{k_1}} {{\delta _2}}}} \right. } {{\delta _2}}}$ and $\tilde \alpha = {\alpha \mathord{\left/ {\vphantom {\alpha {{\delta _2}}}} \right. } {{\delta _2}}}$, function${}_1{F_1}\left({a, b; z} \right)$. is a confluent hypergeometric function ^[52], defined as ${}_1{F_1}\left({a, b; z} \right) = \sum\limits_{n = 0}^\infty {\frac{{{{\left(a \right)}_n}}}{{{{\left(b \right)}_n}}}\frac{{{z^n}}}{{n!}}}$, and ${(c)_n}$ is the Pochhammer symbol, defined as ${(c)_n} = {{\Gamma \left({c + n} \right)} \mathord{\left/ {\vphantom {{\Gamma \left({c + n} \right)} {\Gamma \left(n \right)}}} \right. } {\Gamma \left(n \right)}}$ with $\Gamma (\cdot)$ being the common Gamma function. Equation 4 indicates that the stationary eRNA distribution is independent of mRNA. Thus, the eRNA noise intensity is given by

where

Then, consider the stationary mRNA distribution. Under some assumed conditions, we can derive the approximate expression of this distribution (seeing Appendix A for derivation)

where $\left\langle n \right\rangle$ is given by Eq 5a. This analytical expression shows how eRNA quantitatively impacts the stationary mRNA distribution. Apparently, this impact is nonlinear. We point out that this approximate distribution has been verified by the Gillespie stochastic simulation algorithm ^[53], referring to Figure 3A where the E-P loop length is set as $1500bp.$ From the panels of this figure, we observe that in all the three cases of eRNA regulation, the absolute differences between theoretical and numerical results do not exceed 0.01, indicating that our approximate method (seeing Appendix A for detail) is effective.

The mean mRNA is given by

whereas the mRNA noise intensity by

Both show how the mean of eRNA molecules affects the mean and noise intensity of mRNA. Moreover, $\left\langle m \right\rangle$ is a monotonically increasing function of $\left\langle n \right\rangle$ whereas $\eta _m^2$ is a monotonically decreasing function of $\left\langle n \right\rangle$. Similar to the case of mRNA distribution, the approximate mean and noise intensity of mRNA have also been verified by the Gillespie stochastic simulation algorithm ^[53], referring to Figure 3B (for the mRNA mean) and 3C (for the mRNA noise intensity). In fact, the maximum error between the two mRNA means does not exceed 0.4, whereas the maximum error between the two mRNA noise intensities does not exceed 0.002. Therefore, except that the error is slightly larger at $d = 2000bp$ (no more than 0.5), the two methods are still in good agreement in other E-P looping length. These results further indicate that our approximate method is reasonable.

3.2.   Impact of E-P loop length on mRNA probability distribution

Experimental evidences have indicated that an enhancer may be located in any position of the promoter's upstream or downstream. The enhancer can not only prompt the transcription of the core promoter but can also facilitate that of any promoter near it. Therefore, it is necessary to investigate how the E-P loop length impacts gene expression and further cell-to-cell variability in isogenic cells. Here we address this issue in terms of mRNA distribution, with reasons stated as follows. First, in contrast to a single stochastic trajectory, the distribution can provide more complete information on stochastic dynamics of the underlying system. Second, distributions provide an intuitive understanding for cellular phenotypes and phenotypic diversity as well as phenotypic switching ^[54,55,56].

In order to show the effect of the E-P loop length on the mRNA distribution, we consider four cases: Type-I regulation, Type-II regulation, mixed regulation and non-regulation. Numerical results are demonstrated in Figure 4, where the lines represent the results obtained by the Gillespie algorithm ^[53] whereas the circles represent the results obtained by theoretical prediction. We observe that with the increase of the E-P loop length, the mRNA distribution can change from unimodality to bimodality and then be back to unimodality in all the cases, independent of eRNA regulatory patterns.

In Figure 4A, the first, second, third and fourth columns correspond respectively to Type-I regulation, mixed regulation, Type-II regulation and non-regulation. In all the cases, the E-P loop length can result in unimodal or bimodal mRNA distribution where one peak is close to zero and the other peak is away from zero. Specifically, the mRNA distribution is unimodal for a small $d$(for example, $d$is set as 200$bp$) in all the four cases of eRNA regulation, referring to the first row in Figure 4A. At least $d$is about 600$bp$ is needed to generate bimodality for Type-I regulation, at least $d$is about 800$bp$ for mixed regulation, at least $d$is about 1000$bp$ for Type-II regulation, and at least $d$is about 1200$bp$for non-regulation, referring to the second row in Figure 4A. With the increase of $d$, the bimodal mRNA distribution gradually become unimodal, but the E-P loop length needed for this change is different among the four cases. Specifically, at least $d$is about 3500$bp$ is needed for Type-I regulation, at least $d$is about 4600$bp$for mixed regulation, at least $d$is about 5500$bp$ for Type-II regulation, and at least $d$is about 7000$bp$ for non-regulation, referring to the fourth row in Figure 4A. The panels from the second row to the third row in Figure 4A show how the E-P loop length impacts the shape of the bimodal mRNA distribution in all the four cases of eRNA regulation. Figure 4B demonstrates the time series of the mRNA number in the four cases of eRNA regulation: from top to bottom, the fifth row for bimodality and the sixth row for unimodality. In order to help understand the mechanism of how the E-P loop length induces the transition between ON and OFF states in all the four cases of eRNA regulation, we also plot the time evolution of promoter activity, referring to Figure 4C.

In a word, we numerically find that the shortest E-P loop length is needed to generate bimodality or to switch from bimodality to unimodality in the case of Type-I regulation, but the longest $d$ is needed in the case of non-regulation. In the other two regulation cases, the $d$ needed for bimodality or switching is in between the shortest and longest E-P loop lengths. The results shown in Figure 4 would have important biological implication. Perhaps cells would make use of the E-P loop for better survival in complex environments.

Finally in this subsection, we give a simple interpretation for the appearance of bimodality. First, while the common on-off model of gene expression itself has the potential to generate bimodality, our model has the similar one ^[40,41]. Second, eRNA in our gene model is designed to regulate the expression of the gene in an activating manner. Due to the low production rate and high degradation rate of eRNA, the number of eRNA is generally small, but this activating effect on gene expression indeed exists and is unneglectable, and even can induce phenotypic switching, as demonstrated in Figure 4. In addition, we point out that previous studies showed that a simple genetic auto-activation loop can lead to bistability (meaning that the underlying deterministic system has two stable equilibria) but can also lead to bimodality in the stochastic case ^[57]. In contrast, we showed that the E-P loop can induce bimodality. Therefore, our mechanism is different from the previous one.

3.3.   Effects of the E-P loop length on the mRNA mean

Figure 5 shows how the E-P loop length affects the mean mRNA level. Specifically, Figure 5A demonstrates the dependence of the mRNA mean on the E-P loop length in all the four cases of eRNA regulation, Figure 5B demonstrates the dependence of the ratio of the mRNA mean in all the three cases of eRNA regulation on the E-P loop length, Figure 5C demonstrates the impact of the E-P loop length on transition rate ${k_0}$, and Figure 5D demonstrates the dependence of mRNA mean on both E-P loop length and feedback strength.

From Figure 5A, we observe that the curve in the case of non-regulation is always below the curves in three cases of eRNA regulation for all possible E-P loop lengths that are assumed to change in the interval from 200$bp$ to 8000$bp$ in our simulation, and the curve in the case of Type-II regulation is above that in the case of mixed regulation, which is above that in the case of Type-I regulation. Moreover, the changing range of the mean mRNA level is larger in the case of Type-II regulation than that in the cases of Type-I and mixed regulations, indicating that the Type-II regulation can amplify the mRNA expression level in contrast to the other two ways of eRNA regulation (i.e., Type-I and mixed regulations). We also observe that the mRNA mean first decreases for smaller E-P loop lengths and then increases with the increase of the E-P loop length, that is, the mean mRNA level reaches the minimum in all the four cases of eRNA regulation.

Figure 5B demonstrates the dependence of the mean mRNA ratio (defined as the ratio of the mean mRNA in the case of eRNA regulation over the mean mRNA in the case of non-regulation) on the E-P loop length for three ways of eRNA regulation, where by "relative" we mean that the mean mRNA level in the case of non-regulation is set as unit (or the mean mRNA level in the case of non-regulation is taken as a reference). We observe that these relative changes exhibit different characteristics, depending on the way of eRNA regulation. Specifically, the mean mRNA ratio reaches the minimum in the case of Type-II regulation but the maximum in both cases of Type-I and mixed regulations. In addition, for all possible E-P loop lengths, the relative curve for Type-II regulation is always higher than the ones for two other regulations, and the change trend of the former is opposite to that of the latter two.

In order to verify the above qualitative results, we further analyze the dependence of the mean mRNA level on both E-P loop length and feedback strength, referring to Figure 5D. Interestingly, we observe that for any feedback strength, the mean mRNA level can reach a minimum at about 2000$bp$ E-P loop length. Particularly in the cases of Type-II and mixed regulations, the larger feedback strength leads to a higher mean mRNA level, but the E-P loop length needs to be limited in the interval from 200 to 4000$bp$.

It is worth pointing out that the above qualitatively different relative change tendencies may in turn help us infer the way and the qualitative effect of eRNA regulation.

3.4.   Effects of the E-P loop length on the mRNA noise

Figure 6 shows how the E-P loop length affects the mRNA noise. Specifically, Figure 6A demonstrates dependence of the mRNA noise on the E-P loop length in the four cases of eRNA regulation as specified above; Figure 6B, an enlarged diagram for the red dotted line shown in Figure 6A, demonstrates the dependence of the mRNA noise on the E-P loop length in the case of Type-I regulation; Figure 6C demonstrates the dependence of the ratio of the mRNA noise (defined as the ratio of the mRNA noise in the case of eRNA regulation over the mRNA noise in the case of non-regulation) on the E-P loop length in the three cases of eRNA regulation; and Figure 6D demonstrates the dependence of mRNA noise on both E-P loop length and feedback strength.

From Figure 6A, we first observe that the curve in the case of non-regulation is above the curve in the case of Type-I regulation, which is above the curve in the case of mixed regulation, which is above the curve in the case of Type-II regulation. This implies that each of three eRNA regulation ways can reduce the mRNA noise in contrast to non-regulation. In addition, whatever the E-P loop length, the mRNA noise in the case of Type-II regulation is the lowest, the mRNA noise in the case of Type-I regulation is the highest, and the mRNA noise in the case of mixed regulation is moderate. Then, we observe that in all the three cases of eRNA regulation, the mRNA noise intensity first increases and then decreases with the increase of the E-P loop length, implying that the mRNA noise has the highest level even for Type-I regulation (referring to Figure 6B).

From Figure 6C, we observe that the noise-ratio curves for Type-I regulation (red) and mixed regulation (purple) have the same change trend, that is, they first decrease and then increase until they become flat. The noise ratios reach the minimum around the 2000$bp$E-P loop length. If the E-P loop is long enough (e.g., more than 6000$bp$), the noise level in the case of Type-I regulation is very close to that in the case of non-regulation, and the former is much smaller than the latter otherwise. In addition, the curve for Type-II regulation is always below the ones for the other two regulation ways, and the change trend of the former is opposite to that of the latter two. Furthermore, we analyze the simultaneous dependence of mRNA noise intensity on E-P loop length and feedback strength, referring to Figure 6D. From this diagram, we observe that for any feedback strength, the mRNA noise intensity reaches the maximum around the 2000$bp$ E-P loop length. Particularly in the cases of Type-I and mixed regulations, the greater the feedback strength is, the smaller is the mRNA noise, and this relationship holds for the E-P loop length in the interval from 200 to 4000$bp.$

The combination of Figures 5 and 6 implies that the Type-II regulation is advantageous over Type-I or mixed regulation in controlling gene expression and tuning gene expression noise, and this advantage is independent of the E-P loop length.

3.5.   Effects of feedback strength and eRNA synthesis rate on mRNA noise

Here, we investigate the effects of feedback strength parameters ($\alpha$and $\beta$) and eRNA production rate (${\mu _2}$) on the mRNA noise, for a fixed E-P loop length. Numerical results are demonstrated in Figure 7.

From Figure 7A and B, we observe that the mRNA noise is a monotonously decreasing function of feedback strength $\alpha$ or $\beta$, where red dotted lines correspond to Type-I regulation, blue circled lines to Type-II regulation, and purple cross lines to mixed regulation. Figure 7A shows that with the increase of $\alpha$, the mRNA noise for Type-I regulation is a bit greater than that for mixed regulation. Similarly, Figure 7B shows that with the increase of $\beta$, the mRNA noise for Type-II regulation is apparently greater than that for mixed regulation. In addition, the mRNA noise in Figure 7A is a concave function of $\alpha$, whereas that in Figure 7B has an inflecting point nearly at $\beta = 0.45$: the mRNA noise is a concave function if $\beta < 0.45$ but a convex function if $\beta > 0.45$.

Figure 7C demonstrates how the mRNA noise intensity changes with the eRNA synthesis rate ${\mu _2}$. We observe that in all the three cases of eRNA regulation, with the increasing of synthesis rate ${\mu _2}$, the mRNA noise first decreases rapidly and then maintains at a stable level after the synthesis rate is beyond a certain value. In the three cases of eRNA regulation, although the change trend of the mRNA noise is fundamentally similar, the mRNA noise level displays an apparent difference. Specifically, if ${\mu _2}$ is kept the same, the line for Type-I regulation is above that for mixed regulation, which is above that for Type-II regulation. This indicates that if the eRNA generation rate is set as unit, the mRNA noise for Type-II regulation is the lowest, and the mRNA noise for Type-I regulation is higher than that for mixed regulation. Therefore, from the viewpoint of reducing the mRNA noise, the best way of eRNA regulation is the Type-II regulation, which is in accordance with the conclusion for the mean mRNA.

4.   Conclusions

More and more experimental evidences indicate that non-coding RNAs play important role in controlling gene expression. As a kind of mysterious non-coding RNA and a new regulation factor, eRNA may affect gene expression and its anomalous expression may lead to many diseases ^{[58,59,60,61,34]}. Exploratory interrogation should be conducted to unravel the association of functional eRNAs with specific chromatin regions and other func­tional RNAs ^[62], as well as to understand possible functional motifs, secondary structures and chemical modifications ^[63,64]. On the other hand, mathematical models are a powerful tool for exploring possible functions of eRNAs.

In this work, based on biologically experimental evidence, we have proposed a mechanistic model of stochastic gene expression with eRNA regulation. In this toy model, the eRNA regulatory role is primarily reflected in one of three ways: (1) Type-I regulation by which eRNA accelerates the formation of the E-P loop, (2) Type-II regulation by which eRNA promotes the production rate of its target gene product, (3) mixed regulation which is a combination of the former two regulation ways. We have also provided a theoretical analysis framework to investigate the eRNA's contribution to the target gene expression, and derived analytical formulae for the mRNA's distribution, mean and noise. Interestingly, we have found that with the increase of the E-P loop length, the mRNA distributions can transition to bimodality from unimodality and then be back to unimodality, independent of eRNA regulation ways. We have also shown that in contrast to Type-I and mixed regulations, the Type-II regulation can result in the higher mean mRNA level and the lower mRNA noise, independent of the E-P loop length. These results are independent of the choice of the model's parameter values and are hence qualitative, although we only demonstrated relevant numerical results for particular sets of parameter values. Our results indicate that the E-P loop can induce phenotypic switching, a mechanism different from the traditional one (i.e., feedback can induce phenotypic switching). Both would be exploited by live cells surviving in complex environments.

In our model, eRNA is taken as an enhancer's product, which inter-regulates the expression of its target gene. The made assumption may not be exactly consistent with biological reality of eRNA regulation under some specific environments. In other situations, it would be needed to construct other more biologically reasonable models to investigate the role of eRNA in regulating its target genes, e.g., the constructed models should consider dynamic regulation where eRNA is taken as an input signal of the model systems ^[65]. For this, one may assume that eRNA is a special spatiotemporal regulation factor that wanders in cellular space, and its concentration is a time-dependent function, eRNA(t). Based on this assumption, one can also establish theoretical models as done in this paper. This will be the next task in our future research.

It should be pointed out that in the introduction, we have mentioned three hypotheses on the "functional theory": eRNA facilitates the formation of the E-P loop, eRNA remodels chromatin, and eRNA improves the extension of RNAPII. In this paper, we assumed that the most possible biological function of eRNA is to enhance the mRNA production rate. Although this rate integrates the effects of chromatin remodeling and promoting RNAPII elongation, but we did not specify molecular mechanisms behind these effects. Therefore, the future study needs to construct a more realistic model to investigate whether chromatin remodeling or the extension of RNAIIP can increase the mRNA transcriptional rate, and how this rate affects gene expression and further cell-to-cell variability across a genetically identical population of cells.

Acknowledgments

This paper was supported by grant 11775314 and 11931019 (TS), 11801482(HL), 11701591(TQ) from Natural Science Foundation of P. R. China; grant 202007030004 from Key-Area Research and Development Program of Guangzhou.

Conflict of interest

The authors declare that there is no conflict of interest.

Appendix A: Derivation of analytical results

Let ${Q_i}\left(n \right) = \sum\limits_m {{P_i}\left({m, n} \right)\quad } (i = I, A)$ be the factorial marginal probability density function of eRNA. It follows from Eq 3 in the main text that

This is a solvable model where mRNA does not impact eRNA. To solve this equation group, we introduce probability-generating functions ${G_i}\left(z \right) = \sum\limits_n {{{\left({z + 1} \right)}^n}{Q_i}\left(n \right)}, \; i = I, A,$ and $G\left(z \right) = {G_I}\left(z \right) +$${G_A}\left(z \right)$. Then

where the parameters are normalized by ${\delta _2}$. By adding the above two equations, we get ${G_A}\left(z \right) = \frac{1}{{{{\tilde \mu }_2}}}\partial G\left(z \right)$ and ${G_I}\left(z \right) = G\left(z \right) - {G_A}\left(z \right) = G\left(z \right) - \frac{1}{{{{\tilde \mu }_2}}}\partial G\left(z \right)$. From the first equation of (A2), we have

That is,

Setting $x = z + \frac{{\tilde \alpha }}{{\tilde \alpha + 1}}$, we have

Making Taylor expansion: $G\left(x \right) = \sum\limits_n {{a_n}{x^n}}$, yields the following iterative equation

which gives

Thus, we obtain the following expression of function $G\left(x \right)$

Being back to the original variable yields

Note that $G\left(0 \right) = 1$, which can give

Using the relationship between probability distribution and probability-generating function, we thus obtain the analytical expression for the marginal probability density function of eRNA

Furthermore, the mean of eRNA is given by

Next, let ${R_i}\left(m \right) = \sum\limits_n {{P_i}\left({m, n} \right)}$ be the marginal probability density functions of mRNA. It follows from Eq 3 in the main text that

Now, we make the approximation: replacing ${P_I}(m, n)$ with ${R_I}(m){Q_I}(n)$. In addition, we make the approximation: replacing $\sum\limits_{n = 0}^\infty n {\kern 1pt} {Q_I}(n)$ with $\sum\limits_{n = 0}^\infty {n{\kern 1pt} {\kern 1pt} } Q(n).$ The effectiveness of these approximations has been verified by numerical simulation (seeing the main text). Then, we have the approximate equations

where $\left\langle n \right\rangle$ is given by Eq A12. Since $\left\langle n \right\rangle$ is unrelated to mRNA, it can be taken as a parameter in Eq A14. The other parameters in Eq A14 have been normalized by ${\delta _1}$. Thus, Eq A14 is also a solvable model. In fact, it is a common ON-OFF model where the transition rate from ON to OFF state is ${\tilde k_1}$, the transition rate from OFF to ON state is ${\tilde k_0} + \tilde \alpha \left\langle n \right\rangle$ and mRNA transcriptional rate is ${\tilde \mu _1} + \tilde \beta \left\langle n \right\rangle$. By solving this model, we obtain the following analytical expression of stationary mRNA probability distribution

where ${(c)_n}$ is the Pochhammer symbol, defined as ${(c)_n} = {{\Gamma (c + n)} \mathord{\left/ {\vphantom {{\Gamma (c + n)} {\Gamma (n)}}} \right. } {\Gamma (n)}},$ and ${}_1{F_1}\left({a, b; z} \right)$ is a confluent hypergeometric function, defined as ${}_1{F_1}\left({a, b; z} \right) = \sum\limits_{n = 0}^\infty {\frac{{{{\left(a \right)}_n}}}{{{{\left(b \right)}_n}}}\frac{{{z^n}}}{{n!}}}$.

Appendix B: Supplemental figure

In the main text, we have investigated the effect of eRNA on the mRNA expression in the case that eRNA regulates the transition rate from OFF to ON. However, it is possible that eRNA regulates the transition rate from ON to OFF. In this case, whether the qualitative effect of eRNA affects on the mRNA expression is that in the case that eRNA regulates the transition rate from OFF to ON. We numerically find that both are indeed the same. Here we only demonstrate an example, referring to Figure S1.