Buckling a Semiflexible Polymer Chain under Compression

Instability and structural transitions arise in many important problems involving dynamics at molecular length scales. Buckling of an elastic rod under a compressive load offers a useful general picture of such a transition. However, the existing theoretical description of buckling is applicable in the load response of macroscopic structures, only when fluctuations can be neglected, whereas membranes, polymer brushes, filaments, and macromolecular chains undergo considerable Brownian fluctuations. We analyze here the buckling of a fluctuating semiflexible polymer experiencing a compressive load. Previous works rely on approximations to the polymer statistics, resulting in a range of predictions for the buckling transition that disagree on whether fluctuations elevate or depress the critical buckling force. In contrast, our theory exploits exact results for the statistical behavior of the worm-like chain model yielding unambiguous predictions about the buckling conditions and nature of the buckling transition. We find that a fluctuating polymer under compressive load requires a larger force to buckle than an elastic rod in the absence of fluctuations. The nature of the buckling transition exhibits a marked change from being distinctly second order in the absence of fluctuations to being a more gradual, compliant transition in the presence of fluctuations. We analyze the thermodynamic contributions throughout the buckling transition to demonstrate that the chain entropy favors the extended state over the buckled state, providing a thermodynamic justification of the elevated buckling force.


Introduction
The fundamental physics of instability and transition is common to a broad range of physical phenomena. Important examples include phase transitions in thermodynamic systems, structural transitions in elastic objects under deformation, and instabilities in non-linear dynamical systems. Structural transitions in elastic objects are typically analyzed with an emphasis on geometry, (e.g., object size and shape), material properties, and external perturbations (e.g., applied force and torque), with less focus on fluctuations (either mechanical or thermal). The conditions and nature of a thermodynamic phase transition are fundamentally altered by thermal fluctuations, but the thermodynamic limit typically neglects the system size and geometry. Applied and fundamental studies of physical systems are increasingly focused on molecular structures whose size scale is sufficiently small where both geometrical considerations and thermal fluctuations are central to their manipulation and stability. Here, we consider a physical problem that is both practically important and marries both the mechanical contributions and thermal fluctuations that are central to instability in elastic objects and thermodynamic systems at the microscope.
The buckling of elastic objects arises in many practical applications and has a long history of theoretical analysis. The Euler buckling instability of a rigid rod under a compressive load is the simplest example of elastic instability, where the rod buckles if the compression force exceeds a critical value. Filament buckling under compressive force is widely observed in biopolymers, such as actin filaments, microtubules or DNA [1,2], as well as in synthetic systems like carbon nanotubes polymeric fluids [3] and polymer brushes [4,5]. While fluctuations are crucial in these examples, other physics like excluded volume, neglected here, are also important.
Whereas filament extension is a relatively common experiment, filament buckling experiments are more difficult to conduct and, therefore, are more rare [1,2,[6][7][8]. There, the critical buckling force is usually analyzed using classical elasticity theory [9], which ignores fluctuations or is valid only in the zero-temperature limit. However, thermal fluctuations should be taken into account, since the forces measured in single-molecule experiments lie in the pico-Newton range, which are comparable to thermal forces. It has been suggested [10] that the sharp (second-order) transition to the buckled state predicted by Euler does not occur when fluctuations are present, resulting in a smooth (higher-order) transition to the buckled state.
Polymers are usually modeled as random walks: one assumes that the polymer is composed of a chain of statistically independent segments, whether the chain is undiluted or diluted [11][12][13][14]. This model is, in a sense, the opposite extreme of the elastic rod. The conformations of the random-walk model are determined completely by entropy, whereas the conformations of the rigid rod are determined by energy. However, a semiflexible polymer chain is an important physical system where both entropy and energy are significant for chain conformations and dynamics. In such systems, one needs a statistical mechanical model that incorporates both elasticity and fluctuations.
There are four possible approaches to study buckling: approximations made with simplifying assumptions [15,16], asymptotic expansions [17][18][19], numerical simulations, or exact calculations. Without exact results, it can be difficult to assess the reliability of approaches relying on simplifying assumptions. Not surprisingly, then, contradictory results are obtained [10,16,20]. Asymptotic results are reliable but give useful results over a limited domain only, since one must guess the appropriate base case for expansion.
We consider here exact predictions of buckling and of conformations of fluctuating semiflexible polymers-accompanied by Monte Carlo simulations, which provide a helpful visualization. The simulations also provide a strong check of our analytic results, since the two approaches are very different. Surprisingly, we find that fluctuations can be ignored only for very short chains-when the chain length is less than approximately one percent of its persistence length, l p (defined below). We also show that fluctuations stabilize semiflexible polymer to buckling, in contradiction to approximate approaches. Our results can be explained by analyzing thermodynamic contributions to critical buckling force.

Results
Our starting point is the worm-like chain model [21,22], where the chain is represented as an inextensible space-curve r(s) with arclength variable s that runs from 0 to L. The deformation of the chain is opposed by a quadratic bending energy where κ is the elastic modulus, L is the contour length of the chain, u(s) = ∂r ∂s is the tangent vector of r(s), and inextensibility is enforced by the constraint |u(s)| = 1 at all s. In this work, we consider the response of a chain to a compressive force f , which is captured by a force energy E f = f |r(L) − r(0)| = f R that favors having the ends in close spatial proximity.
We define the persistence length l p = κ k B T as the length of chain over which a curved conformation (with local curvature 1/l p ) has an energy comparable to k B T. Predictions using the worm-like chain model represent the thermal behavior of a bendable but unstretchable elastic thread, and the behavior reduces to that of a classical elastic chain in the limit of zero temperature (or large persistence length l p ). At finite temperature, the persistence length l p defines the length scale of correlation between tangent vectors along the chain.
The bending energy E B provides the necessary input to predict the minimum-energy shape of a chain. Appendix A provides an overview of the buckling of an elastic filament [9], resulting in a theoretical prediction of the critical buckling force f E (or the Euler buckling force) and the shape of the minimum-energy conformation. The critical Euler force f E depends on the rod length L and persistence length l p as f E = k B Tl p π 2 /L 2 . For forces below f E , the minimum-energy conformation is a straight chain with R = L, and at larger forces, the chain adopts a bent (or buckled) conformation.
A thermally fluctuating elastic filament (i.e., a worm-like chain) is governed by the Green function G(R; L), which gives the probability that the end-to-end distance is given by R and is formally found by path integration over all conformations u(s) to be In this work, we exploit a method that calculates the exact Green function G(K; p) for the worm-like chain model, which is the Fourier-Laplace transform of G(R; L) (Fourier transform R → K and Laplace transform L → p). The Green function G(K; p) is written in the continued fraction representation [23,24] G(K; p) = 1 where P l = p + l(l + 1), and a l = l/ (2l + 1)(2l − 1). To use Equation (3), we must inverse Laplace transform from p to the dimensionless chain length N = L/(2l p ) and inverse Fourier transform from K to dimensionless end separation r = R/(2l p ). Our previous work [23,24] provides efficient methods to perform these inversions, and we provide numerical realizations of the Green function G(R; L) on our research group websites [25,26]. The Helmholtz free energy of a worm-like chain subjected to a compressive force f is given by Figure 1 shows the free energy F versus the end extension R for a chain length of L/(2l p ) = 1/4 at four values of the force f , ranging from f / f E = 0 to f / f E = 2. This range of forces corresponds to two forces below the zero-temperature critical force f E and two forces above. The solid curves are the total free energy F, and the dashed curves indicate the force energy E f = f R at each force. On each curve, the circle indicates the location of the free-energy minimum. Figure 1 provides a clear indication of the non-linear response of a worm-like chain to a compressive force. Between f / f E = 0 and f / f E = 0.67, there is very little change in the minimum free-energy end extension R min , indicating a large resistance to compressive force. However, between f / f E = 0.67 and f / f E = 1.33, there is considerable end retraction associated with a dramatic compliance in this force range. We emphasize that all of the conformation images exhibit dramatic conformational fluctuations, despite the fact that the chain length L/(2l p ) = 0.25 is significantly shorter than a persistence length.   Figure 2 shows R min versus f for chain lengths ranging from L/(2l p ) = 0.25 (blue curve) to L/(2l p ) = 2 (red curve). The compressive force ranges from where negative values correspond to a worm-like chain under tension. The dashed curve in Figure 2 shows the zero-temperature (i.e., no fluctuations) end extension [9], and the dotted vertical line gives the zero-temperature critical force f = f E . The inset images show realizations from Monte Carlo simulations at the same four force values as in Figure 1. We note that the Monte Carlo images in Figure 1 have a fixed end extension given by R min , whereas the Monte Carlo images in Figure 2 are from simulations with a fixed force f , and the ends undergo thermal fluctuations as dictated by the shape of the free energy curve. The large fluctuations of chain conformations seen at points C and D in Figure 2 are not surprising, given how flat the free energy curves are for C and D in Figure 1.
The zero-temperature end extension (dashed line in Figure 2) indicates a second-order transition at the critical force f E . In contrast, the chain lengths with finite temperature exhibit a smooth transition for forces near the zero-temperature critical force f E . We interpret buckling of a thermally fluctuating chain as having a nonlinear compliance or a "softening" of the response of the end extension to compressive force. We define a response function α F ≡ − dR min d f , and we define the critical force f c as the force with the largest response in R min versus force f (i.e., a maximum value in the response function α F ). In other words, the critical force f c corresponds to an inflection point on the R min versus f curve.
The left side of Figure 2 provides constant-L curves of R min versus f , so this information can be represented by a surface plot of R min versus f and L/(2l p ). The right side of Figure 2 shows a bird's-eye view of this surface plot with red indicating R min = 1 and blue corresponding to R min = 0. The colored vertical lines indicate the L/(2l p ) values plotted on the left side of Figure 2. The black curve shows the critical force f c for each L/(2l p ) value. One sees that fluctuations actually enhance the stability of the chain to buckling by as much as 80%. As the length of the chain goes to zero, the force naturally converges to the zero-temperature or the Euler-buckling limit. The counterintuitive behavior of the critical force f c versus the length L/(2l p ) can be explained by parsing out the thermodynamic contributions to the free energy. To determine each thermodynamic contribution, we exploit our Monte Carlo simulations with finite degrees of freedom (Appendix B). The average polymer energy E poly − E 0 converges in the limit of large bead number, and the entropy is found using S = ( E poly − E 0 − F)/T. The reference energy E 0 defines a thermodynamic energy that arises from the discretization of the polymer chain. The value of the reference energy is where N b is the number of beads in the discretized chain. This value offsets the thermodynamic contribution from the energy associated with the equipartition theorem for each active mode in the discrete chain (see Appendix B for more details). Figure 3 gives results from a Monte Carlo simulation for N b = 51 beads, and Appendix B demonstrates that these results are invariant to the discretization. . Thermodynamic contributions to the free energy F for a worm-like chain versus the end-to-end distance R. The free energy F (black dots) is determined using our exact analytical treatment (Equation (4)). The average polymer energy E poly − E 0 (blue dots) is found from discretized Monte Carlo simulations, and the entropy S (red dots) is calculated from S = ( E poly − E 0 − F)/T. The dashed curve is the zero-temperature bending energy, and the four Monte Carlo conformations span the range of end retractions. Figure 3 shows the free energy F (black dots), average bending energy E poly − E 0 (blue dots), and entropy TS (red dots) versus end retraction for a worm-like chain with length L/(2l p ) = 1/4 and zero compressive force f . The four Monte Carlo snapshots indicate typical conformations over the range of end retraction. The dashed curve gives the polymer energy E poly in the zero-temperature limit. For this relatively short chain (L = l p /2), the free energy F closely tracks with the zero-temperature bending energy, indicating that the response is largely dictated by the bending energy. However, there is a systematic trend of faster increase of the free energy F with end retraction than the zero-temperature behavior. This observation is consistent with the elevation of the critical buckling force with length L or, equivalently, with the degree of thermal fluctuations.
This subtle trend is dictated by the entropic contribution. The free-energy minimum coincides with the maximum-entropy state, and the entropy S exhibits a mild decrease with the decreasing R. Thus, the extended state is both energetically favored and entropically favored, and the critical force is, therefore, larger than the zero-temperature Euler critical force f E . We note that the entropic favorability arises from the conformational fluctuations that are accessible at each end retraction. This should not be confused with the entropic preference for a flexible Gaussian chain to have zero end-to-end distance, since a flexible chain does not exclude conformations based on bending deformation.

Discussion
In this work, we determine that fluctuations can be ignored only for chains much shorter than their persistence length. In fact, we find that the buckling phenomenon is never purely elastic for macromolecules, except on length scales where the granularity of monomers becomes important, and a continuous rod is a poor assumption. Our exact results show unambiguously that fluctuations indeed stabilize semiflexible chains to buckling in 3D. The counterintuitive behavior of the critical force occurs because the extended state is thermodynamically-both energetically and entropically-favored at the critical region, and the critical force needed to overcome both thermodynamic contributions is elevated above the zero-temperature critical force f E . Our exact calculations provide a strong check for numerical simulations and testable hypotheses for experimental phenomena.
Many critical biological processes involve deformation of semiflexible biopolymers such as DNA and filamentous proteins. The formation of gene regulatory complexes involve looping segments of DNA to bind to protein complexes. The buckling process involved in looping distal segments of DNA is controlled by the elastic mechanics of the DNA strand and the entropic change upon docking the DNA strand into the regulatory complex [27,28]. Single-molecular measurements of the kinetics of DNA binding to Lac repressor for DNA lengths of approximately 100 basepairs (i.e., L/(2l p ) ≈ 0.67) demonstrate that the time scale associated with loop formation is highly sensitive to DNA mechanics (both bending and twisting), leading to time scales around 100 seconds for DNA binding [27]. Our theoretical predictions based on the worm-like chain model for the kinetics of loop formation and loop dissociation exhibit quantitative agreement with the measured in vitro binding and unbinding time scales [27,28]. The in vivo forces involved in loop formation in this regulatory complex arise from protein physical interactions (e.g., binding proteins including Lac repressor and Catabolite Activator protein), motor protein activity, and thermal fluctuations. Our current work provides a physical basis for interpreting the impact of conformational fluctuations in modulating the forces necessary to induce a biopolymer to adopt such a buckled conformation.
Single-molecule measurements have been instrumental in elucidating the response of individual biopolymers to external forces and torques [29][30][31][32][33]. A closely related problem to the Euler buckling phenomenon considered in this manuscript is the response of an elastic filament to externally applied torque [9,34]. For an elastic filament under tension and subjected to torque, elastic strain associated with twist deformation leads in a buckling transition at a critical torque, resulting in the formation of a looped "hackle" in the filament. The hackle causes an abrupt decrease in the end-to-end distance of the filament. This phenomenon is observed in DNA strands subjected to external torque using an angular optical optical trap [32]. However, the experimentally observed end-to-end retraction [32] exceeds predictions from classical elasticity theory [9,34]. This effect is modulated by salt concentration [33], indicating that self repulsion plays a non-trivial role in the twist induced buckling transition. In contrast to Euler buckling, which exhibits a distributed buckle over the length of the polymer, twist-induced buckling localizes the instability within a small loop that would have considerable self-interaction. However, the impact of conformational fluctuations is an important contribution in the twist buckling phenomenon, and our work on the statistical behavior of the worm-like chain model with twist [35] would serve as a valuable starting point for performing analogous predictions as presented in this work.
Current single-molecule experimental methods primarily apply tension to measure mechanical properties of biopolymers and mechanochemical behavior in enzymes [30,31]. In order to experimentally test the predictions in this manuscript, it would be necessary to develop a single-molecule apparatus or a DNA construct that either directly applies a compressive force or is capable of transducing a tensile force to a compressive force. One approach that could be developed is to construct a DNA or RNA nanostructure that can be pulled in one direction, resulting in a compression in another direction. For example, a cube structure can be pulled at opposite vertices, resulting in elongation along the inter-vertex axis and compression along the orthogonal axes. More complicated architectures resulting in tensegrity structures [36] could be exploited in determining biopolymer response using a tensile apparatus. Furthermore, such structures are proposed as a structural motif in cellular mechanics [36] with complex elastic properties, and the study of such structures would be valuable for their intrinsic nonlinear elastic properties. The impact of thermal fluctuations on the response of biopolymer structures and networks is crucial in a variety of biological phenomena, and our work serves to address the buckling phenomenon within the individual segments of a more complex molecular architecture.
It follows that the role of thermal fluctuations for buckling phenomena cannot be neglected during the analysis of collective compression-resistant behavior of, for example, actin networks or nanotubes [37]. Understanding of collective behavior or network response to deformations starts with the simpler problem of how a single semiflexible polymer responds to compressive force. We hope that our findings will encourage experimentalists to seek reliable data on critical buckling force for a single semiflexible filament. It would be interesting to study whether fluctuations generally stabilize transitions in other microscopic systems. The work presented here presents a roadmap for attacking such problems.

Appendix A. Euler Buckling in the Zero-Temperature Limit
At zero temperature, the polymer shape is given by the minimum-energy conformation, which is found by setting the first variation of the energy to zero. The energy of an elastic filament with orientation given by tangent vector u(s), or, equivalently, by polar coordinates in a given coordinate system θ(s), φ(s) is written as The first variation (or differential change) of the end-to-end distance R is written as Now, we align the end-to-end vector to the z-axis, such that R/R =ẑ. The first variation in the energy is now written as The solution to the 2nd Euler-Lagrange equation is φ = φ 0 , which simplifies the 1st Euler-Lagrange equation The Euler-Lagrange equation is solved using a symmetric formulation. Define σ = 2s/L − 1, resulting in the dimensionless Euler-Lagrange equation where am(x|m) is the Jacobi amplitude, which is the inverse function of the elliptic integral The final orientation θ f is the angle at σ = 1, which is found from the transcendental equation There is always a solution θ = 0. However, θ = 0 arises for forces above a critical value of f c = π 2 . Numerical solution of Equation (A12) is used to find θ f , and the filament shape is determine by subsequent evaluation of Equation (A10).

Appendix B. Monte Carlo Simulations
We perform Monte Carlo simulations of a worm-like chain under compressive force. The polymer chain is modeled as a thermally fluctuating isotropic elastic filament whose conformations are dictated by a polymer chain energy E poly . The inextensible worm-like chain model is defined to be a continuous chain, and implementation of this model into a simulation requires discretization of the chain to N b beads. In our previous work [38,39], we demonstrate that the process of discretizing the chain results in quadratic-order deformation modes that go beyond the bare bending energy of the continuous chain. A simple example of this effect is that the continuous inextensible worm-like chain model has an infinite chain stretching modulus; however, a discretized worm-like chain would manifest a stretching mode due to the entropic elasticity from the implicit chain between the beads. At the quadratic level, the discrete polymer chain energy contains a bending, stretching, shearing, and bend-shear coupling that accounts for the underlying thermal fluctuations of the continuous chain [38,39].
The resulting discrete, coarse-grained model for an inextensible worm-like chain utilizes our discrete shearable stretchable worm-like chain (dssWLC) [38,39] with N b discrete beads. The dssWLC defines a universal coarse-graining procedure that is able to accurately capture chain statistics at all length-scales above the discretization ∆ [38,39]. In contrast to the bare discrete worm-like chain, the polymer configuration is additionally represented by an auxiliary set of orientation vectors u i , which capture the coarse-grained tangent of the chain configuration in the vicinity of r i . At a given length-scale of coarse-graining ∆, the conformational free energy of the chain is dictated by six effective elastic parameters that capture the conformational statistics at the level of coarse-graining. The bending modulus B elastically couples adjacent local chain orientation u i , u i−1 . The stretch modulus characterizes the entropic cost of stretching the chain away from its ground-state extension γ. The shear modulus ⊥ penalizes shearing of the local chain orientation away from the inter-bead displacement vector R i = r i − r i−1 . The bend-shear coupling modulus η characterizes the tendency of the chain to shear in the same direction as bending. The chain conformational free energy is given by where R ⊥ i = R i − R i · u i−1 is the component of the inter-bead displacement vector perpendicular to u i−1 . Our research group website [25] contains Matlab (R2017a, MathWorks, Natick, MA, USA) scripts to determine the elastic parameters for the dssWLC model and for performing Monte Carlo and Brownian dynamics simulations of this model.
The average polymer chain energy E poly cannot be determined directly from the continuous worm-like chain free energy F, given by Equation (4). The formulation of the worm-like chain represents the polymer as a continuous elastic thread, having infinite degrees of freedom. Equipartition theorem dictates that each thermally active degree of freedom would contribute k B T/2 to the average energy and k B /2 to the entropy. Thus, the worm-like chain model predicts infinite average energy and infinite average entropy, which results in several unphysical effects such as diverging polymer stress in the continuous-chain limit [40]. Although each contribution independently diverges, the free energy F is finite and well defined, and our exact results provide this thermodynamic contribution.
Our Monte Carlo simulations provide a framework to evaluate the change in the average polymer chain energy E poly with end extension, provided we can reliably translate the discrete chain behavior to the continuous limit. We perform several simulations with an increasing number of beads to determine the average polymer chain energy E poly . The left plot of Figure A1 shows results from Monte Carlo simulations of a polymer chain of length L = l p /2 with bead number N b = 11 (blue), N b = 31 (purple), and N b = 51 (red). These results demonstrate the increase in the average polymer chain energy with the number of beads, reflecting the divergence in the energy as the number of degrees of freedom increases. However, the thermodynamic contribution from these active modes is easily interpreted and eliminated from these results. Since each bead has a position and orientation, there are five degrees of freedom per bead. Thus, a free chain has an average polymer chain energy predicted from the equipartition theorem to be 5 2 k B T(N b − 1). We define a reference energy E 0 = 5 2 k B T(N b − 4) that includes the free-chain energy in order to eliminate the contribution from discretization; the factor of four is a convenient choice for visualization of the results in Figure 3. The right plot of Figure A1 shows our simulation results for average polymer chain energy E poly − E 0 (i.e., subtracting the reference energy E 0 ) for bead number N b = 11 (blue), N b = 31 (purple), and N b = 51 (red). The coarsest discretization (N b = 11, blue) exhibits slight numerical deviations from the results from simulations with N b = 31 (purple) and N b = 51 (red). With increasing N b , the average polymer energy E poly − E 0 is invariant to the number of beads N b , and we can interpret the average polymer energy E poly − E 0 as the thermalized deformation energy of a buckled polymer chain. The entropy S is then found from our exact results for the free energy F and the average polymer energy E poly − E 0 using the thermodynamic expression S = ( E poly − E 0 − F)/T. Our results for these thermodynamic contributions based on Monte Carlo simulations with N b = 51 are shown in Figure 3.  Figure A1. The left plot shows the average polymer energy E poly versus the end extension from Monte Carlo simulations for a chain of length L = l p /2 and number of beads N b = 11 (blue), N b = 31 (purple), and N b = 51 (red). The right plot demonstrates that subtracting the references' energy E 0 = 5 2 k B T(N b − 4) from each average polymer energy E poly gives an average polymer energy that is invariant to the discretization and reflects the deformation energy of the fluctuating polymer chain under compression.