Kinetics of Interior Loop Formation in Semiflexible Chains 6 Changbong Hyeon and D. Thirumalai 0 0 Biophysics Physics Program, 2 n Institute for Physical Science and Technology a J University of Maryland, 9 1 College Park, MD 20742 ] t f o s . Loop formation between monomers in the interior of semiflexible chains describes t a m elementary events in biomolecular folding and DNA bending. We calculate ana- - d lytically the interior distance distribution function for semiflexible chains using a n o mean-field approach. Using the potential of mean force derived from the distance c [ distribution function we present a simple expression for the kinetics of interior loop- 1 v ing by adopting Kramers theory. For the parameters, that are appropriate for DNA, 6 2 the theoretical predictions in comparison to the case are in excellent agreement with 4 1 explicit Brownian dynamics simulations of worm-like chain (WLC) model. Theinte- 0 6 rior looping times (τ ) can be greatly altered in cases when the stiffness of the loop 0 IC / t differs from that of the dangling ends. If the dangling end is stiffer than the loop a m then τ increases for the case of the WLC with uniform persistence length. In con- - IC d n trast, attachment of flexible dangling ends enhances rate of interior loop formation. o c The theory also shows that if the monomers are charged and interact via screened : v i Coulomb potential then both the cyclization (τc) and interior looping (τIC) times X r greatly increase at low ionic concentration. Because both τ and τ are determined c IC a essentially by the effective persistence length (l(R)) we computed l(R) by varying the p p range of the repulsive interaction between the monomers. For short range inter- actions l(R) nearly coincides with the bare persistence length which is determined p largely by the backbone chain connectivity. This finding rationalizes the efficacy of describing a number of experimental observations (response of biopolymers to force and cyclization kinetics) in biomolecules using WLC model with an effective persistence length. 2 I. INTRODUCTION The kinetics of formation of contact between the ends of a polymer chain has a rich history.1,2 Both experiments,1,2 theory,3,4,5,6,7,8,9 and simulations10,11,12 have been used to address the ele- mentary event of the dynamics of end-to-end contact formation (or cyclization kinetics) (Fig.1- A). Contact formation between two reactive groups separated by a certain distance along the chain is a basic intramolecular rate process in a polymer. Recently, there has been renewed interest in understanding the looping dynamics that has been studied both theoretically3,4,5,6,7 and experimentally12,13,14,15,16 because of its fundamental importance in a number of biological processes. The hairpin loop formation is the elementary step in RNA folding,17 structure for- mation in ssDNA,18,19 and protein folding.13,20,21,22,23 Cyclization in DNA has recently drawn renewed attention not only because of its importance in gene expression24,25 but also it provides a way to assess DNA’s flexibility. The promise of using single molecule technique to probe the real time dynamics of polymer chains has also spurred theories and simulations of cyclization kinetics. Using loop formation times between residues that are in the interior as the most ele- mentary event in protein folding, it has been argued, using experimental data and theoretical expression for probability for loop formation in stiff chains, that the speed limit for folding is on the order of a 1 µs.26 These examples illustrate the need to understand quantitatively the elementary event of contact formation between segments of a polymer chain. Even without taking hydrodynamic interactions into account theoretical treatment of cy- clization kinetics in polymer chains is difficult because several relaxation times and length and energy scales are interwined. At the minimum the variation of time scale for cyclization (τ ) c with polymer length is dependent on polymer relaxation time (τ ). In biopolymers additional R considerations due to chain stiffness and heterogeneity of interactions between monomer (amino acid residue or nucleotides) must be also taken into account. Majority of the cyclization kinetics studies on synthetic polymers2 have considered examples in which the contour length (L) of the polymer is much greater than its persistence length (l ). In contrast, loop formation dynamics p in biopolymers have focused on systems in which L/l is relatively small. In disordered polypep- p tide chains L/l can be as small as 3,15,27 while in DNA L/l < 1.16,28 Thus, it is important to p p develop theoretical tools for the difficult problem of loop formation dynamics for arbitrary L 3 and l . Despite the inherent complexities in treating loop formation in biopolymers it has been p found that the use of polymer-based approach is reasonable in analyzing experimental data on cyclization kinetics in proteins8,15 and DNA.16 Inthispaperweareprimarilyconcernedwiththeloopingdynamicsbetween interiorsegments of a semiflexible chain. While a lot of theoretical and experimental works (mentioned above) have been done on the end-to-end looping (Fig.1-(a)), only a few studies have been reported on the contact formation between monomers in the interior of a chain (interior looping) (Fig.1- (b)).29,30,31,32,33 There areafewreasons toconsider kinetics ofinterior looping. (1) Thebiological events such as hairpin formation and DNA looping often involve contact formation between monomers that are not at the ends of the chain. For example, it is thought that the initiation of nucleation in protein folding occurs at residues that are near the loop regions.34 The residues that connect these loops are in the interior of the polypeptide chain. Similar processes are also relevant in RNA folding.35 (2) It is known that for flexible chains with excluded volume interactions (polymer in a good solvent) the probability of loop formation is strongly dependent on the location of the two segments. For large loop length (S) the loop formation probability, P(S), in three dimensions for chain ends Sθ1 where θ 1.9 while P(S) Sθ2 with θ 2.1 1 2 ∼ ≈ ∼ ∼ for monomer in the interior.36 Although the values of θ and θ are similar it could lead to 1 2 measurable differences in loop formation times.32 The rest of the paper is organized as follows. In section II we present the physical consider- ations that give rise to the well-established results for τ for flexible chains. The extension of c the arguments for flexible chains to semiflexible polymers suggests that the local equilibrium approximation can be profitably used to analyze both cyclization kinetics and interior looping dynamics. The basic theory for the equilibrium distance distribution between two interior segments s and s (Fig.1-(b)) is presented in section III. Using the equilibrium distribution 1 2 function and adopting Kramers theory and following the suggestion by Jun et. al.,37 we obtain an analytical expression for time scale τ for interior contact formation in section IV. Explicit IC results of simulations of worm-like chain (WLC), which validate the theory, are presented in section V. In section VI we consider the kinetics of interior loop formation in WLC in which the stiffness of the loop is different from that of the dangling ends. Section VII describes the consequences of screened Coulomb interaction between monomer segments on cyclization 4 kinetics and interior looping dynamics. Because the results in section VIII are expressed in (R) (R) terms of a renormalized persistence length (l ) of WLC we present simulation results for l p p variation for a number of potentials that describe interactions between monomers in section VII. The conclusions of the article are summarized in section IX. II. PRELIMINARY CONSIDERATIONS The pioneering treatment of loop formation dynamics due to Wilemski and Fixman (WF)3,38 has formed the basis for treating cyclization kinetics in flexible polymer chains. Using a generalized diffusion equation for the probability density, ∂P({∂rtN},t) = LFPP({rN},t) − k ( rN )P( rN ,t) ( is a generalized diffusion operator, k is a sink term) for a N-segment FP S { } { } L S polymer,andlocalequilibriumapproximationwithinthesink,WFexpressedthecyclizationtime τ in terms of an integral involving a sink-sink correlation function. From the WF formalism c and related studies it is known that even in the simplest cases (ideal chains or polymers with excluded volume interactions) the validity of the local equilibrium approximation depends on the interplay between τ and the chain relaxation time, τ . If τ τ then the local equilibrium c R c R ≫ approximation is expected to hold because the polymer chain effectively explores the available volume before the monomers at the end (reactive groups) form a contact. In this situation, τ can be computed by considering mutual diffusion of the chain ends in a potential of mean c 2 force (F(R )). For ideal chains, F(R ) = k T logP(R ) 3k TR2/2R where R is the e e − B e ∼ B e e end-to-end distance, R aN1/2 is the mean end-to-end distance, a is the size of the monomer, ∼ T is the temperature, and k is the Boltzmann constant. By solving such an equation sub- B ject to the absorbing boundary condition, Szabo, Schulten, and Schulten (SSS)6 showed that τ = τ N3/2. Simulations4 and theory5 show that if the capture radius for contact formation SSS o is non-zero, and is on the order of a monomer size then τc ∼ hDR2ci ∼ τ1N2ν+1 where ν = 1/2 for Rouse chains and ν 3/5 for polymers with excluded volume, and D is a mutual diffusion c ≈ coefficient. The use of these theories to analyze the dependence of τ onN in polypeptides shows c that the physics of cyclization kinetics is reasonably well described by diffusion in a potential of mean force F(R ) which only requires accurate calculation of P(R ) the end-to-end distri- e e 5 bution function.15,37 For describing interior looping times τ for contact between two interior IC monomers s and s we need to compute P(R , s s ) where R is the distance between 1 2 12 1 2 12 | − | s and s . With P(R , s s ) in hand τ can be computed by solving a suitable diffusion 1 2 12 2 1 c | − | equation. Because the use of F(R ) in computing τ and τ is intimately related to chain relaxation e c IC times it is useful to survey the conditions which satisfy the local equilibrium approximation. By comparingtheconformationalspaceexploredbythechainendscomparedtotheavailablevolume prior to cyclization39 the validity of the local equilibrium approximation in flexible chains can be expressed in terms of an exponent θ = d+g.7 Here d is the spatial dimension, the correlation z hole exponent (des Cloizeaux exponent)40 g describes the probability of the chain ends coming z close together, and z is the dynamical scaling exponent (τ R ). If θ > 1 the local equilibrium R ∼ approximation is expected to hold and τ is determined essentially by the equilibrium P(R ) as c e R a the capture radius. Using the scaling form of P(R ) for small R P(R ) 1 Re g e → e e e ∼ Rd R and R Nν (ν is the Flory exponent) we find τ Nν(d+g). For Gaussian chains ν = 1/2(cid:0) an(cid:1)d c ∼ ∼ g = 0 and hence τ τ N3/2. This result was obtained fifty years ago by Jacobsen and SSS c ∼ ∼ Stockmayer.41 However, in the free-draining case (z = 4, g = 0, ν = 1/2, d = 3), θ < 1 and hence the condition τ τ is not satisfied. In this case τ τ Nzν N2. Thus, for ideal c R c R ≫ ∼ ∼ ∼ Gaussian chains it is likely that τ < τ < τ .42 Indeed, recent simulations show that if the SSS c WF number of statistical segments is large (& 20) then for ideal chains τ N2 which signals the c ∼ breakdownoftheconditionτ τ . Experiments oncyclizationofpolypeptidechainsshowthat c R ≫ τ N3/2a is obeyed for N in the range 10 < N < 20 (see Fig.(5) in Ref.15). Deviations from c ∼ ideal chain results are found for N < 10, either due to chain stiffness15 or sequence variations.43 For polymer chains in good solvents with hydrodynamic interactions (d = 3, g = 5/18, z = 3, and ν = 3/5), θ = 59/54 > 1. Thus, in real chains the local equilibrium approximation may be accurate. For stiff chains bending rigidity severely restricts the allowed conformations especially when the contour length (L) is on the order of the persistence length (l ). Because of high bending p rigidity the available volume is restricted by thermal fluctuations. Clearly in this situation, the chain is close to equilibrium. This may be the case for short DNA segments. In effect these chains satisfy the τ > τ condition which enables us to calculate τ or τ by solving an c R c IC 6 appropriate one dimensional diffusion equation (see below) in a suitable potential of mean force. Effect of chain Stiffness : Many biopolymers are intrinsically stiff and are better described by worm-like chain (WLC) models. The persistence length, which is a measure of stiffness, varies considerably. It ranges from (3 7) ˚A (proteins),15 (10 25) ˚A (ss-DNA19,44 and RNA,45,46) 50 − − nm for ds-DNA. Typically, loops of only a few persistence length form, which underscores the importance of chain stiffness. In order to correctly estimate the loop closure time, consideration of the stiffness in the loop closure dynamics is necessary unless the polymer looping takes place between the reactive groups that are well separated and the chain length L is long. If L l p ≫ (persistencelength),theloopingdynamicswillfollowthescalinglawforflexiblechains. However, at short length scales loop dynamics can be dominated by chain stiffness.16 If the chain is stiff then WLC conformations are limited to those allowed by thermal fluctuations. In this situation, the time for exploring the chain conformations is expected to be less than τ . Thus, we expect c local equilibrium to be a better approximation for WLC than for long flexible chains. Recently, Dua et. al.47 have studied the effect of stiffness on the polymer dynamics based on Wilemski-Fixman formalism and showed, that for free-draining semiflexible chain with- out excluded volume τ N2.2 2.4 at moderate values of stiffness. However, the proce- c ∼ ∼ dure used to obtain this result is not complete, as recognized by the authors, because they 3/2 use a Gaussian propagator G(r,t|r′,0) = 2π r2 (31 φ(t)) exp −23(rr2−(φ1(t)φr(′t))2) which is not h i − h i − valid for WLC. The end-to-end distance d(cid:16)istribution b(cid:17)ecomes(cid:16)a Gaussian (cid:17)at equilibrium, 3/2 lim (r,t r,0) = P (r) = 3 exp 3r2 , which is incorrect for semiflexible chain t→∞G | ′ eq 2πhr2i −2hr2i especially when l L (see Fig(cid:16).2 and(cid:17)Refs.48,4(cid:16)9). (cid:17) p ∼ As an alternative method we include the effect of chain stiffness assuming that local equilib- rium approximation is valid. This is tantamount to assuming that τ > τ which, for reasons c R given above, may be an excellent approximation for WLC.37 In this case we can compute τ c by solving the diffusion equation in a one dimensional potential F(R ) = k T lnP(R ) where e B e − P(R ) is the probability of end-to-end distance distribution for WLC. For the problem of in- e terest, namely, the computation of τ , we generalize the approach of Jun et. al.37 who used IC Kramers theory in the effective potential F(R ) to obtain τ . In general, the time for cyclization e c 7 can be calculated using r 1 L τ = dyeβF(y) dze βF(z) (1) c − D Za Zy where a is the capture (contact) radius of the two reactive groups. We show that Eq.(1) provides accurate estimates of τ , thus suggesting the local equilibrium approximation is guaranteed. c Here, we address the following specific questions: What is the loop formation time between the interior segments in a semiflexible chain? Does the dangling ends (Fig.1-(b)) affect the dynamics of loop formation? How does the effect of interaction between monomer segments (e.g. excluded volume, electrostatic interaction) affect loop closure kinetics in WLC models? III. DISTANCE DISTRIBUTION FUNCTION BETWEEN TWO INTERIOR POINTS A key ingredient in the calculation of the potential of mean force is appropriate distribution functionbetween thetwomonomers thatformacontact. InRefs.49,50 theequilibrium end-to-end (R ) radial distribution function of a semiflexible chain P(R ) was obtained in terms of the per- e e sistence length (l ) and the contour length (L). Despite the mean field approximation employed p in Refs.49,50 the distribution function P(R ) is in very good agreement with simulations.48 The e simplicity of the final expression has served as a basis for analyzing a number of experiments on proteins,51 RNA45 and DNA.52 In this section, we use the same procedure to calculate the distribution function P(R ;l ,s ,s ,L) where 0 < s ,s < L, and R is the spatial distance 12 p 1 2 1 2 12 between s and s . 1 2 For the semiflexible chain in equilibrium we write the distribution function of the distance R between s and s (Fig.1) along the chain contour as 12 1 2 s2 G(R ;s ,s ) = δ(R u(s)ds) 12 1 2 12 MF h − i Zs1 D[u(s)]δ(R s2u(s)ds)Ψ [u(s)] = 12 − s1 MF (2) D[u(s)]Ψ [u(s)] R R MF where u(s) is a unit tangent vector at positionRs. The exact weight for the semiflexible chain is Ψ[u(s)] ∝ exp[ lp Lds ∂u 2] δ(u2(s) 1). The nonlinearity, that arises due to the restric- −2 0 ∂s − tion u2(s) = 1, maRkes th(cid:0)e co(cid:1)mpQutation of the path integral in Eq.(2) difficult. To circumvent 8 the problem we replace Ψ[u(s)] by the mean field weight Ψ [u(s)],53 MF l L ∂u(s) 2 L Ψ [u(s)] ∝ exp[ p ds λ (u2(s) 1)ds δ[(u2 1)+(u2 1)]]. (3) MF −2 ∂s − − − 0 − L − Z0 (cid:18) (cid:19) Z0 The Lagrange multipliers λ and δ, which are used to enforce the constraint u2(s) = 1,54 will be determined using stationary phase approximation (see below). The path integral associated with the weight Ψ [u(s)] is equivalent to a kicked quantum mechanical harmonic oscillator MF with “mass” l and angular frequency Ω 2λ/l . Using the propagator for the harmonic p p ≡ oscillator p 3 Z(u ,u ,s) = πsinh(Ωs) −2 exp( Ω (u2s +u20)cosh(Ωs)−2us ·u0) (4) s 0 p Ω − sinh(Ωs) (cid:18) p (cid:19) and defining Ω Ωlp the isotropic distribution function becomes p ≡ 2 d3k G(R12,s1,s2) = N−1eλL+2δ (2π)3 du0dus1dus2duLe−δu20Z(u0,us1;s1) Z Z × eik·R12−k4λ2|s1−s2|Z(us1 + 2ikλ,us2 + 2ikλ;s2 −s1) × e−δu2LZ(us2,uL;L−s2). (5) i i By writing the distribution function as G(R12,s1,s2) = ∞i dλ ∞i dδexp(−F[λ,δ]) it is clear −∞ −∞ that the major contribution to G (in the thermodynamRic limit RL ) comes from the saddle → ∞ points of the free energy functional [λ,δ], i.e., ∂ = ∂ = 0. The functional [λ,δ] is (see F F F ∂λ ∂δ F Appendix A for details of the derivation) [λ,δ] = (Lλ+2δ) F − 3 sinhΩL 3 λ2 + ln (δ2 +Ω2 +2δΩ cothΩL) ln 2 Ω p p − 2 Q(s ,s ;λ,δ) (cid:18) p (cid:19) 1 2 λ2R2 + 12 (6) Q(s ,s ;λ,δ) 1 2 To obtain the optimal values of λ and δ we first take the L limit and then solve → ∞ stationarity conditions ∂ = ∂ = 0. Technically, the optimal value of δ and λ should be F F ∂λ ∂δ calculated for a given L and then it is proper to examine the L limit. The consequences of → ∞ reversing the order of operation are discussed in Appendix B. Using the first procedure (taking L first) we obtain Q(s ,s ;λ,δ) s s λ in the limit L s s 1, and thus 1 2 2 1 2 1 → ∞ → | − | ≫ ≫ ≫ 9 [λ,δ] becomes F [λ,δ] (Lλ+2δ) F ≈ − 3 eΩL δ 3 3 s s R2 λ + ln[ ( +1)2]+ lnΩ2 + ln | 2 − 1| + 12 2 Ω Ω 2 p 2 λ s s p p 2 1 | − | 3 s s R2 = L Ω λ(1 | 2 − 1| 12 ) 2 − − L s s 2 (cid:18) | 2 − 1| (cid:19) 3 Ω δ 3 + ln[ p( +1)2]+ ln s s 2δ (7) 2 1 2 λ Ω 2 | − |− p where we have omitted numerical constants. The major contribution to the integral over λ and δ comes from the sets of λ and δ which pass the saddle point of a stationary phase contour on the Re plane. Since the term linear in L dominates the logarithmic term in l even when p {F} L/l (1), the stationary condition for λ can be found by taking the derivative with respect p ∼ O to λ by considering only the leading term in L (cf. see Appendix B for details). The stationarity condition leads to λl 3 1 p Ω = = (8) p r 2 41 |s2−s1|r2 − L R where r = 12 with 0 < r < 1. Similarly, the condition for δ can be obtained as s2 s1 | − | 3 δ = Ω . (9) p 2 − Determination of the parameters λ and δ by the stationary phase approximation amounts to replacing the local constraint u2(s) = 1 by a global constraint u2(s) = 1.54 Finally, the h i stationary values of λ and δ in the large L limit give the interior distance distribution function: N 9 s s 2 1 G(R ,s s ) = exp( | − | ). (10) 12 2 − 1 (1− |s2−Ls1|r2)9/2 −8lp(|s2−Ls1|)(1− |s2−Ls1|r2) Themean-fieldapproximationallowsustoobtainasimple expression fortheinternalsegment distance distribution function. The previously computed P(R )49 can be retrieved by setting e s s = L. The radial probability density, for the interior segments, in three dimensions, for 2 1 | − | semiflexible chains is r2 3t P(r;s s ,t) = 4πC exp( ) (11) 2− 1 (1 |s2−s1|r2)9/2 −4(|s2−s1|)(1 |s2−s1|r2) − L L − L 10 where r = R / s s R/ s s and t = s s /l with l = 2l . The normalization 12 | 2 − 1| ≡ | 2 − 1| | 2 − 1| 0 0 3 p 1 constant C is determined using P(r,s;t)dr = 1. The integral is evaluated by the substitution 0 |s2−s1|r = x to yield R L √1+x2 q C = 1 (|s2 −s1|)3/2 x0dxx2(1+x2)e α(1+x2) −1 − 4π L (cid:18)Z0 (cid:19) 4 s s = πα 7/2(| 2 −L 1|)3/2[−2√αx0e−α(1+x20)(15+2α(6 +5x20 +2α(1+x20)2)) − 15 + α2e α√πerf[√αx ](1+3α 1 + α 2)] 1 (12) − 0 − − − 4 whereα = 4|s23−ts1|,x0 = L|s2s−2s1s|1 ,anderf(x)istheerrorfunction. Thepeakinthedistribution L −| − | q function is at η + η2 +14 r = (13) max s 7 s s /L | 2p− 1| where η = 5 3t . For s s = L, r 0 as t and r 1 as t 0. 2 − 4|s2−s1| | 2 − 1| max → → ∞ max → → L In Fig.2 we compare the distribution functions P(R ) and P(R ,s s ). When e 12 2 1 − s s /L = 1 Eq.(11) gives the end-to-end distribution for semiflexible chains. By adjusting 1 2 | − | the value of t (or equivalently l ) we can go fromflexible to intrinsically stiff chains. As the chain 0 gets stiff there is a dramatic difference between the P(r; s s ,t) and P(r; s s = L,t) 1 2 1 2 | − | | − | (see Fig.2-(b)). Contact formation between interior segments are much less probable than cyclization process (compare the green and red curves with the black in Fig.2-(b)). Physically, this is because stiffness on shorter length scales ( s s /L < 1) is more severe than when 1 2 | − | s s /L (1). However, when the chain is flexible (large t) the difference between the 1 2 | − | ∼ O probability of contact between the interior segments and cyclization is small (Fig.2-(a)). In the limit of large t( L/l ) the Hamiltonian in Eq.(3) describes a Gaussian chain for which the p ∝ distance distribution between interior points remains a Gaussian. However, if excluded volume interactions are taken into account there can be substantial difference between P(x, s s ,t) 1 2 | − | and end-to-end segment distribution even when t is moderately large.