Comparison of Conformational Phase Behavior for Flexible and Semiflexible Polymers

We employ the recently introduced generalized microcanonical inflection point method for the statistical analysis of phase transitions in flexible and semiflexible polymers and study the impact of the bending stiffness upon the character and order of transitions between random-coil, globules, and pseudocrystalline conformations. The high-accuracy estimates of the microcanonical entropy and its derivatives required for this study were obtained by extensive replica-exchange Monte Carlo simulations. We observe that the transition behavior into the compact phases changes qualitatively with increasing bending stiffness. Whereas the Θ collapse transition is less affected, the first-order liquid-solid transition characteristic for flexible polymers ceases to exist once bending effects dominate over attractive monomer-monomer interactions.


Introduction
Statistical analysis methods have long been used for the identification, characterization, and classification of phase transitions in complex systems. Advanced computer simulation techniques, mostly based on Markov chain Monte Carlo algorithms, helped advance the field beyond the deadlocks mathematical approaches had run into. Nonetheless, the general idea remained the same. Phase transitions were generally only considered in the thermodynamic limit, which was historically introduced to make studies of macroscopic systems mathematically tractable. In computational statistical physics, finite-size scaling analysis provided the tool for extrapolating the results obtained in simulations of systems of finite size toward this hypothetical limit. These methods proved to be extremely successful, and for decades the scientific community did not see any reason to change paradigms and consider a more general perspective despite the fact that all systems in nature are finite and systems on nanoscales moved into the focus of substantial interest in microbiology and technology. Ignoring surface effects was also convenient.
Consequently, conventional canonical statistical analysis of phase transitions rests on the search for catastrophic behavior of thermodynamic state variables and response functions in thermodynamic parameter spaces. Hence, significant effort has been dedicated to the evaluation of critical exponents in continuous phase transitions and the grouping of systems in universality classes.
In recent years, however, the growing interest in systems on microscopic and mesoscopic scales, which do not satisfy the criteria for the thermodynamic limit, has shifted the perspective. Alternative approaches that enable studies of phase transitions in systems of any size have turned out to be promising. Modern interdisciplinary research of systems on nanometer scales cannot ignore the fact that the surface of the system is at least as important for the structural behavior of the system as the bulk effects. For example, the functionality of biomacromolecules like proteins depends on the folding process into stable geometric conformations. These processes, which occur in a thermal environment, resemble transitions between disordered and ordered phases known from macroscopic systems. In fact, for very long polymers, these structure formation processes are phase transitions, even in the strict conventional sense. However, many types of heterogeneous polymers like proteins cannot be scaled up, but still show clear transition features. This makes the extension of the established theory of phase transitions a necessity.
However, canonical statistical methods create ambiguities, which render the unique characterization of phase transitions problematic. A prominent example is the peak in the specific heat curve of the one-dimensional Ising model. As it is also finite in the thermodynamic limit, one does not associate a phase transition with it. But how to interpret a peak in response quantities, if the thermodynamic limit for this system does not exist? In recent years, it has nonetheless become common to consider peaks or "shoulders" transition signals in curves of response quantities for finite systems. This appears to be an inconsistent approach and renders, in fact, the choice of the temperature as the basic thermodynamic state variable questionable.
Combining results of previous studies of the microcanonical Boltzmann entropy [1] and the principle of minimal sensitivity [2,3], we recently developed a systematic approach to phase transitions in systems of any size by identifying least-sensitive inflection points in the entropy and its derivatives [4]. This analysis method even allows for the classification of phase transitions in analogy to Ehrenfest's classification scheme in thermodynamics, which is based on derivatives of thermodynamic potentials [5].
In this paper, we extend previous Monte Carlo computer simulation studies of a coarse-grained model for flexible polymers [6] and analyze structural transitions in semiflexible polymers by means of the generalized microcanonical least-sensitive inflection point method. We thoroughly investigate the changes of phase behavior of a coarse-grained model for flexible and semiflexible polymers.
The paper is organized as follows. In Section 2, we describe the flexible and semiflexible variants of the coarse-grained polymer model used in this study, as well as the numerical methods for the stochastic simulation of the model and for the statistical analysis of the data obtained in the simulations. The results of the canonical and microcanonical analysis are discussed and the phase behavior of the polymers is interpreted in Section 3. This section also contains a detailed comparison of geometric properties of the lowest-energy conformations found for each system. Eventually, the paper is concluded by the summary in Section 4.

