format: text

text:
# Energetics of Protein Thermodynamic Cooperativity: Contributions of Local and Nonlocal Interactions

## Abstract

The respective roles of local and nonlocal interactions in the thermodynamic cooperativity of proteins are investigated using continuum (off-lattice) native-centric Gō-like models with a coarse-grained Cα chain representation. We study a series of models in which the (local) bond- and torsion-angle terms have different strengths relative to the (nonlocal) pairwise contact energy terms. Conformational distributions in these models are sampled by Langevin dynamics. Thermodynamic cooperativity is characterized by the experimental criteria requiring the van’t Hoff to calorimetric enthalpy ratio $`\mathrm{\Delta }H_{\mathrm{vH}}/\mathrm{\Delta }H_{\mathrm{cal}}1`$ (the calorimetric criterion), as well as a two-state-like variation of the average radius of gyration upon denaturation. We find that both local and nonlocal interactions are critical for thermodynamic cooperativity. Chain models with either much weakened local conformational propensities or much weakened favorable nonlocal interactions are significantly less cooperative than chain models with both strong local propensities and strong favorable nonlocal interactions. These findings are compared with results from a recently proposed lattice model with a local-nonlocal coupling mechanism; their relationship with experimental measurements of protein cooperativity and chain compactness is discussed.

## 1 Introduction

How a globular protein can fold reliably into a particular three dimensional conformation *in vitro*, without the participation of molecular chaperones, is a central puzzle in biophysics. If we wish not only to predict the folded state of a protein, but also to understand the folding phenomenon in terms of physical processes, we need to use physics-based methods: we run computer simulations of self-contained polymer models  that attempt to mimic the behavior of real protein molecules. A complete quantum mechanical simulation, which would include the solvent molecules in addition to all the atoms in the protein molecule, is not yet possible. But in attempting to design simplified models, we face the problem of how to simplify: which characteristics are essential and which can be neglected? What effective energy functions does this imply for the simplified system?

Part of this general question is addressed in this article. Folding experiments on small globular proteins have long shown evidence of thermodynamic and kinetic cooperativity , which indicates a phenomenon similar to a first order phase transition between native and denatured states. As our group has argued recently , this observation can be exploited to constrain the set of possible simplified models and interaction schemes: for a particular simplified model to be a quantitatively accurate representation of protein thermodynamics and kinetics, it is essential that, when appropriately applied to a small globular protein, it can produce the experimentally observed generic cooperative behavior.

Such constraints turn out to be rather stringent. It is nontrivial to construct model interaction schemes that can produce proteinlike cooperativities . A case in point is a class of common Gō-like  models . Their potential functions are native-centric, in that they are explicitly biased to favor a given native structure. Gō-like modeling of proteins has provided important physical insights . These include an increasing number of elegant elucidations of functional protein dynamics under native conditions . As for global folding and unfolding of proteins (in contrast to their near-native dynamics), a detailed discussion of the merits and limitations of Gō-like approaches can be found in Ref. . Notably, common Gō-like models do not appear capable of producing simple two-state folding/unfolding kinetics. Instead, their chevron plots exhibit severe rollovers, which are typical of the class of folding kinetics that is customarily referred to as non-two-state . Nonetheless, common three-dimensional Gō-like protein models seem sufficient to produce apparent two-state thermodynamic behavior , although their two-dimensional counterparts fail to do so .

In the present investigation, we limit our scope to thermodynamic cooperativity. Specifically, we aim to explore how protein thermodynamic two-state-like behavior is affected by the relative strengths of local interactions (between residues close together along the chain sequence) and nonlocal interactions (between residues far apart along the chain sequence). The respective roles of local and nonlocal interactions are an issue of long standing interest in the study of protein energetics , and the effect of analogous interactions on the phase diagram of lattice polymers has also been investigated . Here the issue is addressed by varying the potential function in a series of coarse-grained Gō-like models, which represent the protein as a string of Cα positions in continuum space and which are simulated using Langevin dynamics. In view of the limitations of common Gō models , the present study should be viewed as a first step in tackling the issue of local vs. nonlocal interactions in cooperative continuum protein models. To assess the robustness of our conclusions, results from the continuum Langevin models are also compared with results from lattice model simulations.

We begin in section 2 by providing details of the models. An outline of the thermodynamics involved in interpreting the simulations is given in section 3. Our findings are presented in section 4, and we conclude in section 5 with a discussion of the implications of our results.

## 2 Models and simulation details

### 2.1 Continuum models

For the present continuum Gō-like models we use a representation, introduced by Clementi *et al.* , of the 64-residue truncated form of chymotrypsin inhibitor 2 (CI2). The native contact set corresponds to NCS2 in Ref. .

We use an energy function that is similar to one used previously . The potential energy function $`V`$, from which the conformational force is derived, is given by

$$V=V_{\mathrm{stretching}}+V_{\mathrm{bending}}+V_{\mathrm{torsion}}+V_{\mathrm{native}}+V_{\mathrm{nonnative}},$$
(1)
where

$$V_{\mathrm{stretching}}=\underset{i=1}{\overset{N1}{}}k_l\left(l^il_0^i\right)^2$$
(2)
contains a summation over the virtual bonds between pairs of residues,

$$V_{\mathrm{bending}}=\underset{i=1}{\overset{N2}{}}\epsilon _\theta (\theta ^i\theta _0^i)^2$$
(3)
involves a summation over the virtual bond angles between triplets of residues, and

