Abstract

Here we developed a theoretical approach to predict the effect of mutations on protein stability from non-equilibrium unfolding simulations. We establish a relative measure based on apparent simulated melting temperatures that is independent of simulation length and, under certain assumptions, proportional to equilibrium stability, and we justify this theoretical development with extensive simulations and experimental data. Using our new method based on all-atom Monte-Carlo unfolding simulations, we carried out a saturating mutagenesis of Dihydrofolate Reductase (DHFR), a key target of antibiotics and chemotherapeutic drugs. The method predicted more than 500 stabilizing mutations, several of which were selected for detailed computational and experimental analysis. We find a highly significant correlation of r = 0.65—0.68 between predicted and experimentally determined melting temperatures and unfolding denaturant concentrations for WT DHFR and 42 mutants. The correlation between energy of the native state and experimental denatu ration temperature was much weaker, indicating the important role of entropy in protein stability. The most stabilizing point mutation was D27F, which is located in the active site of the protein, rendering it inactive. However for the rest of mutations outside of the active site we observed a weak yet statistically significant positive correlation between thermal stability and catalytic activity indicating the lack of a stability-activity tradeoff for DHFR. By combining stabilizing mutations predicted by our method, we created a highly stable catalytically active E. coli DHFR mutant with measured denaturation temperature 72°C higher than WT. Prediction results for DHFR and several other proteins indicate that computational approaches based on unfolding simulations are useful as a general technique to discover stabilizing mutations.

Author Summary

However, commonly employed molecular dynamics simulations suffer from a limitation in accessible time scale, making it difficult to model large-scale unfolding events in a realistic amount of simulation time without employing unrealistically high temperatures. Here, we describe a rapid all-atom Monte Carlo simulation approach to simulate unfolding of the essential bacterial enzyme Dihydrofolate Reductase (DHFR) and all possible single point-mutants. We use these simulations to predict which mutants will be more thermodynamically stable (i.e., reside more often in the native folded state vs. the unfolded state) than the wild-type protein, and we confirm our predictions experimentally, creating several highly stable and catalytically active mutants. Thermally stable active engineered proteins can be used as a starting point in directed evolution experiments to evolve new functions on the background of this additional “reservoir of stability.” The stabilized enzyme may be able to accumulate a greater number of destabilizing yet functionally important mutations before unfolding, protease digestion, and aggregation abolish its activity.

Introduction

Most proteins must be folded to carry out their functions in vitro or in vivo. In addition, nonfunctional aggregation of unfolded or par-tially-unfolded proteins can have a deleterious effect on the fitness of an organism and can lead to protein aggregation diseases, which include Alzheimer’s and Huntington’s, in humans [4—6]. Aggregation of poorly folded proteins can also hamper protein production for research and technological purposes [7].

This stabilization can be achieved by either slowing the rate of unfolding or speeding the rate of folding, depending on the role of the mutated residue in the folding nucleation process [13,14]. The unfolding temperature, Tm, at which the free energy of the folded and unfolded states coincide (AG 2 0) serves as a common measure of protein stability. Tm is obtainable by experiment and, in theory, from simulation, although current molecular dynamics simulations are limited in their ability to capture full folding or unfolding trajectories of most proteins (except very small fast folding domains [15]) in a tractable amount of simulation time [16].

However, the performance of these popular methods is still relatively weak [20—22]. Other existing techniques to rationally design proteins with improved stability have involved optimization of charge-charge interactions [23], saturation mutagenesis of residues with high crystallographic B-factors [24], methods based on protein simulation and calculation of free energies [25—27] and comparison to homologous proteins including the ultra-stable proteins of thermophiles [28,29]. We reasoned that better predictions of mutant stability might be obtained by evaluating the unfolding temperature Tm in realistic yet computationally tractable simulations of protein unfolding. Here, we use a Monte Carlo protein unfolding approach (MCPU) with an all-atom simulation method and knowledge-based potential developed earlier in our lab [16,30,31] to simulate unfolding and predict melting temperatures for all possible single point mutants of E. coli Dihydrofolate Reductase (DHFR). DHFR is an essential enzyme in bacteria and higher organisms, and it is an important target of antibiotics [32] and anticancer drugs [33,34]. Its moderate size (18 kDa) makes it amenable to both simulation and experiment. As described in the Materials and Methods section, the Monte Carlo move set consists of rotations about torsional angles. At high temperature, the higher entropy of unfolded states overcomes the increase in energy due to loss of favorable contacts and torsional preferences, leading to unfolding. We eXperimentally determine melting temperatures and catalytic activities for several predicted stabilizing mutants, and for mutants combining multiple stabilizing mutations. Our approach allows us to identify several stabilized mutants of DHFR, and our prediction method marks an improvement over existing stability predictors such as Eris [19], FoldX [17] , and PopMusic [18]. Simulations of non-DHFR proteins likewise indicated that our method is useful as a general approach to simulate protein unfolding and select stabilizing mutations.