Models and Methods
For this study we employed a generic coarse-grained model for polymers in implicit solvent that allows us to investigate and compare the phase behavior of flexible and semiflexible polymer chains by means of advanced Monte Carlo computer simulation methodologies. Model, simulation methods, and the microcanonical inflection-point technique we used for the statistical analysis of the conformational transitions are described in the following.

Generic Model for Flexible and Semiflexible Polymers
The total energy of a polymer chain with N monomers is composed of contributions from bonded and non-bonded interactions between monomers. For the interactions between non-bonded monomers, we employ the standard 12-6 Lennard-Jones (LJ) potential: Here, r is the monomer-monomer distance, σ = 2 −1/6 r 0 is the van der Waals distance associated with the potential minimum at r 0 , and LJ is the energy scale. For computational efficiency, we introduce a cutoff at r c = 2.5σ. Shifting the potential by the constant V shift = V LJ (r c ) avoids a discontinuity at r c . Thus the potential energy of non-bonded monomers is given by The elastic bond between two neighboring monomers is modeled by a potential which combines a Lennard-Jones potential and the finitely extensible nonlinear elastic (FENE) potential [7][8][9]: where the FENE parameters are fixed to standard values R = (3/7)r 0 and K = (98/5)r 2 0 [10]. The bond length r is restricted to fluctuations within the range [r 0 − R, r 0 + R]. With these parameters the minimum of V B is located at r 0 .
For semiflexible polymers, an additional potential related to the bending stiffness is introduced. Any deviation from the reference angle θ 0 between neighboring bonds is subject to an energy penalty of the form: If we represent a conformation of a polymer chain with N monomers by X = (r 1 , ..., r N ), where r i is the position vector of the ith monomer, the total energy of the polymer is given by where r i,j = |r i − r j | is the distance between monomers i and j, and θ l is the bond angle between two adjacent bonds. The parameter κ ≥ 0 controls the stiffness of the polymer chain. For κ = 0, the model describes flexible polymers. In simulations and statistical analysis of the results, we set the basic scales to the following values: k B = 1 (Boltzmann constant), LJ = 1, and r 0 = 1. For the reference bending angle, we chose θ 0 = 0. The flexible chain with N = 55 monomers has already been studied extensively in the past and serves as the reference for the comparison with the semiflexible model. This chain length is sufficiently short to recognize finite-size effects, but long enough for the polymer to form a stable solid phase at low temperatures [11].

Replica-Exchange Monte Carlo Method and Multiple-Histogram Reweighting
For the simulations of the polymer, we employed the replica-exchange (parallel tempering) Monte Carlo method [12][13][14][15][16]. Replicas of the systems were simulated at different temperatures T k ∈ [0.1, 5.0] with k = 1, 2, ..., K. The total number of temperature threads K ranged from 40 to 50 in individual simulations. Up to ten independent runs were performed for all κ values studied, which allowed for the estimation of numerical errors in all microcanonical quantities.
At each temperature T k , Metropolis sampling was performed, which is based on the acceptance probability: Here k B T k is the the thermal energy at the simulation temperature T k , E k is the energy of the proposed state, and E k is the energy of the current conformation.
Local structural changes were achieved by displacement moves. In this Monte Carlo update, a monomer i is randomly picked and its position is modified by a random shift ∆r i within a cubic box with edge lengths r d surrounding the monomer. Before measurements can be performed, r d is determined adaptively for each temperature thread to achieve a Metropolis acceptance rate of about 50%.
In order to explore the conformation space more efficiently, rotational pivot updates were also used. After 20 to 70 sweeps of displacement updates (depending on the temperature), we performed a sweep of pivot rotational updates. In a pivot update, a monomer i is randomly chosen as the pivot point and a rotation axis vector s i . Then the section of the chain following the pivot monomer was rotated about s i by a randomly chosen angle taken from the interval [0, 2π].
Replicas were swapped between neighboring threads k and k + 1 after every 1500 sweeps with the exchange probability: where E k and E k+1 are the energies of the replicas before the swap at temperatures T k and T k+1 , respectively. In order to estimate the density of states of the system, we use the multi-histogram reweighting method [17,18]. For this purpose, we measured the canonical histograms in each thread. By utilizing the canonical distribution function P can (E; T k ) ∼ g(E) exp{−E/k B T k }, estimates for the density of states can be obtained by reweighting the energy histograms h(E; T k ) measured in the different simulation threads k, g(E) = h(E; T k ) exp{E/k B T k }. Using a single histogram only covers a narrow range of energies effectively. Combining the histograms obtained at different temperatures by employing the multi-histogram reweighting method yields an estimator for the density of states that covers the entire energy range [17,18]: Here M k is the total number of sweeps in the thread k and Z k is an estimator for the partition function: Equations (8) and (9) are solved iteratively until reasonable convergence has been achieved.

