Dynamic Formation and Breaking of Disulfide Bonds in Molecular Dynamics Simulations with the UNRES Force Field.

Many proteins contain disulfide bonds that are usually essential for maintaining function and a stable structure. Several algorithms attempt to predict the arrangement of disulfide bonds in the context of protein structure prediction, but none can simulate the entire process of oxidative folding, including dynamic formation and breaking of disulfide bonds. In this work, a potential function developed to model disulfide bonds is coupled with the united-residue (UNRES) force field, and used in both canonical and replica exchange molecular dynamics simulations to produce complete oxidative folding pathways. The potential function is obtained by introducing a transition barrier that separates the bonded and nonbonded states of the half-cystine residues. Tests on several helical proteins show that improved predictions are obtained when dynamic disulfide-bond formation and breaking are considered. The effect of the disulfide bonds on the folding kinetics is also investigated, particularly their role in stabilizing folding intermediates, resulting in slower folding.

Resolution of the excitation-emission spectra of FMN in rigid poly(vinyl alcohol) matrices.

The excitation-emission spectra of flavin mononucleotide (FMN) were measured in rigid PVA films for concentrations ranging from 6.92 x 10(-4)M to 1.03 M. The theoretical three-linear decomposition of the excitation-emission spectra indicated the presence of two absorption and emission centers corresponding to FMN monomer and dimer, respectively. The component of the fluorescence profile corresponding to the FMN monomer has a large negative part which is the mirror image of the emission band profile of the dimer. The elimination of this part by taking a linear combination of the emission components of the monomer and of the dimer resulted in emission spectrum, which is in a very good agreement with the monomer spectrum measured directly. The appearance of a negative part of the monomer emission profile obtained by trilinear decomposition of the emission-absorption spectra of FMN can be explained in terms of the non-radiative reverse energy transfer from the FMN dimers to the FMN monomers. The presented results confirm that the FMN molecules in rigid PVA form dimers but not higher order aggregates. Moreover, they enable to obtain fluorescence spectra of dimers and suggest that FMN dimers may take part in the process of non-radiative energy transfer occurring in photoreception phenomena.

Physics-based protein-structure prediction using a hierarchical protocol based on the UNRES force field: assessment in two blind tests.

Recent improvements in the protein-structure prediction method developed in our laboratory, based on the thermodynamic hypothesis, are described. The conformational space is searched extensively at the united-residue level by using our physics-based UNRES energy function and the conformational space annealing method of global optimization. The lowest-energy coarse-grained structures are then converted to an all-atom representation and energy-minimized with the ECEPP/3 force field. The procedure was assessed in two recent blind tests of protein-structure prediction. During the first blind test, we predicted large fragments of alpha and alpha+beta proteins [60-70 residues with C(alpha) rms deviation (rmsd) <6 A]. However, for alpha+beta proteins, significant topological errors occurred despite low rmsd values. In the second exercise, we predicted whole structures of five proteins (two alpha and three alpha+beta, with sizes of 53-235 residues) with remarkably good accuracy. In particular, for the genomic target TM0487 (a 102-residue alpha+beta protein from Thermotoga maritima), we predicted the complete, topologically correct structure with 7.3-A C(alpha) rmsd. So far this protein is the largest alpha+beta protein predicted based solely on the amino acid sequence and a physics-based potential-energy function and search procedure. For target T0198, a phosphate transport system regulator PhoU from T. maritima (a 235-residue mainly alpha-helical protein), we predicted the topology of the whole six-helix bundle correctly within 8 A rmsd, except the 32 C-terminal residues, most of which form a beta-hairpin. These and other examples described in this work demonstrate significant progress in physics-based protein-structure prediction.

The protein folding problem: global optimization of the force fields.

The evolutionary development of a theoretical approach to the protein folding problem, in our laboratory, is traced. The theoretical foundations and the development of a suitable empirical all-atom potential energy function and a global optimization search are examined. Whereas the all-atom approach has thus far succeeded for relatively small molecules and for alpha-helical proteins containing up to 46 residues, it has been necessary to develop a hierarchical approach to treat larger proteins. In the hierarchical approach to single- and multiple-chain proteins, global optimization is carried out for a simplified united residue (UNRES) description of a polypeptide chain to locate the region in which the global minimum lies. Conversion of the UNRES structures in this region to all-atom structures is followed by a local search in this region. The performance of this approach in successive CASP blind tests for predicting protein structure by an ab initio physics-based method is described. Finally, a recent attempt to compute a folding pathway is discussed.

