Revised Backbone-Virtual-Bond-Angle Potentials to Treat the l- and d-Amino Acid Residues in the Coarse-Grained United Residue (UNRES) Force Field.

Continuing our effort to introduce d-amino-acid residues in the united residue (UNRES) force field developed in our laboratory, in this work the Cα ··· Cα ··· Cα backbone-virtual-bond-valence-angle (θ) potentials for systems containing d-amino-acid residues have been developed. The potentials were determined by integrating the combined energy surfaces of all possible triplets of terminally blocked glycine, alanine, and proline obtained with ab initio molecular quantum mechanics at the MP2/6-31G(d,p) level to calculate the corresponding potentials of mean force (PMFs). Subsequently, analytical expressions were fitted to the PMFs to give the virtual-bond-valence potentials to be used in UNRES. Alanine represented all types of amino-acid residues except glycine and proline. The blocking groups were either the N-acetyl and N',N'-dimethyl or N-acetyl and pyrrolidyl group, depending on whether the residue next in sequence was an alanine-type or a proline residue. A total of 126 potentials (63 symmetry-unrelated potentials for each set of terminally blocking groups) were determined. Together with the torsional, double-torsional, and side-chain-rotamer potentials for polypeptide chains containing d-amino-acid residues determined in our earlier work (Sieradzan et al. J. Chem. Theory Comput., 2012, 8, 4746), the new virtual-bond-angle (θ) potentials now constitute the complete set of physics-based potentials with which to run coarse-grained simulations of systems containing d-amino-acid residues. The ability of the extended UNRES force field to reproduce thermodynamics of polypeptide systems with d-amino-acid residues was tested by comparing the experimentally measured and the calculated free energies of helix formation of model KLALKLALxxLKLALKLA peptides, where x denotes any d- or l- amino-acid residue. The obtained results demonstrate that the UNRES force field with the new potentials reproduce the changes of free energies of helix formation upon d-substitution but overestimate the free energies of helix formation. To test the ability of UNRES with the new potentials to reproduce the structures of polypeptides with d-amino-acid residues, an ab initio replica-exchange folding simulation of thurincin H from Bacillus thuringiensis, which has d-amino-acid residues in the sequence, was carried out. UNRES was able to locate the native α-helical hairpin structure as the dominant structure even though no native sulfide-carbon bonds were present in the simulation.

Introduction

d-amino acids are currently of great interest in drug design.[1−3]d-amino acid substitution can increase protein resistance to proteolysis[4] and can enhance the thermodynamic stability[5] of proteins. This idea has been used in designing CD4 analogues as potential HIV binding domain inhibitors.[1] Moreover, d-amino acids can increase specificity of, e.g., d-tyrosine conjugated with naproxen.[2]

Despite increases in computing power, all-atom molecular dynamic simulations are still limited to either small system size[6] or short time scale.[7] Therefore, coarse-grained models of proteins continue to play a great role in protein simulations,[8,9] including computer-aided drug design.[10] The coarse-grained treatment speeds up the time scale of simulations by 3–4 orders of magnitude compared to all-atom treatment with explicit solvent.[11,12] However, most of the coarse-grained force fields that are capable of predicting the structures of proteins are knowledge-based and, naturally, treat l-amino-acid residues rather than d-amino acid residues, because of an insufficient number of d-amino-acid residues in protein structural databases to derive the respective statistical potentials. By contrast, the united residue (UNRES) model and force field of polypeptide chains,[13−25] which we have been developing for over two decades in our laboratory, is physics-based and can easily be extended to include d-amino-acid residues, because it does not depend on structural databases. The capability of UNRES to predict protein structures is being continuously assessed, since 1998, in the consecutive community wide experiments of the Critical Assessment of Techniques for Protein Structure Prediction (CASP) with very good results.[26−31] In particular, in the recent CASP10 exercise, UNRES, as one of only two methods, predicted correct domain packing of target T0663.[31]

In our previous work,[32] we introduced a general treatment of d-amino-acid residues in UNRES. We demonstrated that the potentials that depend on single amino-acid residue types as far as local interactions are concerned, namely, the side-chain rotamer and third-order backbone-local–backbone-electrostatic correlation potentials can be obtained from the respective alanine-type l-residue potentials by applying the inversion operation. Conversely, we found that only approximate symmetry relations exist between the local potentials that depend on pairs (the torsional potentials) or triplets (double-torsional and virtual-bond-angle-potentials) of consecutive amino-acid residues. In our previous work,[32] we determined the physics-based torsional and double-torsional potentials that involve l- and d-amino-acid residues and applied approximate symmetry relations to estimate the virtual-bond-angle potentials. With these modifications, the UNRES model enabled us to simulate folding pathways of the BBAT1 protein[33] and stability of gramicidin D.[32] The determination of the exact virtual-bond-angle potentials for polypeptides and proteins containing d-amino-acid residues had then been postponed to further work because of substantial effort required to accomplish this task.