Generalized Microcanonical Inflection-Point Analysis Method
Historically, phase transitions have been identified and classified by means of discontinuities or divergences in thermodynamic state variables or response functions. However, systems of finite size do not exhibit such obvious signals, as these only occur in the thermodynamic limit. In this study, we use the recently introduced generalized microcanonical inflection-point analysis method [4] for the systematic identification and classification of transitions in systems of any size. This method, which combines microcanonical thermodynamics [1] and the principle of minimal sensitivity [2,3], has already led to novel insights into the nature of phase transitions. Even the Ising model, which has been excessively studied in almost a century, possesses a more complex phase structure than previously known, as our most recent analysis showed [19].
Adopting Boltzmann's formula, the microcanonical entropy can be written as where g(E) is the density of states with energy E. When not experiencing phase transitions, the curves of entropy S(E) and its derivatives exhibit well-defined concave or convex monotony [4]. From canonical statistical analysis of first-and second-order transitions, it is known that entropy and/or internal energy rapidly change, if the temperature is varied near the transition point. This behavior corresponds to least-sensitive dependencies of microcanonical quantities in the space of system energies. Thus, a phase transition causes a least-sensitive inflection point in the entropy or its higher derivatives and impacts their monotonic behavior. By systematically analyzing these alterations, different types of transitions can be identified and classified. In this scheme, a first-order transition is signaled by a least-sensitive inflection point with energy E tr in S(E). Therefore, the first derivative, which is the inverse microcanonical temperature β(E), possesses a positive valued minimum at E tr , Consequently, if there is a least-sensitive inflection point in β(E), the phase transition is classified as a second-order transition. The derivative of β(E) has a negative-valued peak at the transition energy E tr , More specifically, we call this an independent second-order transition. This implies that there is another transition type: dependent transitions. As the name suggests, dependent transitions are associated with an independent transition of lower-order rank. It is important to note that not all independent transitions have a dependent companion. However, if it exists-and it can only exist at a higher transition energy than its independent partner-it servers as a precursor of the independent transition in the disordered phase. For example, if this inflection point is located within the convex region of the entropy, which is associated with the first-order transition, the transition is a dependent second-order transition. Consequently, there is a positive-valued minimum in the derivative of β(E) at the transition energy E tr , Generally, for an independent transition of odd order (2k − 1) (k is a positive integer), the (2k − 1)th derivative of S(E) possesses a positive-valued minimum, and an independent transition of even order 2k is characterized by a negative-valued peak in 2kth derivative, A dependent transition of even order 2k is identified by a positive-valued minimum in the 2kth derivative, while for odd order 2k + 1: Dependent transitions are less common than independent transitions. In this polymer study we did not identify any dependent transitions.