$$V_{\mathrm{torsion}}=\underset{i=1}{\overset{N3}{}}\left\{\epsilon _\varphi ^{(1)}\left[1\mathrm{cos}\left(\varphi ^i\varphi _0^i\right)\right]+\epsilon _\varphi ^{(3)}\left[1\mathrm{cos}3\left(\varphi ^i\varphi _0^i\right)\right]\right\}$$
(4)
represents the virtual torsional potential between quadruplets of residues. The latter contains a term with a single minimum as well as the traditional three-minimum term . \[We note that there is an apparent typographical error in the corresponding $`V_{\mathrm{torsion}}`$ in , which effectively lists these terms as $`1+\mathrm{cos}\left(\varphi ^i\varphi _0^i\right)`$ and $`1+\mathrm{cos}3\left(\varphi ^i\varphi _0^i\right)`$. But such terms would fold the chain into the mirror image of the PDB structure.\] $`V_{\mathrm{stretching}}`$, $`V_{\mathrm{bending}}`$ and $`V_{\mathrm{torsion}}`$ together account for the local interactions (between residues that are separated by no more than three places along the chain), which include local conformational propensities for the native structure. The local interactions are expressed in this way because it biases the local geometry of the chain. The fourth term,

$$V_{\mathrm{native}}=\underset{\left|ij\right|4}{}\epsilon _{\mathrm{native}}\left[5\left(\frac{r_0^{ij}}{r^{ij}}\right)^{12}6\left(\frac{r_0^{ij}}{r^{ij}}\right)^{10}\right],$$
(5)
sums over the pairwise interactions between residues that are regarded as being in contact in the native structure; this accounts for the nonlocal interactions (between residues that are separated by four or more places along the chain). Finally,

$$V_{\mathrm{nonnative}}=\underset{\left|ij\right|q}{}\epsilon \left(\frac{r_{\mathrm{rep}}}{r^{ij}}\right)^{12}$$
(6)
contains repulsive pairwise interactions between other pairs of residues, in order to ensure the self-avoidance of the chain.

The chain contains $`N`$ residues. $`l^i`$ is the length of virtual bond $`i`$, $`\theta ^i`$ is a bond angle, $`\varphi ^i`$ is a dihedral angle and $`r^{ij}`$ is the distance between two residues $`i`$ and $`j`$. The corresponding values in the native structure are $`l_0^i`$, $`\theta _0^i`$, $`\varphi _0^i`$ and $`r_0^{ij}`$. The range $`r_{\mathrm{rep}}`$ of the repulsive potential between pairs of residues that are not bonded and do not interact via a native contact interaction is set to 4 Å. Length is expressed in units of Å, and energy in units of $`\epsilon `$, the energy parameter of the nonnative repulsive interaction, so that $`\epsilon `$ itself is unity. $`k_l`$, $`\epsilon _\theta `$, $`\epsilon _\varphi ^{(1)}`$, $`\epsilon _\varphi ^{(3)}`$ and $`\epsilon _{\mathrm{native}}`$ are also parameters of the potential energy function. $`k_l`$ is fixed at $`100\epsilon /\text{Å}`$, but the other parameters can be varied; $`\epsilon _\theta =20\epsilon _\varphi ^{(1)}`$ and $`\epsilon _\varphi ^{(3)}=0.5\epsilon _\varphi ^{(1)}`$ are defined in terms of $`\epsilon _\varphi ^{(1)}`$, which can be varied to test the effect of changing the strength of the local interactions, while $`\epsilon _{\mathrm{native}}`$ can also be varied (see below) in order to test the effect of changing the strength of the non-local interactions. The energy of the system is thus controlled by three parameters $`\epsilon `$, $`\epsilon _{\mathrm{native}}`$ and $`\epsilon _\varphi ^{(1)}`$. All interaction parameters are taken to be temperature independent in the present study.

Apart from the variable parameters, this energy function differs from the similar function used in  in two further important ways, as follows. (1) For $`r^{ij}/r_0^{ij}<\sqrt{5/6}`$, we set $`\epsilon _{\mathrm{native}}=\epsilon `$, while for $`r^{ij}/r_0^{ij}\sqrt{5/6}`$, we set $`\epsilon _{\mathrm{native}}=\epsilon _a`$. Then we can vary the native interaction parameter $`\epsilon _a`$, in order to test the effect of changing the strength of the non-bonded attractive interactions between residues, while the short-range repulsive part of $`V_{\mathrm{native}}`$ maintains the self-avoidance of the chain. (2) The value $`q`$, which is the smallest number of places along the chain by which two residues can be separated if they are to interact by $`V_{\mathrm{nonnative}}`$, can be set either to $`q=4`$ (in order to eliminate any double counting of local interactions, in situations where $`\epsilon _\varphi ^{(1)}`$ is not being varied) or to $`q=2`$ (in order to allow $`\epsilon _\varphi ^{(1)}`$ to decrease without compromising the self-avoidance of the chain).

The equation of motion of each residue is

$$m\frac{v^i(t)}{t}=F_{\mathrm{conf}}^i(t)m\gamma v^i(t)+\eta ^i(t),$$
(7)
where $`m`$ is the mass of a residue (set to unity), $`\gamma `$ is the coefficient of friction, $`t`$ is time, and $`v^i(t)`$, $`F_{\mathrm{conf}}^i(t)`$ and $`\eta ^i(t)`$ represent each of the three components of the velocity, conformational force and random force, respectively . The random force is given by

$$\eta ^i(t)=\sqrt{\frac{2m\gamma k_BT}{\delta t}}\xi ^i,$$
(8)
where $`\delta t`$ is the integration time step and $`\xi ^i`$ is a random variable taken from a Gaussian distribution with zero mean and unit variance. The most appropriate time scale can be estimated  by $`\tau =\sqrt{m_0a_0^2/\epsilon _0}`$, where $`m_0`$, $`a_0`$ and $`\epsilon _0`$ are the mass, length and energy scales, respectively. We set $`m_0=m=1`$, $`\epsilon _0=\epsilon =1`$ and $`a_0=4`$ Å(the latter is approximately the length of a virtual bond between two residues and is also the range $`r_{\mathrm{rep}}`$ of the repulsive interaction), and so $`\tau =4`$. We define the integration time step $`\delta t=0.005\tau `$ and the coefficient of friction $`\gamma =0.05\tau ^1`$ in terms of this time scale. The velocity-verlet algorithm  is used to integrate the equations of motion.

