Thermal Stabilization of Dihydrofolate Reductase Using Monte Carlo Unfolding Simulations and Its Functional Consequences

<div><p>Design of proteins with desired thermal properties is important for scientific and biotechnological applications. 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 <i>r</i> = 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 denaturation 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 <i>positive</i> 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 <i>E</i>. <i>coli</i> DHFR mutant with measured denaturation temperature 7.2°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.</p></div

A sample WT DHFR unfolding trajectory at simulation temperature 1.5 (arbitrary simulation units).

<p>In MC simulations, separation of the C-terminal beta hairpin from the rest of the protein (steps 1,000,000 through 1,200,000) is an early event in the unfolding process.</p

Simulated <i>T</i>_m values, based on RMSD, Total Energy and Contact number.

<p>(A) Scatter plot of <i>T</i>_m (RMSD) vs. <i>T</i>_m (Total energy), with <i>T</i>_m (contact number) represented by color (see color bar to right of plot). The green ball denotes WT and the gold ball denotes the destabilized mutant I155A. The correlation coefficients of simulated <i>T</i>_m between RMSD and total energy, RMSD and Contact number, and Contact number and total energy were 0.68, 0.79 and 0.84, respectively. (B) Histogram of <i>T</i>_m values, determined by averaging the values obtained from RMSD, energy, and contact number. The vertical red line denotes WT <i>T</i>_m.</p

WT DHFR unfolding curves from MC simulations, averaged over 2,000,000 simulation steps, with 50 replications.

<p>The <i>T</i>_m value was calculated based on the sigmoidal fit (solid blue line). (A) RMSD vs. simulation temperature. (B) Number of contacts vs. simulation temperature.</p

The simulated and experimental results of the selected single point mutants and WT.

<p>Note: The data were averaged over 50 replications. 2,000,000 MC steps were simulated in total, and the last 1,000,000 steps were used to calculate <i>T</i>_m.</p><p>The units: <i>T</i>_m: °C, <i>C</i>_m: M, <i>k</i>_cat: s⁻¹, <i>k</i>_cat∕<i>K</i>_M: s⁻¹ μM⁻¹</p><p>The simulated and experimental results of the selected single point mutants and WT.</p

Maximum stabilization and destabilization induced by mutations at each residue position.

<p>(A) Plot comparing the minimum and maximum simulated <i>T</i>_m values, for each residue across all 19 simulated mutants. <i>T</i>_m is normalized to WT, by dividing each <i>T</i>_m value by the simulated WT <i>T</i>_m = 1.489 simulation units. Outliers are circled in purple (left), green (middle) and orange (right). (B) DHFR with outlier residues colored according to the color scheme from (A). Purple: residues F153, W30, Y111, L156, L110. Green: residues A107, I155, L112, H114. Orange: residues A6, E154. Excluding outlier residues, the C-terminal beta hairpin is colored yellow, and the rest of the protein is colored cyan.</p

Correlation between the relative simulated and experimental <i>T</i>_m values.

<p>(A) Plot of simulated <i>T</i>_m vs. experimental <i>T</i>_m. The relative <i>T</i>_m values were calculated by normalizing to WT: (<i>T</i>_m(mutant)-<i>T</i>_m(wild type))/ <i>T</i>_m(wild type). Experimental values from this study and from Bershtein <i>et al</i>. [<a href="http://www.ploscompbiol.org/article/info:doi/10.1371/journal.pcbi.1004207#pcbi.1004207.ref048" target="_blank">48</a>] are included. WT is shown as a blue triangle. <i>r</i> = 0.65, <i>p</i> = 3 x 10⁻⁶. (B) Plot of simulated <i>T</i>_m vs. experimental C_m. <i>r</i> = 0.68. <i>p</i> = 6 x 10⁻⁷.</p

Correlation between DHFR activity and stability.

<p>WT is shown as a blue triangle; D27F is shown as a red diamond at zero activity. (A) Plot of <i>k</i>_cat vs. experimental relative <i>T</i>_m. <i>r</i> = 0.46, <i>p</i> = 0.02 (excluding outlier D27F). (B) Plot of <i>k</i>_cat/<i>K</i>_m vs. experimental relative <i>T</i>_m. <i>r</i> = 0.41, <i>p</i> = 0.03 (excluding outlier D27F).</p

Simulation results on non-DHFR proteins.

<p>Error number and error rate describe the number and fraction of mutations not predicted in the correct direction (stabilizing vs. destabilizing)</p><p>Simulation results on non-DHFR proteins.</p

The effect of replication number and number of MC steps on simulation predictive power.

<p>(A) Correlation between simulated <i>T</i>_m and experimental <i>T</i>_m, 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 <i>T</i>_m and experimental <i>T</i>_m with different numbers of MC steps and 50 replications, for the 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 <i>T</i>_m.</p