Results

Predicting the effects of mutations on protein stability from non-equilibrium unfolding simulations

However, despite recent progress in ab initio simulations of protein folding [15] this goal is not attainable for proteins of realistic size and biological relevance. Currently, non-equilibrium unfolding simulations are within reach for sufficiently large proteins and the question arises whether such simulations can be used to assess mutational effects on protein stability, which is an equilibrium property. The following analysis provides an affirmative answer to this question, under certain assumptions. Although the idea of obtaining equilibrium free energy differences from non-equilibrium measurements is not new [35], and protein stabilities have been calculated from molecular dynamics simulations using the Iarzynski equality, e.g., [36—38] , such simulations require application of an external steering force; in the present paper we report the use of multi-temperature Monte-Carlo unfolding simulations in obtaining protein stabilities. Assuming two-state unfolding kinetics [39—42] we can estimate the characteristic time required to cross the unfolding free energy barrier (in fact it is the time spent in the native state waiting for sufficient thermal fluctuation to cross the barrier) as:

When simulation time Tsim approaches If unfolding events are observed in simulation. The apparent “melting temperature”, i.e., the temperature at which unfolding events occur in simulations, therefore depends on the simulation time Tsim according to

However the effect of mutations on stability can be predicted from unfolding simulations. In order to see this we note that the mutational effect on protein stability AAG is related to the change in the unfolding free energy barrier AAG#, the difference between the WT barrier height and the mutant barrier height, shown in Fig. 1. Where i denotes the mutated amino acid and (p1- is the (p-Value for residue 1' which determines the fraction of interactions that this residue forms in the folding/unfolding transition state [40,43,44]. We therefore obtain where ATfif’P 2 T3? — Tip? (WT) is the shift in apparent unfolding temperature upon a specific mutation in the i-th residue. Introducing the relative (to WT) unfolding temperature ATffli) = AT;PP(i)/T;PP(WT) we get i.e. the mutational shift in observed unfolding temperature, normalized to the observed unfolding temperature of the wild-type at the same simulation condition does not depend on the simulation length, provided that the simulation is sufficiently equilibrated in the native basin so that the rules of transition state theory apply. The analysis of extensive kinetic and equilibrium data for multiple proteins shows that for the majority of mutations (except for a small fraction of residues that participate in the folding nucleus) (p1- % 0.24 with remarkable accuracy and con-i.e. AT;f’(i) is independent of simulation time and proportional to the equilibrium free energy effect of mutations, provided that simulations have equilibrated in the native basin of attraction.

Monte Carlo protein unfolding simulation

Unfolding steps of a sample trajectory are shown in Fig. 2, and a flowchart of the simulation method is shown in 81 Fig. The protein was subject to a brief MD energy minimization, beginning from the WT crystallographic native state, followed by unfolding simulations at each of 32 different temperatures using all-atom Monte-Carlo (see Materials and Methods section). As shown in figures 82 Fig—S4 Fig, the RMSD and total energy increased and the number of contacts decreased as each simulation proceeded, and with increasing temperature. (Here, temperature is given in arbitrary simulation units.) Plots of RMSD and contact number vs. temperature showed sigmoidal behavior, with a clearly identifiable transition temperature, and the melting temperature (Tm) could be determined by fitting to a sigmoidal function (Fig. 3). Plots of energy vs. temperature (SS Fig) were sigmoid-like, but with an additional rise in energy at low to intermediate temperatures, perhaps indicating pre-melting to a dry-molten globule state with loosened side chains but native-like topology [46,47]. This deviation from sigmoidal behavior becomes clearer as the simulation length is increased (S6 Fig).

Computational identification of stabilizing single point mutations