### 2.2 Lattice models

The lattice models considered here are 27mers with a maximally compact native (ground-state) conformation. Details of the models have been described elsewhere . We compare three native-centric interaction scenarios which have varying degrees, and different mechanisms, of thermodynamic cooperativity. As an example of a particular native conformation to which these three scenarios can be applied, we choose the one in Ref.  with relative contact order $`0.410`$. In scenario (i), which corresponds to the common Gō model, the native contact interactions are pairwise additive. In scenario (ii), we add an extra favorable energy $`E_{\mathrm{gs}}`$ for the native structure as a whole (as defined by equation 5 in Ref. ). Scenario (iii) introduces, in place of the extra favorable energy, a coupling between the strength of the contact interaction and the local geometry: two residues which are in contact in the native state will interact strongly only when the local geometries of the protein chain around the residues are the same as those in the native state, as described in Ref. . We characterize this mechanism as local-nonlocal coupling or “a cooperative interplay between favorable nonlocal interactions and local conformational preferences” . In this scenario, the strength of the native contact interaction is reduced by an attenuation factor $`a`$ when the local geometry is nonnative. The common (uncoupled) Gō model is equivalent to $`a=1`$ (no attenuation), while $`a=0`$ implies complete coupling; $`a=0`$ is used here. Under scenarios (i) and (iii), the native state has an energy of $`28`$ units, while the extra favorable native energy in scenario (ii) changes the energy of the native state to $`42`$ units. Standard Monte Carlo methods are used for conformational sampling . The permitted chain moves are end flips, corner flips, crankshafts and rigid rotations. Each attempted move is counted as one simulation time step, irrespective of whether the move is accepted by the Metropolis criterion.

## 3 Thermodynamics

All simulations are performed at constant temperature, with no explicit consideration of pressure. This is because the focus of the present study is protein behavior under atmospheric pressure, and the contribution of a $`PV`$ term to protein energetics is small under these conditions . Therefore, for our present purposes, we can consider the Helmholtz and Gibbs free energies to be equivalent.

### 3.1 Calculation of the heat capacity

The specific heat capacity $`C_V(T)`$ of the model protein is given by the standard relation

$$C_V(T)=\frac{1}{k_BT^2}\left[E^2(T)E(T)^2\right],$$
(9)
where $`k_BT`$ is the Boltzmann constant multiplied by the absolute temperature, $`X(T)`$ denotes the Boltzmann average of quantity $`X`$ at temperature $`T`$, and the total energy $`E`$ is the sum $`V+E_K`$ of potential and kinetic energies. We compute the averages by standard histogram sampling techniques .

In lattice studies, the kinetic energy $`E_K`$ is not treated. Therefore, $`E`$ in Eq. (9) has traditionally been taken, in protein modeling, to be the potential energy term $`V`$. This procedure has often been extended to continuum model studies, in Ref.  for example, although $`E_K`$ is accessible and well-defined in off-lattice models. However, $`E^2`$ $``$ $`E^2`$ $``$ $`V^2`$ $``$ $`V^2`$ in general. The two quantities would be equal if $`E_K`$ were a constant, but that would be unphysical. In this study, we have calculated $`C_V(T)`$ using Eq. (9) both with $`E=V+E_K`$ and with the substitution $`EV`$. The results are not identical: an example is given in Fig. 1.

Fig. 1 shows that the difference between heat capacity values obtained using the two methods is small around the transition midpoint $`T_m`$. This is because any energy added to the system during the unfolding transition contributes mostly to the potential rather than to the kinetic energy. The difference is less negligible for the “shoulders” on either side of the heat capacity peak. At very low temperatures, including the kinetic energy contribution can lead to a smaller heat capacity, because the molecule at this temperature is in a relatively fixed state: nearly all of the kinetic energy is accounted for by the oscillation of pairs of residues about the minima of their mutual (bonded or non-bonded) interaction energies. The potential energy and the kinetic energy associated with these oscillations both fluctuate, but their sum fluctuates much less, and so the fluctuations in total energy are smaller than the fluctuations in potential energy, with the result that the calculated heat capacity is smaller when $`E_K`$ is taken into account. Overall, Fig. 1 indicates that while the difference between the heat capacities calculated using the two different methods is not negligible, it is not drastic. Probably this is because $`E_K`$, while not invariant, fluctuates much less than $`V`$. For this reason, we do not expect conclusions drawn from previous calculations of heat capacities , which used $`V`$, to be changed greatly by calculations using $`E`$. Nonetheless, we do expect a proper account of the kinetic contributions to protein heat capacities to be important in addressing the contribution of bond vector motions to the heat capacity .

Since the $`PV`$ term is neglected in the present formulation, $`C_V(T)`$ is effectively equal to $`C_P(T)`$, which is generally measured by calorimetry (and which can be expressed in a form similar to Eq. (9), but with the enthalpy $`H`$ taking the place of the energy $`E`$ ). Therefore, we may refer to the quantity computed using Eq. (9) simply as heat capacity. All subsequent heat capacity curves shown in this article for the continuum models are obtained using the total energy $`E=V+E_K`$.

### 3.2 The free energy

The Helmholtz free energy of the model system is $`F(T)=k_BT\mathrm{ln}Z(T)`$, where $`Z(T)`$ is the partition function at temperature $`T`$. It follows that, in the vicinity of the simulation temperature $`T_{\mathrm{sim}}`$, the Helmholtz free energy of the model system at temperature $`T`$, relative to its value at the simulation temperature, may be approximated using the formula:

$$\frac{\mathrm{\Delta }F(T)}{k_BT}=\frac{F(T)}{k_BT}\frac{F(T_{\mathrm{sim}})}{k_BT_{\mathrm{sim}}}=\mathrm{ln}\left\{\underset{i}{}p(E_i;T_{\mathrm{sim}})\mathrm{exp}\left(E_i\left[\frac{1}{k_BT_{\mathrm{sim}}}\frac{1}{k_BT}\right]\right)\right\},$$
(10)
where the sum is performed over sets of microstates in different energy ranges $`E_i`$. $`p(E_i;T_{\mathrm{sim}})`$ is the probability density at the simulation temperature, and is estimated directly from the Langevin dynamics simulations.

The inset of Fig. 1 provides an example of $`\mathrm{\Delta }F(T)`$, showing that the gradient of the free energy with respect to $`T`$ changes rather abruptly around the transition temperature $`T_m`$ (vertical dotted lines). The transition temperature $`T_m`$ corresponds to the temperature at the peak of the heat capacity curve, which was denoted by $`T_{\mathrm{max}}`$ in Ref. . Apparently, below $`T_m`$, the system spends most time in states in the vicinity of the bottom of the native basin, and so the changes in $`F`$ with respect to temperature are dominated by the behavior of these states. However, as the temperature increases past $`T_m`$, the system, and therefore the rate of change of $`F(T)`$, starts to be dominated by states near the bottom of the denatured basin. As a result, the gradient of $`F(T)`$ changes rather suddenly at $`T_m`$. While the $`F(T)`$ gradient could never be discontinuous because the model system is finite, the kink at $`T_m`$ does indicate that the transition is two-state-like, and therefore that it is similar to a first order phase transition.

### 3.3 Thermodynamic cooperativity

The presence of a peak in the heat capacity at a transition temperature $`T_m`$, as in Fig. 1, indicates that the folding/unfolding transition possesses a degree of thermodynamic cooperativity. As our group has argued, the degree of thermodynamic cooperativity in protein models can be quantified by the ratio $`\kappa _2=\mathrm{\Delta }H_{\mathrm{vH}}/\mathrm{\Delta }H_{\mathrm{cal}}`$ of the van’t Hoff enthalpy $`\mathrm{\Delta }H_{\mathrm{vH}}`$ to the calorimetric enthalpy $`\mathrm{\Delta }H_{\mathrm{cal}}`$ of the transition. This ratio is closely related to that determined experimentally by differential scanning calorimetry . In model studies, the calorimetric enthalpy $`\mathrm{\Delta }H_{\mathrm{cal}}`$ may be determined from an integral of the heat capacity across the transition region,

$$\mathrm{\Delta }H_{\mathrm{cal}}=𝑑TC_P(T),$$
(11)
while the van’t Hoff enthalpy is equal to twice the maximum standard deviation of the enthalpy distribution at the transition midpoint,

$$\mathrm{\Delta }H_{\mathrm{vH}}=2\sqrt{k_BT_m^2C_P^{\mathrm{max}}},$$
(12)
where $`C_P^{\mathrm{max}}`$ is the peak value of the heat capacity. We have followed standard usage in this section by expressing $`\kappa _2`$ in terms of $`H`$ and $`C_P`$. However, as mentioned in the previous section, simulations produce values for $`E`$ and $`C_V`$, which for the present application are essentially equivalent to $`H`$ and $`C_P`$.

As has been pointed out , comparison of simulation heat capacity scans to experiment is often complicated by the fact that the heat capacity tails which we observe in simulations would, if they occurred in a real system, be swamped by the solvent contribution and ignored by the common procedure of using empirical baseline subtraction to calculate $`\mathrm{\Delta }H_{\mathrm{vH}}/\mathrm{\Delta }H_{\mathrm{cal}}`$. In other words, tail contributions that arise from conformational transitions may be masked by solvent contributions in real data analysis . Therefore, for completeness, we also perform empirical baseline subtractions on our simulated heat capacity scans, producing a revised ratio $`\kappa _2^{(\mathrm{s})}`$ (defined in Ref. ) to facilitate comparison with experiment.

### 3.4 Radius of gyration

The radius of gyration $`R_g`$ of a particular conformation of the protein is an indicator of its compactness. It is defined by

$$R_g^2=\frac{1}{N}\underset{i=1}{\overset{N}{}}\left|r_ir\right|^2,$$
(13)
where $`N`$ is the number of residues, $`r_i`$ is the position of the $`i`$th residue, and $`r`$ is the average position (centroid) of the given conformation. The Boltzmann average $`R_g=\sqrt{R_g^2}`$ over a given conformational ensemble is obtained by standard histogram techniques. Two-state-like behavior requires a steplike sigmoidal change in $`R_g`$ upon denaturation at $`T_m`$, with little postdenaturational expansion of the chain .

## 4 Results and discussion

To study the effect of local vs. nonlocal interactions, we first vary the strength $`\epsilon _\varphi ^{(1)}`$ of the local interactions while keeping the strength of the nonlocal interactions fixed in the continuum CI2 construct (Figs. 2, 3). The heat capacity scans, for four scenarios (four models) with different values for $`\epsilon _\varphi ^{(1)}`$, are shown in Fig. 2. They all exhibit a fairly sharp peak except for the model with $`\epsilon _\varphi ^{(1)}=0.25`$. The heat capacity peak signifies substantial heat absorption within a narrow temperature range at the folding/unfolding transition. The absorbed energy propels the chain from its low-energy folded conformations (native ensemble) to its high-energy unfolded conformations (denatured ensemble). However, the mere existence of a relatively sharp peak in the heat capacity function does not necessarily mean that the transition is as cooperative as those observed in small single-domain proteins. Coil-globule transitions in homopolymers are not two-state-like, but their calorimetric heat capacity scans can have very sharp peaks . A more quantitative measure of thermodynamic cooperativity is the traditional calorimetric two-state criterion (see above), which has emerged recently as a powerful modeling tool . Ratios of van’t Hoff to calorimetric enthalpy were calculated for the four models as described above; the ranges of the temperature integrations used in the determination of $`\mathrm{\Delta }H_{\mathrm{cal}}`$ \[Eq. (11)\] were taken to be equal to the ranges shown in Fig. 2.