Canonical Analysis
In the past, the common approach to the characterization of phase transitions has been canonical statistical analysis. For infinitely large systems, one or more derivatives of an appropriate thermodynamic potential such as the free enthalpy or the free energy would exhibit nonanalytic behavior, if one of the potential's natural variables was altered at the transition point. Because it is easy to control in experiment, the analysis has typically been done in the space of the canonical (or heat bath) temperature T can . Discontinuities in the canonical entropy and specific heat help identify and characterize discontinuous and continuous phase transitions, respectively, in the thermodynamic limit. Later, this approach was simply extended to systems of finite size that do not allow for the hypothetical extrapolation toward the thermodynamic limit (for example, heterogeneous polymers such as proteins).
However, transition features in finite systems cannot be linked to nonanalyticities. Hence, peaks and shoulders in response quantities like heat capacity or order parameter fluctuations have usually been considered signals of pseudotransitions in finite systems. The problem with this approach is that extremal points in fluctuating quantities are not safe indicators of phase transitions. The probably most prominent example is the specific-heat curve for the one-dimensional Ising model, which exhibits a pronounced maximum at a certain temperature, but there is no obvious feature that links this extremal thermal activity to a phase transition. In the thermodynamic limit, this peak does not develop into a divergence.
Interestingly, if one plots the heat capacity for the 1D Ising system as a function of the inverse temperature β = 1/T, though, the peak vanishes and there is no pseudotransition signal. Note that in the microcanonical analysis the inverse temperature β(E) = dS(E)/dE is the more obvious thermodynamic variable than the temperature. Consequently, microcanonical inflection-point analysis does not show any transition signal for the 1D Ising system in equilibrium [19].
In fact, historically the (canonical) temperature was introduced when the need arose to assign a reliable scale to the thermoscope to quantify the environmental heat content, which ultimately lead Galilei to the introduction of the thermometer at the end of the 16th century. It is ironic that the nowadays most commonly used centigrade scale, the Celsius scale, was introduced in 1742 by Celsius in inverted form, i.e., Celsius chose 0 • as the boiling point of water under norm conditions and 100 • as the reference point for the melting of ice. At about the same time, French scientist Christin independently used a centigrade scale that essentially inverted the original scale by Celsius, which is what we nowadays refer to as the Celsius Scale [23].
This brief historic recollection shows that the temperature and its scale were introduced in an era, when it was still believed to be a quantitative measure for heat, and the scale was ultimately linked to the instrument and substances used for its measurement. Thus, there is no physical reason that favors the historic quantity temperature over the inverse temperature, which is the more natural variable in the microcanonical analysis. To distinguish the microcanonical temperature from the heat bath temperature, we denote the latter by T can in the following canonical analysis of structural transitions.
Consequently, we plot the canonical mean energy E as a function of β can = 1/T can in Figure 1 for the three models studied: κ = 0 (flexible chain) and κ = 1, 2 (semiflexible variants). All curves show a sharp drop at β can ≈ 0.7. Around β can ≈ 3.2, the flexible chain (κ = 0) experiences another significant drop in energy, which is less pronounced for the semiflexible polymer with bending stiffness κ = 1. No obvious signal is visible on this level for the stiffer chain with κ = 2.
The plots of fluctuating quantities, such as the heat capacity C V = d E /dT can and the fluctuations of the mean square radius of gyration, d R 2 gyr /dT can , allow for a more detailed analysis. Both quantities are shown as functions of β can in Figure 2a,b, respectively. Whereas we certainly see changes in the curvature of the C V curves at β can ≈ 0.7, only the plots of the structural fluctuations show sharp peaks. These indicate the well-known Θ collapse transition between extended, random-coil structures and compact globular conformations. The transition signal is more pronounced for the stiffer chains. Since this transition is more entropy than energy driven, it is not surprising that it shows up more prominently in the structural rather than the energetic fluctuations. This is different for the freezing transition at about β can ≈ 3.2, which is strongest for the flexible chain (κ = 0). Large error bars at very small temperatures (β can ≈ 8.5) prevent a conclusion about another separate transition for the semiflexible polymer with κ = 2, which, if it exists, marks the transition into the global energy minimum basin.