In the present work, using the ab initio potential-energy surfaces of terminally blocked glycine, alanine, and proline determined in our previous work[32] to calculate the respective potentials of mean force, we determine the physics-based virtual-bond-angle potentials for backbone sections that include d-amino-acid residues. This work concludes the extension of UNRES to polypeptides and proteins with d-amino-acid residues, which have not been determined in our earlier work.[32] We use the harmonic-extrapolation method developed in our first paper on the determination of virtual-bond-angle potentials,[20] where we treated only l-amino-acid residues. Moreover, because of limited computational resources available at that time, in that work, energy minimization was carried out with the RHF scheme with single-point MP2 correction, and the N,N′-dimethyl group was used to mimic the link of a residue to the proline residue at the C-terminus. Therefore, in this work, we also revised the virtual-bond-valence angle potentials for systems composed of l-amino-acid residues that were determined in ref (20). We demonstrate the performance of the modified UNRES force field with the examples of thermodynamics of helix formation in model peptides and folding of thurincin H, a natural peptide that contains d-amino-acid residues.

Methods

UNRES Model of the Polypeptide Chain

In the UNRES model, a polypeptide chain is represented as a sequence of α-carbon (Cα) atoms with attached united side chains (SC’s) and united peptide groups (p’s) positioned halfway between two consecutive Cα’s. Only the united side chains and united peptide groups act as interaction sites, while the Cα atoms assist only in the definition of geometry (Figure 1). The effective energy function is expressed by eq 1. This effective energy function originates from the restricted free energy (RFE) or the potential of mean force (PMF) of the chain constrained to a given coarse-grained conformation along with the surrounding solvent.[16,17,21] The PMF is then expanded into factors,[17] which are the cluster-cumulant functions introduced by Kubo.[34] Factors of order m are PMF components that arise from m groups of atomic-detailed interactions, each of which comprises interactions between the atoms that belong to a given coarse-grained site or to two coarse-grained sites. For example, the side-chain–side-chain interactions are factors of order 1, whereas the torsional potentials are factors of order 2.where the U′s are energy terms, θ is the backbone virtual-bond angle, γ is the backbone virtual-bond-dihedral angle, α and β are the angles defining the location of the center of the united side chain of residue i (Figure 1), and d is the length of the ith virtual bond, which is either a Cα ··· Cα virtual bond or Cα ··· SC virtual bond; τ(1) is the SC ··· Cα ··· Cα ··· Cα torsional angle, τ(2) is the Cα ··· Cα ··· Cα ··· SC torsional angle, and τ(3) is the SC ··· Cα ··· Cα ··· SC torsional angle.[35] Each energy term is multiplied by an appropriate weight, w (where x stands for the interactions contained in the respective energy term); these quantities are determined by force-field calibration[19,36] (see the end of this section for more information). The terms corresponding to factors of order higher than 1 in the cluster-cumulant expansion of the PMF are additionally multiplied by the respective temperature factors that were introduced in our earlier work[19] and which reflect the temperature dependence of the first generalized-cumulant term in those factors, as discussed in refs (19) and (37). The factors f are defined by eq 2.where To = 300 K.

Figure 1

UNRES model of polypeptide chains. The interaction sites are peptide-group centers (p), and side-chain centers (SC) attached to the corresponding α-carbons with different Cα ··· SC bond lengths, dSC. The peptide groups are represented as dark gray circles and the side chains are represented as light gray ellipsoids of different size. The α-carbon atoms are represented by small open circles. The geometry of the chain can be described either by the virtual-bond vectors dC (from Cα to Cα), i = 1,2,...,n – 1, and dX (from Cα to SC), i = 2,...,n – 1, represented by thick lines, where n is the number of residues, or in terms of virtual-bond lengths, backbone virtual-bond angles θ, i = 1,2,...,n – 2, backbone virtual-bond-dihedral angles γ, i = 1,2,...,n – 3, and the angles α and β,i = 2,3,...,n – 1 that describe the location of a side chain with respect to the coordinate frame defined by Cα, Cα, and Cα, Cα.