The inset of Fig. 2 shows that the $`\mathrm{\Delta }H_{\mathrm{vH}}/\mathrm{\Delta }H_{\mathrm{cal}}`$ ratio (diamonds) of these models is only weakly dependent on $`\epsilon _\varphi ^{(1)}`$ over an extended range of $`\epsilon _\varphi ^{(1)}`$ values, but that the ratio is significantly smaller when the local interactions are substantially weaker, at $`\epsilon _\varphi ^{(1)}=0.25`$, than the nonlocal interactions. As discussed above, quantitative comparisons between simulated and experimental $`\mathrm{\Delta }H_{\mathrm{vH}}/\mathrm{\Delta }H_{\mathrm{cal}}`$ values require the introduction of model calorimetric baselines  similar to those employed in the interpretation of experimental data. Traditionally, experimental baselines are designed to remove solvation contributions (temperature-dependent effective interactions), in order to extract the heat capacity effects associated with the folding/unfolding transition itself . The present models do not contain temperature-dependent interactions. Therefore, the heat capacity contributions eliminated by the model calorimetric baselines in Fig. 2 can only originate from vibrational motions and conformational transitions. Increasingly, it is being recognized  that similar heat capacity contributions from bond vector motions and more collective conformational transitions might also be “hidden” below traditional baselines constructed for analyzing experimental calorimetric data, although the magnitude of such contributions needs to be elucidated. The three models in Fig. 2 that are relatively more cooperative (with higher $`\kappa _2`$ values) all have modified $`\mathrm{\Delta }H_{\mathrm{vH}}/\mathrm{\Delta }H_{\mathrm{cal}}`$ values ($`\kappa _2^{(\mathrm{s})}`$, circles in the inset), after empirical baseline subtractions, that are very close to unity. (We note that the recent determination of $`\mathrm{\Delta }H_{\mathrm{vH}}/\mathrm{\Delta }H_{\mathrm{cal}}`$ values in an all-atom Gō model  involved baseline subtractions as well: c.f. Fig. 8 of Ref. .)

However, a protein chain model’s ability to attain a near-unity $`\mathrm{\Delta }H_{\mathrm{vH}}/\mathrm{\Delta }H_{\mathrm{cal}}`$ ratio after baseline subtractions does not by itself imply that its thermodynamic behavior is similar to that of real, small single-domain proteins . This is because the heat capacity contributions discarded by certain baselines can actually be symptoms of significant deviations from two-state-like behavior. It has been recognized  that, to clarify this situation, we can use the behavior of the average radius of gyration $`R_g`$ of a protein chain model as an additional evaluation criterion for the model’s thermodynamic cooperativity. Small angle X-ray scattering (SAXS) experiments have demonstrated that the average radius of gyration $`R_g`$ of several small single-domain proteins behaves in an apparently two-state manner , showing very little postdenaturational ($`T>T_m`$) expansion of the chain outside the transition regime that corresponds to the region of the heat capacity peak. We require chain models of small single-domain proteins to exhibit similar behavior . Now, to further assess the four models in Fig. 2 with different local interaction strengths, we calculate their average radii of gyration as a function of temperature (Fig. 3). To ensure adequate sampling, $`R_g`$ for each model is obtained from three different simulation temperatures; the results are thus displayed as three discontinuous curves. Despite some minor discrepancies (owing to sampling uncertainties) between parts of the $`R_g`$ function deduced from different simulation temperatures for the $`\epsilon _\varphi ^{(1)}=0.25`$ case, the general trend in Fig. 3 is very clear. Models with weaker local interactions are less cooperative in that their $`R_g`$ curves show more postdenaturational increase than do those of models having stronger local interactions. For instance, the $`R_g`$ of the $`\epsilon _\varphi ^{(1)}=0.25`$ model increases by $`3.0`$ Å between $`T0.82`$ (the end of the transition region) and $`T1.11`$. In contrast, a similar temperature increase for the $`\epsilon _\varphi ^{(1)}=1.00`$ model from $`T1.11`$ (the end of the transition region) to $`T1.40`$ leads to an increase of only $`1.6`$Å in $`R_g`$. These observations indicate that two-state-like thermodynamic cooperativity cannot be achieved if the local conformational propensities of a protein are much weaker than the favorable nonlocal interactions. This confirms a similar conclusion which was derived recently from a more limited study of a “contact dominant model” .

