Structural transitions in biomolecules - a numerical comparison of two approaches for the study of phase transitions in small systems

We compare two recently proposed methods for the characterization of phase transitions in small systems. The usefulness of these techniques is evaluated for the case of structural transition in alanine-based peptides.


Introduction
When phase transitions are studied in statistical physics, the basic assumption is usually that the dimensions of the macroscopic system are very large compared with that of the constituting elements. However, there are many phenomena in finite systems which resemble phase transitions, for instance, in the physics of clusters of atoms [1]. The main question is how the observed effects in small systems can be related to true phase transitions in macroscopic (or infinite) systems. Attempts in this direction include the exploration of the density of complex zeros of the canonical partition function for finite and small systems by Borrmann et al. [2] and the linear behavior for this limiting density [3,4]. The approach suggested by Janke and Kenna (JK) presents also a new scaling relation to identify the order and strength of a transition from the behavior of small systems [5].
In this paper we try to evaluate the usefulness and validity of the Borrmann et al. and the JK approaches for the investigation of structural transitions in biomolecules, another important example for "phase transitions" in small systems. For this purpose we consider the helix-coil transition in polyalanine and test whether the two approaches allow a characterization of that transition. In addition, we also investigate a second molecule Ala 10 -Gly 5 -Ala 10 in the context of Borrmann et al. approach.