The simulated Tm values were calculated as described above. Of the 3,021 mutations, 523 mutations (17.3%) were predicted to have a stabilizing effect according to all three metrics (energy, contacts, and RMSD), while 421% of Contact number

WT DHFR unfolding curves from MC simulations, averaged over 2,000,000 simulation steps, with 50 replications. The Tm value was calculated based on the sigmoidal fit (solid blue line). (A) RMSD vs. simulation temperature. (B) Number of contacts vs. simulation temperature.

These predictions are in good agreement with statistical analysis of published experimental data and FoldX predictions [8,12]. The simulated Tm values evaluated using RMSD, total energy, and number of contacts are strongly correlated, as shown in Fig. 4A. The distribution of predicted melting temperatures (averaged over the 3 metrics) for all 3021 point mutants is shown in Fig. 4B. Next, we selected a subset of predicted stabilizing mutations for subsequent in depth computational and experimental analysis. To that end we selected the loci where multiple mutations were consistently predicted as stabilizing. Out of this set we selected one mutation at each loci which were predicted as most stabilizing. As a result we arrived at 23 single predicted stabilizing point mutants shown in 81 Table, which we deemed most promising for subsequent in depth computational and experimental analysis. Furthermore, five stabilizing mutations at different sites within DHFR, shown in Fig. 5, were combined to form the multiple mutants listed in Table 1, with the rationale that the combination of individual stabilizing mutants often yields more stable proteins, and these mutants were likewise subjected to computational and experimental analysis.

Computational test of the theoretical analysis

The first prediction is that the apparent unfolding temperature decreases as the length of the unfolding simulation increases (Equation 4). Secondly and most importantly the mutational change in relative (normalized to WT) apparent unfolding temperature is a) robust with respect to simulation time provided that simulations have equilibrated in the native basin and b) directly proportional to the effect of mutations on equilibrium protein stability (Equation 6). We test these predictions using MCPU simulations and experiment.

6). Indeed both predictions of our theoretical analysis are confirmed, i.e., the apparent unfolding temperature decreases with simulation time (Fig. 6A) while the relative unfolding temperature ATfjl is remarkably independent of simulation time (Fig. 6B). We note that due to the nature of the energy function used in our simulations, there is no obvious mapping of simulation temperature to real absolute temperature (i.e., in Celsius or Kelvin). Conversion of simulation temperature to physical temperature would require use of experimental data (e.g., WT unfolding temperature and deviation of temperatures over all mutants) and therefore would not provide a completely simulation or theo-ry-based prediction. Furthermore, as noted above, the apparent absolute value of the transition temperature in the Monte-Carlo unfolding approach depends on simulation time. Therefore, we used relative melting temperature, ATgf’ 2 ATS? / T31” ( WT), when comparing simulation results with experimental results.

In what follows we define the computational unfolding temperature Tm as averaged over Tm values determined using these three criteria.

Experimental characterization of predicted mutants

We cloned, expressed, and purified the 23 single point mutants of DHFR listed in 81 Table, as well as the multiple mutants listed in Table 1 (see Materials and Methods). The biophysical properties of the mutants were measured and compared with WT DHFR, as shown in $2

As many studies have shown that oligomerization can alter protein stability [23,48,49], we first tested whether mutations induce oligomerization and/or aggregation using the gel filtration method [48,50] and light scattering. The results indicated that all of the 23 mutants were monomeric at studied concentrations except for E154V, which appeared aggregation-prone. We excluded E154V from the subsequent analysis. As shown in S2 Table, all single mutants are catalytically active except for D27F. D27 is known to be a key catalytic residue of E. coli DHFR [51].

Both measures of stability were highly correlated, despite the fact that thermal unfolding was irreversible (S7 Fig). Of the selected 22 single point mutations, 10 mutations were stabilizing, according to their Tm or Cm values (82 Table). Given that statistically most random mutations are destabiliz-ing with only a small fraction (less than 18%) stabilizing [8,12], this statistically significant result (p = 0.002 under the null hypothesis that mutations are random) indicates that MCPU is an effective method for selecting stability-enhancing mutations.