The term USC represents the mean free energy of the hydrophobic (hydrophilic) interactions between the side chains, which implicitly contains the contributions from the interactions of the side chain with the solvent. The term USC denotes the excluded-volume potential of the side-chain–peptide-group interactions. The peptide-group interaction potential is split into two parts: the Lennard-Jones interaction energy between peptide-group centers (UpVDW) and the average electrostatic energy between peptide-group dipoles (Upel); the second of these terms accounts for the tendency to form backbone hydrogen bonds between peptide groups p and p. The terms Utor, Utord, Ub, Urot, and Ubond are the virtual-bond-dihedral angle torsional terms, virtual-bond dihedral angle double-torsional terms, virtual-bond angle bending terms, side-chain rotamer, and virtual-bond-deformation terms; these terms account for the local properties of the polypeptide chain. For Utor, Utord, Ub, three types of residues are distinguished: alanine, glycine, and proline where alanine describes all amino acid types except glycine and proline. The terms Ucorr( represent correlation or multibody contributions from the coupling between backbone-local and backbone-electrostatic interactions, and the terms Uturn( are correlation contributions involving m consecutive peptide groups; they are, therefore, termed turn contributions. The multibody terms are indispensable for reproduction of regular α-helical and β-sheet structures.[16,17,38]USC–corr are correlation potentials between local conformations of side-chain and backbone. These potentials significantly improve the description of loop structures.[35]

The energy-term weights are determined by force-field calibration to reproduce the structure and folding thermodynamics of selected training proteins.[19,36] In this work, we used the energy-term weights from the force field calibrated in our earlier work[36] with the tryptophan cage (PDB code: 1L2Y)[39] and the tryptophan zipper 2 (PDB code: 1LE1).[40] It should be noted that the force field of our previous work[36] was calibrated with peptides that contain only l-amino-acid residues. The hierarchical optimization procedure was used[19,36] whose goal is to achieve a free-energy landscape in which the free energy decreases with increasing the native likeness below the folding-trasition temperature and free energy increases with increasing the native likeness above the folding-transition temperature.

Determination of the New Ub Potentials

The procedure developed in our earlier work[20] was used to determine new virtual Cα ··· Cα ··· Cα valence angle θ (Figure 2) bending potentials (Ub). In that procedure, the respective potentials of mean force (PMFs) were first calculated from nonadiabatic energy maps of terminally blocked amino-acid residues, and then a three-dimensional Fourier series was fitted to them. Energy maps had been computed[20] in the angles of rotation λ(1) and λ(2) about the virtual-bond Cα ··· Cα axes; these angles have been defined in ref (41) (Figure 3). Energy was minimized with respect to all degrees of freedom except for λ(1) and λ(2). Therefore, energy maps are minimized in the θ angle.

Figure 2

Illustration of the model terminally blocked tripeptides constructed to compute the integrals of eqs 3. Each X, Y, and Z denotes side-chains of l-Ala, d-Ala, Gly, l-Pro, or d-Pro.

Figure 3

Definition of the dihedral angles λ(1) and λ(2) for rotation of a peptide unit about the Cα ··· Cα virtual bonds of a peptide unit.[41]

To calculate the PMF on a θ grid, harmonic extrapolation has to be used,[20] as expressed by eq 3 (which has the same form as eq 15 in ref (20)); it should be noted that eq 3 includes a harmonic-entropy contribution.where Ub,XYZ denotes the bending potentials of mean force, X, Y, and Z denote l-alanine l-proline, d-alanine, d-proline, or glycine, where alanine represents all amino acids except proline or glycine, β = (RT)−1, T is the absolute temperature, R is the universal gas constant, x′, y′, and z′ are the vectors of variables other than θ, λ(1), and λ(2) of residues X, Y, and Z, respectively. V is the volume of the space spanned by the secondary variables pertaining to the second residue (Y), contained in the y′ vector, H denotes the Hessian matrix (matrix of the second derivatives of energy in geometrical variables); the terms with H contain the harmonic-entropy contributions, and the asterisks indicate the values corresponding to the points on the nonadiabatic energy maps. λ1 – λ4, shown in Figure 3, were determined from eqs 4a–f:where the angles λ(1) and λ(2) are shown in Figure 3.

The second multiple integral in eq 3 corresponds to a sum of torsional and double-torsional potentials and must be subtracted from the PMF due to local interactions in tripeptide systems sketched in Figure 2, consistent with the definition of the PMF factors in UNRES,[17] because the torsional and double-torsional potentials constitute separate UNRES energy terms. These energy terms were determined for systems that include d-amino-acid residues in our earlier work,[32] and here in eq 3.

After determination of Ub,XYZ, a three-dimensional Fourier series (eq 5) was used to fit the free energy surfaces. As variables, γ1, γ2, and θ/2 were used. The use of θ/2 instead of θ is justified by the fact that the squares of interatomic distances within a peptide unit are naturally expressed in terms of powers of cos(θ/2) and sin(θ/2) as shown in our previous work.[20]

In our previous work,[20] we computed the energy maps of terminally blocked amino-acid residues by using ab initio molecular quantum mechanics with the 6-31G(d,p) basis set. Because less powerful resources were available at that time, the energy was minimized in the Restricted Hartree–Fock (RHF) scheme and, subsequently, the single-point Møller–Plesset (MP2) energy correction was computed for the energy-minimized conformation. Previously, the limited resources enabled us to calculate Hessian matrices only in the RHF scheme. Moreover, interfacing a residue to the following proline residue was modeled by blocking it with an N′,N′-dimethyl group. In this work, all energy maps were recalculated by minimizing the energy in the MP2/6-31G(d,p) scheme; additionally, we used the C-terminal pyrrolidine group and not the N′,N′-dimethylamine blocking group for residues preceding proline. Consequently, all diabatic energy maps of trans-N-acetyl-alanyl-N′-methylamide, trans-N-acetyl-alanyl-N′-pyrrolidylamide, trans-N-acetyl-glycyl-N′-methylamide, trans-N-acetyl-glycyl-N′-pyrrolidylamide trans-N-acetyl-prolyl-N′-methylamide, and trans-N-acetyl-prolyl-N′-pyrrolidylamide were computed. The angles λ(1) and λ(2) (Figure 3) were used as grid variables; both varied from −180° to +180° with a 15° step. The energy was minimized with respect to all other degrees of freedom. The Hessians (second-derivative matrices or force-constant matrices) were also calculated in the MP2 scheme at every grid point to determine the harmonic-entropy contribution and harmonic approximation to the PMF value (eq 3).

Thermodynamic Stability of α-Helix Formation in KLALKLALxxLKLALKLA Peptides

In our previous work[32,33] and also in this work (section 3.3), we assessed the ability of the modified UNRES force field to reproduce the structure of proteins with d-amino-acid residues. However, it is equally important to assess if the force field can distinguish different chirality by energy. To assess whether the new potentials can reproduce the change of the energetics of local interactions upon d-substitution, we studied a set of 39 peptides KLALKLALxxLKLALKLA where xx denote two consecutive natural amino-acid residues.[42] All coded l-amino-acid residues and their d-amino-acid counterparts as well as glycine were considered. For these peptides, the free-energy differences related to the formation of α-helical structure were determined from molar ellipticities at λ = 222 nm and reversed-phase high pressure liquid chromatography.[42] The experimental values of free energies could thus be compared with the results of our simulations. To carry out simulations, multiplexed replica-exchange molecular dynamics (MREMD), as implemented in version 3.2 of UNRES (available at www.unres.pl),[12,43] was used. The Berendsen thermostat[44] was used to control temperature. The equations of motions were integrated by using the variable time step (VTS) algorithm[11] as implemented in the UNRES program. Each simulation was performed at the following 16 temperatures with 2 copies at each temperature: 250, 260, 270, 280, 290, 300, 310, 320, 330, 340, 350, 360, 370, 380, 390, and 400 K. Each simulation lasted about 50 ns (approximately 10 million steps, with a step length of 4.89 fs). The weighted histogram analysis method (WHAM)[19,45] was used to evaluate the thermodynamic quantities and conformation-dependent ensemble averages, as implemented in our earlier work.[19]

The Gibbs free-energy changes (ΔΔG) at T = 300 K were calculated from the fractions of amino-acid residue in the α-helical state at this temperature, averaged (using WHAM) over conformations. For each conformation, the fraction of residues in the α-helical state was evaluated using the following two methods. In the first approach, the coarse-grained UNRES structure was converted to all-atom representation with the use of our energy-based method, based on the optimal alignment of peptide-group dipoles to reconstruct the backbone[46] and optimal side-chain packing to reconstruct all-atom side chains from the positions of side-chain centroids.[47] A residue was considered to be in the α-helical state if its ϕ and ψ angles were in the A region of the Ramachandran map, as defined by Zimmerman et al.[48] (ϕ = −75° ± 40° and ψ = −40° ± 35°). In the second approach, a residue was considered to be in the α-helical state if the peptide group preceding it formed a hydrogen-bonding contact with the third succeeding peptide group; the presence of a hydrogen-bonding contact was assessed based on the mean-field energy of interactions, which depends on the distance between the centers of the two peptide groups; the details of this method are described in ref (49). The free energies of α-helical-structure formation, ΔGhel(x), were calculated from eq 6. For comparison with experimental values,[42] we calculated the Gibbs free energies of α-helix formation expressed relative to glycine, ΔΔGhel(Gly → x), where x = d-x or l-x, and the differences between the free energies of the d- and l-amino acid substituted host–guest peptides, ΔΔGhel(l-x → d-x). These quantities are expressed by eqs 7 and 8, respectively.where R is the universal gas constant, T is the absolute temperature, and fhel(x) is the ensemble-averaged fraction of α-helical structures in the ensemble for the host–guest peptide containing a pair of specific residues x. The temperature T was set at 298 K, as in the experiment.[42]

Simulation of Folding Pathway of Thurincin H

To evaluate the predictive power of the UNRES force field extended to treat d-amino-acid residues, a small peptide thurincin H[50] (PDB code: 2LBZ), which contains d-allo-threonine residues, was selected. The sequence of this peptide is l-Asp-l-Trp-l-Thr-l-Cys-l-Trp-l-Ser-l-Cys-l-Leu-l-Val-l-Cys-l-Ala-l-Ala-l-Cys-l-Ser-l-Val-l-Glu-l-Leu-l-Leu-d-Asn-l-Leu-l-Val-d-allo-Thr-l-Ala-l-Ala-d-allo-Thr-Gly-l-Ala-d-Ser-l-Thr-l-Ala-l-Ser. Thurincin H contains the following four sulfur–Cα bonds S(Cys4)–Cα(d-Ser28), S(Cys7)–Cα(d-allo-Thr25), S(Cys10)–Cα(d-allo-Thr22), and S(Cys13)–Cα(d-Asn19), where the sulfur is bonded to the Cβ atom of Cys and to the Cα atom of the d-residue. Such sulfide bonds are rare in natural peptides and proteins, as opposed to the frequently occurring disulfide bonds, which makes this target even more interesting to study. The simulation was started from the all-extended conformation, with all virtual-bond-dihedral angles γ set at 180°. The UNRES version and technical details were those described in section 2.3 except that the duration of the simulation was 300 ns (about 60 million steps) to ensure the convergence of the heat-capacity curve.

Because the difference between allo-threonine and threonine is too small at the coarse-grained level to warrant the development of separate side-chain-interaction parameters for allo-threonine for UNRES, d-threonine was used in the simulation instead of d-allo-threonine. Although thurincin H contains four sulfur–Cα bonds, no restraints were imposed on the S–Cα distances in the simulation and, consequently, no information about the structure of the peptide was input. 2000 last snapshots from each trajectory (snapshot frequency being 20 000 steps) were analyzed as described in our earlier work.[19] First, WHAM was run to compute the heat-capacity curve and the weights of conformations at the selected temperature, Ta, (10 K below the heat-capacity maximum; in this study Ta = 300 K), then cluster analysis with the use of Ward’s minimum variance method[51] was carried out at Ta to partition the set of conformations into families and the families were ranked according to decreasing fraction in the ensemble.

Results and discussion

Virtual-Bond-Angle Potentials

A total of 125 virtual-valence-angle-bending potentials, of which 63 are not related by symmetry to each other, were determined for systems with each of the C-terminally blocking group (NHMe or Pir). These numbers are obtained as follows. In each position of the XYZ triplet, Gly, l-Ala, d-Ala, l-Pro or d-Pro can appear, which provides 5 possibilities. This gives 5 × 5 × 5 = 125 possible PMF surfaces. Of these, one corresponds to the Gly-Gly-Gly triplet, the other 124 involve at least one nonglycine residue. Those 124 consist of symmetry-related pairs obtained by flipping the chirality of all residues involved and, consequently, the corresponding pairs of PMF surfaces are related by the change of the signs of both the γ1 and γ2 angles. Thus, there are 124/2 + 1 = 63 symmetry-unrelated PMF surfaces.

Maps of the virtual-bond-angle potentials in the γ2 and θ angles (see Figure 2 for the definition of the θ, γ1, and γ2 angles) are displayed in Figures 4A–L and 5A–L, respectively, for the following 8 selected systems: l-Ala-l-Ala-l-Ala-NHMe (Figure 4A–C), l-Ala-l-Ala-d-Ala-NHMe (Figure 4D–F), l-Ala-d-Ala-d-Ala-NHMe (Figure 4G–I), l-Ala-d-Ala-l-Ala-NHMe (Figure 4J–L), l-Ala-l-Ala-l-Ala-Pir (Figure 5A–C), l-Ala-l-Ala-d-Ala-Pir (Figure 5D–F), l-Ala-d-Ala-d-Ala-Pir (Figure 5G–I), and l-Ala-d-Ala-l-Ala-Pir (Figure 5J–L), and for the following selected values of γ1: γ1 = 60° (Figures 4A,D,G,J and 5A,D,G,J), γ1 = 180° (Figures 4B,E,H,K and 5B,E,H,K) and γ1 = 60° (Figures 4C,F,I,L and 5C,F,I,L). The angle γ2 and not γ1 was chosen as the second variable because the U′s depend on this angle more strongly than on γ1, as already concluded by Levitt[52] based on a statistical analysis of protein structures. Figure 4 shows all symmetry-unrelated PMF surfaces for alanine-type residues. As already mentioned in this section, the PMF surfaces for the other alanine-type triplets can be obtained from those by inversion in the γ1 and γ2 angles. It can be seen that two broad regions of minima (blue color) appear in the PMF surfaces of all 8 model systems shown in Figures 4A–L and 5A–L. The first region occurs around θ = 90°, a value characteristic of α-helical or turn structures, and the second one is centered around θ = 140°, a value characteristic of extended or β-sheet structures. The sizes of the θ = 90° region are always comparable regardless of the values of γ1, whereas the θ = 140° region is the largest for γ1 = 60° (or γ1 = −60° for the symmetry-related surfaces). For γ1 = 180°, the two regions nearly merge and there is only a small free-energy barrier between them. For γ1 = 60°, there is a noticeable barrier of ≈2 kcal/mol. For example, for a potential of mean force for U (Figure 4A–C and 5A–C), when γ1 = 180° and θ = 140°, the free-energy difference between γ2 = 0° and γ2 = 75° is ≈8 kcal/mol.

Figure 4

Sample contour plots of the valence bond bending potentials U(θ, γ(1), γ(2)) as functions of θ and γ(2) angles for alanine-type tripeptides with NHMe terminal groups, where (d,l)-Ala indicates a d- or an l-Ala residue, for three selected values of the γ(1) angle =60°, 180°, −60°. γ1 is always fixed and its value is printed at the top of each panel; γ2 is a variable. Finally, the energy scale (kcal/mol) is on the top of the figure.

Figure 5

Same as Figure 4, but for (d,l)-Ala-Pir.

Because of a weak dependence of U on γ1 and a strong dependence on γ2, U (Figures 4A–C and 5A–C) and U (Figures 4J–L and 5J–L) are similar to each other whereas U (Figures 4A–C and 5A–C) is almost a mirror image of U (Figures 4D–F and 5D–F) and U (Figures 4G–I and 5G–I). This observation is universal for each amino acid type for all γ1,γ2. However, this pseudosymmetry is not a full symmetry and does not follow the symmetry eq 9 as:where X, Y, and Z denote amino-acid-residue types of a given chirality (l or d), and a bar over a symbol denotes the change of enantiomer chirality; this operation does not affect glycine.

Comparison of the PMF surfaces corresponding to systems with the NHMe C-terminal blocking groups (panels A–L of Figure 4) with their counterparts with the Pir C-terminal blocking groups (panels A–L of Figure 5) reveals that the potentials are similar to each other; however, the regions of minima are narrower for the Pir C-terminal blocking group. Therefore, if the residue next in sequence is a proline residue, this restricts the range of the available θ angle slightly. This tendency is the same as for the torsional and double-torsional potentials.[32] It is, however, less noticeable because the virtual-bond-angle potentials are potentials corresponding to higher-order terms of the cluster-cumulant expansion of the PMF.

Thermodynamic Stability of α-Helix Formation in KLALKLALxxLKLALKLA Peptides

The calculated and experimental[42] free-energy differences of α-helix formation of the host–guest peptides with general sequence KLALKLALxxLKLALKLA (where x denotes a natural l-amino-acid residue or its d-counterpart), expressed relative to glycine (eq 7), as well as the free-energy differences of helix formation in the l- and d-amino-acid-substituted peptides (eq 8), are summarized in Table 1. It can be seen from Table 1 that the two methods used to estimate the fraction of α-helical structure described in section 2.3 resulted in similar ΔΔGhel(Gly → x) and ΔΔGhel(l-x → d-x), in the respective “contact” and “all-atom” columns, which suggests that the calculated values are not critically sensitive to the method. For example, for lysine, the ΔΔGhel(Gly → d-Lys) determined by the contact-based method is −0.69 kcal/mol and that determined based on the backbone φ and ψ angles is −0.95 kcal/mol. The value of ΔΔGhel(Gly → l-Lys) determined by the contact-based method is −3.33 kcal/mol and that determined based on the backbone ϕ and ψ angles is −3.30 kcal/mol. These values give ΔΔGhel(d-Lys → l-Lys) = 2.64 and 2.35 kcal/mol, as determined by the contact-based method and from the backbone ϕ and ψ angles, respectively. These differences are even smaller than those arising from use of different experimental methods to estimate ΔΔGhel (ref (42)).

Table 1

Calculated and Experimental Gibbs Free Energy Differences of α-Helical Structure Formation in the KLALKLALxxLKLALKLA Host–Guest Peptides for the d- [ΔΔGhel(Gly → d-x)] and the l-Amino-Acid Residues [ΔΔGhel(Gly → l-x)] Expressed Relative to Glycine, And the Free-Energy Differences of α-Helical Structure Formation between the d- and l-Amino-Acid Substituted Host–Guest Peptides [ΔΔGhel(l-x → d-x)] at 300 K

residue	ΔΔGhel(Gly → d-x)a			ΔΔGhel(Gly → l-x)a		ΔΔGhel(l-x → d-x)b
name	contactsc	all-atomd	expe	contactsc	all-atomd	contactsc	all-atomd	expe	expf
Pro	1.44	0.45	6.04	0.48	–0.25	0.96	0.70		–0.6
Lys	–0.69	–0.95	3.44	–3.33	–3.30	2.64	2.35		4.1
Arg	–1.01	–1.21	2.87	–3.48	–3.39	2.48	2.17		4.2
His	–1.00	–1.24	–0.63	–3.55	–3.42	2.54	2.18	1.76	1.6
Asp	–0.90	–1.12	–0.32	–3.49	–3.28	2.58	2.16	0.66	1.8
Glu	–1.12	–1.19	1.48	–3.48	–3.41	2.36	2.23		2.6
Asn	–0.83	–1.19	1.76	–3.36	–3.34	2.54	2.15		3.2
Gln	–1.00	–1.20	1.61	–3.50	–3.43	2.49	2.23		2.2
Ser	–0.99	–1.23	1.97	–3.31	–3.42	2.31	2.19		3.4
Thr	–0.89	–1.21	4.74	–3.65	–3.40	2.76	2.18		6.9
Ala	–1.00	–1.14	4.38	–3.30	–3.38	2.30	2.23		4.2
Tyr	–1.37	–1.39	5.68	–3.75	–3.44	2.37	2.06	4.52	3.8
Trp	–1.02	–1.19	3.84	–3.58	–3.47	2.56	2.28	4.69	4.1
Val	–0.77	–0.82	4.86	–3.61	–3.36	2.84	2.54	5.95	7.1
Leu	–0.99	–1.13	2.77	–3.61	–3.47	2.62	2.33		3.0
Ile	0.06	–0.43	5.11	–4.01	–3.46	4.07	3.03	5.84	6.0
Phe	–1.08	–1.16	4.82	–3.84	–3.48	2.76	2.32		3.1
Cys	0.08	–0.69	1.58	–1.50	–3.42	1.58	2.72		1.0
Met	–0.89	–1.18	2.95	–3.73	–3.55	2.84	2.37		2.8

See eq 7 for definition.

See eq 8 for definition.

Calculated based on peptide-group contacts; see section 2.3.

Calculated based on the backbone ϕ and ψ angles obtained by conversion of united-residue polypeptide chains to all-atom chains by using our energy-based procedure;[46,47] see section 2.3.

Based on molar ellipticity; Tables 1 and 2 of ref (42).

Based on chromatographic capacity factors determined by RP-HPLC (the bar heights in Figure 5 of ref (42)).

As can be seen from Table 1, the calculated values of ΔΔGhel(l-x → d-x) > 0, indicating that inserting a pair of d-residues into the host sequence destabilizes the α-helical structure compared to the insertion of a pair of the corresponding l-amino-acid residues, consistent with the experimental data except for proline.[42] For proline, the experimental value of ΔΔGhel(l-Pro → d-Pro) is slightly negative (−0.6 kcal/mol[42]), which means that insertion of two consecutive l-proline residues destabilizes α-helical structure slightly more than the insertion of two consecutive d-prolines. Conversely, the calculated values of ΔΔGhel(l-Pro → d-Pro) are slightly positive, independent of the method (0.96 kcal/mol for the contact-based method and 0.70 kcal/mol for the method based on the backbone ϕ and ψ angles, respectively). This observation does not indicate any gross disagreement between the UNRES and the experimental results but suggests that proline affects α-helix stability to a similar extent independent of chirality. Thus, although the extended UNRES force field was not specifically recalibrated to reproduce the thermodynamical properties of d-substituted polypeptides and proteins, it reproduces the helix-destabilizing trend of d-amino acid substitution.

It can be seen from Table 1 that the calculated values of ΔΔGhel(Gly → x) are mostly negative, which indicates that UNRES predicts that glycine destabilizes α-helix formation more than most d-amino acids, in contrast to experiment. Two or more glycine residues in a row act as a clear helix breaker when inserted into clearly helix-forming l-amino-acid-residue sequences, such as poly(hydroxybutylglutamine)[53] but the experimental data of ref (42) demonstrate that a pair of d-residues except d-His and d-Asp is a stronger helix breaker than the Gly-Gly sequence (because the experimental values of ΔΔGhel(Gly → x) are negative only for d-His and d-Asp; see column 3 of Table 1). This observation indicates that parametrization of UNRES must be enhanced to include more thermodynamic data, as is planned for further work.

Simulations of the Folding of Thurincin H

The experimental Model 1 structure of thurincin H is shown in Figure 6A. Its structure is an α-helical hairpin (Figure 6A), which is stabilized by four S–Cα bonds, which were not present in the simulations (see section 2.4) because this would impose distance restraints, while we intended to assess whether the extended UNRES can reproduce the topology of the fold in unrestricted conformational search. Nevertheless, the UNRES force field was able to fold this peptide (Figure 6). The lowest-root-mean-square deviation over the Cα atoms (Cα-RMSD between structures in Figure 6A,B) obtained after equilibration at T = 310 K was 3.86 Å. This value is comparable to the value of 3.04 Å of the mean RMSD of the structures from the experimental conformational ensemble of thurincin H determined by NMR (PDB code: 2LBZ) from the average structure of the ensemble. Clustering results reveal that the best folded structure belongs to the most probable cluster (28%) at T = 310 K. The average RMSD of the top cluster at 310 K is 5.72 Å (Figure 6C). The results suggest that thurincin H is able to fold without formation of the S–Cα bonds but, on the other hand, that the formation of these bonds can provide additional stability of the structure. When the S–Cα bonds are included, the lowest RMSD drops to 2.40 Å (Niadzvedtski, A; Sieradzan, A. K., unpublished). Because the resulting structure is not very different from the lowest-RMSD structure shown in Figure 6B, it is not shown.

Figure 6

(A) Model 1 from the ensemble of NMR structures of thurincin H (PDB code: 2LBZ);[50] black lines mark the sulfide (S–Cα) bonds, (B) the calculated structure of thurincin H with the lowest CαRMSD (3.86 Å) obtained from the MREMD simulation, and (C) average structure of the top cluster of the calculated structures of thurincin H at 310 K; Cα-RMSD from the NMR Model 1 structure is 5.72 Å.

Summary

In this work, we determined the physics-based Cα ··· Cα ··· Cα backbone-virtual-bond-angle potentials for systems that include d-amino-acid residues, and we also revised the backbone-virtual-bond-angle potentials for l-amino-acid residues that were determined in our earlier work[20] at a less advanced level of theory and are used with the present version of UNRES. Together with the virtual-bond-torsional, virtual-bond-double-torsional, side-chain-local, and backbone-electrostatic-correlation potentials determined in our recent work,[32] the new virtual-bond-angle potentials now enable us to simulate peptides and proteins, containing d-amino-acid residues, with the UNRES force field. To derive the virtual-bond-angle potentials, we used the energy surfaces of terminally blocked glycine, l-alanine, and l-proline calculated in our previous work[32] at the MP2/6-31G(d,p) ab initio level.

Because the change of the chirality of all residues is equivalent to the inversion of the sign of the virtual-bond dihedral angles γ1 and γ2, there are only 63 independent virtual-bond-angle bending potentials for each of the two types (NHMe or pyrrolidine) of C-terminal blocking groups. For l-amino-acid residues and the NHMe C-terminal blocking group, the potentials are very similar to the less accurate potentials determined in our previous work.[20] In that work,[20] the potential-energy surfaces used to determine the PMF surfaces were calculated by energy minimization with the RHF scheme with single-point MP2 correction; the Hessians were also calculated in the RHF/6-31G(d,p) scheme. In this work, the potential-energy surfaces were obtained by full minimizations of the MP2/6-31G(d,p) energy and the Hessians were also calculated with the MP2/6-31G(d,p) scheme. Moreover, in our previous work,[20] the N′,N′-dimethylamine group and not the pyrrolidine group was used to represent a succeeding proline residue.

The extension of UNRES to d-amino-acid residues, which includes the torsional, double-torsional, local, and correlation potentials determined in our previous work[32] and the virtual-bond-angle (θ) potentials determined in this work, were tested with regard to reproducing the change in the stability of the α-helical structure of the KLALKLALxxLKLALKLA peptides upon substitution of d-amino-acid residues. The force field was able to reproduce the experimentally measured[42] decrease of ΔG values of α-helical-structure formation upon d-substitution of their l-counterparts in a qualitative manner. However, to reproduce the thermodynamics more quantitatively, enhancement of terms corresponding to local interactions is required, presumably deriving potentials that are more specific to sequence. Recently[35] we started to introduce sequence-specific torsional potentials into UNRES.

We also tested the enhanced UNRES force field by studying the simulated folding of thurincin H, which has d-amino-acid residues. The UNRES force field was able to fold the thurincin H peptide to native-like conformations in unrestricted MREMD simulations, without the S–Cα covalent links that occur in the experimental structure. This suggests that thurincin H first folds into an appropriate structure even without the presence of the S–Cα bonds and then its conformational stability is increased by the formation of these bonds. The influence of S–Cα bond formation on folding and stability of thurincin H will be investigated.