Analysis of Partition Function Zeros in Small Systems
In the canonical ensemble a system is completely described by its partition function Introducing the variable u = exp(−kβ) with conveniently defined constant k allows to write the partition function for discrete energy models as a polynominal: We label the complex zeros in order of increasing imaginary part as u j = exp(−kβ j ), j = 1, 2, .... As established by Yang and Lee [6] and later by Fisher [7], the statistical theory of phase transitions can be described by the distribution of complex partition function zeros. In the case of a temperature driven phase transition, we expect the zeros (or at least the ones close to the real axis) condense for large enough system size L on a single line u j = u c + r j exp(iϕ) . (3) As the system size L increases, these zeros will move towards the positive real u-axis and for large L the corresponding value is β c , the inverse of the physical critical temperature T c . Crucial information on phase transitions can be obtained from the way in which the first zero approaches the real u-axis. However, such an analysis depends on the extrapolation towards the infinite large system and does not allow characterization of the situation in small systems. One possible extension of the above ideas to "phase transitions" in small systems is the classification scheme by Borrmann et al. [2]. This classification scheme computes density of zeros taking into account numerical estimates of the first four complex zeros. Now, writing the complex zeros as u k = Re(u k )+i τ k , where τ k stands for Im(u k ), the assumed distribution of zeros on a straight line allows to define two parameters α u and γ u : with k labelling the first zeros, and Note, that our notation differs from that in Ref. [2] in that we define the discrete line density φ in function of the u-zeros instead of the inverse temperature β. Following the classification scheme by Grossmann and Rosenhauer [3,4], phase transitions can now be classified according to the values of these two parameters: for α u = 0 and γ u = 0 one has a phase transition of first order, it is of second order if 0 < α u < 1 and arbitrary γ u , and for α u > 1 and arbitrary γ u one has a higher order transition. Of course, as expected for a phase transition, the imaginary part of the leading zero (τ 1 ) should move towards the real axis to obtain in the thermodynamic limit a real temperature. However, for small systems with finite τ 1 , we may have α u < 0 [2], which translates the condition of first order transition to α u ≤ 0 and γ u = 0. This classification of phase transitions has been tested for the finite Bose-Einstein condensates in a harmonic trap, small magnetic clusters and nuclear multifragmentation [2]. Another extension of partition function zero analysis to small systems is the approach by Janke and Kenna (JK) [5], which uses that the theoretical average cumulative density of zeros [5] can be written in the thermodynamic limit and for a first order transition as Here, the slope at the origin is related to the latent heat, ∆e ∝ g ∞ (0). Equations (3) and (7) imply that the distance r j of a zero from its critical point can be written for large enough lattice sizes as Im u j (L) since Re u j (L) ∼ u c . Hence, in this limit Eq. (6) and (7) lead to the following scaling relation for the cumulative density of zeros as an equation in j and L, A necessary condition for the existence of a phase transition is that a 3 is compatible with zero, else it would indicate that the system is in a well-defined phase. The values of the constants a 1 and a 2 then characterize the phase transition. For instance, for first order transitions should the constant a 2 take values a 2 ∼ 1 for small r, and in this case the slope of this equation is related to the latent heat through the relation [5] ∆e = k u c 2πa 1 , with u c = exp(−kβ c ). On the other hand, a value of a 2 larger than 1 indicates a second order transition whose specific heat exponent is given by α = 2 − a 2 . The above approach was originally developed and tested for systems with well defined first order phase transitions such as the 2D 10-state Potts and 3D 3-state Potts model. The obtained results agree with previous work from numerical simulations and partition function zeros analysis of systems up to L = 64 [8,9] (L = 36 for the 3D case [10]). However, some difficulties in the identification of the order of the transition appear when applied to a model with a very weak first order transition [11].
Here we explore further these two approaches for a non-trivial example of structural transitions in biomolecules, namely the helix-coil transition in polyalanine. The characteristics of this so-called helixcoil transition have been studied extensively [12]. For instance, evidence was presented [13,14,15] that polyalanine exhibits a phase transition between the ordered helical state and the disordered random-coil state when interactions between all atoms in the molecule are taken into account. These previous results, obtained independently by different methods, are used here as benchmarks for the test of the usefulness and to point out the difficulties in applying the classification schemes of Borrmann et al. and of Janke and Kenna for research of structural transitions in biomolecules.
Our investigation of the helix-coil transition for polyalanine is based on a detailed, all-atom representation of that homopolymer. Since one can avoid the complications of electrostatic and hydrogenbond interactions of side chains with the solvent for alanine (a nonpolar amino acid), explicit solvent molecules were neglected. The interaction between the atoms was described by a standard force field, ECEPP/2, [16]. Chains of up to N = 30 monomers were considered, and our results rely on multicanonical simulations [17] of N sw Monte Carlo sweeps starting from a random initial conformation, i.e. without introducing any bias. Our statistics consists of N sw =400,000, 500,000, 1,000,000, and 3,000,000 sweeps for N = 10, 15, 20, and 30, respectively. Measurements were taken every fourth Monte Carlo sweep. Additional 40,000 sweeps (N = 10) to 500,000 sweeps (N = 30) were needed for the weight factor calculations by the iterative procedure described first in Ref. [17]. In contrast to our first calculation of complex zeros presented in Ref. [14], where we divided the energy range into intervals of lengths 0.5 kcal/mol in order to make Eq. (2) a polynomial in the variable u = e −β/2 , we avoided any approximation scheme in the present work. This because the above approximation works very well for the first zero, but not for the next ones. Since we need high precision estimates also for the next zeros we applied the scan method (see Ref. [18] and references therein). In Table 1 we present our first four partition function zeros, although the fourth one is less reliable due to the presence of fake zeros [19]. Using these zeros we first calculated the parameters α u and γ u which characterize in the Borrmann et al. approach phase transitions in small systems. As one can see in Table  2 the so obtained values for polyalanine are characterized by large fluctuations. It seems that the median of the α u values is α u = 0 which would indicate a first order transition. However, our data are not good enough to draw such a conclusion. Our error estimates were obtained by means of the Bootstrap method [20,11]. For this reason, we tried instead the JK scaling relations. Table 3 lists the parameter a 3 (N) of  Eq. (8). Here, the average cumulative density of zeros is replaced by where we have translated the linear length L as N 1/d [14]. Therefore all finite size scaling relations can be written in terms of the number of monomers N. In Fig. 1 we show the cumulative distribution of zeros in a log-log plot for chain lengths N = 10, 15, 20 and 30. The values a 3 are compatible with zero for chains of all length indicating that we have indeed a phase transition. In order to evaluate the kind of transition we also calculate the parameters a 1 (N) and a 2 (N) which we also summarize in Table 3. Our analysis procedure differs from JK in that we are calculating a 2 for a fixed N before performing a finite-size scaling analysis of that quantity, i.e. we account for the dependence of the average cumulative density on the system size. The parameter a 2 (N) decreases with system size and a log-log plot of this quantity as a function of chain length, as shown in Fig. 2, suggests a scaling relation  A numerical fit of our data to this function leads to a value of a 2 = 1.31(4) with goodness of fit Q = 0.95. Using α = 2 − a 2 we find α = 0.70(4) which is barely compatible with our previous value of α = 0.86 (10) in Ref. [14], obtained from the maximum of specific heat. A fit of all four chain lengths can also not exclude a value a 2 = 1 (i.e. a first order phase transition) since we can find acceptable fits with Q > 0.55 in the range 0.92 < a 2 < 1.44. Especially, a close examination of Fig. 1 and Fig. 2 shows that the N = 30 data point exhibits a considerable deviation from the trend suggested by the smaller chain lengths. Since the N = 30 data are the least reliable, we also evaluated Eq. (11) omitting the N = 30 chain. This leads to a value of a 2 = 1.16(1) and a critical exponent α = 0.84(1) which is now compatible with our previous result α = 0.86 (10).
While the JK approach is able to reproduce results for polyalanine obtained in previous work [14], it does not allow to establish the order of the helix-coil transition from simulation of small chains. Our results for the parameter a 2 seem to favor a second order transition, but are hampered by large errors (such that a first order transition cannot be excluded) and disputed in Ref. [15] were indications for a finite latent heat were found.
The Borrmann et al. approach led to even less decisive results for polyalanine since the fluctuation were too large to draw firm conclusions. However, the advantage of this approach is that it does not require a finite-size scaling. We have therefore used this method to investigate a second molecule, the slightly more complicated Ala 10 -Gly 5 -Ala 10 . A detailed study of this molecule will be published elsewhere [21].
We describe the interactions between the atoms again by the ECEPP/2 [16] force field (as implemented in the program package SMMP [22]). Our results rely on multicanonical simulations [17] of 4, 000, 000 Monte Carlo sweeps starting from a random initial conformation, i.e. without introducing any bias. Measurements were taken every tenth Monte Carlo sweep. Additional 500,000 sweeps were needed for the calculations of the weight factor.
Analyzing the partition function zeros for this peptide, we find a "phase transition" at two temperatures, each being characterized by a line of complex zeros. The corresponding first four zeros for each characteristic line are listed in Table 4. Both temperatures also correspond with peaks in the specific heat which is displayed in Fig. 3 as a function of temperature. In Ref. [21] it is shown that the higher temperature T 1 = 480 K marks the helix-coil transition, i.e. the temperature where secondary structure elements (α-helices) are formed. At the lower temperature T 2 = 260 K, the peptide then folds into its native structure, a U-turn-like bundle of two (antiparallel) α-helices connected by a turn. Here, we use the Borrmann et al. approach to study these two transitions in more detail. Using the zeros from Table 4 we calculated the parameters α u and γ u which characterize in the Borrmann et al. approach phase transitions in small systems. For the first transition, at T = 480 K, which corresponds to the above studied helix-coil transition in polyalanine, we find α u = −1.46 and γ u = 0.49. While the negative value of α u indicates a first order transition, our results do again not allow us to draw a firm conclusion on the nature  of the transition since the values of the two parameters α u and γ u are subject to large fluctuations when we vary the number of sweeps considered in the calculation of the zeros. This is not surprising since we were also not able to establish clearly the order of the helix-coil transition in the previous example of polyalanine. However, our data are more decisive in the case of the second transition, at T = 260 K, which marks the compactification and folding of the peptide. Here we find α u = 0.38 and γ u = −1.07. These values indicate a second-order transition which is consistent with what one would expect for a transition between extended and compact structures. We hope that by doubling the number of Monte Carlo sweeps we will be able to verify this results and also obtain reliable values for the Borrmann et al. parameter characterizing the helix-coil transition at T = 480 K. These simulations are now under way.

Conclusion
Let us summarize our results. We have evaluated two recently proposed schemes for characterizing phase transitions in small systems. Simulating polyalanine molecules of chain lengths up N = 30 residues by multicanonical Monte Carlo, we calculated the partition function zeros of these molecules. Analyzing these zeros by the JK approach we were able to reproduce for polyalanine results obtained in previous work [14] while the Borrmann et al. approach led here to inconclusive results. However, the later method seems more appropriated to study structural transitions in proteins since it does not require a finite size scaling. When applied to Ala 10 -Gly 5 -Ala 10 , the Borrmann et al. approach allowed us to characterize structural transitions in this molecule. However, our results are limited again by the precision with which one can calculate the first four partition function zeros. While these numerical limitations restrict the applicability of the Borrmann et al. and the JK approaches, our results demonstrate that both methods can be useful tools for the study of "phase transitions" in biomolecules.