In particular, the stability of the quintuple mutant (T68N,Q108D,T113V,E120P,S138Y) was found to be substantially higher than that of the wild type protein (Table 1), with Tm 7.2°C higher than WT, and Cm, the urea concentration at the mid-unfolding point, was 0.43M higher than WT. All multiple mutants were catalytically active, and the quintuple mutant and triple mutant (T113V,E120P,S138Y) were found to be more catalytically active than WT. We note that while combination of stabilizing mutations generally increases stability, the effect is less than additive (88 Fig); for instance, the quintuple mutant is about 4°C less stable than predicted under the assumption of additive ATm (a 7.2°C stability increase vs. predicted 9.6°C).

The correlation coefficient between experimental relative Tm and simulated relative Tm for the 42 mutants was about 0.65, as shown in Fig. 7A. To address the issue that both simulated Tm and DSC measurements are not strictly at equilibrium, we plotted the relation between simulated Tm and equilibrium measurement of stability in chemical denaturation by urea. The denaturation mid-transition urea concentration Cm and computationally determined unfolding temperature exhibit even a slightly higher correlation of r = 0.68 (Fig. 7B), demonstrating that our non-equilibrium simulation method shows good agreement with the equilibrium measurement of urea denaturation, as predicted by Equation 6.

As shown in Fig. 8A, the prediction accuracy is sensitive to the number of replications. To achieve reliable Tm predictions, at least 20 replications should be used. However, the number of MC steps did not greatly affect prediction accuracy, provided simulations were run for at least ~ 200,000 steps (see Fig. 8B). In the context of the theory developed in the earlier section: “Predicting the effects of mutations on protein stability from non-equilibrium unfolding simulations,” this initial equilibration period may allow time for equilibration within the native basin, after which simulation length does not appreciably affect the consistency of results with equilibrium stability measurements.

Stability and activity do not tradeoff for DHFR

Our data, however, paints a different picture for DHFR—of a weak positive correlation between Tm and kcat or kcat/KM (r = 0.46, p = 0.02 and r = 0.41, p = 0.03 respectively) with one notable outlier D27F, where the stabilizing mutation is made right in the active site (Fig. 9). The D27F mutant has high thermal stability but, as noted above, is not catalytically active, indicating that there is in fact a stability-activity tradeoff for this active-site residue.

Evolutionary analysis

Mutation of a non-consensus residue to the consensus residue generally resulted in protein stabilization [29]. In 4/ 16 of the experimentally stabilizing mutations, a residue was changed to the consensus residue, while only 2/29 destabilizing mutations resulted from a change to consensus. Likewise, in 18/29 destabilizing mutations, a residue was changed away from the consensus residue, while this was true for only 5/ 16 of stabilizing mutations.

Simulated melting temperatures by residue

10A). There is a weak positive correlation between minimum and maximum melting temperatures (r = 0.30, p 2 10‘4). Apparently, protein loci where mutations can cause significant stabilization are statistically less susceptible to destabilizing mutations and vice versa, which may be expected: once a residue is already at its most stabilizing amino acid variant, the protein cannot be stabilized further by mutation. Distinct outliers correspond to the loci with the strongest stabilizing or destabilizing effects of mutations. Interestingly, these outliers, which may represent structural weak spots in DHFR, tend to fall on the interface connecting the C-terminal beta hairpin and the rest of the protein

10B). This is in fact the interface that is the first to dissociate in the Monte Carlo simulations (see Fig. 2).

Comparison with other methods

(S3 Table). The MCPU performs better than these methods on DHFR mutants. PopMusic shows also strong performance with highly statistically significant 1’ = 0.55 between theory and experiment, however the limitation of this method is that it can consider only single point mutations. To further evaluate MCPU performance we tested it on four additional proteins from four different SCOP structural classes: the Cro repressor protein from bacteriophage lambda

m

Simulated T Replication Number

The effect of replication number and number of MC steps on simulation predictive power. (A) Correlation between simulated Tm and experimental Tm, averaging over different numbers of replications, for the DHFR wild type and mutants. Each protein was simulated for 2,000,000 MC steps, following MD minimization and equilibration at low temperature. (B) Correlation between the simulated Tm and experimental Tm with different numbers of MC steps and 50 replications, forthe DHFR wild type and mutants. Each protein was first simulated for the number of steps given on the x-axis, and the next 100,000 steps were averaged in determining the simulated Tm.

cat