Microcanonical Analysis of Phase Transitions
In this section, we perform a full-scale microcanonical inflection-point analysis of the three models that aims at the identification of all structural transitions in these systems up to third order. As in the canonical analysis, we first discuss the low-β (or high-temperature) regime in the relevant energy space. The entropies plotted in Figure 3a for all cases do not possess least-sensitive inflection points and thus there is no first-order transition, as expected. However, least-sensitive inflection points in the first derivative (β) signal second-order phase transitions in all three systems, which reflect the mostly entropic Θ collapse from random-coil to globular polymer structures. The corresponding peak locations in the second entropy derivative (γ), shown in Figure 3c, allow for a unique determination of the transition points in energy space and thus also in β, thereby rendering the ambiguous canonical analysis of response functions in the previous section obsolete. Since there are no least-sensitive inflection points in the γ plots, none of the systems undergoes a third-order transition in this energy region.  Figure 4 shows the same microcanonical quantities as plotted in Figure 3, but for an intermediate energy range that covers the inverse temperatures in the interval β = (2, 4). As expected for flexible polymers, the entropy curve for κ = 0 does exhibit a least-sensitive inflection point, which corresponds to the minimum in the backbending region found in the β plot (Figure 4b) at about ∆E (κ) = 36.5. The inverse temperature associated with it is approximately β = 2.95, which confirms earlier results for flexible polymers [6]. The polymer undergoes a freezing transition from the liquid-globular states into the solid phase, which for this model is dominated by icosahedral structures. Remarkably, for the semiflexible polymers with κ = 1, 2, the situation changes. For κ = 1, there is still a weak inflection point in the entropy, but the β curve has already almost plateaued at about the same energy value, where we found the first-order transition for the flexible polymer. Consequently, the peak value of the next derivative (γ) is virtually zero at this energy (Figure 4c). Therefore, the freezing transition seems to turn from first order for the flexible polymer (κ = 0) to second order for the semiflexible case with κ = 1. However, even more surprisingly, the transition behavior changes again for the stiffer chain with κ = 2: In fact, the transition signal has completely vanished. The subsequent structural analysis of ground-state properties will lend deeper insight into the reasons for these changes.
The second derivative γ in Figure 4c shows signs of independent third-order transitions for κ = 0 and κ = 1 in the ordered phase, but not for κ = 2. These results suggest that freezing into a unique and characteristic global energy minimum state is not possible if bending effects overwhelm non-bonded interactions that are commonly responsible for tertiary structure formation in polymeric systems. Local bending effects turn into restraints that inhibit or at least suppress the formation of symmetries. Previous studies have shown that bending restraints play an important role in the formation of stable conformations in tertiary assemblies of helical segments [24].
Eventually, the sequence of figures shown in Figure 5 covers the lowest-energy regime for all models compared in this study. In the flexible case, no additional transitions are expected and, correspondingly, we do not see signals in any of the derivatives up to third order that would indicate a transition. The situation is potentially different for the systems with bending restraint. Whereas we do not find indications of transition signals in the entropy and inverse temperature curves for κ = 1 and κ = 2 in Figure 5a,b, respectively, the γ results shown in Figure 5c might tell a different story. The error bars are too large for an ultimate conclusion, but it looks like a third-order transition develops close to the ground state for the semiflexible polymer models.
The quantitative results of the microcanonical inflection-point analysis obtained for the polymer systems studied here are listed in Table 1.

Structural Analysis of Lowest-Energy Conformations
The lowest energies found in the simulations, which we consider best estimates of the ground state energies, are E Since the bonded interactions (bond vibrations) are negligible near the ground state, the differences in the ground-state energy estimates for the different models must be attributed to the penalty paid for chain bending. The ground-state conformation for the flexible polymer that does not experience these restraints is formed by optimizing the non-bonded interaction and results, for this model, in an ideally icosahedral structure. Turning on the bending restraint and setting κ = 1, the competition between nonlocal attraction and local repulsion increases the tension in the icosahedral structure, although it still stays in shape. For κ = 2, however, it cannot maintain the optimal icosahedral arrangement of monomers and rather forms an entangled structure of longer, less bent, segments. As we have seen in the thermodynamic analysis of the transitions, the phase behavior changes significantly with increasing bending stiffness.
Representative lowest-energy conformations are shown in Figure 6. For a more quantitative analysis let us take a look at the pair distribution functions and the contact maps. We introduce the pair distribution function as where ∆(r − r i,j ) = 1, |r − r i,j | < r t , 0, otherwise.
As a robust threshold for the necessary binning of the r space, we chose r t = 0.01. The histograms for the lowest-energy conformations shown in Figure 6 are plotted in Figure 7. The perfect icosahedral structure is only found for κ = 0, whereas its decay is already visible for the weaker semiflexible polymer (κ = 1). The maximum number of nearest-neighbor contacts (r i,j ≈ 1) found for κ = 0 is not reached in the semiflexible cases. The broadening of the peaks and additional spikes not present for κ = 0, 1 are clear indicators that the ground-state structure for the semiflexible polymer with κ = 2 is not icosahedral. Bending restraints prevent the formation of perfect symmetries and thus a characterization of the κ = 2 ground-state conformation as a distinct crystalline or quasicrystalline structure is not possible. It rather resembles tertiary folds of protein conformations, where effective bending restraints and the local stable secondary segments purposefully prevent symmetric arrangements of monomers. This enables different heteropolymers of similar size to form distinct and functional individual conformations. This is obvious from the plots of the contact maps for the lowest-energy states shown in Figure 8. In the grid spanned by the monomer labels, we mark pairs of monomers with distances r i,j < 1.2. Bonded monomer pairs are not included. For the flexible polymer, the contact map does not exhibit particularly remarkable structural features. Since there are no bond angle restraints, each monomer tries to maximize its number of nearest neighbors for the energetic benefit. Icosahedral conformations with one of the end monomers in the center are optimal. For the semiflexible polymer with κ = 1, we already see that changes occur. Monomers do not only try to maximize the number of contacts, but three monomers connected by two adjacent bonds now have to cooperate to minimize bending. Close to the diagonal, we find that short anti-diagonal streaks form. These are clear indicators of turns with two linear strands in contact with each other (hairpins). In protein folds, these would be referred to as building blocks for β-sheets. For stronger bending rigidity (κ = 2), there are fewer, but longer such segments. We also observe the formation of streaks parallel to the diagonal, which can be associated with helical alignments. Both types of these secondary structure elements can be seen in the geometric representation in Figure 6. The hairpin-like sections are located in the interior and the helical segments wrap around it. Of course, under these conditions, energetically favored structures are not icosahedral. In fact, no global symmetries can be identified in this tertiary fold. Although it is tempting to think of these structures in analogy to protein folds, the lack of dihedral (torsion) constraints, which separates secondary structure elements from each other in the tertiary fold, prevents a more direct comparison.

