Real sequence effects on the search dynamics of transcription factors on DNA

Recent experiments show that transcription factors (TFs) indeed use the facilitated diffusion mechanism to locate their target sequences on DNA in living bacteria cells: TFs alternate between sliding motion along DNA and relocation events through the cytoplasm. From simulations and theoretical analysis we study the TF-sliding motion for a large section of the DNA-sequence of a common E. coli strain, based on the two-state TF-model with a fast-sliding search state and a recognition state enabling target detection. For the probability to detect the target before dissociating from DNA the TF-search times self-consistently depend heavily on whether or not an auxiliary operator (an accessible sequence similar to the main operator) is present in the genome section. Importantly, within our model the extent to which the interconversion rates between search and recognition states depend on the underlying nucleotide sequence is varied. A moderate dependence maximises the capability to distinguish between the main operator and similar sequences. Moreover, these auxiliary operators serve as starting points for DNA looping with the main operator, yielding a spectrum of target detection times spanning several orders of magnitude. Auxiliary operators are shown to act as funnels facilitating target detection by TFs.

and theoretically [29][30][31][32][33] . Finally, the fact that genes, that interact via local TFs, are statistically proximate along the prokaryotic genome (colocalisation) was argued to be due to the increased interaction rates (rapid search hypothesis) [34][35][36] . In line with the increasing knowledge of the microscopic details of gene regulation many computational studies appeared that go beyond the typical idealisations 19,37,38 .
Motivated by recent experiments showing that on encounter the target operator is not detected with certainty by a TF sliding along the DNA 15 , we here combine theoretical and simulations analyses to quantify the sliding motion of a TF along the real nucleotide sequence of a common E. coli strain in the presence of crowding proteins on the DNA. We establish a model including search and recognition states of the TF in combination with the barrier discrimination model 10,24 with a position weight matrix (PWM) based binding energy approach 39 . We also include looping effects-as often studied in thermodynamic models 40 -in the present model: the TF, for instance, the lac repressor dimer, can simultaneously bind to two operators, mimicking the intersegmental transfer mechanism 5,9,10 .

Blockers and movers, and the role of auxiliary operators
We describe the sliding motion of a TF for its target operator along DNA, on which N block other proteins are bound, so-called blockers or roadblocks 18 . We focus on immobile blockers, keeping in mind that mobile blockers may add another layer of complexity 41 . The N block non-overlapping blockers are positioned randomly and partition the DNA into N block + 1 intervals. We assume that the TF cannot by-pass the blockers, see Fig. 1. Where the DNA is not occupied by a blocker, the TF can bind to the DNA in two orientations. In the case of palindromic sequences the binding energies in both orientations are equal (see also the score values in Methods).
We first focus on the processes in the target region carrying possible binding positions between the two nearest roadblocks to the left and to the right of the main operator O1. Such roadblocks could be proteins like H-NS or HU 42 . We only consider configurations in which the main operator is accessible. From both simulations and an approximate analytical approach we determine the probability p t that the TF detects the target in the correct orientation before dissociation. The TF starts from a random position in this target region.

Simulation scheme
We focus on base pairs 359,990 to 370,010 of E. coli strain K-12 MG1655 from ecocyc.org 43 , comprising the genes lacA, lacY, and lacZ as well as the three operators O1, O2, and O3, to which the lac repressor (LacI) can bind 44 . The sequence length is 10,021 base pairs (bps). Since the binding motif of LacI covers w 21 = bps we obtain 10,001 possible binding positions in two orientations. We choose N 71 block = blockers of size w to match the occupation fraction of Tabaka et al. 21 .
The general simulation scheme is depicted in Fig. 1. At each position the TF can be either in the loosely bound search state or in the tightly bound recognition state. In the search state the TF has four possible actions: the particle can move to the left or to the right, it can dissociate, or it can change to the recognition mode at its position. If the latter occurs at the position of the main operator O1, the corresponding time is saved as a first target detection. We later deal with dissociation from the DNA. Once in the recognition mode, we assume that the binding is so tight that the TF cannot move to neighbouring positions. As looping is neglected in this first, linear version of the model, its only option is to return to the search state at this position. The rates at which these transitions occur depend on the energetic barriers that need to be crossed during the internal protein dynamics. These are determined by the standard Gillespie algorithm 6,45 . Methods contains a detailed description of the simulations. Times are measured in units of the inverse attempt rate 0 λ from Eq. (6) in Methods. The energy E s in the TF search state and the barrier E bs for sliding to a neighbouring base pair are assumed to be independent of the DNA sequence 10 is the difference between the score at position i and the average score in the data set. γ = − .
1 3378 is a proportionality factor (Methods). The volatility parameter α tunes the sensitivity of E bc i , to the DNA sequence. If 0 α= the barrier height does not change with the sequence and therefore this corresponds to blind testing of the sequence. If 1 α= , an induced fit mechanism is at work. The closer the probed sequence is to that of the target, the faster the TF switches to the target-sensitive recognition mode since the barrier height changes exactly as much as the energy in the recognition mode. To obtain the target detection probability before dissociation shown in Fig. 2  = × independent simulations starting from random positions in the target region were performed and it was counted in how many cases the target was reached. As we show here our model (1) for the energy score relation together with the additional element of the volatility α elucidate the role of the sequence sensitivity in the speed stability tradeoff of TF search processes.

