Thermal breakage of a semiflexible polymer: Breakage profile and rate

Understanding fluctuation-induced breakages in polymers has important implications for basic and applied sciences. Here I present for the first time an analytical treatment of the thermal breakage problem of a semi-flexible polymer model that is asymptotically exact in the low temperature and high friction limits. Specifically, I provide analytical expressions for the breakage propensity and rate, and discuss the generalities of the results and their relevance to biopolymers.


Introduction
From man-made materials to biopolymers, semi-flexible polymers are ubiquitous in science and engineering, and better understanding of their stability is of high importance. A semi-flexible polymer can naturally be broken by stretching or bending via external forces, but thermal fluctuations alone will also induce breakage. Fluctuation-induced breakage is particularly relevant in the biological world as many biopolymers are stabilized by hydrophobic interactions and hydrogen bonds between proteins, which are relatively weak compared to covalently bonded synthetic polymers. Indeed, for amyloid fibrils, a kind of polymer implicated in numerous human diseases including Parkinson's and Alzheimer's [1], it has been advocated that thermal breakage could be a key mechanism underlying amyloid fibril proliferation [2][3][4]. Despite the importance of understanding how polymers break, it is surprising that the basic physics remains to be elucidated. For instance, there is currently no concensus on the breakage profile [5][6][7]: some have advocated that breakage happens predominantly in the middle of the polymer [8] while others have assumed uniform breakage propensity [2,3,[9][10][11]. This confusion is due in part to the fact that most existing results rely Content from this work may be used under the terms of the Creative Commons Attribution 3.0 licence. Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI. on considering simplified polymer models in one dimension [12][13][14] or numerical simulations [15]. Here I will provide for the first time analytical expressions of the breakage profile and rate of a highly-rigid polymer in the high friction (highly damped) and low temperature regimes.

Minimal model
I will employ a bead-and-stick representation of a semiflexible polymer in which the sticks are massless and the beads experience isotropic thermal fluctuations (see figure 1(a)). In the highly damped regime, the equations of motion (EOM) for the ith bead is where ζ is the drag coefficient for the beads and η i denotes Gaussian noise terms with zero means and unit variance. Also, the tensile and bending rigidities are enforced by two energy potentials U and V : where r k,k−1 ≡ | r k,k−1 | ≡ | r k − r k−1 | and θ k ≡ arccos( r k,k−1 · r k+1,k r k,k−1 r k+1,k ).
Thermal fluctuations induce tensile and bending strains on the polymer and the polymer is broken if the bond is stretched or the angle is bent beyond certain thresholds. Since I have previously analysed thermal breakage by fluctuations-induced tensile strain for a polymer on a one dimensional track [13], I will focus here on breakage by bending, and I will comment on how the results are modified if both breakage by stretching and bending are possible in Discussion. With regard to the breakage criterion, I assume that the polymer is broken if one of the angles is greater than some material-dependent threshold angle c (figures 1(b) and (c)), which is expected to be much less than 1 for a highly rigid polymer. With this definition, we are ready to pose the two questions that will be answered in this paper: (i) What is the rate for one of the angles to be pushed beyond c due to fluctuations? (ii) What is the breakage propensity as a function of the monomer position in the polymer?

Coordinate transformation
Since the breakage criterion concerns the set of angles θ i , it is more natural to consider the EOM of the θ i instead of r i . Using the chain rule, where the second equality comes from equation (1), and the summation index is now restricted to [i − 1, i + 1] since variations of beads outside this range do not affect θ i . Focusing on the highly rigid polymer so that c is expected to be small, the EOM can be rewritten as (appendix A): where is the energetically optimal distance between neighbouring beads, L is the mobility matrix given in equation (A.4), and the noise terms are now defined by ξ i (t) = 0 and ξ i (t)ξ j (t ) = 2(L −1 ) ij δ(t − t ). In other words, the set of angles may be viewed as under the influence of the potential energy U tot ({θ }) ≡ k U(θ k ) and a thermal heat bath with the nondiagonal mobility matrix L [16].

Two dimensions
Equipped with the EOM for θ i , let us now calculate the escape rate and the breakage profile by employing the quasi-static approximation. While the arguments below will be heuristic, the results can be shown to be exact in the asymptotic limit of β → ∞ with β ≡ 1/k B T [17,18]. In this quasi-static limit, we assume that the probability distribution is normalized to one and the distribution is equilibrated to be at the Gibbs state: (c)). Since at low T , the distribution is highly centred around θ k = 0, P ({θ }) is well approximated as N M e −βA k θ 2 k /2 where A = U (0). By integrating this multi-dimensional Gaussian distribution, we find that for a (M + 2)-bead polymer, N M = (βA/2π) M/2 . As T → 0, breakage is highly improbable and is dictated by the configurations where one of the angle is at c (say θ 1 ∼ c ) while the rest remain close to 0. Around this breakage boundary, we can expand the energy landscape in a Taylor series. As a result, the probability distribution at this breakage boundary P 1 ≡ P ({θ 1 ∼ c ∩ θ k>1 ∼ 0}) can be written as where E = U( c ) − U(0) is the 'Arrhenius' activation energy term, and b ≡ |U ( c )| is proportional to the gradient of the potential energy at the breakage point ( figure 1(b)).
To calculate the breakage rate R 1 , we integrate the probability flux across the breakage boundary θ 1 = c . The probability flux at this boundary is given by the expression and hence where the second equation is obtained by substituting into equation (8) the expression for P 1 in equation (6). Since the diagonal elements in the mobility matrix L are identical and the probability distribution at distinct breakage boundaries P m , 1 m M, are the same, all angles have identical breakage rate. In the asymptotic limit of T → 0, the likelihood of reaching a particular breakage boundary is proportional to the corresponding breakage rate [17,18]. We can therefore conclude that the breakage propensity is uniform along the polymer.
We can also now calculate the average breakage rate of the polymer. Since there are M angles and there are two breakage positions per angle in 2D (θ m = ± c ), the average breakage rate of the polymer is where once again A = U (0) and b = |U ( c )|. The predictions on the breakage propensity and rate are supported by Brownian dynamics simulations as shown in figure 2(a).

