article: According to recent experimental evidence, the interaction between chromatin loops, which can be characterized by three factors—connection pattern, distance between regulatory elements, and communication form, play an important role in determining the level of cell-to-cell variability in gene expression. These quantitative experiments call for a corresponding modeling effect that addresses the question of how changes in these factors affect variability at the expression level in a systematic rather than case-by-case fashion. Here we make such an effort, based on a mechanic model that maps three fundamental patterns for two interacting DNA loops into a 4–state model of stochastic transcription. We first show that in contrast to side-by-side loops, nested loops enhance mRNA expression and reduce expression noise whereas alternating loops have just opposite effects. Then, we compare effects of facilitated tracking and direct looping on gene expression. We find that the former performs better than the latter in controlling mean expression and in tuning expression noise, but this control or tuning is distance–dependent, remarkable for moderate loop lengths, and there is a limit loop length such that the difference in effect between two communication forms almost disappears. Our analysis and results justify the facilitated chromatin–looping hypothesis.
  Regulatory elements and their interactions play a critical role in the spatial, temporal, and physiological control of gene expression [1,2]. These cis-acting elements can be divided into two distinct families: the one is composed of both the promoter and the regulatory elements with the distance from the transcription start site being less than 1kb [3,4], and the other contains the regulatory elements with the distance of much more than 1kb, which may be enhancers, silencers, insulators, or locus control regions [3]. Complexity of transcription activation is mainly in the ordered interactions among regulatory elements of the underlying gene, also including formation of DNA loops and interactions between DNA loops. These interactions are in essential biochemical, leading to stochastic gene transcription and further cell-to-cell variability in gene expression. Such stochasticity is essential for many cellular functions [5,6], and has been identified as a key factor underlying the observed phenotypic variability of genetically identical cells in homogeneous environments [7]. In general, transcription of genes is regulated by promoter–proximal DNA elements and distal DNA elements that together determine expression patterns of these genes. Remarkably in eukaryotic genomes, enhancers can be several hundreds of kilo–bases away from the promoter they regulate [8–10], and the intervening DNA can contain other promoters and other enhancers [11–14]. Enhancers activate promoters by directly contacting binding sites for transcription factors (TFs) and chromatin–modifying enzymes via DNA looping [15–19]. In theory, specific interactions between DNA elements can either assist enhancer–promoter looping by bringing the enhancer and the promoter closer together or are thought to interfere with enhancer–promoter looping by placing them in separate loop domains [20]. Moreover, these DNA elements may form chromatin loops in different manners, e. g. , promoter–tethering elements in Drosophila that allow activation by specific enhancers over long distances are proposed to form DNA loops between sequences near the enhancer and the promoter [21,22]. Other examples include that in bacteriophage λ, the CI protein forms a 2. 3-kb DNA loop that brings a distal stimulatory site close to RNA polymerase at the PRM promoter [23], and in the mouse β–globin locus, the Ldb1 protein appears to form a bridge necessary for efficient enhancer–promoter contact by binding to proteins at the locus control region and at the promoter [24]. To understand how a chromatin loop forms in a nucleus or how it is involved in gene regulation, Li, et al. [2] proposed that the major feature determining loop formation is the flexibility of chromatin, which may be modulated by histone acetylation and other modifications. Single DNA looping and its functions have been intensively and extensively studied. For example, it has been shown that DNA looping can maintain a stable lysogenic state in the face of a range of challenges including noisy gene expression, nonspecific DNA binding, and operator site mutations [25], speed up the search process by bypassing proteins that block the sliding track close to the target [26], enhance lysogenic CI transcription in phage lambda [27], increase the range of bistability in a stochastic model of the lac genetic switch [28], and enhance or suppress transcriptional noise depending on conditions [29]. In addition, Boedicker, eta al. provided a quantitative characterization of the way that critical regulatory parameters modulate the output of transcriptional circuits involving DNA looping [30], Vilar, et al. , showed that regulation based on DNA looping, in addition to increasing the repression level, can reduce the fluctuations of transcription and, at the same time, decrease the sensitivity to changes in the number of regulatory proteins [31], and Choudhary, et al. , found that DNA loop formation is so fast that small bursts are averaged out, making it impossible to extract their size and frequency from the data [32]. In a word, previous studies that are restricted to analysis of single looping treated the looping mechanism as a local role, but according to the recent experimental evidence that reveals DNA–looping interactions, we currently recognize that the looping role is global (e. g. , one DNA loop may interact with another DNA loop) [33]. The classic models of gene expression assume direct spatial contact between a distal enhancer and a proximal promoter [34–39]. However, recent chromatin capture technologies such as Chromosome Conformation Capture [40–41] have shown that enhancers and promoters are connected in a highly complex network of DNA–looping interactions [13,42,43]. For example, Wouter and his colleagues [44,45] showed that long–range (over tens or even hundreds of kilobases) gene regulation in eukaryotic cells involves spatial interactions between transcriptional elements, with intervening chromatin looping out. More interestingly, Priest, et al. [20], used two well–characterized DNA–looping proteins: Lac repressor and phage λ CI, to measure interactions between pairs of long DNA loops in E. coli cells in three possible topological arrangements of two pairs of interacting sites on DNA, namely, side-by-side loops, nested loops, and alternating loops. They found that the first loop pair does not affect each other; the second pair assists each other’s formation consistent with the distance–shortening effect; and the third pair, where one looping element is placed within the other DNA loop, inhibits each other’s formation, thus providing clear support for the loop domain model for insulation. They also argued that combining loop assistance and loop interference can provide strong specificity in long–range interactions. Another related yet important work is that Savitskaya, et al. , experimentally observed [46] that when a pair of repressors and their binding sites are in between the enhancer and the promoter, the gene expression level is not lowered but is raised. For this counter–intuitive phenomenon, Mirny, et al. , conjectured that this repressor pair shortens the distance between the enhancer and promoter pair, leading to the rise of the expression level [33]. Another interesting conjecture is on communication between enhancers and promoters. The related questions are what form the enhancer–promoter communication takes (for example, what are the molecules that are transmitted between regulatory element and promoter), when this takes place, and whether this is the same for all classes of enhancers. This conjecture was put forward as early as in 1988 [47] and was later specified by Blackwood and Kadonaga [48] but has not been solved until now. Studies on specific loci, together with genome–wide approaches, suggest that there may be many common mechanisms involved in enhancer–promoter communication [1,49] (particularly see [49] wherein the author summarized 4 different kinds of communication models). In recent years, two views on communication between enhancer and promoter, namely the direct looping model and the facilitated–tracking model (referring to Fig 1), have become the main stream or popular. The former model assumes a direct interaction between two chromosomal regions by looping out the intervening DNA sequence. For such communication mechanisms, various proteins have been proposed to bridge enhancers and promoters together as looped chromatin structures [50–52], and enhancer RNAs have also been proposed to be physically involved in establishing enhancer–promoter ‘looping’, and involving the integrator [53]. The latter model assumes that enhancer–bound proteins move progressively in a unidirectional manner towards the promoter, sometimes without leaving the enhancer sequence, thus resulting in the formation of a progressive loop that increases its size until it reaches the promoter to form a stable conformation. By modeling and analysis, we will justify the facilitated chromatin–looping hypothesis [1,2]. In a word, experiments have provided evidence for the influence of loop connection pattern, loop distance and communication form on gene expression, but how these factors impact cell-to-cell variability in gene expression remains not fully understood. We address this issue by developing on a stochastic model that maps connection patterns between interacting DNA loops into a multistate model of gene expression. In the case that communication form is not considered, we show that nested loops promote gene expression but reduce expression noise whereas alternating loops reduce the mean expression level but enlarge the noise, compared to side-by-side loops. In the case that communication form is considered, however, we find that the facilitated–tracking mechanism performs better than the direct looping mechanism in enhancing (or reducing) expression and reducing (or enlarging) noise, depending on connection patterns. Moreover, the effects of controlling expression and tuning noise are distance–dependent, remarkable for moderate loop lengths. In addition, we find that there is a limit loop length such that the difference in effect between two communication forms almost disappears. Our results imply that living organisms or cells would use the facilitated–tracking mechanism to deal with information about their environments.
  A pair of insulators Su and Hw found in the gypsy retrotransposon are the most potent enhancer blockers in the Drosophila melanogaster, but they do not prevent the enhancer–promoter communication, mainly because of their pairing interaction that results in mutual neutralization. Savitskaya, et al. [46], experimentally showed that long–distance functional interactions between Su and Hw can regulate communication between the enhancer and the promoter. Specifically, this insulator pair can interact with each other over considerable distances, across interposed enhancers or promoters and coding sequences, whereby enhancer blocking may be attenuated, cancelled, or restored. They also specified the role of the distance between these insulators in blocking enhancer–promoter communication. Recall that Priest, et al. [20] measured interactions between pairs of long DNA loops in E. coli cells in three possible topological arrangements, using a pair of DNA–looping proteins (i. e. , Lac repressor and phage λ CI). Thus, we can accordingly assume that the Su-Hw loop (called as the green loop in this paper) and the enhancer–promoter loop (the blue loop) interact with each other also in three connection patterns: alternating loops (called as the cross–type structure), nested loops (the inline–type structure), and side-by-side loops (the independence–type structure), although other connection patterns are possible in cases that space factors associated with DNA looping are considered. Furthermore, we assume that gene expression is enhanced if and only if the enhancer and the promoter form a loop, although it is possible that the enhancer–promoter loop represses gene expression. In addition, we assume that the fundamental production rate of mRNA is zero. In this paper, we do not consider the regulatory roles of TFs although they may influence chromatin looping and enhancer–promoter communication [1]. Apart from connection patterns between interacting DNA loops, the rates of chromatin looping are also important factors affecting gene expression. In general, the looping rate for a pair of regulatory elements (e. g. , Su and Hw) depends on not only the distance between this element pair but also the distances between other pairs of regulatory elements. In our case, related experimental results that support our model settings are stated as follows. In the alternating structure, the Su and Hw forms a blocking structure, which interferes with (referring to the red X in Fig 2) the enhancer–promoter looping, thus reducing the looping rate and further decreasing gene expression. For example, if λ is the looping rate of enhancer and promoter, then this rate can be reduced to 0. 3λ after its DNA loop is repressed by the other loop [20,46]. In the nested structure, the formation of the Su and Hw loop decreases the length of the enhancer–promoter loop, thus increasing the enhancer–promoter looping rate. For example, if λ is the looping rate of enhancer and promoter, then this rate can be increased to 8λ after its DNA loop is enhanced by the other loop [20,46]. In the side-by-side structure, the Su-Hw looping and the enhancer–promoter looping do not affect each other, so each has its looping rate. See Subsection 2. 3 for more details. Finally, although communication form between regulatory elements may be diverse, three representative mechanisms [1] have been proposed based on experimental evidence: the first mechanism is the linking model in which an activator protein first binds the promoter at a proximal sequence and facilitates the recruitment of a second TF at a site located just downstream of the former [54–56]. The second mechanism is the tracking model and/or a facilitated–tracking model in which histone acetylation and TF complexes are transiently detected in the intervening sequence and precedes transcription [52,57–59]. The third mechanism is the looping model that suggests a direct interaction between two chromosomal regions by looping out the intervening DNA sequence [53,60]. Currently, the latter two views have been extensively received, but it is unclear which communication mechanism is more possibly used by living organisms or cells. We will address this issue. First, the green loop formed by a pair of insulators (Su and Hw) and the blue loop formed by a pair of expression elements (enhancer and promoter) interact with each other in the one of three possible connection patterns: cross–type structure (due to alternating loops), inline–type structure (due to nested loops), and independence–type structure (due to side-by-side loops). Refer to the first column of schematic Fig 2. For convenience, we denote by d1 the length of the blue loop along the DNA line, and by d2 the length of the green loop along the underlying DNA line. Experimental evidence supports that alternating loops give loop interference, nested loops give loop assistance, and side-by-side loops do not interact [20,46]. In spite of this, an unexplored question is how these interacting DNA loops affect gene expression (including mean outcome and expression noise). In order to address this issue, we first introduce physical models for three fundamental patterns (alternating loops, nested loops and side-by-side loops), referring to the second column of Fig 2. Note that in theory, the enhancer and promoter pair may form a loop but also may not form any loop, i. e. , there are two possibilities. Similar cases are for the Su and Hw pair. Thus, there are in total four possibilities for each of three patterns. In order to help understanding, some details are listed below. It is needed to point out that to derive the physical models shown in Fig 2, we make simplifications, e. g. , we do not consider binding of large protein complexes to promoters and enhancers, sequential recruitment of TFs to enhancers and promoters, and histone acetylation (and other modifications). All these factors or other complex processes would affect formation of DNA loops and interactions between them. Then, we further map three physical models into a common multistate model of gene expression at the transcription level (the third column in Fig 2). This mapping can bring us great conveniences for analysis and simulation. After mapping, the looping rates, which are functions of loop lengths, currently become transition rates between promoter activity states. Note that once two DNA loops are formed, any one of them can affect the length of the other, often in a nonlinear manner. This impact can lead to changes in transition rates and further in the gene outcome. Also note that the ON1 or ON2 state indicated in the theoretical model of Fig 2 means that the enhancer and the promoter form a loop (i. e. , the blue loop) whereas the Su and the Hw may form a loop (i. e. , the green loop) but also may not form any loop. In contrast, the OFF1 or OFF2 state means that the enhancer and the promoter do not form a loop whereas the Su and Hw pair may form a loop (i. e. , the green loop) but also may not form any loop. It is worth pointing out that the proposed–above mapping method is easily extended to cases of complex interactions among arbitrarily many pairs of chromatin loops. Note that looping rates are actually transition rates between on and off states shown in the third column of Fig 2 according to the above mapping relationships. In order to quantify how two interacting DNA loops affect the level of cell-to-cell variability in gene expression, it is needed to know the dependence of transition rates between active and inactive states in the mapped gene model on the lengths of these two loops (along the DNA) since different loop lengths would lead to different expression levels. For cases of single DNA loops, previous works studied the influence of the loop length on the looping rate, and even gave some experiential formulae between them [37,61]. In our case, these formulae read λ14=kloop (1) =kRonexp (−ud1−vlog (d1) +wd1+z) λ12=kloop (2) =kRonexp (−ud2−vlog (d2) +wd2+z) (1) where kloop represents the rate of DNA looping, ad d is the loop length along the DNA line (i. e. , the inter–operator distance). Some parameter values are set as kRon=1, u = 140. 6, v = 2. 52, w = 0. 0014, and z = 19. 9, which are obtained by fitting experimental data [37,61]. In general, for two interacting DNA loops, the rate of each DNA looping depends on not only its own loop length but also the length of the other loop, and is therefore a function of two variables. Note that parameter λ23 represents the looping rate of the blue loop after the green loop has been formed, implying that it is influenced by the former looping more than by the latter looping rate. Also note that the looping rate of the green loop relies on its own length d2. Therefore, λ23 is in principle a function of d1 and d2. Similarly, λ43 is also a function of d1 and d2. However, the existing experimental data only supported the quantitative relationship between the looping rate and the length of the blue loop [33]. Based on the above analysis and without loss of generality, we may set λ23=k1λ14, λ43=k2λ12 (2) where the size of parameter k can determine the connection patterns of interacting DNA loops, according to Refs. [20,33]. In particular, according to experimental data [33] with small modifications, we may set k={4e−0. 5d+1>1nested loops1side-by-side loops0. 5<1alternating loops (3) where d represents the length of the underlying loop. First, this setting is reasonable since the formation of one loop in the nested, side-by-side, or alternating pattern enlarges, does not change, or reduces the length of the other loop, respectively. Second, the setting is based on observations from Fig 3A in Ref. [33] that for the nested pattern, k depends on the DNA loop distance in an exponential manner, whereas for the alternating pattern, the effect of reducing the loop length is remarkable when the blue loop distance is small but almost disappears when this distance is large. Next, we consider the setting of λ14 and λ12 in the case that the so–called tracking (more precisely facilitated–tracking) mechanism between looping elements is considered. To help the reader understand this mechanism, we imagine a DNA loop as a string with some fixed length, two ends of which represent looping elements (e. g. , Su and Hw). If one element slides along this string (here we only consider the sliding of one element in the green loop since we have assumed that the blue loop promotes gene expression), then this will affect the range that the enhancer and the promoter form the blue loop. Thus, the tracking mechanism leads to the increase in looping rates. Specifically, if we denote respectively by λ˜14 and λ˜12 the looping rates of the blue and green loops in the case that the facilitated–tracking mechanism is considered, and respectively by λ14 and λ12 the natural looping rates of these two loops, then we have λ˜14=λ14+Δ1, λ˜12=λ12+Δ2 (4) (correspondingly, λ23 and λ34 need to be modified). Note that for the facilitated–tracking mechanism, the longer the DNA loop length is, the more is the range that one regulatory element tracks another regulatory element, implying Δ∼d. Thus, it is reasonable to set Δ = rd, where r is a nonnegative parameter. Also note that no tracking or direct looping corresponds to r = 0 whereas tracking corresponds to r ≠ 0, so r characterizes the difference between the two mechanisms. Therefore, we call this parameter as the tracking ratio, which can be also understood as the probability that the enhancer and the promoter track to each other along the DNA line. In order to study qualitative and quantitative effects of the above three pairs of interacting DNA loops on gene expression (including the mean mRNA expression level and the mRNA noise intensity), we establish a mathematical model for the schematic diagram shown in the third column of Fig 2. Assume that the gene has transcription only at ON1 to ON2 states with transcription rates denoted respectively by μ1 to μ2, and the mRNA degrades in a linear manner with the degradation rate denoted by δ. Let λ14 and λ41 stand respectively for transition rates from OFF1 to ON1 and vice versa, λ12 and λ21 respectively for transition rates from OFF1 to OFF2 and vice versa, λ23 and λ32 respectively for transition rates from OFF2 to ON2 and vice versa, and λ34 and λ43 respectively for transition rates from ON1 to ON2 and vice versa. Note that these transition rates are all functions of DNA loop lengths (along DNA lines). Denote by P1 (m; t), P2 (m; t), P3 (m; t) and P4 (m; t) the probabilities that the mRNA has m copies at OFF1, OFF2, ON1 and ON2 states respectively and at time t. Then, the chemical master equation for the full reaction system takes the following form ∂P1 (m; t) ∂t=λ21P2 (m; t) +λ41P3 (m; t) − (λ14+λ12) P1 (m; t) +δ (E−I) [mP1 (m; t) ]∂P2 (m; t) ∂t=λ12P1 (m; t) +λ32P4 (m; t) − (λ21+λ23) P2 (m; t) +δ (E−I) [mP2 (m; t) ]∂P3 (m; t) ∂t=λ14P1 (m; t) +λ34P4 (m; t) − (λ41+λ43) P3 (m; t) +μ1 (E−1−I) [P3 (m; t) ]+δ (E−I) [mP3 (m; t) ]∂P4 (m; t) ∂t=λ23P2 (m; t) +λ43P3 (m; t) − (λ34+λ32) P4 (m; t) +μ2 (E−1−I) [P4 (m; t) ]+δ (E−I) [mP4 (m; t) ] (5) where E with the reverse E−1 is the step operator and I is the unit operator. In Eq (5), the quantitative dependences of all λ on loop lengths have been given in the previous subsection. It is worth pointing out that if more complex connection patterns for interactions between DNA loops are considered, then the corresponding chemical master equation can be generalized as [62] dP (m; t) dt=AP (m; t) +Λ (E−1−I) [P (m; t) ]+δ (E−I) [mP (m; t) ] (6) where E and E−1 are vectors of step operators and I is the vector of identity operators. In Eq (6), matrix A describes the interactions between chromatin loops, with the entries representing looping rates that are functions of loop lengths, diagonal matrix Λ describes the exits of transcription, with the diagonal elements representing transcription rates, and diagonal matrix δ depicts the degradation of mRNA, with the diagonal elements representing degradation rates. It is in general difficult to solve a chemical master equation. However, the binomial moment approach that we developed previously [39] can well solve this question since it can conveniently give numerical solutions and even can give analytical solutions in some cases. In the following, we simply introduce this approach for the general gene model described by Eq (6). For this, we first introduce factorial binomial moments for factorial distributions Pi (m; t), which are defined as bk (i) =∑m=k∞ (mk) Pi (m; t), k=0,1, 2, ⋯ (7) Then, we can derive the following linear ordinary differential equations with constant coefficients from Eq (6) ddtbk (t) =Abk (t) +Λbk−1 (t) −δkbk (t) (8) Our interest is in finding steady–state distribution. For convenience but without loss of generality, we assume that all the degradation rates are the same and the common degradation rate is 1. At steady state, factorial binomial moments satisfy the following algebraic equation (−A+kI) bk=Λbk−1 (9) from which we obtain bk=1∏i=1kdet (iI−A) ∏i=k1[uN (iI−A) *Λ]b0, k=1,2, ⋯ (10) where bk=∑i=1Nbk (i) represents the total binomial moment with k called the order of this binomial moment, uN = (1,1, ⋯, 1) is an N–dimensional row vector, and (iI − A) * and det (iI − A) are the adjacency matrix and the determinant of matrix (iI − A), respectively. In Eq (10), the column vector b0 is analytically given in S1 Text. Furthermore, we can use the following formula P (m) =∑k≥m (−1) m−k (km) bk, m=0,1, 2, ⋯ (11) to reconstruct the steady–state mRNA distribution, where P (m) =∑i=1NPi (m) represents the total distribution. In particular, we can conveniently use binomial moments to express the noise intensity defined as the ratio of variance over the square of mean. In fact, if this intensity is denoted by ηm, then we have ηm=2b2+b1−b12b12 (12) It is worth pointing out that formulae (11) and (12) are exact and can be directly used to numerical calculation and even derivation of analytical results. Note that the higher the orders of binomial moments are, the exacter are the resulting distributions calculated according to Eq (11) (S1 Text). In addition, the above formulae hold still in the dynamic case. In fact, we show numerical results on time evolutions and distributions of the mRNA number in Figs C–H in S1 Text, implying that our binomial moment approach can be also conveniently used to analysis of stochastic dynamics.
  Here, we apply the above analysis to the model considered in this paper, i. e. , to Eq (5). In our case, we have bk (t) = (bk (1) (t) bk (2) (t) bk (3) (t) bk (4) (t) ), A= (− (λ14+λ12) λ21λ410λ12− (λ21+λ23) 0λ32λ140− (λ41+λ43) λ340λ23λ43− (λ32+λ34) ), Λ= (0000000000μ10000μ2) Assume that two transcription rates are the same, i. e. , μ1 = μ2 = μ. Then, we can show from Eq (10) that the binomial moments have the following explicit expressions (S1 Text) bk=b˜0 (3) μkk! ∏i=1k (i+γ1 (1) ) (i+γ2 (1) ) (i+γ3 (1) ) (i+α1) (i+α2) (i+α3) +b˜0 (4) μkk! ∏i=1k (i+γ1 (2) ) (i+γ2 (2) ) (i+γ3 (2) ) (i+α1) (i+α2) (i+α3) (13) where all γ are constants depending on system parameters (we omit their concrete expressions due to complexity), b˜0 (i) =b0 (i) / (b0 (3) +b0 (4) ), i=3,4, −α1, −α2, −α3 are nonzero characteristic values of M–matrix A, and b0 (k) =∏i=13 (βi (k) /αi), 1≤k≤4 with −β1 (k), −β2 (k), ⋯, −βN−1 (k) being nonzero characteristic values of matrix Mk (the minor of the element akk of matrix A). To that end, using Eq (11), we can further arrive at the following analytical expression of mRNA distribution P (m) =b˜0 (3) P1 (m) +b˜0 (4) P2 (m) (14) where Pi (m) =μmm! (γ1 (i) ) m (γ2 (i) ) m (γ3 (i) ) m (α1) m (α2) m (α3) mF33 (m+γ1 (i) m+γ2 (i) m+γ3 (i) m+α1m+α2m+α3|; −μ), i=1,2 (15) in which Fnn (a1⋯anb1⋯bn|; z) is a confluent hypergeometric function [62,63], and (c) n is the Pochhammer symbol and defined as (c) n = Γ (c+n) /Γ (c) with Γ (c) being the common Gamma function. Eq (14) indicates that the mRNA distribution is a linear combination of two distributions with each similar to that in the common two–state gene model at the transcription level [62]. Correspondingly, the mRNA noise intensity is given by Eq (12) with b1=μ[C (1) b0 (3) +D (1) b0 (4) ] (1+α1) (1+α2) (1+α3), b2=μ22[C (2) W33 (1) +D (2) W34 (1) ]b0 (3) +[C (2) W43 (1) +D (2) W44 (1) ]b0 (4) (1+α1) (i+α2) (1+α3) (2+α1) (2+α2) (2+α3) (16) where C (i), D (i) and all W are given in S1 Text. Note that Eqs (14) – (16) are exact for any values of the system parameters, and hence can be directly used to numerical calculations. The above general yet exact results can be simplified in some cases. For example, according to experimental evidence [20,33], we can set λ14 = λ23, λ12 = λ34, λ41 = λ32, and λ21 = λ43 for the side-by-side pattern of interacting DNA loops; 0. 5λ14 = λ23,0. 5λ12 = λ43, λ41 = λ32, and λ21 = λ34 for the alternating pattern; and (4e−0. 5d + 1) λ14 = λ23, (4e−0. 5d + 1) λ12 = λ43, λ41 = λ32, and λ21 = λ34 for the nested pattern. In the following, we consider two particular cases. First, if λ23 = λ43 = kλ (k is a positive constant but possibly depends on DNA loop lengths along the DNA lines) and λ12 = λ14 = λ21 = λ34 = λ32 = λ41 = λ, then it is interesting that we find that for the side-by-side pattern, the steady–state mRNA number follows a distribution of the form (S1 Text) P (m) =μ˜mm! (a1) m (a2) m (b1) m (b2) mF22 (m+a1m+a2m+b1m+b2|; −μ˜) (17) where μ˜, ai and bi are constants depending on system parameters including parameter k, and in particular, they are functions of two transition rates λ12 and λ14, each that is further functions of DNA loop lengths given by Eq (1) but may take an arbitrarily positive value. This distribution is in accord with our intuition since two DNA loops do not interact to each other, which leads to the bursting expression of the gene in a manner similar to that in the common two–state gene model at the transcription level. The corresponding noise intensity is given by ηm=[ (1+a1) (1+a2) (1+b1) (1+b2) +1μ˜]b1b2μ˜a1a2−1 (18) Next, we consider an extreme case, i. e. , all the transition rates are the same (this is possible only for the side-by-side pattern). In this case, we find that the mRNA number follows a simpler distribution of the form (also S1 Text) P (m) =μmm! (λ) m (2λ) mF11 (m+λ, m+2λ; −μ) (19) The corresponding noise intensity is reduced to ηm=14 (1+2λ) +2μ (20) The above analytical results in principle show how three connection patterns and two loop lengths as well as tracking probability or tracking ratio (characterized by parameter r) impact gene expression (including mRNA distribution, mean expression and noise intensity), but the results are implicit and not intuitive. Here, we perform numerical calculation to give intuitive results for this impact. Since the output level is directly related to the length of the blue loop (along the DNA line), we first let this length change but keep the green loop length fixed. Also since our focus is on effects of connection pattern and loop length on gene expression, we do not consider the impact of the communication mechanism (i. e. , we do not consider tracking between loop elements) in this subsection, and will leave it to the next subsection for investigation. Figs 3C and 2D show how the mean expression level and the mRNA noise intensity depend on the blue loop length for three connection patterns: alternating loops, nested loops and side-by-side loops, where the green loop length is fixed at 1500 whereas the blue loop length changes in the interval (0,1000). We observe that the mean level in each of the three structures is fundamentally a monotonically decreasing function of the blue loop length whereas the noise intensity is fundamentally a monotonically increasing function of this length, except for a small range of the loop length where both have local extreme values (more specifically, the mean expression level has a local maximum at a certain value of the DNA loop length whereas the noise intensity has a local minimum at another value of the length). Moreover, the mean level is generally lower in the case of alternating loops or higher in the case of nested loops than in the case of side-by-side loops, implying that the cross–type structure reduces the mean expression whereas the inline–type structure enlarges the mean level. On contrary, the noise intensity is generally higher in the case of alternating loops or lower in the case of nested loops than in the case of side-by-side loops, implying that the cross–type structure enlarges the expression noise whereas the inline–type structure reduces the expression noise. From these results, however, one cannot clearly see qualitative differences in the effect of the blue loop length on gene expression between alternating and the nested structures. Since two loops in the side-by-side structure are independent of each other or since the green loop does not influence the blue loop length that is directly related to the expression efficiency, we then take the side-by-side structure as a reference or control to show the effects of alternating loops and nested loops on gene expression, respectively. Fig 3A and 3B show how the mean expression level and the noise depend on the blue loop length. From these two diagrams, we observe that the nested structure increases the mean expression level but reduces the expression noise. In contrast, the alternating structure decreases the mean expression level but enlarges the expression noise. We also observe that the effects are most remarkable in the cases of small blue loop lengths. All these observations are in good accord with experimental observations [20,46] and also with our intuition. This is because in the alternating structure, the looping of Su and Hw intervenes the formation of the enhancer–promoter loop, thus reducing the latter looping rate and further weakening gene expression. In the nested structure, the Su and the Hw form a loop, which decreases the length of the enhancer–promoter loop but increases the looping rate of enhancer and promoter. In addition, we demonstrate numerical results on time evolutions and distributions of the mRNA number in S1 Text, referring to Figs C–F in S1 Text wherein more numerical results are shown. This demonstration has two aims: the one is to show that the results predicted by theory are in accord with those obtained by the Gillespie stochastic simulation algorithm; the other is to show that our binomial moment approach can be easily used to analysis of transient dynamics. After having understood the fundamental functions of alternating and nested structures, we now turn to considering the effects of the forms that the enhancer–promoter communication takes on gene expression. By comparing differences in effects induced by direct looping and facilitated tracking between these two connection patterns, we try to answer the question of which of these two representative mechanisms is more reasonable or more possibly used by living organisms. In the previous subsection, we considered the case that the blue loop length changes but the green loop length is fixed. In this subsection, we will consider the case that the green loop length changes but the blue loop length is fixed. This consideration is mainly for comparing qualitatively different effects of direct looping (i. e. , no tracking, characterized by r = 0) and facilitated tracking (characterized by 0 < r ≤ 1) mechanisms on gene expression. First, we consider the case of the alternating structure. Note that in this structure, the Su and Hw pair indirectly represses gene expression by negatively influencing the looping of enhancer and promoter. Since parameter r can represent the probability that one DNA element tracks the other DNA element, we call it as the tracking ratio. Also note that for this cross–type structure, the larger the size of r is, the more difficult is the enhancer–promoter looping, implying that gene expression is weakened. In addition, note that for the alternating loop pattern, the tracking of looping elements is limited, and this limitation implies that the role that the Su-Hw looping weakens indirectly gene expression is also limited. In order to show the effect of the communication mechanism on gene expression in the case of the cross–type structure, we first compare the difference in effect between the cases: tracking exits (r = 0. 1) and no tracking exits (r = 0), referring to Fig 4A and 4B. We observe that the mean level in the case of tracking is in general lower than that in the case of no tracking but the noise intensity is in general higher in the former case than in the latter case. Then, we show the dependences of the mean level and the noise intensity on the tracking ratio for a fixed green loop length, referring to Fig 4C and 4D. From these two diagrams, we observe that the mean level is a monotonically decreasing function whereas the noise intensity is a monotonically increasing function. The monotonicity becomes more apparent for small tracking ratios but disappears when the tracking ratio is beyond a threshold. Fig 4E or 4F shows a panoramic view for how the mean expression level or the noise intensity depends both on the green loop distance and on the tracking ratio. We observe that (1) only for moderate loop distances, does the mean expression level or the noise intensity have apparent differences between tracking and no–tracking cases; (2) the dependence of the mean level on both the loop length and the tracking ratio is fundamentally opposite to that of the noise intensity on both the loop length and the tracking ratio. These results are in accord with those shown in Fig 4A–4D. Then, we consider the case of the inline–type structure. Note that in this structure, the Su and Hw pair indirectly represses gene expression by positively influencing the looping of enhancer and promoter. Also note that for this structure, the larger the size of r is, the enhancer and the promoter form a loop more easily (just opposite to the case of the cross–type structure), implying that the gene expression is promoted since we have assumed that the enhancer–promoter looping promotes gene expression. In addition, note that for the nested loop pattern, the tracking between looping elements has no limitation, implying that the role of what the Su-Hw looping enhances directly gene expression is further promoted with increasing the size of r. In addition, we use Fig 5E or 5F to show a panoramic view for how the mean expression level or the noise intensity depends on both the green loop distance and the tracking ratio. Similar to the case of the cross–type structure, we still observe the following two points: (1) only for moderate loop distances, does the mean expression level or the noise intensity have apparent differences between tracking and no–tracking cases; (2) the dependence of the mean level on both the loop length and the tracking ratio is fundamentally opposite to that of the noise intensity on both the loop length and the tracking ratio. These results are in accord with those shown in Fig 5A–5D. In order to show the effect of communication mechanisms on gene expression in the case of the inline–type structure, we first compare the difference in effect between the cases that tracking exits (r = 0. 1) and no tracking exits (r = 0), referring to Fig 5A and 5B. We observe that the mean level in the case of tracking is lower than that in the case of no tracking but the noise intensity is higher in the former case than in the latter case. Then, we show the dependences of the mean level and the noise intensity on the tracking ratio for a fixed green loop length, referring to Fig 5C and 5D. From these two diagrams, we observe that the mean level is a monotonically increasing function whereas the noise intensity is a monotonically decreasing function. This monotonicity becomes more apparent for small tracking ratios, but finally disappears when the tracking ratio is beyond a threshold. In addition, we use Fig 6 to show the quantitative dependence of the relative change ratio that is defined as the ratio of the difference between the mean expression level (or the noise intensity) in the tracking case and that in the non–tracking case over the mean expression level (or the noise intensity) in the non–tracking case, on the green loop length. We observe that for the inline–type structure, the relative change ratio of the mean expression is more than zero whereas that of the noise intensity is less than zero, remarkably for moderate lengths of the Su-Hw loop. For the cross–type structure, however, this ratio of the mean expression is less than or equal to zero whereas that of the noise intensity is equal to or greater than zero, remarkably for moderate lengths of the green loop. In particular, the largest relative change ratio for the mean expression is 2% lower in the non–tracking case than in the tracking case for the alternating structure, whereas it is 4% higher in the former case than in the latter case for the nested structure. In a word, the results shown in Fig 6 indicate that nested loops play a role of enhancing gene expression and reducing noise while alternating loops play a role of repressing gene expression and enlarging expression noise. Moreover, there is a limit length of the green loop such that the effect of this control almost disappears. The combination of Figs 4–6 implies that from the viewpoint of controlling gene expression, the facilitated–tracking communication (i. e. , r ≠ 0) is advantageous over the direct–looping communication (i. e. , r = 0), independent of connection patterns between interacting DNA loops. In addition, there is an optimal loop length such that the mean expression level or the expression noise intensity is maximal or minimal. This is an interesting phenomenon. From the viewpoint of mathematics, the occurrence of this phenomenon is generated because the looping rates are nonlinear functions of loop lengths along DNA lines. From the viewpoint of biology, however, this implies that if two loop elements are very close or very far away, then the expression level is low, and that there is a suitable location of these two loop elements such that the expression level reaches the highest value. In order to show effects of the number of interacting DNA loops on the obtained-above qualitative results, we also investigate three interacting chromatin loops in Fig M in S1 Text, and find that the qualitative results are invariant.
  We have studied the influence of interactions between DNA loops on gene expression, using a quantitative model established by mapping fundamental structures of DNA–looping interactions into a multistate model of stochastic gene expression. In contrast to previous gene models that assume direct spatial contact between a distal enhancer and a proximal promoter [34–38], our model has the advantage in many aspects, e. g. , three main factors (connection pattern, loop length and communication form) characterizing DNA–looping interactions are easily incorporated in the model (even including those in more complex cases). Also for example, it can be analytically solved using the binomial moment method that we previously developed [39], with results that can well show how key yet experimentally testable parameters associated with the interactions between DNA loops such as looping rates, loop lengths and the tracking ratio affect gene expression levels. Thus, our model provides a possible platform for experimentally testing effects of specific connection patterns, loop lengths and communication forms on gene expression. In fact, we have found good agreement between results obtained by our model and experimental measurements [20,46]. We have analyzed effects of three fundamental structures of DNA–looping interactions (alternating loops, nested loops, and side-by-side loops) on gene expression including the mean level and the output noise, and shown that different structures have different effects. However, as pointed out in the introduction, enhancers and promoters may be connected in a highly complex network of DNA–looping interactions [13,42,43], remarkably in eukaryotic cells. One main question that needs to be addressed is at which step during gene activation, various nucleoprotein complexes assemble at distant enhancers, and how these complexes then contribute to promoter accessibility, the preinitiation complex recruitment and/or assembly, and transcription initiation and elongation. Enhancers have been shown to have a role in the preinitiation complex recruitment at target promoters [64–68], the removal of proteasome complexes at promoters [69], the generation of intra-chromosomal loops between regulatory regions [70], and the regulation of elongation [71–75]. Enhancers are also involved in the removal of repressive histone modifications [76–79], suggesting that they also contribute to the delivery of enzymes that regulate histone modifications [80]. In a word, enhancers in eukaryotic genomes can be many hundreds of kilobases away from the promoter they regulate [8–10], and the intervening DNA can contain other promoters and other enhancers [11–14]. For these complex cases, how the regulatory influence of distal elements is exerted efficiently and specifically at the correct promoters is worth further investigation. In this paper, we have also shown the influences of two representative mechanisms of loop element tracking (namely direct looping and facilitated tracking) on gene expression. By comparing their respective performance in controlling mean expression and in tuning expression noise, we conclude that from the perspective of controlling gene expression and tuning expression noise, the facilitated tracking is superior to the direct looping, remarkably for moderate loop lengths. In spite of this, a question is whether increasing the enhancer–promoter distance in higher organisms favors particular mechanisms of interactions between regulatory elements, for example, looping rather than tracking or linking. In fact, given a larger distance (>10 kb and up to 100s of kb) separating many enhancers from their target promoter, it is difficult to envisage a mechanism in which the intervening chromatin is directly involved in a mechanism of enhancer–promoter communication. Therefore, tracking mechanisms are likely limited to enhancers that are close (1–10 kb) to their target promoters. In addition, the stiffness of the chromatin fibre might restrict short–range enhancer–promoter interactions, with a minimal estimated length of 10 kb for uninterrupted 30-nm chromatin fibres and 0. 5 kb for naked DNA [71,81,82]. In a more general sense, genes transcribed by RNA polymerase have two distinct families of cis-acting elements: the promoter and more remote cis-regulatory elements. The former, the length of which is generally less than 1kb, is composed of a core promoter [83,84], whereas the latter, the length of which is in general more than 1kb, can be enhancers, silencers, insulators, or locus control regions [3]. The exact composition of core promoter elements may be a key determinant of enhancer–promoter specificity [85,86]. These regulatory elements may communicate to each another in different manners [1]. When various possible regulatory elements and the distances among them are simultaneously considered, how the corresponding communication mechanisms affect gene expression still remains poorly understood and further studies are required. Finally, although numerous studies in several multi–gene clusters have shown that gene activation by a remote enhancer is associated with chromatin loop formation [87,88], the nature of its formation and involvement in gene regulation is not fully understood. Li, el al. [2], proposed that the major feature that determines loop formation is the chromatin flexibility, which is modulated by histone acetylation (and other modifications). Based on this, they also concluded that histone modifications lead to the determination of the probability of interaction between a remote enhancer and its cognate genes, and these features constitute the main elements of the facilitated chromatin–looping hypothesis. Our model analysis has partially verified reasonability of this hypothesis, but relevant experimental tests are required.
summary: The classic model of gene expression assumes direct spatial contact between a distal enhancer and a proximal promoter. Recent experimental evidence supported that enhancers and promoters are linked in a highly complex network of DNA–looping interactions, but it is unclear how these interactions including communication forms affect cell-to-cell variability in gene expression. Here we establish a quantitative model by mapping three possible patterns for two interacting DNA loops into a multistate model of gene expression. We first show that different connection patterns have different effects on gene expression, remarkable for short loop lengths. Then, we compare the results on the effects of two representative communication forms and find that the facilitated tracking performs better than the direct looping in controlling mean expression and in tuning expression noise, remarkable for moderate loop lengths. As such, we conclude that the former communication form is superior to the latter, rationalizing the facilitated DNA–looping hypothesis.
keywords: engineering and technology
  gene regulation
  signal processing
  binomials
  messenger rna
  noise reduction
  dna transcription
  mathematics
  algebra
  dna
  epigenetics
  chromatin
  polynomials
  chromosome biology
  gene expression
  biochemistry
  rna
  enhancer elements
  cell biology
  nucleic acids
  genetics
  biology and life sciences
  physical sciences
title: Effect of Interaction between Chromatin Loops on Cell-to-Cell Variability in Gene Expression