Theoretical approach
We compare the simulations results of Fig. 2 to a theoretical model based on a target region with N possible binding positions. For mathematical details see Methods.
The fundamental parameters are the sliding rate Γ to neighbouring positions, the rate k st of a conformational switch to the recognition mode at the target site resulting in direct target detection, and the dissociation rate k off from any site. At all non-target positions we assume constant rates for the changes between recognition and search modes, denoted by k sr and k rs . We place the target at bp m and the TF starts at a random position. As detailed in Methods these quantities determine the mean target detection time N m τ , (see below) and the probability to reach the target before dissociation p N m t ( , ), written as   there is no longer an energy barrier to be crossed at the target site and thus no more changes are observed. Lines of matching colour in Fig. 2 are results of the analytical model, Eq. (2). The target is either centred (full lines) or located at the boundary of the target region (dashed).
The simulated data scatter nicely between the two limiting theoretical lines for centred and boundary target positions over three orders of magnitude in the size N of the target region. p t decreases monotonically with N , as large target regions on average imply longer paths which have to be traversed en route to the target, implying a higher risk to dissociate. Larger α values, corresponding to a searcher which checks more often for the target, lead to a higher detection probability. Another effect of α concerns the influence of the target position. For small values of α the corresponding curves nearly coincide, i.e., there is no significant target position dependence. For higher α values, centred targets effect a substantially higher detection probability as the full lines lie above the corresponding dashed ones. Thus, only when the target detection probability on an individual encounter reaches substantial values, a suitable position of the target pays off.
We see that for the target detection probability the theoretical model, in which all energies on non-target sites are replaced by average values, nicely reproduces the results of the simulations based on sequence specific binding energy values. Fig. 3 the mean detection times N m τ , to the target are shown for the same α values used in Fig. 2. Since the particles can dissociate, N m τ , is a conditional time: given that the particle detects the target with the probability shown in Fig. 2, at what time will this occur on average. The symbols in Fig. 3 show the simulations results, the lines correspond to the theoretical model with a centred target (full lines) and a target at the boundary (dashed).

