Computational Modeling of the Thermodynamics of the Mesophilic and Thermophilic Mutants of Trp-Cage Miniprotein.

We characterize the folding-unfolding thermodynamics of two mutants of the miniprotein Trp-cage by combining extended molecular dynamics simulations and an advanced statistical-mechanical-based approach. From a set of molecular dynamics simulations in an explicit solvent performed along a reference isobar, we evaluated the structural and thermodynamic behaviors of a mesophilic and a thermophilic mutant of the Trp-cage and their temperature dependence. In the case of the thermophilic mutant, computational data confirm that our theoretical-computational approach is able to reproduce the available experimental estimate with rather good accuracy. On the other hand, the mesophilic mutant does not show a clear two-state (folded and unfolded) behavior, preventing us from reconstructing its thermodynamics; thus, an analysis of its structural behavior along a reference isobar is presented. Our results show that an extended sampling of these kinds of systems coupled to an advanced statistical-mechanical-based treatment of the data can provide an accurate description of the folding-unfolding thermodynamics along a reference isobar, rationalizing the discrepancies between the simulated and experimental systems.

Introduction

The quantitative characterization of a folding–unfolding process via theoretical–computational approaches represents one of the main challenges in the field. In fact, although there has been a continuous increase in computational power, the possibility of observing several folding–unfolding events during a molecular dynamics simulation is still limited by the time scale of these events, which usually ranges from milliseconds to seconds at physiological temperature. To this end, proteins formed by a few residues, i.e. miniproteins, have been synthesized with the objective of accelerating the folding kinetics, making the folding process faster: i.e., on the scale of microseconds. Among others, the Trp-cage miniprotein is one of the most studied and characterized by both experimental and computational approaches.[1−16] NMR experiments show that this miniprotein is characterized by a stable fold in solution,[13] formed by helical structures in the N-terminal region and a polyproline II helix at the C-terminus. The feature of a (sort of) tertiary structure is due to the Trp-cage fold determined by the interactions among Tyr3, Trp6, Gly11, Pro12, Pro18, and Pro19.

Due to the very limited size of Trp-cage, it was also possible to carefully investigate the effect of Trp-cage point mutations on its folding–unfolding thermodynamics. Very interestingly, these experiments have shown that single mutations are able to significantly alter the melting point with respect to the wild-type construct.[15,17] For example, alanine insertions within the N-terminal α-helix have been shown to increase the Trp-cage helicity as well to stabilize the folded state.[17] Replacement of Leu, Ile, Lys, or Ser residues by Ala resulted in a folding free energy decrease of 1.5 kJ/mol for each substitution, without any structural changes.[15] In addition, selected mutations on the N-terminal or C-terminal regions have been shown to significantly change the melting temperature, providing mesophilic or thermophilic constructs.[11] Therefore, we present here a theoretical–computational characterization of two mutants that experimentally show an increase (TC10b, Tm = 328–329 K) and a decrease (S14A, Tm = 294 K) in their melting temperature with respect to the wild-type construct (TC5b, Tm = 316 K). By means of an advanced statistical–mechanical based approach, we are able to characterize the structural–dynamic behavior of these two peptides (represented in Figure ) in water and, when a clear two-state folding–unfolding behavior was observed by means of extended molecular dynamics simulations at different temperatures, we were able to reconstruct its thermodynamics along a reference isobar. Such an approach allowed us to rationalize the differences between the experimental data and the computational estimates of the thermodynamic properties.

Figure 1

Cartoon representations and the correspondence sequences of the wild-type (left), TC10b (middle), and S14A (right) constructs.

Theory