Summary
The statistical analysis of structural transitions by the recently developed generalized microcanonical inflection-point method [4] had already yielded promising results in studies of flexible polymers [6]. For this study, we extended the coarse-grained model by incorporating bending restraints of different strengths to investigate the impact of bending stiffness on the microcanonical entropy and its derivatives, which are used as indicator functions for the systematic identification and classification of phase transitions in systems of any size. In this coarse-grained model for semiflexible polymers, the bending restraint is controlled by the bending stiffness parameter κ.
Despite its simplicity, the model required a careful numerical treatment. For this purpose, we performed extensive parallelized replica-exchange computer simulations and verified the structural behavior in the lowest-energy basin by means of stochastic methods optimized for this purpose such as simulated annealing. Microcanonical thermodynamic properties were obtained from the density of states, which is the main output of the multiple-histogram reweighting method applied to the raw simulation data.
Based on the microcanonical results obtained in the simulations, we compared the phase behavior of three model variants: purely flexible polymers (κ = 0) and semiflexible polymers with κ = 1 (identical energy scales of attractive nonbonded monomer-monomer interactions and repulsive bending strength) and κ = 2 (bending restraints dominate over monomer-monomer attraction). Our simulations reproduced previous results for the flexible reference system very well, creating sufficient confidence for the subsequent studies of the semiflexible polymers. Because of the additional restraints, it was significantly more challenging to achieve in simulations of semiflexible polymers the data quality necessary to enable an accurate microcanonical analysis.
For the flexible polymer, we identified the known structural transitions, i.e., the characteristic second-order transition associated with the Θ collapse from extended random coils to liquid globules and the first-order liquid-solid transition, which is accompanied by an independent third-order transition. The lowest-energy conformation is icosahedral, as expected. The behavior does not change qualitatively for the weaker semiflexible case with κ = 1, where the effectively repulsive bending effects are competitive to attractive monomer-monomer interaction. However, the larger number of possible attractive monomer-monomer contacts represents an entropic advantage and thus the icosahedral solid phase is still maintained. This is not the case anymore if the chain is further stiffened by choosing κ = 2. Whereas the Θ collapse remains widely unaffected, there is no obvious transition into the solid phase anymore. Signals at very low energies at the third-order transition level are too noisy to provide a clear picture. However, it is obvious that the lowest-energy conformation found is a compromise of satisfying both compactness and bending restraint. Longer, only slightly bent segments form and assemble in a well-organized, coil-like structure. It is an intriguing future task to extend our microcanonical study to stronger bending restraints and to construct the complete hyperphase diagram in the spaces of bending stiffness and inverse temperature and to apply it to other modern applications in polymer science like hyperbranched polymers [25].