The results of the benchmarking calculations are shown in table 1. Energies are reported relative to the energy of the GGG conformation (the minimum on the PES at the G3X level of theory). Mean and RMS deviations from the G3X results are reported for each of the levels of theory investigated, together with the maximum deviation seen. The HF and MP2(full) results are the same as those described by Görbitz [9] and have been included here for comparison. Very little deviation was seen between the geometries of the different conformations at the different levels of theory.

The Amber energies in column 2 show the most significant deviation from the benchmarks, both in terms of the RMS deviation (2.3 kJ mol⁻¹) and in having the largest discrepancy for any one configuration (TGT being predicted to be 5.4 kJ mol⁻¹ too stable relative to GGG). Most importantly, confirming earlier reports [8], the order of the stabilities of the critical points is not consistent with the benchmark G3X results. The comparison between the HF results and the new G3X benchmarks was similarly poor, with an RMS of 2.0 kJ mol⁻¹ and a maximum deviation of 4.5 kJ mol⁻¹. The order of the critical points was better than observed for the Amber force field calculations but was still not correct in some cases.

The density functional results, with both the 6-31G(d) and 6-31G(2df,p) basis sets, also showed surprisingly poor agreement with the benchmark relative energies. Although the order of the minima was now correctly described, both methods predicted the GGG′ and GGT conformations to be roughly equal in energy, and likewise the G′GT and TGT minima. The G3X calculations showed these to be separated by 1.6 and 1.7 kJ mol⁻¹, respectively. In addition, higher energy structures were, in most cases, predicted to be too stable, that is, the PES is predicted to be too flat. This is a significant problem for this work, where the higher energy structures are those of greatest interest and need to be described as accurately as possible. The RMS deviations were, however, significantly smaller than those calculated with either Amber or HF.

The MP2 calculations were again found to give by far the best agreement with the benchmarks. The RMS deviation was only 0.6 kJ mol⁻¹ and the maximum deviation 1.2 kJ mol⁻¹. MP2(full) was, therefore, chosen as the most reliable level of theory with which to calculate the 3D-PES.

The PES is displayed in the form of contour plots for increasing values of χ₃. Figure 2 shows slices through the surface for χ₃ between 60 and 130°. Contour plots for χ₃ values outside this range can be found in Appendix 1. Note that the contour plots for negative χ₃ values are related by symmetry to those for χ₃ > 0. In all contours the energies are shown relative to the energy of the absolute minimum of the 3D-PES, χ₂ = 70°, χ₃ = 90°, χ₂′ = 70°. The colour scheme has been chosen so that blue regions are low in energy (<10 kJ mol⁻¹). Disulfides with energies in this region are expected to be stable. Cyan indicates higher energy regions (energy between 10 and 15 kJ mol⁻¹) and green represents very high energy regions (energy between 15 and 20 kJ mol⁻¹). Disulfides in these regions are expected to be less stable and more likely to be involved in redox processes. All other colours correspond to extremely high energy regions (>20 kJ mol⁻¹). Disulfides are not expected to be found with such high torsional strain.