In this section, the essential equations derived by the use of quasi-Gaussian entropy theory are reported.[18] By such an approach, it is possible to express the free energy of a complex molecular system by its potential energy distribution function. Using a rather general and accurate model to describe the energy fluctuation distribution as provided by the Gamma distribution, the thermodynamic properties can be then expressed as a function of the temperature. In an isothermal–isobaric ensemble, the Gibbs free energy (G) of a system composed of one solute molecule embedded in ns solvent molecules iswhere Δ is the partition functionIn eq , 1/β = kT (k is the Boltzmann constant), Q is the canonical partition function, p is the equilibrium pressure, and the summation is over all the possible instant volumes V of the system (the volume fluctuations are considered as discrete variations with the difference between two consecutive volumes virtually corresponding to a differential). Using the derivations recently published for a similar system,[18] where the energy–enthalpy fluctuation distribution was modeled by means of a γ distribution, the chemical potential of the solute (at a given reference low concentration) in the native (N) or denatured (D) conformational state can be expressed by the equationswhere T0 is the reference temperature used with the zero subscript indicating that the property is obtained at T0, the prime indicates reduced properties as obtained with no quantum vibrational contribution, Qv is the solute quantum vibrational partition function ,and as usual μ, H, S, and Cp are the chemical potential, the (molecular) enthalpy, entropy, and heat capacity, respectively. With the folding–unfolding equilibrium temperature Tm chosen as the reference temperature, the (standard) unfolding partial molecular free energy, enthalpy, entropy and heat capacity along the isobar with Tm as reference temperature can be expressed asMore details on the mathematical treatment used to derive the above expressions can be found in our previous work.[18]

Methods

Thermophilic and mesophilic mutants of the Trp-cage miniprotein were studied by means of molecular dynamics simulations. The structure of the thermophilic Trp-cage (named TC10b), characterized by a melting temperature of Tm = 329 K, was taken from the Protein Data Bank for TC10b (pdb code: 2JOF), whereas the structure of the mesophilic construct (named S14A) has been modeled by manually mutating the thermophilic structure by means of the Pymol software package.[19] Microseconds-long MD simulations at different temperatures were performed for both the systems in the range 320–390 K (see Table for details). The set of temperatures was chosen according to our previous work[18] as well as by a few short preliminary MD simulations.

Table 1

Temperatures and Lengths of the MD Simulations of the TC10b and S14A Systems Using the Amber99sb-ildn Force Field

system	T (K)	t (μs)
TC10b	325 (folded)	1.1
	325 (unfolded)	3.0
	340	26.3
	364	10.0
	390	8.8
S14A	320	26.3
	340	10.0
	360	8.8

All of the MD simulations have been performed using the Gromacs 2019 software package.[20] The structures have been solvated with the SPC water model[21] in a triclinic box large enough to maintain at least a 1.3 nm distance between the protein and the box faces. We used the Amber99sb-ildn force field[22] and the leapfrog algorithm with 2 fs time-step integration. An additional subset of MD simulations has been repeated using the CHARMM 36m force field,[23] which is particularly suited to investigate unstructured systems.[24] The V-rescale algorithm has been used to keep the temperature fixed.[25] All of the bonds were constrained using the LINCS algorithm.[26] We used the periodic boundary conditions and the PME to compute long-range interactions;[27] the cutoff radius for short-range interactions was set to 1.1 nm with the Verlet cutoff scheme.[28] After minimization and thermalization, the volume of the box at each temperature has been tuned to reproduce the pressure of 560 bar, the value at which the SPC model has a density corresponding to the experimental water density.[29]

Results

The combination of the theoretical model briefly described in Theory (see our previous work[18] for more details) and a set of MD simulations performed at different temperatures (see Table ) allows us to model the folding–unfolding thermodynamics of these kinds of systems along a reference isobar. The RMSD of the C-α carbon has been used to discriminate between folded and unfolded states. Such a metric has been successfully used in our previous work[18] on the wild type (wt) of the Trp-cage miniprotein. Note that, to reconstruct the whole thermodynamics along the isobar using the eqs –8, the equilibrium populations of the folded and the unfolded states have to be known (at least) at one temperature value. Due to the increase of the folding and unfolding kinetic constants with an increase in temperature, a better estimate of the equilibrium populations was possible at high temperature in the case of TC10b. On the other hand, the slow kinetics of the S14A construct as well as the appearance of metastable states prevented a reliable estimate of the equilibrium populations even at high temperatures, and thus, only a structural characterization of this construct was possible.

Trp-Cage Thermostable Mutant