We next extend our analysis by applying the same computational procedure to varying the strength $`\epsilon _a`$ of the favorable nonlocal interactions while keeping the strength of the local interactions fixed. Consistent with the seminal study of Gō and Taketomi , Figs. 4 and 5 show that variations in nonlocal $`\epsilon _a`$ have a more prominent effect on thermodynamic cooperativity than variations in local $`\epsilon _\varphi ^{(1)}`$. While the peak heat capacity values for the three models in Fig. 2 with $`\epsilon _\varphi ^{(1)}0.5`$ are similar, the peak heat capacity values for the three models in Fig. 4 with $`\epsilon _a0.5`$ show a significant monotonic increase with $`\epsilon _a`$. In addition, for the $`\epsilon _a=0.5`$ model in Fig. 4, the difference between unity and the $`\mathrm{\Delta }H_{\mathrm{vH}}/\mathrm{\Delta }H_{\mathrm{cal}}`$ ratio after baseline subtraction is not negligible ($`\kappa _2^{(\mathrm{s})}=0.91`$). Despite these differences, the trends in Figs. 4, 5 are in large measure similar to those in Figs. 2, 3. In particular, Fig. 4 shows that the model with $`\epsilon _a=0.25`$, like the $`\epsilon _\varphi ^{(1)}=0.25`$ case in Fig. 2, has a significantly lower $`\mathrm{\Delta }H_{\mathrm{vH}}/\mathrm{\Delta }H_{\mathrm{cal}}`$ ratio than the other three models considered in the same figure. The $`R_g`$ data in Fig. 5 shows that thermodynamic cooperativity increases with $`\epsilon _a`$, as manifested in a smaller amount of postdenaturational conformational expansion with increasing $`\epsilon _a`$; this is comparable to the effect of increasing $`\epsilon _\varphi ^{(1)}`$ in Fig. 3. Taken together, the results in Figs. 2–5 suggest that a high degree of thermodynamic cooperativity, similar to that in real, small single-domain proteins, requires both strong local and strong nonlocal interactions. Apparently, a high degree of thermodynamic cooperativity is incompatible with either a much weakened local conformational preference relative to the favorable nonlocal interactions ($`\epsilon _\varphi ^{(1)}\epsilon _a`$) or much weakened favorable nonlocal interactions relative to the local conformational preference ($`\epsilon _a\epsilon _\varphi ^{(1)}`$).

Although three-dimensional Gō-like models with strong local and nonlocal interactions appear to satisfy the thermodynamic criterion of calorimetric two-state cooperativity, it has recently been noted that they are unable to produce simple two-state folding/unfolding kinetics . This is because the thermodynamic cooperativity of these models is not sufficiently high. As a result, and in spite of the native-centric nature of the common pairwise additive Gō-like interactions, kinetic trapping becomes significant under strongly native conditions, leading to folding rate slow-downs and chevron rollovers . More recent lattice model investigations indicate that simple two-state folding/unfolding kinetics require a high degree of thermodynamic cooperativity that may be characterized as “near-Levinthal” , necessitating many-body interactions beyond those postulated by the common Gō model .

In view of this recent development, and to facilitate the construction and investigation of continuum models that incorporate these new ideas, it is instructive to compare in more detail the thermodynamic behavior of the common lattice Gō construct (with only pairwise additive contact energies) with the behavior of models which have many-body interactions and enhanced cooperativity. We also wish to investigate whether results obtained from lattice models supply additional support for the conclusions which we have derived from our continuum model results. To this end, Figs. 6–8 compare three 27mer lattice models.

Because of their intrinsic restrictions on conformational possibilities, it is more straightforward to construct cooperative lattice models than to construct off-lattice continuum models that are similarly cooperative. Recently, using evidence from kinetic simulations of chevron plots, our group has proposed that a $`\mathrm{\Delta }H_{\mathrm{vH}}/\mathrm{\Delta }H_{\mathrm{cal}}`$ ratio of $`\kappa _2>0.9`$ before baseline subtractions \[as for models (ii) and (iii) in Fig. 6a\] is likely to be required in order for a lattice protein chain model to produce chevron plots with linear regimes similar in extent to those observed for real, small single-domain proteins . However, this numerical criterion is not readily generalizable to off-lattice continuum models. This is because the heat capacity effects of bond vibrations and kinetic energy have to be taken into account in continuum models, whereas these effects are absent in lattice models. Thus, in the characterization of a model’s thermodynamic cooperativity, more detailed information concerning, for example, the behavior of the average radius of gyration, has to be relied upon more heavily for continuous models (see above) than for lattice models.

For the three lattice models studied here, the $`R_g`$ plots in Fig. 6b show that, while the postdenaturational conformational expansion of the common lattice Gō construct (solid curve in Fig. 6b) is considerably milder than that of its continuum counterpart (solid curves in Figs. 3, 5), the more cooperative lattice models with many-body interactions exhibit much less (dotted curve in Fig. 6b) or nearly non-existent (dashed curve in Fig. 6b) postdenaturational conformational expansion. The fluctuations in $`E`$ and $`R_g`$ near the transition midpoint, shown in Figs. 7 and 8, indicate further that the transitions between the native and denatured ensembles are sharper and more two-state-like for the more cooperative models (ii) and (iii) with many-body interactions \[parts (b) and (c) of Figs. 7, 8\] than for the common Gō model \[parts (a) of Figs. 7, 8\]. The corresponding fluctuations in the fractional number of native contacts $`Q`$  (data not shown) were also found to exhibit a trend very similar to that of the energy fluctuations in Fig. 7.

The lattice model results, shown in Figs. 6–8, are compatible with the conclusion, reached above on the basis of continuum model results, that both local and nonlocal interactions are important for thermodynamic cooperativity. The common lattice Gō model of scenario (i) includes only nonlocal interactions, analogous to the interactions encoded by $`V_{\mathrm{native}}`$ in the off-lattice model. Scenario (iii), which takes account also of the local geometry of the chain, displays greater cooperativity than scenario (i).

