Abstract:We study intrachain dynamics of intrinsically disordered proteins, as man-ifested by the time scales of loop formation, using atomistic simulations, experiment-parametrized coarse-grained models, and one-dimensional theories assuming Markov or Markov dynamics along the reaction coordinate. Despite the generally non-Markov character of monomer dynamics in polymers, we find that the simplest model of one-dimensional diffusion along the reaction coordinate (equated to the distance between the loop-forming monomers) well captures the mean first passage times to loop closure measured in coarse-grained and atomistic simulations, which, in turn, agree with the experimental values. This justifies use of the one-dimensional diffu-sion model in interpretation of experimental data. At the same time, the transition path times for loop closure in longer polypeptide chains show significant non-Markov effects; at intermediate times, these effects are better captured by the generalized Langevin equation model. At long times, however, atomistic simulations predict long tails in the distributions of transition path times, which are at odds with both the one-dimensional diffusion model and the generalized Langevin equation model.
This project was published in theJournal of Physical Chemistry Bon April 7th 2020 [171]. The research for it was done by Rohit Satija and Dmitrii Makarov. My contribution was to run and make available end-to-end trajectories extracted from the polymer model developed in this thesis, and to help in revising the manuscript (as did Jörg Enderlein and Atanu Das). Atanu Das has further made available end-to-end trajectories from MD simulations of GS-repeats.
In this project, the end-to-end dynamics in GS-repeats are compared from different viewpoints:
Experimentally measured loop formation rates and trajectories from Brownian dynamics (BD) and molecular dynamics (MD) simulations are compared with (i) one-dimensional Markovian diffusion and (ii) one-dimensional non-Markovian diffusion. Trajectories from our Brownian dynamics model are particularly valuable in this context because – unlike other coarse-grained polymer models – it is specifically parametrized to match experimental data for GS-repeats and – unlike experimental results – yields complete trajectories and not just contact rate values.
From the perspective of this thesis, the most interesting result of this project is that one-dimensional descriptions fail in reproducing certain dynamical features of the end-to-end dy-namics described by the BD and MD simulations. Nevertheless, the mean first passage time (inverse contact rate) is reasonably well captured by the one-dimensional picture. This match is in agreement with our finding that SSS-theory (which assumes one-dimensional Markovian end-to-end dynamics, see page 142) correctly predicts the scaling behavior of the contact rate as a function of chain length. For further reading, we refer to the published manuscript [171].
APPEN
A
C ONVERGENCE OF DISCRETE CURVATURES TO THE WORM - LIKE CHAIN
The context of this appendix is explained in section 2.3.2 on page 60. Here, we consider the convergence of the end-to-end distance’s probability distribution in discretized worm-like chains towards the continuous limit. This convergence is being compared between three definitions of discrete curvature, which are
(A.1)
κ1=2 tanθ/2
∆s , κ2=2 sinθ/2
∆s and κ3= θ
∆s
whereθis the bending angle at the considered chain vertex and∆sis the arc length associated with it. For all chains considered here, the contour length equals 1, and the persistence length equals 1/3. The corresponding bending stiffnesses of the discrete chains were chosen according to equations 2.181 and 2.182 (page 68) for all discretizations and curvature definitions respectively.
The probability distributions were obtained as histograms using MC sampling. An analytical approximation for the “true” end-to-end probability distribution in a continuous worm-like chain was taken from reference [172] whose author’s made available aMathematicanotebook, which I have used. The following three figures show this analytical approximation and the numerically sampled probability distributions for chains consisting ofN∈{5, 10, 20, 40} vertices.
It turns out that the tangent-based curvature definition yields the fastest convergence towards the continuous worm-like chain limit.
0.0 0.2 0.4 0.6 0.8 1.0 0.0
0.5 1.0 1.5 2.0 2.5 3.0
end - to - end distance r
probability density
analytic N=40 N=20 N=10 N=5
Figure A.1: The calculated end-to-end probability distribution for the tangent-based curva-ture for chains consisting of 5, 10, 20 and 40 vertices. The black curve shows an analytical approximation of the continuous worm-like chain limit.
0.0 0.2 0.4 0.6 0.8 1.0 0.0
0.5 1.0 1.5 2.0 2.5 3.0 3.5
end - to - end distance r
probability density
sine
analytic N=40 N=20 N=10 N=5
Figure A.2: The calculated end-to-end probability distribution for the sine-based curvature for chains consisting of 5, 10, 20 and 40 vertices. The black curve shows an analytical approximation of the continuous worm-like chain limit.
0.0 0.2 0.4 0.6 0.8 1.0 0.0
0.5 1.0 1.5 2.0 2.5 3.0
end - to - end distance r
probability density
analytic N=40 N=20 N=10 N=5
Figure A.3: The calculated end-to-end probability distribution for theθ-based curvature for chains consisting of 5, 10, 20 and 40 vertices. The black curve shows an analytical approximation of the continuous worm-like chain limit.
APPEN
B
P AWULA ’ S THEOREM
Proofs of Pawula’s theorem, similar to the one presented here, can be found in references [97]
and [104]. The context of this appendix is explained in section 2.2.3 on page 36. There, we found that the jump momentsαk(x) uniquely define a Markov process from the Kramers-Moyal expan-sion (in orders ofk) of its Master equation. Now we proof that the Kramer-Moyal expansion either ends after one term, two terms, or it never ends. This fact is known as Pawula’s theorem [104].
The generalized Schwartz inequality reads
(B.1) ³Z
f(z)g(z)ρ(z) dz´2
≤
³Z
f2(z)ρ(z) dz´
·
³Z
g2(z)ρ(z) dz´
whereρ(z) is some nonnegative function and f(x) andg(x) are arbitrary functions. In the context of section 2.2.3, we chooseρ(z)=p(z,t+τ|x,t) which is by definition nonnegative, f =(z−x)n and g(z)=(z−x)n+mfor an arbitrary, fixed xandn,m∈N0. Inserting these choices into the Schwartz inequality yields
(B.2) (Ψ(2n+m))2≤Ψ(2n)·Ψ(2n+2m) where
(B.3) Ψ(k)=
Z
(z−x)np(z,t+τ|x,t) dz
=α(k)τ+O(τ2) .
As we found in the main text, this relation between Ψ(k) and α(k) only holds for k≥1. For n=m=0, as well as form=0 whilen≥1, the inequality above is trivially fulfilled.
(B.4)
Ψ ≤ ·Ψ
⇔ [α(m)τ+O(τ2)]2≤α(2m)τ+O(τ2) |:τ, lim
τ→0
0≤α(2m)
which is interesting but also clear from the definition of the jump moments. In conclusion, we only need to considern,m≥1, for which the last line of equation B.3 holds true. Inserting it into the inequality, then dividing byτ2and considering the limitτ→0, we arrive at
(B.5) (α(2n+m))2≤α(2n)·α(2n+2m) for n,m≥1 . From this inequality follow two conclusions:
(B.6) α
(2n)=0 ⇒ α(2n+m)=0
and α(2n+2m)=0 ⇒ α(2n+m)=0 .
This means that if one finds some even moment α(2r)(x)=0,r∈N, all higher moments also vanish (first line).
Furthermore, since we may choose any combination (n,m)∈{(1,r−1), (2,r−2), . . . , (r−1, 1)} in the second line such that n+m=r, we learn thatαr+1(x),αr+2(x) . . . ,α2r−1(x) also vanish. Among these latter moments, there might then be more even moments which we now know to equal zero as well, and hence the whole argument can be repeated. The reader is highly encouraged to e.g. assumeα8=0 and see what follows. It turns out that only three scenarios are possible: The Kramers-Moyal expansion (ii) ends after the first term (deterministic dynamics of the Liouville equation), (ii) ends after the second term (diffusive dynamics) or (iii) has an infinite number of nonzero terms including all even moments (discontinuous jump processes).
A
PPENC
E NTROPIC PSEUDOFORCE
The physical context of this appendix is explained in section 2.2.6 on page 48. Here, we consider the statistical properties of a stochastic process whose dynamics is governed by an unconstrained Langevin equation in 3N dimensions which reads
(C.1) dt~r=T−1D·£
−~∇E¤ +p
2B·~ξ .
The Boltzmann distribution p∝exp(−E/T) is sampled as long as the relation D=B·BT is fulfilled. For convenience, we choose the unity matrix D=T13N×3N (absorbing the physical dimension into the time scale). The dynamics shall closely approximate the constraints
(C.2) gk(~r)=0 where k=1, . . . ,M
due to the presence of strong restoring forces−~∇Eacting on the system. Fluctuations around the minimum ofEare very small such that it becomes a quadratic function of the (infinitesimal) deviationsδgk from the constraints. The resulting Langevin equation then reads
(C.3) dt~r= −αXM
j=1
δgj~∇gj+. . .
where the abbreviated terms remain constant as the elastic constantα→ ∞. According to the chain rule, the dynamics of each of thek=1, . . . ,Mdeviations satisfies
(C.4)
dtδgk=~∇gk·dt~r
=~∇gk· Ã
−αXM
j=1
δgj~∇gj+. . .
!
= −αXM
j=1
~∇gk·~∇gj·δgj+. . . .
(C.5) dtδ~g= −αG·δ~g+. . . where Gk j=∂igk·∂igj .
Next we transform this equation into the eigensystem of matrix Gby multiplying both sides (from the left) with the orthonormal matrixS=(ST)−1containing its eigenvectors:
(C.6)
dt~v= −αGˆ ·~v+. . . where ~v=S·δ~g
and Gˆ =S·G·ST .
Since ˆGis diagonal, the dynamics of the components of~voccur independently, all experiencing harmonic potentials
(C.7) Ek=α
2λkv2k
where theλkare the eigenvalues (diagonal entries) of ˆG. The probability (Boltzmann) distribution forvkis hence a Gaussian with zero mean and variance
(C.8) σ2k∝ 1
λk
.
The entropy of a Gaussian with varianceσ2equals 1
2ln(2πeσ2), and therefore the entropySkof the fluctuations ofvk are
(C.9) Sk= −1
2ln(λk)+const. .
It follows that the total entropySC of the fluctuations of the vector~vis
(C.10)
For the constrained Langevin equation (main text) we only require the gradient of this entropy with respect to the Cartesian coordinates~r. For this reason the neglected terms (includingαand T) in the calculation above are indeed inconsequential as they vanish upon differentiation.
The matrixGdoes depend on~r, however, because the constraintsgkthemselves do, and therefore the force term
(C.11) ~∇SC= −1
2~∇ln detG is all we need to report back to the main text.
Inextensibility constraints The constraint entropy’s gradient can be evaluated for specific constraints. In this thesis in particular, we consider inextensibility constraints for a discrete chain consisting ofN beads:
Making use of the fact that only the gradients of the next diagonals ofGare nonzero, the gradient of ln detGcan be calculated using the following identity [66]:
(C.14) The last factor can easily be evaluated analytically, and the first factor (of the final sum) requires knowledge of the next diagonal elements of the inverse matrixG−1, whose numerical evaluation is straightforward (but unfortunately depends on the Cartesian coordinates~r).
Bound-end constraint Binding one end of a chain to a fixed point in space (the origin) can be achieved by implementing the constraint
(C.15) g0(~r)= |~r1| −1=0 .
For a chain consisting ofN beads, this additional constraint turns theN−1×N−1-dimensional matrixGinto an N×N-dimensional matrix. Apart from its size, it is almost identical to the matrix for the inextensibility constraints alone. The only difference is that its first element is G1,1=1, and the first “tangent vector” must be taken to be~r1/|~r1|(assuming the chain end is bound to the origin of our coordinate system).
APPEN
D
A NGULAR CORRELATIONS IN WORM - LIKE CHAINS
The context of this appendix is explained in section 2.4 on page 66. The following calculation can also be found in the textbook by Doi and Edwards [115].
For a discrete WLC with a bending energy present, the average cosine of a bending angleθcan be directly calculated from the Boltzmann distribution via
(D.1) 〈cosθ〉 =
Z 1
−1
y·exp µ
− σ 2d
³
2 tan³arccosy 2
´´2¶ dy Z 1
−1
exp µ
− σ 2d
³
2 tan³arccosy 2
´´2¶ dy
where y:=cosθ. Let us now consider the direction of a tangent vector~ti. Since no energy is associated with rotations of~tiaround~ti−1, all such rotational (torsion) angles are equally probable and we have
(D.2) 〈~ti〉 = 〈cosθ〉~ti−1 .
Multiplying both sides with~tk, and then averaging over its orientations yields (D.3) 〈~ti·~tk〉 = 〈cosθ〉 〈~ti−1·~tk〉 .
This is a recursion rule for the vectorial correlations〈~ti·~tk〉wherekis fixed andiis the variable index.
(D.4) 〈~ti·~tk〉 = 〈cosθ〉|i−k| ,
which is an exponential decay along the chain. Defining the persistence lengthlP= − d ln〈cosθ〉, we may alternatively write
(D.5) 〈~t1·~tk+1〉 =exp
µ
−kd lP
¶
wheredis the bond length of the discrete WLC.
The analogous behavior of the continuous WLC immediately follows from this result by consider-ing the limitd→0 while holdinglP fixed, in which casek·d turns into the continuous variables called arc length.
APPEN
E
S OLVING CONSTRAINTS
The context of this appendix is explained in section 4.2 on page 87. Here, we merely insert the equations of motion
(E.1) d~ri=
N
X
j=1
Di j·h
~Fucj +λj~tj−λj−1~tj−1i dt where
(E.2) ~Fuc= −~∇(E−SC)+D−1·(~∇T·D)+p
2 (BT)−1d~W/dt into the constraints
(E.3)
0=!~tT0·d~r1 0=!~tTµ·¡
d~rµ+1−d~rµ+1¢
, µ=1, . . . ,N−1 , with the goal of rewriting them in the simple form
(E.4)
N−1
X
ν=0
Mµνλν=vµ µ=0, . . . ,N−1 . Before we start, let us define the following 3×3-matrices:
(E.5)
Gµν:=Dµν−Dµ,ν+1, µ=1, . . . ,N, ν=1, . . . ,N−1 Hµν:=Dµ+1,ν−Dµν, µ=1, . . . ,N−1, ν=1, . . . ,N
Cµν:=Hµν−Hµ,ν+1, µ=1, . . . ,N−1, ν=1, . . . ,N−1
and remark thatµcounts the constraints andνcounts the corresponding Lagrangian multipliers λ.
(E.6) Hence we can read off that
(E.8)
The most general case of inextensibility as well as bound-end constraints being present is thus treated. The case of a freely swimming, inextensible chain immediately follows from the line of reasoning above – one must only setλ0=0 and ignore the corresponding constraint (µ=0) such thatµ=1, . . . ,N−1 in equation E.4. Therefore, all entries with an index which is zero can be discarded from equation E.8 and the inextensibility constraints may be written as
(E.9)
N−1
X
ν=1
Mµνλν=vµ whereµ,ν=1, . . . ,N−1,
(E.10)
vµ= −
N
X
j=1
~tTµ·Hµj·~Fucj and Mµν=~tTµ·Cµν·~tν .
APPEN
F
D IFFUSION COEFFICIENT – BENCHMARKS
The context of this appendix is explained in section 4.7 on page 101. Here, we show mean square displacements (MSD) as functions of time for simulated model chains. These chains have two different chain lengths;N=5 andN=16. For both chain lengths, we investigate the three cases of (i) no hydrodynamic interactions (HI) (a=0), (ii) strong HI (a=1), and (iii) strong HI (a=1) with a fluorophore attached. The latter fluorophore mimics the real fluorophore Atto655 and has hydrodynamic radiusa1=6/3.8 and a fixed bond length ofd1=7/3.8. The bending stiffness of all chains is set toσ=1 and the excluded volume radius toREV=p
2 /4. For each simulation, we display the MSD of the chain’s center of mass, as well as that of the first bead (which is the fluorophore if applicable).
We observe again that hydrodynamic interactions make the beads fluctuate more cooperatively.
Here, this fact is reflected by the offset between the first bead’s MSD and the diagonal line being smaller in case of hydrodynamic interactions. In general, the center of mass of a simulated chain also undergoes weak intrachain dynamics, as all beads are weighted equally in the center of mass, but some are more mobile than other (the free ends for instance). Hence the offset between the center of mass’ MSD and the diagonal line is also allowed.
All simulated MSDs were divided by the theoretically predicted diffusion coefficient (see equa-tion 4.22, page 96). Their slope is thus supposed to converge towards 1. “BD” stands for Brownian dynamics, i.e. simulated values, and “MC“ stands for Monte-Carlo sampling, i.e. calculated values.
0 10 20 30 40 0
10 20 30 40
t
< R
2> / 6 D
BD-center of mass BD-one bead MC-theory
0 10 20 30 40
0.0 0.5 1.0 1.5 2.0
t
d / dt < R
2> / 6 D
BD-center of massBD-one bead MC-theory
Figure F.1:N=5,a=0. Top: Simulated MSDs from Brownian dynamics simulations (BD) of the chain’s first bead (blue) and its center of mass (black), as well as the theoretically predicted long-term behavior (red), over time. Bottom: slopes of top plot.
0 10 20 30 40 0
10 20 30 40 50
t
< R
2> / 6 D
BD-center of mass BD-one bead MC-theory
0 10 20 30 40
0.0 0.5 1.0 1.5 2.0
t
d / dt < R
2> / 6 D
BD-center of massBD-one bead MC-theory
Figure F.2:N=5,a=1. Top: Simulated MSDs from Brownian dynamics simulations (BD) of the chain’s first bead (blue) and its center of mass (black), as well as the theoretically predicted long-term behavior (red), over time. Bottom: slopes of top plot.