At the lowest energy point on the surface, χ₃ = 90°, the different minima reported by Görbitz [8] are seen as stable (dark blue) areas, with the saddle point regions between them also, in most cases, being of low energy (light blue). The area of the dark blue circles generally gives a good indication of the relative stability of the minima. Likewise the area of the green regions indicates the relative instability of the maxima in this slice. The plot appears largely as expected but, surprisingly, a minimum is not seen for χ₂ ∼ − 60°, χ₃ ∼ 90°, χ₂′ ∼ − 60°, that is, for the G′GG′ conformation. Detailed checks of the relative energy of each gridpoint in the G′GG′ region confirm that the absence of this minimum is not an artefact of the positioning of the contours. The same situation is found for the contour with χ₃ = 80°; while a very shallow minimum, with a depth of ∼1 kJ mol⁻¹, is seen for χ₃ = 100°. We note here that, whereas the χ₃ values of the fully optimised benchmark calculations were approximately 90° for almost all conformations, for G′GG′ the value had increased to approximately 110° with all levels of theory. This unexpected instability of the G′GG′ region of the PES is due to steric clashes between the terminal methyl groups. In proteins, these correspond to the Hα atoms. As shown in figure 1, the Hα atoms of a G′GG′ (χ₂ ∼ χ₂′ ∼ − 60) disulfide are forced into particularly close proximity when χ₃ = 90°. They, therefore, experience strong repulsive interactions and the system is destabilised. These steric interactions are reduced as χ₃ moves away from 90°, thus allowing minima to appear for χ₃ ≤ 70° and χ₃ ≥ 110°. Interestingly, for χ₃ = 70° the minimum is surprisingly deep (4.4 kJ mol⁻¹), although it becomes shallower for smaller χ₃ values.

The effects of steric interference are particularly important for lower values of χ₃. When the χ₂ and χ₂′ dihedral angles are both small, the terminal methyl hydrogens come into very close contact as χ₃ is reduced. This results in the high energy feature near the origin (actually at χ₂ = χ₂′ ∼ 20°) which grows rapidly as χ₃ is reduced below 90°. The growth of this feature also has a significant adverse effect on the stabilities of the GGG′ and G′GT conformations for χ₃ ≤ 70°. Although these minima are not seen on the contour plots for χ₃ = 60 and 70°,they do exist. For both contours the minima are very shallow (∼2 kJ mol⁻¹), but they become deeper again for χ₃ < 50°. For χ₃ < 40° the steric repulsion is so great that the entire contour plot lies in the extremely high energy region, above 20 kJ mol⁻¹. Disulfides are not expected to occur in these regions.

As χ₃ increases above 90°, the contour plot becomes more symmetrical due to the reduction of the steric interactions between the methyl groups. In particular, the GGG′ conformation drops in energy so that for χ₃ = 100° it is equal in energy with GGG, and for χ₃ = 110 and 120° it is actually the most stable conformation on the PES. When χ₃ is increased to 120° there is no longer steric strain in the G′GG′ region so that the G′GG′ conformation is now of equal stability to GGG. The PES continues to look effectively symmetrical for all higher values of χ₃. For χ₃ values of 150° and above, the entire surface is more than 20 kJ mol⁻¹ above the minimum. Again, disulfides with these large χ₃ values are not expected to exist.

For comparison, contour plots of the Amber force field 3D-PES are shown in Appendix 2. Note that the absolute minimum on the Amber PES is at χ₂ = 180°, χ₃ = 80°, χ₂′ = 180°, rather than at χ₂ = 70°, χ₃ = 90°, χ₂′ = 70° as seen with MP2(full). However, all energies are shown relative to the MP2(full) absolute minimum to give the clearest comparison between the two methods. Areas shown in very dark blue lie below the MP2(full) minimum. We note, however, that the difference between χ₂ = 180°, χ₃ = 80°, χ₂′ = 180° and χ₂ = 70°, χ₃ = 90°, χ₂′ = 70° on the Amber PES is only 1.5 kJ mol⁻¹.

All the significant features seen in the MP2(full) 3D-PES are also found in the Amber force field surface, albeit shifted to slightly lower χ₃ angles. The Amber PES does, however, seem to be rather flatter than the MP2(full) version, with the energy not rising as quickly as χ₃ moves away from 90° (even when the difference in reference energy is taken into account). The most significant discrepancy, is in the prediction of the relative stabilities of the local minima. Comparison of figures 2 and A2 (Appendix 2) clearly shows, amongst other problems, an inversion of the order of the GGG, GGT and TGT conformations for most values of χ₃. This issue has already been discussed in detail in Section 3.1.