The investigation of the effects of counterions in protein dynamics simulations.

Molecular simulations able to exactly represent solvated charged proteins are helpful in understanding protein dynamics, structure and function. In the present study we have used two different starting structures of papain (a typical, stable, globular protein of intermediate net charge) and different modeling procedures to evaluate some effects of counterions in simulations. A number of configurations have been generated and relaxed for each system by various combinations of constrained simulated annealing and molecular dynamics procedures, using the AMBER force field. The analysis of trajectories shows that the simulations of solvated proteins are moderately sensitive to the presence of counterions. However, this sensitivity is highly dependent on the starting model and different procedures of equilibration used. The neutralized systems tend to evince smaller root mean square deviations regardless of the system investigated and the simulation procedure used. The results of parameterized fitting of the simulated structures to the crystallographic data, giving quantitative measure of the total charge influence on the stability of various elements of the secondary structure, revealed a clear scatter of different reactions of various systems' secondary structures to counterions addition: some systems apparently were stabilized when neutralized, while the others were not. Thus, one cannot unequivocally state, despite consideration of specific simulation conditions, whether protein secondary structures are more stable when they have neutralized charges. This suggests that caution should be taken when claiming the stabilizing effect of counterions in simulations other than those involving small, unstable polypeptides or highly charged proteins.

Synthesis, activity on NK-3 tachykinin receptor and conformational solution studies of scyliorhinin II analogs modified at position 16.

Two analogs of a tachykinin family peptides - scyliorhinin II (ScyII): [Aib(16)]ScyII and [Sar(16)]ScyII were synthesized by the solid-phase method using Fmoc chemistry. Conformational studies in water and DMSO-d(6) on these peptides were performed using a combination of two-dimensional NMR and theoretical conformational analysis. The solution structure of the peptides studied is interpreted as an equilibrium of several conformers with different statistical weights. The structure of [Sar(16)]ScyII in water appeared to be more flexible, especially in the C-terminal fragment. A better defined structure for this analog was obtained in DMSO-d(6), in which the analysis resulted in a family of conformers with similar shapes. Some of these conformers were characterized by the presence of a 3(10)-helix in the N-terminal fragment and middle part of the molecule. The introduction of the Aib residue in position 16 significantly rigidifies the structure. For [Aib(16)]ScyII in both solvent systems very similar populations of conformations were obtained which are characterized by the presence of a 3(10)-helix in the 13-18 fragment. A common structural motif was found in conformationally constrained Cys(7)-Cys(13) fragment, which resembles the Greek letter 'omega'. The differences in the solution structure of the C-terminal fragment of the peptides studied are responsible for their specificity. [Aib(16)]ScyII showed 25% the agonistic activity of selective NK-3 agonist - senktide, but it also showed antagonist effect vs. this peptide, whereas [Sar(16)]ScyII appeared to be a full agonist of NK-3 tachykinin receptor.

Determination of conformational equilibrium of peptides in solution by NMR spectroscopy and theoretical conformational analysis: application to the calibration of mean-field solvation models.