Target detection time. In
The features of Fig. 3 fall into two cases. For N  100, as with the detection probabilities above the simulations agree well with the theoretical model for all α values. Again, a clear ordering with α occurs: volatile TFs (large α) find the target quicker than nearly blind TFs with 0 1 α= . (cyan). Moreover, only in the case of large α, when individual encounters with the target have a substantial probability for target detection, the target position comes into play (e.g., for the red lines). This is one of our central results.
For N 100 ⪆ , apart from simulations data consistent with the theoretical lines a second branch of results appears with target detection times nearly two orders of magnitude longer than expected. This effect can be rationalised by the presence of the auxiliary operator O3 in the target region. It resides 92 nucleotides away from the main operator O1 such that only target regions with a size larger than that can contain both operators 1 . If both operators are in the target region, the TF can change to the recognition mode at the auxiliary operator and thus become trapped away from the main operator. Such time consuming checks for the target may occur at any non-target position. However, at O3 this is particularly severe since it has a rather strong binding energy (see Fig. 4). The gapped energy spectrum yields search times which are way above the values of the theoretical model, since the latter assumes all non-main target sites to be energetically equivalent. Inspection of the upper branch of the results in Fig. 3 indicates that it barely contains data obtained with small α values (cyan and green). This can be explained by comparison with Fig. 2: in these cases even the probability to detect O1 is rather small. This effect is even more pronounced for the considerably weaker O3. However, when such TFs change to the recognition state at the auxiliary operator, they will spend more time there than particles with a larger α, since these face a larger barrier to be crossed (Eq. (1)). As not all target regions of size N  100 contain the auxiliary operator, the lower branch of results still coexists. Here the conditional target detection time increases with N but levels off to a plateau.
Conversely, for rather volatile searchers (red data points) in regions comprising both operators, for there is a slight tendency that the mean search time decreases with N . This results from the fact that these regions, which are only marginally longer than the distance between the two operators, by definition have both operators near the boundaries. This yields longer search times, similar to the case of shorter target regions, for which the dashed lines are always above the corresponding full lines in Eq. (3). We consider the influence of the location of the operators with respect to the non-specific blockers in more detail in the following paragraph.
Preference of O1 over O3. In the hypothetical situation of two equally strong operators in the target region, only their relative position in the target region would influence which one of them is more likely to be detected first. The biologically relevant situation considered here with two different operators is more subtle. When both O1 and O3 are in the target region we registered which one was detected first. The preference for O1 shown in Fig. 5 is given by the probability that O1 is detected first. The shift by 1 2 / leads to positive values when the probability is larger to detect O1 first.
To single out geometrical effects, the axis x quantifies which of the two operators is more central in the target region and thus has-from a geometric point of view-higher chances to be hit first. We define As expected, since O1 is the stronger operator, most of the data points are positive. For small α values (cyan and green in Fig. 5) it is more probable to detect O1 first, but the relative positions of the two operators are not significant. Increasing the volatility from 0 1 α= . to 0 3 . leads to a monotonic increase in the accuracy of discrimination between the two operators. For even larger values of α this accuracy decreases, since now the particle checks for the target often enough to detect the auxiliary operator with sufficient probability. Then, geometric effects become more important, as seen from the increasing slope  of the black dotted line for 0 4 α= . and the red dot-dashed line for 0 5 α= . . In the latter case, some negative values of the preference are observed, indicating that a volatile TF is more likely to detect the auxiliary operator first, if its position is much closer to the centre of the target region.
Intermediate α values enable the TF to detect the main operator first, without losing time from binding to the auxiliary operator. In terms of the search model presented so far, the occurrence of O3 appears like a design bug instead of a useful feature, since it delays the detection of the main operator. We now show that auxiliary operators in a more realistic scenario indeed act as funnels for TFs towards the main binding site.
Auxiliary operators make sense in presence of looping. As evident from Figs. 3 and 5 the presence of the auxiliary operator O3 in the target region significantly influences the rate of target detection. In an extension of our model several configurations can be distinguished depending on whether or not the two auxiliary operators are accessible. In a living cell the occupation with non-specific binders and thus the probability for a particular blocker conformation change in time.
To model the complete search process of a TF with two binding motifs such as LacI, we consider what happens after a dissociation from DNA. After dissociation the time spent in 3D is assumed to be exponentially distributed with mean time b τ . For the jump length x jump -like all the following lengths measured in bps-on the DNA effected by 3D excursions we assume the cumulative distribution for the length x stored in a random loop formed by a polymer chain occur due to the equivalence of polymers to random walks. For a random chain in three dimensions 1 5 η= . , while in the presence of excluded volume interactions the exponent increases to 2 2 η ≈ . [50][51][52] . Here we chose the lower exponent 1 2 η= . following the data by Priest et al. 53 . To obtain the cumulative distribution (5), we integrate the power law f x ( ) in between the lower and upper cutoffs x min and x max , and normalise this expression. Note that our results are not overly sensitive to the exact value of the exponent β, as in the free energy it corresponds to a logarithmic dependence on x.
Here we assume that a power law similar to Eq. (5) also applies to the jump statistics. Whenever the particle jumps out of the 10 kbps range that we study, we place it at a random position in our system, mimicking the complete loss of correlation with the dissociation position for long jumps. Unlike during sliding motion, it can change the orientation during a 3D relocation. To simplify matters we coarse-grain events outside the target region, since we are not interested in the sliding motion far away from the target. We then first simulate the mean dissociation times from all N block regions that do not contain the target. To this end, simulations are performed as outlined in the paragraph Simulation scheme, where the code is run f 20 avg = times multiplied by the length of the corresponding interval measured in bps to guarantee reasonable statistics.
Whenever the TF detects and binds to one of the auxiliary operators, apart from returning to the search state at this position there is the possibility to form a DNA loop with O1. For this event to occur, an initiation time is drawn from an exponential distribution with mean init τ , which is assumed to be the time needed to form a non-specific complex with the target region. To keep the number of parameters as low as possible we assume these initiation times to be equal for both auxiliary operators. To the loop initiation time we add a time lag, since after landing with its second half in the target region, the TF has to actually detect the main operator. The latter is obtained from a simulation as defined in Simulation scheme. The same process is possible the other way round: starting from binding to the main operator and, before switching to the search mode, closing a loop with one of the auxiliary operators. To simplify matters we do not model direct looping between the auxiliary operators. The times for releasing a loop are calculated similarly to the mean dissociation times above. We now study the full model with looping for N 180 ≈ .