The TC10b thermostable mutant of the Trp-cage miniprotein, characterized by a melting temperature higher (Tm = 329 K) than that of the wild type (Tm = 316 K),[15] was investigated by a set of extended MD simulations at different temperatures (see Table ). As explained above, the conformational behavior has been characterized by the C-α RMSD with respect to the experimental structure. The RMSD along the MD trajectories at 340, 364, and 390 K show the presence of two well-defined conformational states: one characterized by RMSD values below 0.2 nm and one with values greater than 0.2 nm. Therefore, the distribution minimum at 0.2 nm has been considered as the threshold value to discriminate between the folded (RMSD ≤ 0.2 nm) and unfolded states (RMSD > 0.2 nm). The RMSD and the corresponding probability density for this system as provided by the MD simulation at T = 390 K are shown in Figure . From the RMSDs calculated at different temperatures, it can be seen that our best estimate of Δμunf is at T = 390 K (−7.5 kJ/mol), where several folding–unfolding transitions have been observed (see Figure S1 in Supporting Information for details). Note that a 390 K MD simulation, repeated under the same conditions using the CHARMM 36m force field, provides a very similar value of Δμunf (−8.9 kJ/mol). Therefore, from the Δμunf and the ΔHunf values at 390 K the Δsunf value has been obtained, whereas the ΔCp value has been obtained from the slope of a plot of ΔHunf vs T (Figure a). These values allowed us to reconstruct the entire thermodynamics along the isobar using eqs –8; the result of eq (Δμunf vs T) is shown in Figure b. In the same figure, the Δμunf vs T plot as provided by the experimental parameters reported in Table is shown for comparison. Our results at the melting temperature Tm are shown in Table together with the corresponding experimental estimates. In Table , we also extend the comparison to the wt Trp-cage sequence, which was studied in our previous work.[18] Note that in Table we report two values for each thermodynamic property: one obtained from our best estimate of ΔCp and one obtained using ΔCp = 0. This is due to the small value of this property, which makes an accurate estimate very difficult: i.e., potentially affected by a remarkably large error, as confirmed by the relatively wide range of the experimental values reported in the literature.[4,11,15] Using these two ΔCp values, our computational estimates of the melting temperature are both in agreement with the corresponding experimental values. In fact, our predicted values of Tm between 328 and 332 K are very close to the experimental values (328–329 K). Interestingly, our estimates of both ΔHm and ΔSm are lower than the corresponding experimental measurements. Such an effect, already observed in the wt Trp-cage,[18] is probably due to the force-field inaccuracies.[16] However, the underestimation of Δhm is compensated by the underestimation of ΔSm, resulting in a Tm value very close to the experimental values. From a comparison with the wt Trp-cage, it appears that the higher thermostability of TC10b predicted by our theoretical–computational approach, in agreement with experimental data, is due to the enthalpic gain of the folded conformation of TC10b with respect to the wt miniprotein, overcompensating the entropic effect which favors the unfolded state more in TC10b than in wt. This is partially in agreement with the experimental data, where the stabilizing effect of the mutations is mainly due to an enthalpic effect, the Δsm values between the two constructs being almost identical.[4,15] However, an independent set of experimental values reports that the higher melting temperature of TC10b with respect to the wt is entirely due to entropic effects.[11]

Figure 2

RMSD trajectory (a) and the corresponding normalized distribution (b) of the TC10b mutant as obtained by a MD simulation at T = 390 K. The red line represents the minimum value between the peaks of the folded and unfolded states.

Figure 3

Thermodynamic properties of the TC10b mutant. (a) Unfolding enthalpy as a function of temperature. The line represents the linear fit of the unfolding enthalpy obtained from the MD simulations. The slope of this line corresponds to the unfolding molar heat capacity ΔCp. (b) The Δμ curves as predicted using the thermodynamics model described with MD-derived (our model) and experimental parameters (exp 1, exp 2, and exp 3) reported in Table . The error interval associated with our data is shown as a shaded light blue area.

Table 2

Comparison of the Theoretical–Computational Models for the Mutant and Wild Type Peptides with the Corresponding Experimental Dataa

	Tm (K)	ΔHm (kJ/mol)	ΔSm (J mol–1 K–1)	ΔCp (J mol–1 K–1)