Maximum Tm and Gln-25 ribonuclease T1 from Aspergillus oryzae (1RN1). Our predictions were compared with Eris and SDM. We did not compare MCPU results with FoldX and PopMusic as these mutations were selected in the training dataset for the two methods. As shown in Table 2, the correlation coefficient between MCPU predictions and the experimental Tm values, averaged over all proteins, is about 0.71, which is higher than that provided by Eris (-0.05), for Which predictions were quite poor for both DHFR and other proteins, and SDM (0.63). If we consider only the binary prediction of Whether a mutation is stabilizing or destabilizing, MCPU can correctly predict 11 out of 16 mutations, While Eris and SDM correctly classify 9 and 8 mutations respectively.

Entropy of the native state is an important contributor to stability

It is widely believed that in the low-entropy folded state energetic factors dominate. If so that would imply that we can get an equally good correlation between prediction and experiment by estimating the mutational effect on energy of the native state as is the case for most empirical methods. To that end we evaluated the correlation between the energy of the minimized (after long MC equilibration) native state and the experimental Tm and found only a weak correlation with experimental melting temperatures (Table 2, last column), indicating that protein entropy, which is accounted for in the MCPU, in addition to enthalpy, is important in determining protein stability.

Discussion

However using MCPU we were able to efficiently explore stabilities of all possible point mutants for an essential enzyme of a typical size (159 amino acids) in a manageable amount of computational time (approx. one hour for every 1,000,000 MC steps). Although the use of rapid Monte Carlo simulations reduces simulation time and allows for a greater number of replicates, our method to predict stability effects of mutations based on non-equilibrium unfolding simulations represents a general approach that could be modified for use with conventional MD simulations, especially given the current rate of improvement in simulation speed and accuracy.

