Partition Function Zeros and Heat Capacity Decomposition Reveal HP Protein Foldability

TL;DR

This paper shows how a thermodynamic method can be used to evaluate how well proteins fold based on their heat capacity and transition behavior.

Contribution

The extension of heat capacity decomposition using partition function zeros to assess HP protein foldability is novel.

Findings

Significant crossovers in thermal transitions indicate slow folding behavior.

Absence of crossovers characterizes good protein folders.

Kinetic Monte Carlo simulations validate the foldability criterion.

Abstract

The heat capacity decomposition method, a well-established analytical approach in polymer thermodynamics for elucidating thermal transitions in homogeneous polymers, is extended here to heterogeneous systems. We demonstrate that the decomposition of heat capacity based on partition function zeros allows the identification of transition-like crossovers originating from compact low-energy states, thereby enabling the evaluation of the foldability of HP sequences. The occurrence of significant crossovers between the collapse and folding transitions indicates slow folding behavior, whereas their absence characterizes good folders. This criterion is further validated through kinetic Monte Carlo simulations of two representative sequences.

Figures

Partition function zeros for sequence B. Red and blue dots indicate the first zeros associated with the collapse and folding transitions, while there are intermediate orange dots between them.

Heat capacity component CV,k of sequence B. The dashed line represents the heat capacity CV. The red and blue lines, respectively, represent the heat-capacity components contributed by the red and blue first zeros shown in Figure 3.

The native states for (a) sequence A and (b) sequence B. The dark blue beads represent H amino acids, and the pink beads represent P amino acids.

Full text

1. Introduction

Protein folding is fundamentally a thermodynamic process, at least for small protein [1]. The unique native state is the minimum free energy state. How a protein is able to locate its native state among an astronomical number of excited conformations so rapidly remains a central puzzle [2]. A widely accepted explanation is that evolved proteins fold on a free-energy funnel landscape [3,4]. In practice, however, depicting the funnel is challenging due to the complexity of the vast, high-dimensional conformational space. The theoretical study of protein folding has a long history [5,6,7,8,9]. More recently, AlphaFold [10,11,12] has achieved a groundbreaking breakthrough in protein structure prediction. It can predict the three-dimensional structure of a given amino acid sequence with several hundred residues in minutes or hours—tasks that previously could take years. The unexpected success of AlphaFold relies on innovations in neural network architectures, accelerated by advances in computer hardware, and the availability of extensive sequence and structure data in protein databases. It is essentially a computational, knowledge-based approach, which differs significantly from previous theoretical methods. While AlphaFold has achieved impressive accuracy in structure prediction, a fundamental theory of protein folding based on basic physical principles still remains to be fully established.

In this paper, we utilize partition function zeros, a technique from statistical mechanics, to probe the foldability of protein sequences. The partition function is the central quantity in statistical mechanics, from which all thermodynamic functions can be derived. Partition function zeros are defined as the roots of the equation “partition function = 0”. These roots fully determine the properties of the partition function and are distribute in the complex temperature plane (for the canonical ensemble) or in the complex fugacity or magnetic field plane (for the grand canonical ensemble). Since the pioneering work of Yang and Lee [13,14], which related partition function zeros to singularities in phase transitions, steady theoretical development has been made in understanding the their behaviors. From the loci and density of partition function zeros, one can extract critical properties such as critical points, the order of phase transitions, and critical exponents [15,16,17]. However, when applied to systems with competing interactions—such as the heterogeneous HP model [18], in contrast to homogeneous interacting self-avoiding walks—the zeros typically form scattered patterns without clear loci [19,20]. To address this limitation, we developed the strategy of heat capacity decomposition [21], in which each partition function zero is associated with a corresponding heat capacity component, thereby providing the zero with a physical interpretation. This additional information enables the analysis of complex patterns of partition function zeros to become feasible. A brief introduction of heat capacity decomposition is provided in Section 2.

The HP model is the simplest coarse-grained model of protein folding. In this model, an HP chain is a self-avoiding walk composed of two types of amino acid beads. The 20 natural amino acids are classified as either hydrophobic (H) or polar (P), as hydrophobicity is considered the dominant driving force of protein folding [22,23]. The energy is defined as

