Computational Biology and Chemistry 33 (2009) 404–407
Contents lists available at ScienceDirect
Computational Biology and Chemistry journal homepage: www.elsevier.com/locate/compbiolchem
Brief communication
Modeling the folding and hydrogen production of Clostridium acetobutylicum and Clostridium saccharobutylium mutants using electrostatic potential surfaces and molecular dynamics Mark A. Plummer a,∗ , Scott M. Plummer b a b
MPr&d, LLC, 87 Silver Fox Drive, Greenwood Village, CO 80121, USA H2 OPE Biofuels, LLC, 2433 Vine, Denver, CO 80205, USA
a r t i c l e
i n f o
Article history: Received 8 January 2009 Received in revised form 19 July 2009 Accepted 25 July 2009 Keywords: Electrostatic potential surfaces Modeling protein folding Hydrogen production
a b s t r a c t Electrostatic potential surfaces (EPS) were used with molecular dynamics to model the folding mechanisms and kinetics of hydrogenase mutants from wild types Clostridium acetobutylicum and Clostridium saccharobutylium. The purpose of the EPS approach was to incorporate long range electrostatic forces between widely separated regions of the mutants which contain 575 amino acids. Also, it was demonstrated that the ratio of positive to negative EPS of unfolded mutants could be used to predict the production of molecular hydrogen from the folded mutants. Using the prediction model, mutant compositions were determined that should yield hydrogen of up to 40 times that obtainable from the wild type C. acetobutylicum. It is expected that the developed EPS techniques can be used to study the folding of other proteins and to predict the reactivity of the folded protein structures. © 2009 Elsevier Ltd. All rights reserved.
1. Introduction This work had dual objectives. The first was to see if the folding of a hydrogenase protein could be more efficiently modeled using electrostatic potential surfaces (EPS) along with molecular dynamics (MD) versus using MD only. EPS approaches have been previously used to successfully model adsorption and solvent mediated processes (Plummer et al., 2004; Plummer and Cowley, 2007). In this work, EPS was used to model long range electrostatic forces between widely separated amino acids in a hydrogenase. The second objective was to develop a model based on the EPS of unfolded hydrogenases to predict hydrogen productions from Clostridium acetobutylicum and Clostridium saccharobutylium mutants which contain about 575 amino acids. The model would then be used to find hydrogenases yielding the highest levels of hydrogen production utilizing an experimental data set of 256 mutants (Plummer, 2007). 2. Theoretical approach The modeling of protein folding ranges from quick calculations with low accuracy, e.g. knowledge based systems, to accurate and time consuming Ab Initio computations (Zaki and Bystroff, 2007; Merz and LeGrand, 1994; Kretsinger et al., 2004). In most of these
∗ Corresponding author. Tel.: +1 303 771 0749. E-mail address:
[email protected] (M.A. Plummer). 1476-9271/$ – see front matter © 2009 Elsevier Ltd. All rights reserved. doi:10.1016/j.compbiolchem.2009.07.016
calculations, electrostatic interactions are a key element. However, many computational chemistry programs only make short range electrostatic calculations. Hence, a need exists to incorporate long range electrostatic interactions of widely separated amino acid regions in protein folding computations (Braun-Sand and Warshel, 2005). It is a proposal of this work that a surface of constant electrostatic potential can be used as a measure of the effort needed to share or transfer an electron from one region of a molecule to another region. And, that electron sharing in hydrogenases should affect their folding mechanisms and kinetics. If EPS were all negative or positive, there would be no electron sharing and no folding of widely separated regions. Hence, it is proposed that a mixture of negative and positive EPS regions will result in an optimally folded hydrogenase and the maximum amount of hydrogen production. It is also proposed that positive and negative EPS regions can be quickly folded and matched by variation of bond angles between regions followed by MD structuring and energy minimization. This approach should minimize computational time for modeling folding kinetics and mechanisms over that required if only MD calculations were used. The mechanisms and kinetics of protein folding have been related to the folding energy divided by the protein surface area that is accessible to water (Sternberg, 1996). However, it is proposed in this work that the ratio of water accessible volume to water accessible surface area is an improved parameter. Accessible surface areas and volumes are calculated by mathematically rolling a water probe sphere of 1.4 Å in radius over
M.A. Plummer, S.M. Plummer / Computational Biology and Chemistry 33 (2009) 404–407
405
Table 1 Parameters for electrostatic potential surface, molecular dynamics and energy minimization calculations. Electrostatic potential surface Grid mesh size—fine Contour level—64
Horizontal and vertical grid points—64 Rendering—wire mesh
Molecular dynamics—OPLS Scale factors: dielectric—1 Electrostatic—0.5, van der Waals—0.125 Heat time—0.1 ps Cool time—0.0 ps Step size—0.001 ps Bath relaxation time 0.001 ps Start temperature—5 K below constant operating temperature Energy minimization Optimization algorithm—Polak-Ribiere (conjugate gradient) Terminal RMS gradient—0.2 kcal/[Å mole]
the protein making maximal contact without overlap (Creighton, 1997). 3. Modeling techniques An X-ray crystal structure for the hydrogenase used in the folding studies was not available as the target for error evaluations. Hence, the target structure for this hydrogenase was generated using the SWISS-MODEL program (Guex and Peitsch, 1997) followed by MD structuring and energy minimization. In all other work, the HyperChemTM program (Hypercube, 2002) and the OPLS united atom force field were employed. The OPLS parameters were obtained via liquid simulations (Jorgensen and Tirado-Rivers, 1988). Hence, water is implicitly included in the calculations. Detailed HyperChem parameters are given in Table 1 The initial step was to construct each hydrogenase in an ␣helical structure followed by energy minimization. Then areas for each of the positive and negative EPS regions at 0.1 electrostatic units per Å were calculated. If there was overlap, the smaller area, positive or negative, was subtracted from the larger area. Finally, the positive areas were summed and divided by the sum of the negative areas. Then folding studies were done in two steps. First, MD only calculations were made. Next, the major positive and negative EPS surfaces of the ␣-helical structure were closely located in space by bending atom to atom angles between these regions. Then MD folding and energy minimization calculations were performed. Finally, the initial EPS fold was again ESP folded five more times. After each ESP fold, MD and energy minimization calculations were made. However, each ESP fold was made on structures that had not been MD folded and energy minimized. In MD calculations, folded structures were obtained after each time interval of 25–50 picoseconds (ps). Then, the structure was energy minimized and compared to the target structure to obtain a root mean squared (RMS) matching error. Also at each time interval, a water accessible surface area and volume were calculated along with a folding energy. Folding energy was obtained by subtracting the total molecular energy of the initial ␣-helical structure from that for the folded structure. 4. Results and discussion MD only calculations showed that folding starts at both ends of the entire hydrogenase mutant and works its way toward the middle. The center point for this folding is the amino acid 267. Hence to obtain shorter computation times, the rest of the folding work dealt with the 267–575 portion of this mutant. This region shows the greatest diversity of positive (green) and negative (red) EPS regions in the initial ␣-helical structure as shown in Fig. 1a.
Fig. 1. Molecular structure of the 267–575 amino acid region of the studied hydrogenase mutant. (a) Hydrogenase prior to EPS folding. (b) Hydrogenase after the first EPS folding.
Consequently, the 267–575 portion of the target hydrogenase was used in all folding error evaluations. This target was MD structured over conditions of 323–773 K and 50–150 ps followed by energy minimization at 0 K. The most thermodynamically feasible folding energy was obtained at 573 K and 100 ps. However, the difference between the MD structure and the original SWISS-MODEL structure for this region was minimal. Protein folding kinetics has been shown (Huang, 2005) to exhibit two folding regions. First, the protein quickly folds from a linear state into a globular structure. Next, the globular structure folds slowly to the final tertiary state. The results from this work support these literature findings. That is, a correlation between RMS error and the water accessible volume to surface area ratio was found to exhibit two straight line segments of considerable different slopes. It was concluded that these separate straight line segments represent the fast and slow folding regions. Single correlations of this type were obtained for all folding temperatures and times studied. In the MD only work at 323 K and up to 150 ps, RMS error decreased rapidly from 128 Å for the initial ␣-helical structure to the globular state which exhibits an RMS error of about 21 Å. Over this region, the water accessible volume to surface area ratio increased from 1.33 to 1.40 as depicted in Fig. 2a. The slope in this fast region is 1252 Å/Å. Further increases in MD only folding time from 200 to 300 ps and in folding temperature from 323 to 673 K did not change the globular structure or reduce its RMS error. Hence, MD only folding could not generate the slow folding mechanisms needed to get from the globular state to the tertiary target structure. The probable reason for this is that MD only folding located a globular structure with a folding energy nearly equal to that for the tertiary state as shown in Fig. 2b. Therefore, there was an insufficient folding energy gradient to proceed from the globular to the tertiary structure. The same conclusion was obtained when using folding time as the correlating parameter. The first EPS fold centered on amino acids 315 and 465. As depicted in Fig. 1b, this fold reduced the total area of the negative red regions compared to that for the positive green regions. This EPS folding was then followed by MD structuring at conditions of 50–100 ps and 323–873 K followed by energy minimization at 0 K. The most thermodynamically feasible structure was obtained at 673 K and 100 ps. Both fast and slow folding regions were noted in the correlation of RMS error versus volume to surface area ratio, as demonstrated in Fig. 3a, for the first EPS/MD fold. The correlation slope for the fast region was essentially equivalent at 1213 Å/Å to that obtained with MD only folding. The slow folding region exhibited a correlation slope of 74 Å/Å and minimum RMS error of 21.9 Å at a volume/surface area ratio of 1.45. More importantly, the correlation for the slow folding region extrapolates to a RMS error of 7 Å at the target’s volume to surface area ratio of 1.66 Å. This last result shows that further EPS folding was warranted.
406
M.A. Plummer, S.M. Plummer / Computational Biology and Chemistry 33 (2009) 404–407
Fig. 2. Fit of folded structure to tertiary structure and folding energy versus volume to surface area ratio of folded structure via MD only folding. (a) RMS fit. (b) Folding energy.
Fig. 3. Fit of folded structure to tertiary structure versus volume to surface area ratio of folded structure via EPS/MD folding. (a) For the first EPS/MD fold. (b) For EPS/MD folds one through six.
EPS/MD folds two through six resulted in decreasing the minimum RMS error down to 18.1 Å and increasing the volume to surface area ratio up to 1.51 Å. The RMS error versus volume to surface area ratio correlation for all the six folds exhibits a slope of 103 Å/Å for the slow folding region as shown in Fig. 3b. More importantly, this correlation extrapolates to a RMS error of 1.8 Å at the volume to surface area ratio of the tertiary target structure of 1.66 Å. The MD conditions in this portion of the effort were varied from 323 to 973 K and from 50 to 200 ps. The most thermodynamically feasible folded structure was obtained at 673 K and 150 ps. All of the above RMS errors were obtained with folded structures containing the backbone and side chains. When only the backbones of the folded structures were compared to the backbone for the target structure, no significant improvement in RMS error was noted. The folding energy versus the volume to surface area ratio correlation for all the six EPS/MD folds also exhibits a fast and a slow folding region. And, the correlation for the slow region extrapolates to within 53 kcal/gmole of the folding energy (3504 kcal/gmole) of the target tertiary structure as depicted in Fig. 4. These RMS error and energy results suggested that additional EPS/MD folding should yield the target structure whereas MD only folding would not accomplish this objective. Unfortunately, further EPS folding did not significantly change the ratio of positive to negative EPS regions and these regions remained essentially like those for the sixth fold. Also, the minimum RMS error and volume to surface area ratio could not be changed compared to those for the sixth fold. From this it is concluded that further improvements toward obtaining the tertiary structure would have to include the region containing amino acids 1 through 266.
Based on the above results, the EPS concept was tested for predicting hydrogen productions of hydrogenase mutants using their unfolded ␣-helical structure. Initial -sheet structures were also evaluated. However, ␣-helical structure yielded more definable positive and negative EPS regions than did the -sheet structure. Hence, the -sheet structure was not extensively evaluated. In the prediction modeling, a seven hydrogenase data set (Plummer, 2007) was used which included six mutants and the wild type C. acetobutylicum as the control. The experimental hydrogen
Fig. 4. Folding energy for EPS/MD folds one through six versus volume to surface area ratio of the folded structure.
M.A. Plummer, S.M. Plummer / Computational Biology and Chemistry 33 (2009) 404–407
Fig. 5. Measured and predicted hydrogen productions versus ratio of positive to negative EPS of unfolded structures.
productions for the mutants varied from 0 to 4 times that obtained with the wild type control. Also, several correlating parameters were tested. Since in the first step of hydrogen production an electron carrier must be attracted to the hydrogenase, the total positive EPS area was evaluated. For the same reason, the total positive EPS area minus the total negative EPS area was tested. Neither of these two parameters yielded a reliable correlation with hydrogen production. Then, the ratio of total positive EPS area to total negative EPS area was evaluated. The reason behind this approach is that folding is affected by long range electrostatic forces as demonstrated by the above results. Those mutants exhibiting a certain EPS ratio between zero (no positive EPS regions) and infinity (no negative EPS regions) would fold correctly to yield the maximum hydrogen production. And, hydrogen production should taper off to zero as the EPS ratio approaches zero or infinity. Log Normal and Gaussian equations fit these properties. The best fit between hydrogen production and the EPS ratio was obtained with the Log Normal equation yielding a Chi squared error of 0.2 as shown in Fig. 5 with black squares. Using the Log Normal equation, maximum hydrogen production of 40 times that obtained with the wild type control is predicted for a positive to negative EPS ratio of 41. Based on the above results, fifteen of the 256 mutants produced in the experimental work were evaluated. Five of these mutants have predicted hydrogen productions in the range of 6–40 times that for the wild type control – see the black triangles in Fig. 5. The remaining nine were estimated to yield no hydrogen production at all. It is realized that the maximum predicted hydrogen production is ten times that found experimentally. However, the model is based on sufficient theory to proceed with further experimental work to verify the predicted hydrogen production values. 5. Summary This effort has shown that incorporating long range electrostatic forces in the modeling of protein folding has significant advantages in obtaining the final tertiary structure. Long range forces were incorporated by folding and matching positive to negative EPS over widely separated regions followed by standard MD modeling.
407
The EPS/MD approach yielded results supporting literature data showing that protein folding exhibits a fast and a slow folding region. The two folding regions were demonstrated using RMS fitting error and the water accessible volume to surface area ratio exhibited by the folded hydrogenase. When RMS error was correlated to the volume to surface area ratio, two straight line segments of considerably different slopes were obtained. The larger slope segment representing the fast folding region occurs at smaller volume to surface ratios. The smaller slope segment representing the slow folding region occurs at larger ratios. And, the slow folding segment extrapolates to near zero RMS error at the volume to surface area ratio equal to that exhibited by the tertiary structure. A correlation between EPS/MD folding energy and volume to surface area ratio also exhibits two straight line segments representing the fast and slow folding regions, respectively. And again, the slow folding segment accurately extrapolates to the folding energy of the tertiary state at its volume to surface area ratio. With MD only folding, which incorporates just short range electrostatic forces, the slow folding region could not be computationally generated. This was demonstrated by correlating RMS fitting error with either folding time or the volume to surface area ratio parameter. Likewise, MD only folding energies could only be generated for the slow folding region. Lastly, it has been demonstrated that a ratio of positive to negative EPS surfaces of unfolded hydrogenases can be used to accurately predict experimental hydrogen productions of the folded tertiary structures. Using a Log Normal correlation between positive to negative EPS ratio and experimental hydrogen production, five hydrogenases were found that should yield hydrogen productions of 23–40 times that of the wild type control. References Plummer, M.A., Kim, M., Way, J.D., Baldwin, R.M., 2004. Determination of mechanisms via computational chemistry for xylene and hydroxylnaphthalene separations on beta-cyclodextrin. Mol. Phys. 102 (2), 183–189. Plummer, M.A., Cowley, S.W., 2007. Solvent effects in the conversion of anthraquinone to anthrahydroquinone and simultaneously hydrogen sulfide to hydrogen and sulfur. Mol. Phys. 105 (4), 437–443. Plummer, S.M., 2007. Photosynthetic Hydrogen Production from the Green Alga Chlamydomoas reinhardtii. Colorado School of Mines T-6295, Golden, CO. Zaki, M.J., Bystroff, C., 2007. Protein Structure Prediction, second ed. Humama Press, New Jersey. Merz Jr., K.M., LeGrand, S.M., 1994. The Protein Folding Problem and Tertiary Structure Prediction. Birkhauser, Boston. Kretsinger, R.H., Ison, R.E., Hovmoller, S., 2004. Prediction of protein structure. In: Brand, L., Johnson, M.G. (Eds.), Methods in Enzymology 383, Numerical Computer Methods Part D. Elsevier Academic Press, San Diego, pp. 1–27. Braun-Sand, S., Warshel, A., 2005. Electrostatics of proteins: principles, models and applications. In: Buchner, J., Kiefhaber, T. (Eds.), Protein Folding Handbook, vol. 1. Wiley–VCH, Weinheim, pp. 163–199. Sternberg, M.J.E., 1996. Protein structure prediction—principles and approaches. In: Sternberg, M.J.E. (Ed.), Protein Structure Prediction, A Practical Approach. Oxford University Press, New York, pp. 1–27. Creighton, T.E., 1997. Protein Folding. W.H. Freeman & Company, New York. Guex, N., Peitsch, M.C., 1997. SWISS-MODEL and the Swiss-Pdb Viewer: an environment for comparative protein modeling. Electrophoresis, 18. Hypercube, Inc., 2002. Gainesville FL, 32601, USA. Jorgensen, W.L., Tirado-Rivers, J., 1988. The OPLS potential function for proteins. J. Am. Chem. Soc., 113. Huang, K., 2005. Lectures on Statistical Physics and Protein Folding. World Scientific Publishing, New York, pp. 95–101.