Higher resolution data such as that in Figs. 7 and 8 opens up future avenues for the assessment of different mechanisms of cooperativity using comparisons between model predictions and experimental measurements of, for example, conformational sizes and fluctuations . It is noteworthy that in the model (ii) scenario, with an extra favorable energy for the native structure as a whole, the native ensemble does not exhibit much energetic or conformational variation (horizontal line segments at low $`E`$ and low $`R_g`$ values in Figs. 7b and 8b). On the other hand, in the model (iii) scenario with local-nonlocal coupling, there is considerable variation in the native ensemble (c.f. low $`E`$ and low $`R_g`$ fluctuations in Figs. 7c and 8c). Yet the variation in the denatured ensemble is smaller in model (iii) than in model (ii) (c.f. high $`E`$ and high $`R_g`$ fluctuations in parts (b) and (c) of Figs. 7, 8), resulting in a more two-state-like average $`R_g`$ transition for model (iii) than for model (ii), manifested in a near-immediate postdenaturational ($`T>T_m`$) saturation of the dashed curve for model (iii) in Fig. 6b compared to a more gradual postdenaturational saturation of the dotted curve for model (ii) in the same figure. All these differences in conformational properties are in principle detectable through experiments on real proteins. Hence, future experimental efforts along the lines suggested here would help to verify or falsify different proposed scenarios and interaction mechanisms  in the endeavor to decipher the physical origins of cooperativity in real proteins.

## 5 Conclusions

The present study suggests strongly that both local conformational preferences and favorable nonlocal interactions are crucial to protein thermodynamic cooperativity. This result points to a useful constraint on simplified models of protein molecules: they should take account both of local and of nonlocal interactions.

As emphasized above, the scope of the present study is limited. Only native-centric interaction schemes are considered; and here we have elected only one particular physically plausible way to classify energy contributions into “local” and “nonlocal” terms in the continuum models. In addition, by using a native-centric interaction scheme both for local and for nonlocal interactions, we have avoided the possibility of energetic frustration, which might be significant if a more realistic interaction scheme were used . To further elucidate the answers to the questions we have posed, much remains to be investigated.

Nonetheless, our results show clearly that a high degree of thermodynamic cooperativity is compatible neither with a much weakened local interaction nor with a much weakened nonlocal interaction, indicating that both local and nonlocal interactions are important components in protein energetics . This finding is consistent with the notion that a cooperative interplay between local and nonlocal interactions  is a critical ingredient underlying the apparent simple two-state cooperativity of real, small single-domain proteins.

With regard to Gō-like native-centric modeling (see, e.g., discussion in Refs. ), we observe that significant differences in model predictions can result from different Gō-like interaction schemes, even though all of the schemes are designed to bias the chain towards the same native structure. This underscores our point that requiring a consistent account of cooperativity can be a more productive approach to protein modeling than simply designing a model heteropolymer to fold to a target structure . In this context, the present coarse-grained representations constitute only a first step in the understanding of protein cooperativity. Ultimately, atomistic origins of local and nonlocal interactions such as sidechain packing  must be taken into account in an effort to provide the necessary physical underpinning for the cooperative mechanisms proposed here.

## Acknowledgments

The research reported here was partially supported by the Canadian Institutes of Health Research (CIHR grant no. MOP-15323), PENCE, a Premier’s Research Excellence Award from the Province of Ontario, and the Ontario Centre for Genomic Computing at the Hospital for Sick Children in Toronto. H. S. C. holds a Canada Research Chair in Biochemistry.

Figure Captions

FIGURE 1

Heat capacity as a function of temperature, for one particular set of interaction parameters. Solid curve: heat capacity calculated using the total energy; dashed curve: heat capacity calculated using only the potential energy. Parameters are $`\epsilon _\varphi ^{(1)}=1.00`$, $`\epsilon _a=1.00`$, $`q=4`$; simulation temperature $`T_{\mathrm{sim}}=1.02`$. Inset: free energy as a function of temperature, for the same model, showing a sharp change in gradient around $`T_m`$. The vertical dotted lines in the figure and the inset mark the transition midpoint temperature $`T_m`$.

FIGURE 2

Heat capacity as a function of temperature, for varying local interaction energy $`\epsilon _\varphi ^{(1)}`$. Other parameters $`\epsilon _a=1.00`$ and $`q=2`$ are fixed. (From left to right) dotted curve: $`\epsilon _\varphi ^{(1)}=0.25`$, $`T_m=0.74`$; short dashed curve: $`\epsilon _\varphi ^{(1)}=0.50`$, $`T_m=0.84`$; long dashed curve: $`\epsilon _\varphi ^{(1)}=0.75`$, $`T_m=0.94`$; solid curve: $`\epsilon _\varphi ^{(1)}=1.00`$, $`T_m=1.03`$. These scans are obtained by histogram techniques from simulations performed at $`T_{\mathrm{sim}}=`$ $`0.73`$, $`0.84`$, $`0.94`$, and $`1.03`$ respectively. The $`\mathrm{\Delta }H_{\mathrm{vH}}/\mathrm{\Delta }H_{\mathrm{cal}}`$ cooperativity coefficients $`\kappa _2`$ without baseline subtractions are $`0.33`$, $`0.43`$, $`0.44`$, and $`0.44`$ respectively. Modified cooperativity coefficients $`\kappa _2^{(\mathrm{s})}`$ after subtraction of the baselines (indicated by thin lines in the figure) for $`\epsilon _\varphi ^{(1)}=`$ $`0.50`$, $`0.75`$, and $`1.00`$ are $`0.97`$, $`0.98`$, and $`0.99`$ respectively. No value for $`\kappa _2^{(\mathrm{s})}`$ was calculated for $`\epsilon _\varphi ^{(1)}=0.25`$ because the shape of its heat capacity curve does not suggest any clear choice of baselines that are intuitively more reasonable than others. The inset shows $`\kappa _2`$ (diamonds) and $`\kappa _2^{(\mathrm{s})}`$ (circles) as functions of $`\epsilon _\varphi ^{(1)}`$.

FIGURE 3