Three dimensions
Going from 2D to 3D introduces a new degree of freedom as the polymer can rotate around its longitudinal axis in 3D. Using the longitudinal axis of the polymer as the z-axis, θ k corresponds to the polar angle while the azimuth angle is denoted by φ k . Due to this extra degree of freedom, for a (M + 2)-bead polymer, the normalisation factor at the quasistatic distribution is modified to Once again, defining P 1 as P ({θ 1 ∼ c ∩ θ k>1 ∼ 0}) which has the same expression as in equation (6), the breakage rate R 1 can be calculated as Therefore, for a (M + 2)-bead polymer in 3D, the breakage propensity is again uniform along the polymer and the average breakage rate is: These analytical predictions are supported by Brownian dynamics simulation ( figure 2(b)).

Discussion
We have seen that for the minimal semi-flexible polymer model considered here, the answers to the two questions posed earlier are analytical tractable in the T → 0 limit. I will now discuss the generalities and limitations of the results.

Beyond thermal systems
Although the model formulation focuses on polymers at thermal equilibrium, the results remain valid as long as the dynamics of the polymer is well approximated by equation (1). In other words, the results also apply to polymers enclosed in a volume that is at 'local' equilibrium, or under active fluctuations of the form depicted in equation (1) [19,20].

Low T limit versus high energy barrier limit
Mathematically, these two limits are not reversible: In the first limit, there is only one small parameter in the equation of motions (equation (1)), while in the high energy barrier limit, both the fluctuation strength and the distance from the bottom of the well to the escape boundary become effectively small. Physically, since the units of time and length can be set arbitrarily, the high energy barrier and the low temperature limits are in principle interchangeable with the caveat that the high friction assumption has to remain valid. This leads to the next comment.

Drag coefficient
The results presented apply to the high drag regime where the inertia term is neglected. The drag coefficient may be seen as an effective measure of the friction due to both the internal degrees of freedom within the beads (e.g. internal degrees of freedom of a protein), and bead-solvent interactions. These two effects combined seem to lead to the general validity of considering biomolecular kinetics in the highly damped regime [21].   (10) and (13)) to the simulation results for a 9-bead polymer in 2D (a) and in 3D (b) as β E varies. The ratio in both cases approaches one as β E increases, thus showing the convergence of analytical and simulation results as temperature decreases. The inset plots show the histograms of the breakage locations for the case β E = 15. Performing the Pearson's chi-squared test on these histograms with the uniform distribution as the null hypothesis indicates that the fluctuations observed can arise from a uniform distribution with probabilities 89% and 70% for the 2D and 3D cases, respectively. Hence, there is little reason to reject the null hypothesis. The energy functions employed and simulation procedure are detailed in appendix B.

Threshold on bending
The analysis presented here assumes that the threshold bending angle c is small. For biopolymers, I am only aware of one paper on microtubules that provides access to this parameter [22]. Specifically, it was observed that upon bending, a microtubule forms an arc with curvature of around 1 µm −1 before breaking. Since the size of the tubulin dimers making up the microtubule is around 10 nm, one may estimate that for microtubules, c is in the order of 0.01 rad.