Peptides occur in solution as ensembles of conformations rather than in a fixed conformation. The existing energy functions are usually inadequate to predict the conformational equilibrium in solution, because of failure to account properly for solvation, if the solvent is not considered explicitly (which is usually prohibitively expensive). NMR data are therefore widely incorporated into theoretical conformational analysis. Because of conformational flexibility, restrained molecular dynamics (with restraints derived from NMR data), which is usually applied to determine protein conformation is of limited use in the case of peptides. Instead, (a) the restraints are averaged within predefined time windows during molecular dynamics (MD) simulations (time averaging), (b) multiple-copy MD simulations are carried out and the restraints are averaged over the copies (ensemble averaging), or (c) a representative ensemble of sterically feasible conformations is generated and the weights of the conformations are then fitted so that the computed average observables match the experimental data (weight fitting). All these approaches are briefly discussed in this article. If an adequate force field is used, conformations with large statistical weights obtained from the weight-fitting procedure should also have low energies, which can be implemented in force field calibration. Such a procedure is particularly attractive regarding the parameterization of the solvation energy in nonaqueous solvents, e.g., dimethyl sulfoxide, for which thermodynamic solvation data are scarce. A method for calibration of solvation parameters in dimethyl sulfoxide, which is based on this principle was recently proposed by C. Baysal and H. Meirovitch (Journal of the American Chemical Society, 1998, Vol. 120, pp. 800--812), in which the energy gap between the conformations compatible with NMR data and the alternative conformations is maximized. In this work we propose an alternative method based on the principle that the best-fitting statistical weights of conformations should match the Boltzmann weights computed with the force field applied. Preliminary results obtained using three test peptides of varying conformational mobility: H-Ser(1)-Pro(2)-Lys(3)-Leu(4)-OH, Ac-Tyr(1)-D-Phe(2)-Ser(3)-Pro(4)-Lys(5)-Leu(6)-NH(2), and cyclo(Tyr(1)-D-Phe(2)-Ser(3)-Pro(4)-Lys(5)-Leu(6)) are presented.

Fluorescence decay time distribution analysis of cyclic enkephalin analogues; influence of solvent and Leu configuration in position 5 on conformation.

Lifetime distribution analysis were performed to study the influence of Leu configuration in position 5 on changes of the peptide chain of cyclic analogues of enkephalins containing a fluorescence donor and acceptor in different solvents. The configuration change of Leu5 in all the analogues of enkephalins studied which contain donor-acceptor pairs has no apparent influence on Trp lifetime distributions. In contrast, there is a significant solvent effect on the shape of lifetime distribution.

Influence of solvent and configuration of residues at positions 2 and 3 on distance and mobility of pharmacophore groups at positions 1 and 4 in cyclic enkephalin analogues.

The analgesic activity of opioid peptides is mainly connected with their affinity and selectivity for the mu-receptors. The biological activity of cyclic opioid analogues depends on mutual orientation and conformational freedom of aromatic pharmacophore groups at positions 1 and 4. The distance and distance distributions between chromophores at positions 1 [Phe(p-NO(2)), p-nitrophenylalanine] and 4 [Nal, beta-(2-naphthyl)alanine], which constitute an energy donor-acceptor pair, were calculated based on measured fluorescence intensity decays of a donor (Nal). The influence of the solvent and configuration of the residues at position 2 and 3 on donor-acceptor distance distribution and mobility of pharmacophore groups at position 1 and 4 in cyclic enkephalin analogues are discussed.

Influence of solvents and leucine configuration at position 5 on tryptophan fluorescence in cyclic enkephalin analogues.

The fluorescence decay of tryptophan is a sensitive indicator of its local environment within a peptide or protein. In this study we carried out fluorescence measurements of the tryptophan residue of cyclic enkephalin analogues of a general formula X-c[D-Dab(2)-Gly(3)-Trp(4)-Y(5)] where X = Cbz or H and Y = D- or L-Leu, in four solvents [water, methanol, acetonitrile, and dimethyl sulfoxide (DMSO)]. An analysis of the tryptophan fluorescence decays using a discrete-exponential model indicates that tryptophan fluorescence decay can be described by a double exponential function in all solvents studied. Lifetime distribution analysis yields a bimodal distribution in protic solvents (water and methanol), whereas an asymmetric, unimodal distribution in an aprotic solvent (DMSO) and uni- or bimodal distributions in acetonitrile solution, depending on leucine configuration. The data are interpreted in terms of the rotamer model, in which the modality and the relative proportions of the lifetime components are related to the population distribution of tryptophan chi(1) rotamers about the C(alpha)--C(beta) bond. The chirality of the Leu(5) residue and solvent properties affect the local environment of the tryptophan residue and therefore influence the distribution of side-chain rotamers. These results are consistent with the results of theoretical conformational calculations.

Recent improvements in prediction of protein structure by global optimization of a potential energy function.