Results II
Influence of the volatility. Of particular biological interest are the time spans free τ during which the operator is unoccupied, as in these intervals RNA polymerase can bind to the promoter and start transcription. We start with a conformation in which looping is precluded by blocking both auxiliary operators with non-specific binders. In Fig. 6 the distribution of free τ is shown for four values of the volatility parameter α.
In all cases we obtain two distinct peaks separating a short and a long time scale. For increasing values of α the first peak, located at around 10 0 time units, grows relative to the second one, located at around 10 6 time units. Since the total simulation time was fixed, the total number of events grows as well: The peak at short times is due to events when a TF, after switching from the recognition to the search mode, performs just a few sliding steps before returning to the recognition state at the target. Conversely, the long time peak corresponds to events when a TF dissociates, possibly multiple times, from DNA and loses correlation with the unbinding position, and thus leads to long time spans, in which the target operator is vacant. That the first peak gains in importance for larger values of α is due to the fact that, as seen above, the individual target detection probability is higher in that case.
We note that to initiate transcription, RNA polymerase must bind the promoter while the TF is not at the operator. If the repressor rebinds to the operator before an RNA polymerase manages to find the promoter, the cell does not "feel" these quick occupancy fluctuations and experiences only a single effective binding event of the repressor, and no transcription takes place (compare Ref. [54]).

Influence of looping and the average time spent in 3D.
We now choose a configuration in which the auxiliary operator O2 is vacant and we fix α to a value of 0 5 . . The corresponding results are shown by black lines in Fig. 7. Full lines are for the same values of b τ and init τ as in Fig. 6, dashed lines represent the case when both are ten times larger. We observe that both full lines still feature two peaks centred at 1 ≈ and 10 6  τ ≈ . × , these events can be self-consistently interpreted as return events to the target due to looping: the DNA was looped between the main and an auxiliary operator, dissociates from the main operator and reestablishes the loop. That the peak for fast rebinding events has a reduced size is due to the fact that our looping algorithm counts all fast fluctuations of the occupancy during which the loop still exists as a single long-lived event.
In the simulations without looping these events appeared explicitly (Fig. 6). Accordingly, the remaining events in the reduced first peak correspond to target rebinding without an existing loop to an auxiliary operator. Given that fast rebinding has no biological meaning, looping introduces a new intermediate time scale, and typical return times to the operator are greatly reduced, resulting in improved repression. This agrees with the observations of Choi and co-workers according to which DNA looping enables the cell to regulate gene expression on many time scales via distinct forms of dissociation events 55 . Comparing this behaviour to the black dashed line in Fig. 7, when both 3D excursions and looping take around ten times longer, shows that both the looping peak and the rightmost peak are shifted to larger times underlining the physicality of our interpretation.
If both auxiliary operators are accessible and O3 is in the target region (red lines in Fig. 7), the results are similar to the previous ones (Fig. 6). Three peaks are observed, and increase of b τ and init τ shifts the peaks-apart from the b τ -independent fast rebinding peak-to the right. There is one major difference between the two settings: When both auxiliary operators are accessible, the size of the third peak is nearly as large as the second one. Thus, very long return times occur more often when both auxiliary operators are present. The significant changes of the target search times in the presence of the auxiliary operators are our other central result.