Average radius of gyration as a function of temperature, for varying local interaction energy $`\epsilon _\varphi ^{(1)}`$; other parameters $`\epsilon _a=1.00`$ and $`q=2`$ are fixed, as in Fig. 2. The correspondence between line styles and $`\epsilon _\varphi ^{(1)}`$ values is identical to that in Fig. 2. For each value of $`\epsilon _\varphi ^{(1)}`$, simulations were performed at three different values of $`T_{\mathrm{sim}}`$ to ensure adequate sampling across the entire temperature range shown. (From left to right) for $`\epsilon _\varphi ^{(1)}=0.25`$ (dotted curves), $`T_{\mathrm{sim}}=0.73`$, $`0.93`$, $`1.13`$; for $`\epsilon _\varphi ^{(1)}=0.50`$ (short dashed curves), $`T_{\mathrm{sim}}=0.84,`$ $`1.04`$, $`1.24`$; for $`\epsilon _\varphi ^{(1)}=0.75`$ (long dashed curves), $`T_{\mathrm{sim}}=0.94,`$ $`1.14`$, $`1.34`$; and for $`\epsilon _\varphi ^{(1)}=1.00`$ (solid curves), $`T_{\mathrm{sim}}=1.03`$, $`1.23`$, and $`1.43`$.

FIGURE 4

Heat capacity as a function of temperature, for varying nonlocal interaction energy $`\epsilon _a`$. Other parameters $`\epsilon _\varphi ^{(1)}=1.00`$ and $`q=4`$ are fixed. (From left to right) dotted curve: $`\epsilon _a=0.25`$, $`T_m=0.43`$; short dashed curve: $`\epsilon _a=0.50`$, $`T_m=0.65`$; long dashed curve: $`\epsilon _a=0.75`$, $`T_m=0.84`$; solid curve: $`\epsilon _a=1.00`$, $`T_m=1.02`$. These scans are obtained by histogram techniques from simulations performed at $`T_{\mathrm{sim}}=`$ $`0.42`$, $`0.64`$, $`0.84`$, and $`1.02`$ respectively. The $`\mathrm{\Delta }H_{\mathrm{vH}}/\mathrm{\Delta }H_{\mathrm{cal}}`$ cooperativity coefficients $`\kappa _2`$ without baseline subtractions are $`0.28`$, $`0.40`$, $`0.44`$, and $`0.46`$ respectively. Modified cooperativity coefficients $`\kappa _2^{(\mathrm{s})}`$ after subtraction of the baselines (indicated by thin lines in the figure) for $`\epsilon _a=`$ $`0.50`$, $`0.75`$, and $`1.00`$ are $`0.91`$, $`1.00`$, and $`0.99`$ respectively. No value for $`\kappa _2^{(\mathrm{s})}`$ was calculated for $`\epsilon _a=0.25`$ for the same reason that no $`\kappa _2^{(\mathrm{s})}`$ was provided for $`\epsilon _\varphi ^{(1)}=0.25`$ in Fig. 2. The inset shows $`\kappa _2`$ (diamonds) and $`\kappa _2^{(\mathrm{s})}`$ (circles) as functions of $`\epsilon _a`$.

FIGURE 5

Average radius of gyration as a function of temperature, for varying nonlocal interaction energy $`\epsilon _a`$. Other parameters $`\epsilon _\varphi ^{(1)}=1.00`$ and $`q=4`$ are fixed, as in Fig. 4. The correspondence between line styles and $`\epsilon _a`$ values is identical to that in Fig. 4. For each value of $`\epsilon _a`$, simulations were performed at three different values of $`T_{\mathrm{sim}}`$ to ensure adequate sampling across the whole range shown. (From left to right) for $`\epsilon _a=0.25`$ (dotted curves), $`T_{\mathrm{sim}}=0.42`$, $`0.62`$, $`0.82`$; for $`\epsilon _a=0.50`$ (short dashed curves), $`T_{\mathrm{sim}}=0.64`$, $`0.84`$, $`1.04`$; for $`\epsilon _a=0.75`$ (long dashed curves), $`T_{\mathrm{sim}}=0.84,`$ $`1.04`$, $`1.24`$; and for $`\epsilon _a=1.00`$ (solid curves), $`T_{\mathrm{sim}}=1.02`$, $`1.22`$, and $`1.42`$.

FIGURE 6

Thermodynamic cooperativities of the three representative 27mer lattice models described in the text. Model definitions and simulation details are given in Refs. . Heat capacity (a) and average radius of gyration (b) are determined by histogram techniques based upon Monte Carlo sampling performed at $`T_{\mathrm{sim}}=T_m`$. (a) Heat capacity of (i) the common Gō model with pairwise additive contact energy (solid curve, left); (ii) the model that assigns an extra favorable energy to the native structure as a whole (dotted curve, right); and (iii) the model with local-nonlocal coupling (dashed curve, middle). The transition temperatures for the models are $`T_m=`$ $`0.701`$ (i), $`1.13`$ (ii), and $`0.755`$ (iii). Their $`\mathrm{\Delta }H_{\mathrm{vH}}/\mathrm{\Delta }H_{\mathrm{cal}}`$ without baseline subtractions are $`\kappa _2=`$ $`0.86`$, $`0.98`$, and $`0.99`$, and the corresponding ratios after subtracting the baselines shown are $`\kappa _2^{(\mathrm{s})}=`$ $`1.00`$, $`1.00`$, and $`1.00`$ respectively. (b) Average radius of gyration as a function of model temperature for the three models \[represented by the same line styles as in (a)\].

FIGURE 7

Representative trajectories of the three models in Fig. 6 at their respective transition temperatures. Variations in the potential energy of models (i)–(iii) are shown in (a)–(c) respectively.

FIGURE 8

Same as Fig. 7, except that variations in the radius of gyration are shown here.


text:
# Phase engineering of controlled entangled number states in a single component Bose-Einstein condensate in a double well

## References