Recent improvements of a hierarchical ab initio or de novo approach for predicting both alpha and beta structures of proteins are described. The united-residue energy function used in this procedure includes multibody interactions from a cumulant expansion of the free energy of polypeptide chains, with their relative weights determined by Z-score optimization. The critical initial stage of the hierarchical procedure involves a search of conformational space by the conformational space annealing (CSA) method, followed by optimization of an all-atom model. The procedure was assessed in a recent blind test of protein structure prediction (CASP4). The resulting lowest-energy structures of the target proteins (ranging in size from 70 to 244 residues) agreed with the experimental structures in many respects. The entire experimental structure of a cyclic alpha-helical protein of 70 residues was predicted to within 4.3 A alpha-carbon (C(alpha)) rms deviation (rmsd) whereas, for other alpha-helical proteins, fragments of roughly 60 residues were predicted to within 6.0 A C(alpha) rmsd. Whereas beta structures can now be predicted with the new procedure, the success rate for alpha/beta- and beta-proteins is lower than that for alpha-proteins at present. For the beta portions of alpha/beta structures, the C(alpha) rmsd's are less than 6.0 A for contiguous fragments of 30-40 residues; for one target, three fragments (of length 10, 23, and 28 residues, respectively) formed a compact part of the tertiary structure with a C(alpha) rmsd less than 6.0 A. Overall, these results constitute an important step toward the ab initio prediction of protein structure solely from the amino acid sequence.

Solution conformational study of Scyliorhinin I analogues with conformational constraints by two-dimensional NMR and theoretical conformational analysis.

Two analogues of Scyliorhinin I (Scyl), a tachykinin with N-MeLeu in position 8 and a 1,5-disubstituted tetrazole ring between positions 7 and 8, introduced in order to generate local conformational constraints, were synthesized using the solid-phase method. Conformational studies in water and DMSO-d6 were performed on these peptides using a combination of the two-dimensional NMR technique and theoretical conformational analysis. The algorithm of conformational search consisted of the following three stages: (i) extensive global conformational analysis in order to find all low-energy conformations; (ii) calculation of the NOE effects and vicinal coupling constants for each of the low energy conformations; (iii) determining the statistical weights of these conformations by means of a nonlinear least-squares procedure, in order to obtain the best fit of the averaged simulated spectrum to the experimental one. In both solvents the three-dimensional structure of the analogues studied can be interpreted only in terms of an ensemble of multiple conformations. For [MeLeu8]Scyl, the C-terminal 6-10 fragment adopts more rigid structure than the N-terminal one. In the case of the analogue with the tetrazole ring in DMSO-d6 the three-dimensional structure is characterized by two dominant conformers with similar geometry of their backbones. They superimpose especially well (RMSD = 0.28 A) in the 6-9 fragments. All conformers calculated in both solvents superimpose in their C-terminal fragments much better than those of the first analogue. The results obtained indicate that the introduction of the tetrazole ring into the Scyl molecule rigidifies its structure significantly more than that of MeLeu.

Molecular simulation study of cooperativity in hydrophobic association.

To investigate the cooperativity of hydrophobic interactions, the potential of mean force of two- and three-molecule methane clusters in water was determined by molecular dynamics simulations using two methods: umbrella-sampling with the weighted histogram analysis method and thermodynamic integration. Two water models, TIP3P and TIP4P, were used, while each methane molecule was modeled as a united atom. It was found that the three-body potential of mean force is not additive, i.e., it cannot be calculated as a sum of two-body contributions, but requires an additional three-body cooperative term. The cooperative term, which amounts to only about 10% of the total hydrophobic association free energy, was found to increase the strength of hydrophobic association; this finding differs from the results of earlier Monte Carlo studies with the free energy perturbation method of Rank and Baker (1997). As in the work of Rank and Baker, the solvent contribution to the potential of mean force was found to be well approximated by the molecular surface of two methane molecules. Moreover, we also found that the cooperative term is well represented by the difference between the molecular surface of the three-methane cluster and those of all three pairs of methane molecules. In addition, it was found that, while there is a cooperative contribution to the hydrophobic association free energy albeit a small one, the errors associated with the use of pairwise potentials are comparable to or larger than this contribution.

Calculation of protein conformation by global optimization of a potential energy function.