TC10b model	328	28.5	87	314
TC10b model	332	35.4	107	0.0
wild type model	316	18.1	57	350
wild type model	322	30.0	93	0.0
TC10b expt data[11]	329	51.4	156	–566
TC10b expt data[15]	329	58.2	176	–220
TC10b expt data[4]	328	58.2	176	176
wild type expt data[15]	316	56.5	179	97

In the table, the theoretical–computational melting data are obtained by both using the explicit estimate of ΔCp and setting it to zero (in this latter case ΔHm and ΔSm evaluated at 350 K are reported).

Mesostable Trp-Cage Mutant

Similarly to the TC10b construct, the S14A mutant was simulated at three different temperatures (see Table for details). Although the simulation at T = 360 K, shown in Figure , exhibits similar behavior with respect to the thermophilic mutant, the simulations at lower temperatures show an anomalous behavior.

Figure 4

RMSD trajectory (a) and the corresponding normalized distribution (b) of the S14A mutant at T = 360 K. The red line represents the minimum value between the peaks of the folded and unfolded states.

In fact, as shown in Figure , at T = 320 K the RMSD values are characterized by low fluctuations in several intervals along the MD trajectory. These RMSD values are on average higher than those of the folded state and lower than those of the unfolded state (see Figure S2 in the Supporting Information for more details), indicating the presence of “metastable” states differing from both the unfolded and folded ones. We called these states “metastable” because these conformations are stable for several nanoseconds. Such behavior suggests a free energy profile with several local minima, represented by these metastable states, deep enough to trap the system between the folded and unfolded states. This is qualitatively confirmed by a principal component analysis of the S14A MD trajectory. In fact, the trajectory projection on the main conformational subspace, as defined by the first two eigenvectors, shows that each metastable state spans a very limited region of such a subspace (see Figure S3 in the Supporting Information).

Figure 5

RMSDs of the S14A mutant as obtained from the MD simulation at 320 K. Central panel: RMSDs of the S14A mutant with respect to the starting structure. The structures showing low fluctuations of the RMSD (indicated by dots along the RMSD trajectory) are represented in the upper panels (“metastable structures”). Bottom panel: the RMSD matrix as obtained by comparing all the structures sampled during the MD simulations.

An analysis of these conformations characterized by low RMSD fluctuations reveals that the S14A accesses several structurally similar states. These states are well identified by the RMSD matrix (Figure ) by the rectangles near the diagonal of the matrix, which indicate that similar structures are observed consecutively for several nanoseconds. However, such states are sampled at T = 320 K, preventing the applications of our approach on the basis of the knowledge of the thermodynamics properties at different temperatures: i.e., along an isobar in this case. To exclude that such a behavior is due to the specific force field used, we repeated a few MD simulations of the S14A with the CHARMM 36m force field (see Figure S4 in the Supporting Information). The two force fields both show a similar structural–dynamics description of the S14A mutant, confirming the presence of “metastable” states at temperatures between 320 and 340 K. Therefore, the existence of these metastable states kin combination with the lack of an appropriate sampling of the conformational transitions between folding and unfolding states prevent us from reconstructing the thermodynamic properties of such a peptide.

Conclusions

In this work we applied a statistical–mechanics-based approach to an extended set of all-atom molecular dynamics simulations in order to evaluate the thermodynamics of two mutants of the Trp-cage miniprotein. Our results, in line with available experimental data, show that when it is possible to sample several folding–unfolding events by MD simulations at a reference temperature, our theoretical approach is able to accurately model the complete thermodynamics of these systems and its temperature dependence. Our data also confirm that the MD simulations underestimate the ΔHunf value due to the force-field inaccuracies, which are compensated by the underestimation of the ΔSunf value, thus providing a melting temperature very close to the experimental value. Using the same approach, we have also studied the thermodynamics of a mesophilic mutant of the Trp-cage. However, in this case, the lack of well-defined folding–unfolding transitions as well as its slow folding kinetics prevented us from reconstructing the unfolding thermodynamics.

Our results show that it is possible to theoretically describe the thermodynamic behavior of peptides in solution and to understand in detail the discrepancies between the model and the real system on chemical/physical grounds.