Discussion
One-dimensional sliding of a TF along the DNA is a vital ingredient of the facilitated diffusion model. Sliding is indispensable in the final step of the search for the specific binding site by the TF, namely, the recognition of the binding sequence. For a real bacterial DNA sequence we here analyse in detail the dynamics of the TF sliding in a region around the main operator in the presence of roadblocks, e.g., proteins like H-NS or HU 42 . For a minimal set of parameters we unveil the role of the density of the roadblocks and the DNA sequence on the detection speed of the target sequence. Our results underline the special role played by auxiliary TF operators. These auxiliary operators act as a funnel for the TF to facilitate the target search in the nucleotide sequence.
More specifically, combining a simplified theoretical model and simulations we follow a TF moving in a region around the main operator delimited by two non-specifically bound roadblocks while switching between a search mode, in which it shuttles along the DNA while being blind to the target, and a recognition mode, in which it cannot move along DNA but which is essential to detect the target. The interconversion rates depend on the underlying sequence. Motivated by recent experiments showing that not every target encounter leads to detection of the target sequence 15 , we interpolated between the extreme cases of nearly blind switching between the modes and an induced fit situation, in which the energetic barrier to be crossed changes as much as the specific binding energy. Numerical results for the probability to detect the target before dissociation and for the mean detection time demonstrate impressive agreement with our theoretical model (see Fig. 2 and upper branch in Fig. 3), as long as no further binding sites of similar strength are present. If an auxiliary operator is within the target region, an intermediate rate of checking for the target yielded the highest accuracy in discrimination between main and auxiliary operators. However, while auxiliary binding sites act as traps in the simplified model, in the more realistic situation when DNA looping is allowed, they can be seeding points for the formation of loops joining two operators. In the second part we therefore included looping in the simulation. For . In all curves: 0 5 α= . . our parameters, this leads to quick rebinding events to the main operator and thus increases significantly the local effective TF density, in accordance with classical observations. This approach can be easily transferred to other TFs with known binding motif.
Given the fairly large number of the parameters involved and the complexity of the dynamics conveyed by the broad range of apparent time scales, definite quantitative statements of this problem are hard to give. Furthermore, the target search here was modelled for a single TF, while in a living bacteria cell approximately a dozen lac repressors perform this task simultaneously. Additionally, other TFs could partially block the specific binding site of the TF under consideration, and could impede the establishment of the lac specific loop. As recent studies showed (for instance, see Ref. [20]) the effects of additional binders are not always obvious and require careful analysis. Since the binding to the operator(s) is rather strong, it is questionable to assume that the TFs are independent and we face a multi-body problem. However, the concentration effects and the expression output in terms of the occupancy of all three operators were successfully studied in terms of thermodynamic models 40 using similar language. As our simplified theoretical model for events in the target region yields such a good agreement with the numerical simulations and given that more and more quantitative experimental results appear, it seems to be a logical extension to equip thermodynamic models with rates obtained from our model presented herein. Moreover, the accessibility of the three operators could be modulated in time to mimic the mobility of nonspecific binders which can block the operators. In this spirit we believe that the results reported herein represent an important step forward toward the quantitative understanding of gene regulation in living prokaryotic cells, and form the basis for future, more detailed models.

Methods
Here we describe the simulations method and the calculations for the above results.
Numerical simulation of the simplified model. The TF is present in either the loosely bound search state or tightly bound recognition state. In the search state at position i the TF can either slide to the neighbouring sites i 1 − or i 1 + while remaining in the search state, it can dissociate or switch to the recognition state at the same position (Fig. 1). Such a switching event at the target site (the operator O1) corresponds to detecting the target. This differs from the approach of Ref. [56], in which a further target detection step was used after changing to the recognition state at the target site. If the particle is next to a blocker and tries to move onto the excluded site, the move is cancelled. The standard Gillespie algorithm is used to draw the rates for the above events. A central role is played by the energetic barriers which need to be crossed, measured in units of k T B with respect to the unbound state of zero reference energy (similar to Refs. 22,57).
In the recognition state we assume the TF to be immobile. In the first version of the model without looping the TF can only return to the search state. Generally when going from state a with energy E a to state b with energy E b , separated by an energetic barrier E ba , the rate, k ab for this step is given by In absence of a barrier (E E ba a < ) and when the energy of the final state is smaller than that of the initial state (E E b a < ) the reaction is assumed to occur with attempt rate 0 λ , which is the inverse of the elementary time step in which all times are measured. To convert our results to real times, this time step can be related to the known 1D diffusion coefficient of a given TF. We note that our approach differs from the convention of Ref. [58,59], in which the specific binding barrier has to be crossed each time the TF slides to a neighbouring position.
We fix the energy difference between the specific binding energy at the main operator O1 and the energy in the search state as 15 3 . 60 . With the choice E 7 s = − applied in the main text this implies E 22 3 O1 = − . (Fig. 4). The proportionality factor γ can be determined once all values of the score matrix are known via the above mentioned demand E . At the lower end of the energy spectrum the three operators can be recognised. Note that there is an energetic gap to all other binding sites, see the discussion of such a gapped situation in Ref. [63].
, , This set of equations is more conveniently treated in Laplace space with respect to time, where we denote the variable conjugate to t by u and the corresponding functions with a tilde. For convenience we omit the explicit argument u in the following.
If the particle starts its motion in the search state, initially the probabilities p N j ,  vanish. This is due to the simple proportionality between p N j In this limit the existence of the recognition state away from the target simply slows down the mean target detection time via the prefactor q 1 + of the second term. In the limiting case q 0 → , when the recognition state is never entered unless the particle is on the target site, this result reduces to the classical solutions for incoherent exciton hopping,