A novel hierarchical approach to protein folding has been applied to compute the unknown structures of seven target proteins provided by CASP3. The approach is based exclusively on the global optimization of a potential energy function for a united-residue model by conformational space annealing, followed by energy refinement using an all-atom potential. Comparison of the submitted models for five globular proteins with the experimental structures shows that the conformations of large fragments (approximately 60 aa) were predicted with rmsds of 4.2-6.8 A for the C alpha atoms. Our lowest-energy models for targets T0056 and T0061 were particularly successful, producing the correct fold of approximately 52% and 80% of the structures, respectively. These results support the thermodynamic hypothesis that protein structure can be computed solely by global optimization of a potential energy function for a given amino acid sequence.

Protein structure prediction by global optimization of a potential energy function.

An approach based exclusively on finding the global minimum of an appropriate potential energy function has been used to predict the unknown structures of five globular proteins with sizes ranging from 89 to 140 amino acid residues. Comparison of the computed lowest-energy structures of two of them (HDEA and MarA) with the crystal structures, released by the Protein Data Bank after the predictions were made, shows that large fragments (61 residues) of both proteins were predicted with rms deviations of 4.2 and 6.0 A for the Calpha atoms, for HDEA and MarA, respectively. This represents 80% and 53% of the observed structures of HDEA and MarA, respectively. Similar rms deviations were obtained for approximately 60-residue fragments of the other three proteins. These results constitute an important step toward the prediction of protein structure based solely on global optimization of a potential energy function for a given amino acid sequence.

Energy-based de novo protein folding by conformational space annealing and an off-lattice united-residue force field: application to the 10-55 fragment of staphylococcal protein A and to apo calbindin D9K.

The conformational space annealing (CSA) method for global optimization has been applied to the 10-55 fragment of the B-domain of staphylococcal protein A (protein A) and to a 75-residue protein, apo calbindin D9K (PDB ID code), by using the UNRES off-lattice united-residue force field. Although the potential was not calibrated with these two proteins, the native-like structures were found among the low-energy conformations, without the use of threading or secondary-structure predictions. This is because the CSA method can find many distinct families of low-energy conformations. Starting from random conformations, the CSA method found that there are two families of low-energy conformations for each of the two proteins, the native-like fold and its mirror image. The CSA method converged to the same low-energy folds in all cases studied, as opposed to other optimization methods. It appears that the CSA method with the UNRES force field, which is based on the thermodynamic hypothesis, can be used in prediction of protein structures in real time.

Maximum entropy approach to the determination of solution conformation of flexible polypeptides by global conformational analysis and NMR spectroscopy--application to DNS1-c-[D-A2,bu2,Trp4,Leu5]enkephalin and DNS1-c-[D-A2bu2,Trp4,D-Leu5]enkephalin.

A method is proposed to determine the conformational equilibrium of flexible polypeptides in solution, using the data provided by NMR spectroscopy and theoretical conformational calculations. The algorithm consists of the following three steps: (i) search of the conformational space in order to find conformations with reasonably low energy; (ii) simulation of the NOE spectrum and vicinal coupling constants for each of the low energy conformations; and (iii) determining the statistical weights of the conformations, by means of the maximum-entropy method, in order to obtain the best fit of the averaged NOE intensities and coupling constants to the experimental quantities. The method has been applied to two cyclic enkephalin analogs: DNS1-c-[D-A2bu2,Trp4,Leu5]enkephalin (ENKL) and DNS1-c-[D-A2bu2,Trp4,D-Leu5]enkephalin (ENKD). NMR measurements were carried out in deuterated dimethyl sulfoxide. Two techniques were used in conformational search: the electrostatically driven Monte Carlo method (EDMC), which results in extensive search of the conformational space, but gives only energy minima, and the molecular dynamics method (MD), which results in a more accurate, but also more confined search. In the case of EDMC calculations, conformational energy was evaluated using the ECEPP/3 force field augmented with the SRFOPT solvation-shell model, while in the case of MD the AMBER force field was used with explicit solvent molecules. Both searches and subsequent fitting of conformational weights to NMR data resulted in similar conformations of the cyclic part of the peptides studied. For both ENKL and ENKD a common feature of the low-energy solution conformations is the presence of a type II' or type IV beta-turn at residues 3 and 4; the ECEPP/3 force field also gives a remarkable content of type III beta-turn. These beta-turns are tighter in the case of ENKL, which is reflected in different distributions of the D-A2bu(N gamma H)...D-A2bu(CO) and D-A2bu(N gamma H)...Gly3(CO) hydrogen-bonding distances, indicating that the D-A2bu(N gamma H) amide proton is more shielded from the solvent than in the case of ENKD. This finding conforms with the results of temperature coefficient data of the D-A2bu(N gamma H) proton. It has also been found that direct (MD) or Boltzmann (EDMC) averages of the observables do not exactly conform with the measured values, even when explicit solvent molecules are included. This suggests that improving force-field parameters might be necessary in order to obtain reliable conformational ensembles in computer simulations, without the aid of experimental data.