$$[eqn]$$

where the amino acid types $[eqn]$ and $[eqn]$ can be either H (Hydrophobic) or P (Polar), and $[eqn]$ means the ith amino acid and jth amino acid have nonconsecutive nearest-neighbor contacts. There are three types of contact energies, $[eqn]$ , which satisfy the inequalities $[eqn]$ and $[eqn]$ . Here we adopt an integer-scaled version of the interaction parameters: $[eqn]$ , corresponding to the values $[eqn]$ , and $[eqn]$ used in [24]. With integer-valued energies, the partition function takes the form of a polynomial, which makes root-finding more tractable. We test other interaction parameters as well and find they does not lead to significant qualitative differences. The HP model, much like the Ising model, is related to many other problems. For example, determining the ground state of the HP model is an NP-complete optimization problem [25]. It also serves as a testing ground for numerous computational algorithms [26]. Constructing basic principles with such simplified models is often more efficient, and these principles are expected to remain valid for more realistic models as well.

After analyzing the folding properties of HP sequences using partition function zeros and heat capacity decomposition, we verified the results with Gillespie-type kinetic Monte Carlo simulations. Unlike the traditional Metropolis algorithm, this approach focuses on simulating real physical motion. Details of this method are provided in the next section.

2. Method

2.1. Heat Capacity Decomposition

With integer-valued contact energies, the partition function of an HP protein sequence with N amino acids can be expressed as a polynomial in the variable $[eqn]$ :

$$[eqn]$$

where $[eqn]$ denotes the number of structures with energy E, M = $[eqn]$ is the degree of the polynomial, and $[eqn]$ is the k-th zero of the partition function. The heat capacity can be obtained from the partition function as $[eqn]$ and $[eqn]$ , where U is the internal energy and $[eqn]$ . After changes of variables,

$$[eqn]$$

Since $[eqn]$ is a polynomial in x, its roots occur as complex conjugate pairs, with the exception of very few real roots located on the negative real axis. For each complex or real $[eqn]$ ,

$$[eqn]$$

Thus, we may define the heat capacity component as

$$[eqn]$$

and

$$[eqn]$$

In this way, the heat capacity $[eqn]$ is decomposed into individual components $[eqn]$ , each of which is a real function of the real variable x and is contributed by the zero $[eqn]$ .

In the thermodynamic limit, i.e., as the system size $[eqn]$ , the partition function zeros will touch the positive real axis if a phase transition occurs [13,14,27] and $[eqn]$ will diverge. For finite systems, the divergence manifests as a peak, and the peak position can reasonably be defined as the transition point. The roots closest to the real axis, refered to the first zeros, provide the leading contribution of $[eqn]$ due to $[eqn]$ term in the denominator, where $[eqn]$ denotes the first zero. Near the real part of $[eqn]$ , the associated heat capacity component exhibits a pronounced peak, the position of which can also be used to define the transition point of the phase transition in small systems [15,28,29].

Similar to an interacting self-avoiding walk (SAW) [30], a protein-like HP chain undergoes a collapse transition, the most common phase transition in polymers, moving from coil to globule states. A second transition, the folding transition from the globule to the native state, is inherently dependent on the specific amino acid sequence. Using the peak positions of the corresponding heat capacity components, the collapse temperature $[eqn]$ and the folding temperature $[eqn]$ can be defined. Both transitions contribute to the heat capacity, causing the shape of $[eqn]$ to vary markedly across different sequences. The $[eqn]$ peaks of different transitions often mask one another. Through heat capacity decomposition, the overlapping contributions of the two transitions can be separated, providing a clearer view of their respective roles. This heat capacity decomposition technique is specifically designed for small systems that undergo multiple phase transitions, typically indicated by multiple or irregularly shaped heat capacity peaks. Additional details on the method can be found in [21].

2.2. Kinetic Monte Carlo Simulation