Low (p-value residues, which acquire contacts with other residues late in the folding process and lose contacts early in the unfolding process [14] constitute the majority of residues in proteins, with (p-value roughly constant around 0.24 as noted in [45]. Combining this observation with assumptions of transition state theory, we found that for the majority of residues (those not part of the folding nucleus [14,57] exhibiting anomalously high (p-values) the observed simulation Tm relative to WT is proportional to the equilibrium stability change AAG, as verified by simulation and experiment. We establish that relative Tm is independent of simulation length, demonstrating that non-equilibrium simulations can in fact be used to quantify relative protein stability.

It has been shown that the source of ultra-stability in hyperthermophiles generally arises from slowing the unfolding rate, rather than increasing the folding rate [28] , so our method may be particularly suitable for discovering biologically relevant stabilizing mutations. In addition, our results might be particularly applicable to in vivo studies, where protease digestion and/or aggregation proceed from the partially-unfolded state. We note, however, that some stabilizing residues predicted by MCPU lie in the region of the protein that is late to unfold in simulations, including I61V, which raises the experimental melting temperature by 1.7°C. These mutants, along with the destabilized outlier 1155A for which relative Tm depends on simulation length (Fig. 6), are appealing candidates for further study, as they may reflect a breakdown in the simplifying assumptions of 2-state kinetic theory for proteins.

This conclusion was reached in [53,58], based on the exploration of stability effects of mutations in the active site of beta-lactamase [53] and rubisco [58]. Fersht and coauthors also found several stabilizing mutations in the active site of Barnase rendering the protein inactive [59]. While we observe a similar effect with the D27F mutation in DHFR, Fig. 9 shows that exploring only mutations in the active site provides a biased view on the tradeoff between activity and stability. Rather a vast majority of mutations throughout the protein show a qualitatively opposite trend. The likely explanation of the distinction between an apparent tradeoff when mutations are made in the active site and the opposite trend for mutations outside of the active site is that “carving” an active site requires special selection of catalytic amino acids, which could indeed have a destabilizing effect, overall. However our observation of a small positive correlation argues against an obligate relation between global protein dynamics and activity for DHFR, at least for the aspects of dynamics that are correlated with stability. Warshel and colleagues reached a similar conclusion in their theoretical analysis of the role of dynamics for DHFR and other proteins in [60]. This point has likewise been made by Bloom et al. [11], who noted that a number of proteins have been stabilized experimentally without loss of activity, and Taverna and Goldstein argued that marginal stability is an inherent property of proteins due to the high dimensionality of sequence space and not due to a requirement of reduced stability in order to generate sufficient flexibility [61].

It is also important to note that a weak yet statistically significant positive correlation between activity and stability for DHFR can be revealed only when stabilizing mutations are included in the analysis. Our earlier study [48] analyzed a smaller set of primarily destabilizing mutants and did not reveal any statistically significant trend (positive or negative) in the stability-activity relation for DHFR.

It has been postulated that protein stability places a fundamental constraint on the evolutionary pathways available to a protein [29,62] which has particular significance in the development of antibiotic resistance: higher protein stability can provide the microorganism with an increased capacity to evolve to evade antibiotic drugs [63] or, more generally, with capacity to evolve new functions [62]. We plan to use an approach developed in our lab [48] to endogenously introduce stabilized DHFR mutants into the bacterial chromosome and we will evaluate mutant fitness relative to wild-type using growth rates and competition experiments. These experiments will allow us to assess whether an evolutionary tradeoff exists between stability and fitness in vivo, particularly in the presence of antibiotics.

Comprehensive experimental analysis of fitness and/ or stability effects of mutations [64] could be useful in assessing the predictive capabilities of this method. In addition to predicting mutant stabilities, MCPU can provide atomic-detail molecular trajectories to rationalize the stability effects of mutations; such analysis is left to future study.

Materials and Methods

Monte Carlo simulations

Briefly, the energy function is a sum of contact energy, hydrogen-bonding, torsional angle, and sidechain torsional terms, with an additional term describing orientation of nearby aromatic residues. The move set consists of rotations about (I), q], and x dihedral angles, with bonds and angles held fixed. Moves are accepted or rejected according to the Metropolis criterion.

An initial minimization was carried out in NAMD [67] for 5,000 steps, using the default minimization algorithm and par_a1127_prot_lipid.inp parameter file (without waters). An additional minimization step was carried out by running the Monte Carlo simulation program at low temperature (0.100 in simulation units) for 2,000,000 steps. A 2,000,000-step simulation was then run at each of 32 temperatures, averaging over all 2,000,000 steps to obtain Energy, RMSD, and number of contacts. These results were averaged over 50 simulations, for each temperature. Data was then plotted and fit to a sigmoid to obtain the computationally-predicted melting temperature, for each of Energy, RMSD, and number of contacts. To assess dependence of melting temperature on simulation length, longer simulations of 20,000,000 steps were carried out with 30 replications, averaging over the final 2,000,000 steps. For DHFR, 1,000,000 steps took approximately one hour of simulation time, on a single CPU.

DHFR protein sequences from 290 bacterial species were aligned using the program MUSCLE and online server [68]. MATLAB, with the Bioinformatics Toolbox, was used to create sequence logo representations and to determine the consensus sequence.

Effect of number of replications on simulation accuracy

As shown in S9 Fig, the prediction accuracy is sensitive to the number of replications, but converges to a constant value after approximately 20 replications. In addition, we saw that increasing the number of MC steps beyond 2,000,000 steps does not increase prediction accuracy when the protein has been simulated with at least 20 replications, despite the fact that not all simulations have converged by 2,000,000 steps (82 Fig—S4 Fig).

Simulation analysis

Sigmoidal fits were accomplished using the module “Sigmoidal, 4PL” using the software pro

Method availability

The tool is accessible from Shakhnovich lab website http://faculty.chemistry.harvard.edu/ shakhnovich/ software

Site-directed protein mutagenesis of DHFR

Single point mutations of DHFR were constructed based on a two-step PCR-mutagenesis strategy [70], in which the template for the PCR is the plasmid of WT DHFR. The multiple-mutant variants of DHFR were constructed based on the same method with the single point mutation, but the template of PCR was the plasmid of the corresponding dhfr mutant. To verify the mutations of dhfr, DNA sequencing was performed at the GENEWIZ Incorporation (MA, US). The verified plasmids were transformed into competent E. coli BL21(DE3) cells for expression.

Protein expression and purification

WT DHFR and all mutants used in this study were cloned into a pET24 expression vector and overexpressed in the BL21(DE3) pLys E. coli strain.

When the OD600 of the culture reached 0.6, isopropyl B-D-1-thiogalactopyranoside (final concentration, 0.4 mM) was added. Cultures were incubated for an additional 12—16 h at 25°C. The cells were then collected by centrifugation and disrupted by sonication. The recombinant proteins were purified with Ni-NTA Superflow (QIAGEN, U.S.) according to the manufactur-er’s instructions. Then, the collected protein sample was run with SuperdeX 75pg Column and was desalted with the desalting Column in AKTA protein purification system (GE Healthcare, US). The final concentration of the purified protein was determined using the BCA protein assay kit (PIERCE CHEMICAL, USA) or the NanoDrop instrument (GE Healthcare, U.S.).

Enzyme kinetics

A Scientific stopped flow apparatus, RX.2000 Rapid kinetics system (Applied Photophysics, UK) was used with absorbance monitoring at 340 nm, under single-turnover conditions. NADPH was preincubated with DHFR for 5 min in syringe 1 at the temperature 25°C in a thermostated syringe compartment, and then the reaction was initiated by rapidly mixing the contents with dihydropholate (DHF) from syringe 2. The final assay conditions are 25 nM DHFR, 120 M NADPH(D), and 25 [AM DHF in MTEN buffer (50 mM 2-(N-morpho-lino)ethanesulfonic acid, 25 mM tris(hydroxymethyl)aminomethane, 25 mM ethanolamine, and 100 mM sodium chloride, pH 7.6). The kinetics parameters (kcat and KM) were derived from progress-curves analysis using Global Kinetic explorer [72].

Stability measurements

Briefly, DHFR proteins in Buffer A (10 mM potassium-phos-phate buffer pH 7.8 supplemented with 0.2 mM EDTA and 1 mM beta-mercaptoethanol) were subjected to a temperature increase of 1°C/min between 20 to 90°C (nano-DSC, TA instruments, U.S.), and the evolution of heat was recorded as a differential power between reference (buffer A) and sample (120 [AM protein in buffer A) cells. The resulting thermogram (after buffer subtraction) was used to derive apparent thermal transition midpoints (Tm app). Thermal unfolding appeared irreversible for all DHFR proteins tested [48], and the two state scaled model provided in NanoAnalyze software (TA INstruments, U.S.) was used to fit the Tm app value. The mutants constructed in this study and the ones published earlier [48] were determined with different DSC instruments with slightly different calibration leading to a small offset of about 2°C for the WT DHFR for earlier published data[48].

Proteins (25 uM in buffer A) were diluted in urea (0.2 mM increments up to a final urea concentration between 0 and 6 M), preequilibrated overnight at 25°C for 3 hours, and the change in the folded fraction was monitored by a circular dichroism signal at far-uv wavelength (221 nm) at 25°C (I-710 spectropolarimeter, Iasco). Fitting to a two-state model was used to derive the chemical transition midpoint (Cm).

Supporting Information

RMSD is averaged over 50 replications. (TIFF) S3 Fig. The total energy of DHFR wild type and mutants 11 15A and 1155T vs. Monte Carlo step at temperatures from 0.1 to 3.2. Total energies are averaged over 50 replications. (TIFF) S4 Fig. The number of contacts for DHFR wild type and mutants 11 15A and 1155T vs. Monte Carlo step at temperatures from 0.1 to 3.2. Number of contacts were averaged over 50 replications.

Data points are averaged over the 2,000,000 step simulation and 50 separate runs. (TIFF)

Temperature simulated melting curve for WT DHFR and mutants 11 15A and 1155T, averaging over the last 2,000,000 steps of a 20,000,000- step simulation, and 30 separate runs. The mutant melting temperatures show the same trend as in SS Fig, although all melting temperatures are shifted to lower values (since the protein is given more time to unfold at each temperature), and the curve deviates more notably from a sigmoid. (A) Overlaid data points from all three mutants. (B) Fits to a sigmoidal function (blue line), for each of the three mutants. (TIFF) S7 Fig. Correlation between urea denaturation midpoint, Cm, and melting temperature

Change in experimental melting temperature relative to WT is predicted by summing melting temperature changes of individual mutants. This predicted ATm is plotted relative to the observed ATm

Sequence Logo for the alignment, generated using the MATLAB Bioinformatics Toolbox. (TIF)

The simulated Tm values of the selected single point mutations and WT DHFR. (DOCX)

The experimental results for all single point mutations of DHFR. (DOCX)

Acknowledgments

We thank Bharat Adkar and Shimon Bershtein for help with the purification and Biophysical characterization of the DHFR mutants and Muyoung Heo for help and guidance in the use of the Monte Carlo simulations.

Author Contributions

Performed the experiments: IT ICW. Analyzed the data: EIS IT ICW. Contributed reagents/materials/ analysis tools: AW. Wrote the paper: EIS IT ICW.