Fluorescence decay time distribution analysis of cyclic enkephalin analogues. Influence of the solvents and configuration of amino acids in position 2 and 3 on changes in conformation.

The lifetime distribution calculations were applied to study the influence of configuration of amino-acid residues in positions 2 and 3 on changes in conformation of the peptide chain of cyclic analogues of enkephalins containing a fluorescence energy donor and acceptor in different solvents. In all the solvents studied the lifetime distributions were bimodal. This testified to the presence of two families of conformations. In this paper the relationship between the population of each conformation and configuration of the residues in position 2 and 3, and the solvent used is discussed.

New developments of the electrostatically driven Monte Carlo method: test on the membrane-bound portion of melittin.

The electrostatically driven Monte Carlo (EDMC) method has been greatly improved by adding a series of new features, including a procedure for cluster analysis of the accepted conformations. This information is used to guide the search for the global energy minimum. Alternative procedures for generating perturbed conformations to sample the conformational space were also included. These procedures enhance the efficiency of the method by generating a larger number of low-energy conformations. The improved EDMC method has been used to explore the conformational space of a 20-residue polypeptide chain whose sequence corresponds to the membrane-bound portion of melittin. The ECEPP/3 (Empirical Conformational Energy Program for Peptides) algorithm was used to describe the conformational energy of the chain. After an exhaustive search involving 14 independent runs, the lowest energy conformation (LEC) (-91.0 kcal/mol) of the entire study was encountered in four of the runs, while conformations higher in energy by no more than 1.8 kcal/mol were found in the remaining runs with the exception of one of them (run 8). The LEC is identical to the conformation found recently by J. Lee, H. A. Scheraga, and S. Rackovsky [(1988) "Conformational Analysis of the 20-Residue Membrane-Bound Portion of Melittin by Conformational Space Annealing," Biopolymers, Vol. 46, pp. 103-115] as the lowest energy conformation obtained in their study using the conformational space annealing method. These results suggest that this conformation corresponds to the global energy minimum of the ECEPP/3 potential function for this specific sequence: it also appears to be the conformation of lowest free energy.

Design of a knowledge-based force field for off-lattice simulations of protein structure.

Prediction of protein structure from amino-acid sequence still continues to be an unsolved problem of theoretical molecular biology. One approach to solve it is to construct an appropriate (free) energy function that recognizes the native structures of some selected proteins (whose native structures are known) as the ones distinctively lowest in (free) energy and then to carry out a search of the lowest-energy structure of a new protein. In order to reduce the complexity of the problem and the cost of energy evaluation, the so-called united-residue representation of the polypeptide chain is often applied, in which each amino-acid residue is represented by only a few interaction sites. Once the global energy minimum of the simplified chain has been found, the all-atom structure can easily and reliably be constructed. The search of the lowest-energy structure is usually carried out by means of Monte Carlo methods, though use of more efficient global-optimization methods, especially those of deformation of original energy surface is potentially promising. Monte Carlo search of the conformational space can be accelerated greatly, if the chain is superposed on a discrete lattice (the on-lattice approach). On the other hand, the on-lattice approach prohibits the use of many efficient global-optimization methods, because they require both energy and its space derivatives. The on-lattice methods in which the chain is embedded in the continuous 3D space are, therefore, also worth developing. In this paper we summarize the work on the design and implementation of an off-lattice united-residue force field that is underway in our group, in cooperation with Professor HA. Scheraga of Cornell University, U.S.A.