We employ the kinetic Monte Carlo simulation based on the Gillespie algorithm [31,32] to detect the folability of HP sequences. The principles and procedures of this method are briefly described below. Suppose there are K possible moves for a HP chain. The transition probability per unit time for each move i is $[eqn]$ , with $[eqn]$ . Let $[eqn]$ represent the probability that move i occurs at time t. This probability can be expressed as $[eqn]$ , where $[eqn]$ is the probability that no move occurs before t, and $[eqn]$ is the probability that the move i occurs within the interval $[eqn]$ . Deviding t into n equal parts, $[eqn]$ can be writen approximately as $[eqn]$ , where $[eqn]$ . In the limit $[eqn]$ , $[eqn]$ . Thus, we get the final result:

$$[eqn]$$

In other words, the waiting time t for one move to occur follows an exponential distribution, and the probability that the move is type i is given by $[eqn]$ . The value of t can be simulated using the standard transformation method:

$$[eqn]$$

where $[eqn]$ is an uniformly distributed stochastic variable in the range $[eqn]$ . It is worth noting that $[eqn]$ is updated after every move; hence, the exponential function and the probability ratio $[eqn]$ are not fixed. Howerver, $[eqn]$ remains independent of these quantities, and $[eqn]$ = $[eqn]$ . Therefore, $[eqn]$ , which is consistent with expectations.

The move probability $[eqn]$ here is taken to be the Metropolis type:

$$[eqn]$$

i.e., for $[eqn]$ , $[eqn]$ , and for $[eqn]$ , $[eqn]$ , where $[eqn]$ is the energy change associated with move i. When considering actual physical time, the rate $[eqn]$ must be multiplied by an additional constant, which represents the frequency of attempted moves for an amino acid. One step of the kinetic Monte Carlo simulation consists of selecting a move i from all possible moves according to the probability ratio $[eqn]$ , and advancing the simulation time by $[eqn]$ . After this move is performed, the list of possible moves needs to be updated. The procedure is then repeated to select subsequent move and update the system. Through successive moves, the system approaches thermal equilibrium at temperature T. We record the total time when the HP chain first reaches the lowest-energy state, i.e., the native state. This first-passage time is relevant to the kinetic accessibility of the native state at a constant temperature and can serve as a measure of the folding ability of a sequence [33].

3. Result

We consider all HP sequences with 16 HP amino acids on the simple cubic lattice. There are 2^16^ = 65,536 different HP sequences. Using exact enumeration [34], we compute the partition functions of all these sequences and filter out those without a unique ground state. This procedure leaves 5248 sequences as candidates for HP protein sequences. We then calculated the partition-function zeros of all these sequences and performed heat-capacity decomposition for several hundred of them. Table 1 presents 20 sequences randomly selected from the examined set. The table lists their ground-state energy $[eqn]$ , $[eqn]$ and $[eqn]$ . It can be observed that while $[eqn]$ varies in a relatively regular manner, $[eqn]$ exhibits considerable irregularity. This behavior closely resembles that found in the three-dimensional ISAW model [21,30], where the freezing transition temperature is irregular with respect to chain length because it strongly depends on the specific ground-state structures.

We can roughly classify the 20 sequences into two groups based on their heat capacity components between $[eqn]$ and $[eqn]$ . The first group shows no significant heat capacity peaks in this range, while the second group does. This classification should be regarded as qualitative, since the significance of the peaks varies in a continuous manner. The appearance of these peaks in the second group suggests the existence of low-energy states that compete with the native state, forcing the native state to emerge only at lower temperatures and thus reducing $[eqn]$ . Consequently, $[eqn]$ is generally higher in the first group and lower in the second. To further clarify these points, we choose two sequences in Table 1 as representatives of the two groups to present the results in details.

The two sequences are sequence No. 6 from the first group (referred to as sequence A) and sequence No. 18 from the second group (referred to as sequence B). They are not exceptional; similar conclusions would be obtained if other sequences were selected as examples. Sequence A can be regarded as a protein-like sequence that folds more rapidly, whereas Sequence B is a non-protein-like sequence that folds more slowly than Sequence A because it must pass through a crossover to reach the native state. Here, the term crossover refers to a smooth change in the physical state, not as abrupt as a phase transition, but still capable of producing a peak in the heat capacity. Finally, we will perform kinetic Monte Carlo simulations to verify these predictions.