An important test of the usefulness of our PES is to check that the disulfide conformations that are predicted to be the most stable actually correspond to those most commonly seen in proteins. Figure 3 shows the distribution of χ₃ values obtained from our database of high resolution X-ray structures. Superimposed on this distribution is the (one-dimensional) PES for torsion of χ₃ in the most populous GGG conformation. This gives a rough guide as to how the shape of the 3D-PES changes with χ₃. We note that the 1D-PESs for most other conformations have the same general shape, although the minima may be at slightly different dihedral angles. It can be seen that the most stable point on the PES, at about 90°, does indeed correspond to the highest population of disulfides. The population falls rapidly as χ₃ increases or decreases and the energy rises. The slight shift of the histogram towards higher χ₃ values is probably associated with the minima for some conformations being as high as 110° (G′GG′).

It is also interesting to compare how the disulfides in the PDB are distributed amongst the possible conformations. Figure 4 shows a scatter plot of the conformations adopted by disulfides with χ₃ between 85 and 95°. This is superimposed on the PES slice for χ₃ = 90°. Clearly the majority of the disulfides are located in either the low energy (blue) or high energy (cyan) regions, with very few examples in the very high energy zones (green, yellow, etc.). By far the majority of the disulfides adopt the lowest energy GGG conformation (bottom LH corner). However, concentrated patches are seen for each of the different minima. We note that for χ₃ = 90° the saddle points are also in the low energy region of the PES and significant populations of disulfides are found with the corresponding sets of dihedral angles.

What is most interesting is that an appreciable number of disulfides is also seen in the high energy G′GG′ region (top RH corner). Further analysis of the structures which adopt this conformation has revealed that, in almost all cases, the disulfide is fixed in this conformation by the protein secondary structure. In particular, most of these disulfides are found to bridge two neighbouring strands in an antiparallel β-sheet. This secondary structure motif is known as a cross-strand disulfide [4,18,20]. In addition to the strain due to the unfavourable conformation, the associated β-sheets are also significantly distorted by the presence of these disulfides. They are, therefore, highly strained and present very promising candidates for involvement in redox processes. Further analysis of this interesting class of disulfides has been reported elsewhere [5].

Finally, the 3D-PES was used to predict the strain in each of the disulfide bonds found in our database of high resolution structures from the Protein Data Bank. This was done using a simple three-dimensional linear interpolation on the calculated PES. The effects of strain in the χ₁ and χ₁′ dihedral angles were not taken into account in this investigation.

Using our MP2(full) PES, the mean strain energy of the disulfides in our database was found to be 7.1 kJ mol⁻¹, with a standard deviation of 4.8 kJ mol⁻¹. Seventy-nine percent of disulfides were found to have relatively low energy (<10 kJ mol⁻¹ above the minimum), with a further 18% being in the high energy region (between 10 and 15 kJ mol⁻¹). Only 3% had a relative energy higher than 15 kJ mol⁻¹.

A histogram showing the energy distribution of the disulfides in the PDB can be found in figure 5. Energy distributions calculated using the Amber torsional potential and the Amber 3D-PES have been included in addition to the MP2(full) results. The average disulfide energies are also shown along with the standard deviations. As discussed in Methods, only disulfides with χ₃ between 60 and 130° are included, resulting in a slightly different MP2(full) average compared to that for the complete dataset stated above. Using energies calculated with the Amber torsional potential, the histogram is seen to be strongly biased towards lower energies; this is a natural result of the prediction by this potential that all the configurations on the PES are equally stable and thus all have a relative energy of zero. Using the energies from the Amber 3D-PES, the histogram is much more consistent with that from the MP2(full) surface, but is still biased towards lower energies. This is clearly seen by examining the modal energies for each method. The histogram peaks at 2.5 kJ mol⁻¹ if only the Amber torsional energy is considered; at 5.0 kJ mol⁻¹ when the entire Amber potential is considered; and at 7.5 kJ mol⁻¹ using the more exact MP2(full) calculations.

Further analysis of the relative energies associated with each of the different disulfide conformations as well as with various secondary structure elements will be reported elsewhere [18].