Extensile versus bending breakage
I have omitted discussion of breakage by extensile fluctuations here because the results have already been derived previously [13], albeit with the extra assumption that the polymer is on a one-dimensional track. But since the bending fluctuations are orthogonal to the stretching fluctuations. In the low T limit, breakage by stretching of a polymer in higher dimensions is analytically identical to the 1D case. If both extensile and bending breakage events are possible, then the one with a lower Arrhenius activation energy will dominate. If both breakage events have the same activation energy, then the relative proportion of bending versus tensile breakage will be given by the ratio of the corresponding prefactors [17,18].

Other types of potential energy
The analytical results derived in this work apply to any potential function U that has a unique minimum and a nonvanishing gradient at the threshold angle c . For different kinds of potential functions, such as the single-hump potential function usually employed in the discussion of Kramers escape rate, one needs to solve the corresponding first passage time problem of the multidimensional Langevin equation with correlated noise terms (equation (5)). Analytical treatments for such a problem are rare [23] and it would be highly interesting to see whether the breakage propensity will be modified when the bonds are governed by different kinds of potential functions.

Internal structure and end effects of a polymer
The model polymer considered ( figure 1(a)) is certainly a drastically simplified version of a real polymer, but in the spirit of a ghost chain in polymer physics, each bead may be seen as a unit of polymer structure with internal dynamics that are decoupled from the bending dynamics considered here [16]. If such decoupling is valid, the prediction of uniform breakage propensity along the body of the body of the polymer should hold true. However the ends of the polymers may be subject to different kinds of energy function. For instance, in the case of protein amyloid fibrils where each fibril is a bundle of filaments, the filaments close to the ends of the fibril may not be as tightly bound and so the bending rigidity at the ends may differ from that in the middle of the fibril. As a result, the breakage propensity close to the ends of the polymer may differ from that in the main body of the polymer. Another complication is the possibility of having a rugged energy landscape that connects two neighbouring beads [24,25]. In this case, the drag coefficient may need to be modified using the theory of diffusion on rugged landscape [26].

Summary and outlook
I have considered thermal breakage of a semi-flexible polymer model and have obtained analytical results in the highly rigid, highly damped and low temperature limits. My calculations indicate that a semi-flexible polymer satisfying these conditions has uniform breakage propensity and analytical expressions of the breakage rates in 2D and 3D were provided. All analytical results were verified by Brownian dynamics simulations. The generalities of the results and their relevance to biopolymers were discussed. Future work on this problem will include a thorough analysis of the thermalization kinetics of polymerization since both the equilibrium configuration [27,28] and the breakage kinetics are now known.
Other interesting directions are the calculations of higher order corrections, and the investigation of polymer breakage under nonequilibrium and anisotropic fluctuations, such as under shear flows [29,30] and sonications [31,32], which constitute two standard experimental procedures in investigating biopolymer selfassembly.

Acknowledgments
I thank P Sartori (MPI-PKS) for helpful comments.

Appendix A. EOM for the θ k
Let us first focus on the following terms in equation (4): The terms of the form ∇ r k θ i · ∇ r k V (r h,h−1 ) are of order O( c ), To see this, consider the particular term ∇ r k θ k · ∇ r k V (r k,k−1 ) in 2D. Without loss of generality, assume r k,k−1 is along the x-axis, ∇ r k θ k = −1 θ kx + −1ŷ + O(θ 2 k ) and ∇ r k V (r k,k−1 ) =x. Since θ k < c before breakage, ∇ r k θ k · ∇ r k V (r k,k−1 ) = O( c ). Other terms of the same form can similarly be shown to be of the same order, which are small for rigid polymer where c 1. Physically, what it means is that the fluctuations that contribute to bond extension is orthogonal to the fluctuations contributing to bending in the small angle limit. The same arguments apply in 3D.
To O( c ), we are thus left with the following terms in equation (A.2): The only remaining terms in equation (4) that we have to consider are from the coupled noise terms: ∇ r k θ i · η k ≡ G ik .
To O( c ), G is a M × (M + 2) matrix of the form The matrix G is nondiagonal because there are M variables (θ k ) and M + 2 noise terms (η k ). One could therefore reduce the number of the noise terms by 2 by introducing a new set of M noise terms ξ k with statistics depicted below: ξ k (t)ξ h (t ) = 0 for h = k, k ± 1, k ± 2. (A. 10) In other words, the covariance matrix of ξ k is exactly the motility matrix L, this explains the EOM in equation (5).

Appendix B. Simulation procedure
The simulation is based on numerically integrating the set of stochastic differential equations shown in equation (1) by using the following update: where the entries of the vector g are Gaussian distributed random variables with zero mean and unit variance. The time increment t is set to be 2 × 10 −6 , the drag coefficient ζ is one and the energy function U and V in H are of the form: where A = 60 and B = 400. The polymer is always initialized in its lowest energy state and the simulation stops if one of the angles θ is beyond c = 0.1. One thousand sample runs are performed at each distinct k B T .