The partition functions of sequences A and B are provided in the Appendix A. Following the convention described in [35], and considering both rotational and chiral symmetries, the degeneracies of the linear, planar, and three-dimensional structures are 1, 4, and 8, respectively. Consequently, the coefficient of the leading term of all 5248 partiton functions, $[eqn]$ , is always equal to 8. The partition function zeros, i.e., the roots of the two polynomials, are solved with Mathematica.

Figure 1 shows the the partition function zeros of sequence A, and Figure 2 displays the associated heat capacity components $[eqn]$ in the same colors. Note that in the two figures, the variables x and T are inversely related: $[eqn]$ , where $[eqn]$ is set to 1 by adjusting the temperature unit. The locus of the partition function zeros reveals a clear inner ring corresponding to the collapse transition and a diffuse outer ring corresponding to the folding transition. The first zeros are marked in red and blue, respectively. Since two complex-conjugate zeros contribute the same heat capacity component, we plot $[eqn]$ for the pair of $[eqn]$ and its complex conjugate, except in very few cases where the zeros lie on the negative real axis; in those cases, only the heat capacity component contributed by one zero is shown. Individual heat capacity components may take negative values. However, when summed, these negative parts are cancelled out by other positive parts, resulting in an entirely positive total heat capacity. The peak regions of the components contribute most to the total heat capacity, and these are, of course, always positive.

Let us compare the results of sequence B with those of sequence A. Figure 3 and Figure 4 show the the partition function zeros and associated heat capacity components of sequence B. The ring of partition function zeros corresponding to the collapse transition remains evident; however, additional orange zeros appear between the red and blue zeros, such that the locus of zeros resembles three rings rather than two. The orange heat capacity components likewise exhibit significant peaks between the red and blue peaks. A reasonable explanation is that, in this case, more metastable compact globules compete with the native state. As shown in Figure 5, the native state of sequence A is efficiently designed to maximize HH contacts and to greatly reduce the possibility of competing low-energy structures. In contrast, the native state of sequence B is not uniquely favored.The number of structures $[eqn]$ , listed in the Appendix A, further indicates that sequence B possesses a greater number of low-energy states than sequence A. Since an abundance of metastable states increases the likelihood of kinetic trapping, the folding process is slowed down. Consequently, the folding kinetics of sequence B is expected to be slower than that of sequence A. To verify this prediction, we perform kinetic Monte Carlo simulations for both sequences.

Partition function zeros for sequence B. Red and blue dots indicate the first zeros associated with the collapse and folding transitions, while there are intermediate orange dots between them.

Heat capacity component CV,k of sequence B. The dashed line represents the heat capacity CV. The red and blue lines, respectively, represent the heat-capacity components contributed by the red and blue first zeros shown in Figure 3.

The native states for (a) sequence A and (b) sequence B. The dark blue beads represent H amino acids, and the pink beads represent P amino acids.

In our Monte Carlo simulations, trial moves are restricted to corner flips and 90-degree end flips, while 180-degree end flips are excluded. Thus, in each Monte Carlo step, only one bead moves, and it is displaced by $[eqn]$ lattice spacing. Other local moves, such as crankshaft moves and pull moves [36], are not included because they involve the movement of multiple beads. We employ this small-step motion to mimic physical movement. Nonlocal moves, such as pivot moves [37] and bond-rebridging moves [38], are considerably more efficient for searching the native state and computing the density of states; however, since they are not physical, they are also excluded from the simulations.

At each fixed temperature, we conduct 10,000 kinetic Monte Carlo simulations to measure the first-passage time required for a HP sequences to reach its native state from a straight chain. Figure 6 shows the mean first-passage time at various temperatures. Each data point has a relative error of about 1%. The curves exhibit a typical U-shaped profile, a pattern commonly observed in various optimization problems [39]. The underlying reason is that, at high temperatures, entropic effects dominate, reducing the probability of reaching the native state, and at low temperatures, energetic effects dominate, leading to more severe trapping in metastable states.

The minimum of the U-curve indicates the optimal temperature for reaching the ground state. This temperature is approximately 4.5 for sequence A and 5.0 for sequence B. For both sequences, this optimal temperature lies between $[eqn]$ and $[eqn]$ . The specific values are $[eqn]$ and $[eqn]$ for sequence A, and $[eqn]$ and $[eqn]$ for sequence B. Across all temperatures, sequence A reaches its native state more rapidly than sequence B. This result confirms the prediction based on the analysis of partition function zeros and heat capacity decomposition. At higher and lower temperatures, data with reduced statistical precision also support this conclusion.

Klimov and Thirumalai [40] proposed a criterion for the foldability of proteins: a sequence is considered an efficient folder if $[eqn]$ is small. The present results are consistent with this criterion. Furthermore, heat capacity decomposition provides additional insight into the critical region between $[eqn]$ and $[eqn]$ . When the collapse and folding transitions are separated by intermediate crossovers induced by metastable compact globules, the folding speed of an HP chain may be reduced. This scenario is exemplified by sequence B. Being grounded in thermodynamic properties, this new criterion exhibits a high degree of generality and applicability to sequences across different models.

4. Discussion

Heat capacity decomposition based on partition function zeros is an entirely new approach that focuses directly on thermodynamic properties rather than relying on statistical perspectives. It provides an X-ray–like view of the heat capacity and offers new insight into the sequence’s folding behavior. In this paper, we apply heat capacity decomposition to the HP model for simplicity; however, this approach is general and easily extendable to other model types, including off-lattice or all-atom models. For instance, Alves and Hansmann [41] demonstrated the applicability of the partition function zeros approach to systems using the standard all-atom ECEPP/2 force field. By discretizing the energy into intervals of equal lengths, the partition function remains a polynomial. Therefore, the partition function zeros are readily calculable, permitting the application of the heat capacity decomposition method.

Our ongoing follow-up work focuses on two primary tasks: the design of a folding network to explore metastable states and the computation of partition functions for longer chains. In addition to sequences A and B, we have analyzed the partition function zeros of all the remaining 5246 HP sequences. The inner circle of zeros associated with the collapse transition exhibits a high degree of similarity, suggesting that this transition is a universal characteristic of polymers and is not dependent on the specific sequence. In contrast, the zero distribution corresponding to the folding transition is quite divergent, as illustrated by sequences A and B. Its pattern primarily depends on the number of metastable states that the sequence can form. Metastable states compete with the native state, leading not only to kinetic trapping but also to a weakened folding transition. Exploring these metastable states requires a depiction of the energy landscape, which is a technically challenging task. The kinetic Monte Carlo simulations employed in this study can be used to map the energy landscape. The state space can be represented by a folding network, in which only low-energy states serve as nodes. The weighted connections between these nodes are obtained from Monte Carlo trajectories and can be converted into distance metrics. Combining this two-dimensional folding network with one-dimensional energy yields a simplified three-dimensional energy landscape. The folding properties of a sequence can be inferred from whether its energy landscape is funnel-like or glassy. The folding network is currently under construction.

The most time-consuming step of the entire calculation is computing the partition function. We employ exact enumeration, which can generate partition functions for longer chains when not all sequences are considered. In the future work, we plan to select several 27-residue HP sequences with well-designed native structures that fit into a 3 × 3 × 3 cube and compute their partition functions. Longer chains are advantageous because they produce a denser and smoother locus of partition function zeros, thereby enabling more quantitative analyses. For even longer chains, approximation methods are necessary. Wang–Landau sampling [42] is a type of multicanonical algorithm in which the density of states is determined by adaptively updating to achieve flat-histogram sampling of energies. Wang-Landau simulation has been shown to be a practical and effective method to obtain the partition functions in the HP model [26,43]. To ensure comprehensive sampling of the conformational space, this simulation will incorporate a complete set of efficient moves, including pull, pivot, and bond-rebridging moves.

In summary, we use heat capacity decomposition based on partition function zeros to investigate the foldability of HP sequences. This approach provides a novel perspective not available in previous methods. It can be directly applied to various other protein folding models and is a promising approach for future research.