Evaluating the stability of disulfide bridges in proteins: a torsional potential energy surface for diethyl disulfide

Disulfide bonds formed by the oxidation of cysteine residues in proteins are the major form of intra- and inter-molecular covalent linkages in the polypeptide chain. To better understand the conformational energetics of this linkage, we have used the MP2(full)/6-31G(d) method to generate a full potential energy surface (PES) for the torsion of the model compound diethyl disulfide (DEDS) around its three critical dihedral angles (χ2, χ3, χ2′). The use of ten degree increments for each of the parameters resulted in a continuous, fine-grained surface. This allowed us to accurately predict the relative stabilities of disulfide bonds in high resolution structures from the Protein Data Bank. The MP2(full) surface showed significant qualitative differences from the PES calculated using the Amber force field. In particular, a different ordering was seen for the relative energies of the local minima. Thus, Amber energies are not reliable for comparison of the relative stabilities of disulfide bonds. Surprisingly, the surface did not show a minimum associated with χ2 ∼ − 60°, χ3 ∼ 90, χ2′ ∼ − 60°. This is due to steric interference between Hα atoms. Despite this, significant populations of disulfides were found to adopt this conformation. In most cases this conformation is associated with an unusual secondary structure motif, the cross-strand disulfide. The relative instability of cross-strand disulfides is of great interest, as they have the potential to act as functional switches in redox processes.


Introduction
Disulfide bondsb etweeno xidised cysteiner esidues are generally viewed as structurally stabilising elements in proteins. However,anewrole for asubset of disulfides as redox switches is emerging. Redox switching of disulfide bonds has been demonstrated in both reversible and irreversible redox regulation of proteins. Reversible systems include thosei nvolved in redox signalling such as the peroxide sensor, OxyR, where disulfidebondf ormation activates thet ranscriptionf actor in response to oxidative stress [1]. Irreversible redox regulation mediated by disulfide reduction and subsequent irreversible conformationalc hange has also been described. Fore xample, reduction of disulfides and subsequent cleavage of protein chains has been demonstrated in ovotransferrin and plasmin [2,3] and is likelytobeanimportant regulatory mechanism for many otherp roteins.
In principle it shouldb ep ossible to differentiate betweenr edox-active and structurally-stabilising disulfides by analysis of proteins tructuresa nd ultimately proteinsequences. Our previous studies haveinvestigated high disulfide torsional energies as indicatorso fr edox activity as well as identifying structural motifsassociated with redox activity [4,5]. Torsional strain on the disulfide bond is imposedb yt he geometric constraintso ft he protein-fold. As the force constants for torsion around the dihedral angles are much lower than for the stretching and compressing of bond lengths and bond angles, this strain is expected to be mostly accommodated via torsion of the five criticald ihedral angles of the disulfide (figure1). From analysisofthese torsional angles in high resolution X-ray structureso fp roteins, it should be possiblet om akeag ood predictiono ft he strain of a disulfide bond, and thus determine howl ikely it is to undergo redox processes.
In previous work [ 4], we estimated the stability of a disulfide bond usingt he torsionalp otential energy function [6] from the Amber force field [7]: While this function has the right general form, it clearly does nottake into account the steric interactions within the system.Thisresults in apotential energy surface (PES) in which all of the local minima are, incorrectly,predicted to be equally stable. Thefunction is, therefore, not accurate when comparing disulfide stabilities in real systems. Ab etter description of the relative stabilities which includess teric effects can be found using fullA mber calculations, i.e. including non-bonded terms, but these havea lso been reported to give inaccurate results [8], ac laim which will be investigated in this work.
In 1994, Gö rbitz [ 9] attempted to improvet he understandingo ft he conformational preferences of disulfide bridges by performing ab initio calculations on the model compound, diethyl disulfide (DEDS). This model system is too small to allowinvestigation of the x 1 and x 1 0 dihedral angles, however,i ti sp ossible to determine howt he energy of the system changesa s x 2 , x 2 0 and x 3 are varied. This is where the majority of the inaccuracy in the Amber function is expected to lie. Gö rbitz employed both the Hartree -Fock (HF) and MP2(full)l evels of theory,w ith basiss ets up to 6-311G(2d,p), to investigate the relative stabilities of the minima and saddle points on the PES defined by x 2 , x 2 0 and x 3 .The dihedral angles x 2 and x 2 0 were found to prefer values of þ 60, 2 60 and 1808 (which Gö rbitz labelled G, G 0 and T, respectively), while x 3 preferred to be þ 90 or 2 90 (G and G 0 ). This gaver iset os ix symmetrically distinct minima and eightd istinct saddlep oints. Gö rbitz found that methods which involved electron correlation (MP2(full)) weren ecessary to correctlyp redictt he relatives tabilities of the critical points. At the MP2(full) level,t he preferred conformation was found to be GGG ( x 2 , 608 , x 3 , 908 , x 2 0 , 608 ), followed by G 0 GG, TGG, G 0 GT,T GT and G 0 GG 0 .W ith HF,G GG was correctly predicted to be the mosts table,h owever,T GG and TGT were not appreciably higherine nergy.
While Gö rbitz'scalculations were state-of-the-art at the time they were reported, there haveb een significant developments both in quantum chemical methods and in computational power.I np articular,d ensity functional methods (sucha sB 3LYP [10 -14]) haveb ecome widely accepted; and schemes,s uch as G3 [15] and G3X [16], haveb een developed for performing highly accurate calculations at relativelyl ow computational cost.
In order to distinguish between disulfide bridges that are simplyp erforming as tructural role and those which are likelytoberedox active, we need to be able to accurately predict the relative stabilities of disulfide bonds. It is necessary not only to understand the relative stability of the torsional minima but also to haveagood description of the entire PES. Like Gö rbitz, we havechosentof ocus on the three central dihedral angles, x 2 , x 2 0 and x 3 ,t hus reducing av ery large five-dimensionalp roblem to af ar more tractable three dimensions. We expect that the torsiona round the carbon-carbon bonds, x 1 and x 1 0 , shouldb er elativelyw ell described in the Amber force field. Also, x 1 and x 1 0 do not, in general, showsignificant deviation from their optimal values. The goal of this work, therefore, is to create an ew three-dimensional potential energy surface (3D-PES) for the torsionofDEDSaround the x 2 , x 2 0 and x 3 dihedral angles.

Methods
Benchmarking calculations were initially carriedo ut in order to determine the most reliable and cost effective level of theory with which to determinet he 3D-PES. Reference energies for the minima and lowl ying saddle points were calculated usingt he G3X method [16]. G3X involves optimisingt he geometrya tt he B3LYP/ 6-31G(2df,p) level of theory,t hen performing as inglepointe nergyc alculationw ith QCISD(T)/6-31G(d).
Further calculations are then used to correct this energy for the effects of including diffuse and higher polarisation functions in the basis set (at the MP4l evel), including correlation of core electrons and even higher polarisation functions (at the MP2 level)and including gfunctions on the sulfur atom (at the HF level). AG 3X electronic energy is thus obtained. Usually the zero-point energy is also included( calculated using B3LYP/6-31G(2df,p)), however,i nt his case it was omitted as we neededt ou se the reference energies to findt he level of theory which gavet he best possibleprediction of the electronic energy of DEDS. Thea dditionalc alculationo fz ero-point energies at everyp oint of our 3D-PESw ould be far too computationally expensive. G3X is reported to yield an accuracyo f^4kJmol 2 1 fort he calculationo f heat of formationf roma tomisation energies. It is, therefore, expected to be highly accuratefor the prediction of the relative energies of theD EDSm inimaa nd saddle points. The G3X electronicenergies were then compared with the energies from fully optimised calculations for each of thec riticalp oints, calculated usingH F/6-31G(d), B3LYP/6-31G(d), B3LYP/6-31G(2df,p)a nd MP2(full)/6-31G(d). Amber energies for each of these criticalp oints were also calculated for comparison.
The 3D-PES was calculated at the MP2(full)/6-31G(d) level of theory.E nergies were calculated at ten degree increments in x 2 , x 2 0 and x 3 to give the full 3D grid. As Amber calculations for DEDS were relatively cheap to perform,asimilar 3D-PES was created usingA mber for comparison. This was only done for x 3 values between 60 and 1308 as these represent the x 3 values adopted by over 99% of the high resolutionX -ray structures found in the Protein Data Bank (vide infra).
The small increments used to calculate the PES resulted in as urface which was sufficiently fine grained that a simple linear interpolation couldb eu sed to predict the energies of disulfides with ag iven set of x 2 , x 2 0 and x 3 dihedral angles. This methodology was used to predict the relative stabilities of the disulfides in our database of high resolution disulfidesi nt he Protein Data Bank [17]. Please see Ref. [18] for details of howt his database was constructed.
All calculations were performed using the Gaussian 03 suite of programmes [19].T hiss uite includest he current version of Amber,version 9. Calculations wereperformed on the SC and LC computer clusters of the Australian Partnership of Advanced Computing (APAC) National Facility and on the localcomputer clusterofthe Computer ChemieCentrum, Erlangen, Germany.

Benchmark calculations
The results of the benchmarking calculations are shown in table 1. Energies are reported relativetothe energy of the GGG conformation (theminimum 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, togetherw ith the maximum deviation seen. The HF and MP2(full) results are the same as those described by Gö rbitz [9] and havebeen includedhere for comparison. Ve ry little deviation was seen betweent he geometries of the different conformations at the different levelso ft heory.
The Amber energies in column 2s howt he most significant deviation from the benchmarks, both in terms of the RMS deviation (2.3 kJ mol 2 1 )a nd in having the largest discrepancy for any one configuration (TGTbeing predicted to be 5.4 kJ mol 2 1 too stable relativet oG GG). Most importantly,confirming earlier reports [8], the order of the stabilities of the criticalpointsisnot consistent with the benchmark G3X results. The comparison betweenthe HF results and the newG 3X benchmarks was similarly poor,w ith an RMS of 2.0 kJ mol 2 1 and am aximum deviation of 4.5 kJ mol 2 1 .T he order of the critical points was better than observed for the Amber forcefi eld calculations but was still not correct in somecases.
The density functional results, with both the 6-31G(d) and 6-31G(2df,p) basis sets, also showed surprisingly poor agreement with the benchmarkrelative energies. Although Table 1. Relative energies (kJ mol 2 1 )ofthe diethyl disulfide minima and lowenergy saddle points at various levels of theory.Also included are mean, RMS and maximum deviations from the highest levelo ftheory,G 3X.
the order of the minima was nowcorrectly described, both methods predicted the GGG 0 and GGT conformations to be roughly equal in energy,a nd likewise the G 0 GT and TGT minima. The G3X calculations showed these to be separated by 1.6 and1 .7 kJ mol 2 1 ,r espectively. In addition, higher energy structures were, in mostc ases, predicted to be too stable, that is, the PES is predicted to be too flat. This is as ignificant problem for this work, where the highere nergy structuresa re those of greatest interest and need to be described as accurately as possible.

The 3D-PES
The PES is displayed in the form of contour plots for increasing values of x 3 .Figure2shows slices through the surface for x 3 between 60 and 1308 .C ontour plots for x 3 valuesoutside this range can be found in Appendix1.Note that the contour plots for negative x 3 values are related by symmetry to thosefor x 3 . 0. In all contours the energies are shown relative to the energy of the absolute minimum of the 3D-PES, x 2 ¼ 708 , x 3 ¼ 908 , x 2 0 ¼ 708 .The colour scheme has been chosens ot hat blue regions are lowi n energy ( , 10 kJ mol 2 1 ). Disulfides with energies in this region are expected to be stable.C yan indicates higher energy regions (energy between 10 and 15 kJ mol 2 1 )a nd green represents very high energy regions (energy between1 5a nd 20 kJ mol 2 1 ). Disulfides in these regions are expectedt ob el esss table and morel ikely to be involved in redox processes. All othercolours correspond to extremely high energy regions( . 20 kJ mol 2 1 ). Disulfides are not expected to be found with such high torsional strain.
At the lowestenergy point on the surface, x 3 ¼ 908 ,the different minima reported by Gö rbitz [8] are seen as stable (dark blue) areas, with the saddlep oint regions between them also, in most cases, beingoflow energy (light blue). The area of the dark blue circles generally gives ag ood indication of the relative stability of the minima. Likewise the area of the greenr egionsi ndicatest he relative instability of the maximai nt his slice. The plot appears largely as expected but, surprisingly,aminimum is not seen for x 2 , 2 608 , x 3 , 908 , x 2 0 , 2 608 ,t hat is, for the G 0 GG 0 conformation. Detailedchecksofthe relative energy of each gridpoint in the G 0 GG 0 region confirm that the absence of thism inimumi sn ot an artefacto ft he positioning of the contours.T he same situationi sf ound for the contour with x 3 ¼ 808 ;w hile av ery shallow minimum,w ith ad epth of , 1kJmol 2 1 ,i ss een for x 3 ¼ 1008 .W en ote here that, whereas the x 3 values of thef ully optimised benchmarkc alculations were approximately 908 fora lmosta ll conformations,f or G 0 GG 0 the value had increased to approximately 1108 with all levels of theory.T his unexpected instability of the G 0 GG 0 region of the PES is due to steric clashes between the terminal methyl groups.Inproteins, thesecorrespond to the H a atoms. As shown in figure 1, the H a atoms of a G 0 GG 0 ( x 2 , x 2 0 , 2 60)d isulfidea re forcedi nto particularly close proximity when x 3 ¼ 908 .T hey, therefore, experience strong repulsive interactions and the system is destabilised. These steric interactions are reduced as x 3 moves away from 908 ,thus allowing minima to appear for x 3 # 708 and x 3 $ 1108 .Interestingly,for x 3 ¼ 708 the minimum is surprisingly deep (4.4 kJ mol 2 1 ), although it becomes shallower for smaller x 3 values.
Thee ffects of sterici nterference arep articularly important for lower values of x 3 .W hen the x 2 and x 2 0 dihedral angles are both small, the terminal methyl hydrogens come into very closec ontact as x 3 is reduced. This results in the high energy featuren ear the origin (actually at x 2 ¼ x 2 0 , 208 )w hich grows rapidlya s x 3 is reduced below9 0 8 .T he growth of this feature also has a significant adverse effect on the stabilities of the GGG 0 and G 0 GT conformations for x 3 # 708 .A lthough these minima are not seen on the contour plots for x 3 ¼ 60 and 708 ,they do exist. Forboth contours the minima are very shallow(, 2kJmol 2 1 ), but they become deeperagain for x 3 , 508 .For x 3 , 408 the steric repulsion is so greatthat the entire contour plot lies in the extremely high energy region, above2 0kJmol 2 1 .D isulfides are not expected to occur in these regions. As x 3 increases above9 0 8 ,t he contour plot becomes more symmetricald ue to the reduction of the steric interactions betweenthe methyl groups.Inparticular,t he GGG 0 conformation drops in energy so that for x 3 ¼ 1008 it is equal in energy with GGG,and for x 3 ¼ 110 and 1208 it is actually the most stable conformation on the PES. When x 3 is increased to 1208 there is no longer stericstrain in the G 0 GG 0 region so that the G 0 GG 0 conformation is now of equal stability to GGG.T he PES continues to look effectivelysymmetrical for all higher values of x 3 .For x 3 values of 1508 and above, the entire surface is more than 20 kJ mol 2 1 abovet he minimum. Again, disulfides with these large x 3 values are not expected to exist.
All the significant features seen in the MP2(full) 3D-PES are also found in the Amber forcefieldsurface,albeit shiftedtoslightly lower x 3 angles. The Amber PES does, however,s eem to be rather flatter than the MP2(full) version, with the energy not rising as quickly as x 3 moves

Comparison of PES with observed disulfide conformations
An important testofthe usefulness of our PES is to check that the disulfide conformations that are predicted to be the most stable actually correspond to thosem ostc ommonly seen in proteins. Figure3shows the distribution of x 3 valueso btainedf rom our database of high resolution X-ray structures. Superimposed on this distribution is the (one-dimensional) PES for torsiono f x 3 in the most populous GGG conformation. This gives arough guide as to howthe shape of the 3D-PESchanges with x 3 .W enote that the 1D-PESs for mosto therc onformations havet he same general shape,a lthough the minima may be at slightly different dihedral angles. It can be seen that the most stable point on the PES, at about 908 ,d oes indeed correspond to the highest population of disulfides. The population falls rapidlya s x 3 increases or decreases and the energy rises. The slight shift of the histogram towards higher x 3 values is probably associated with the minima for someconformations being as high as 1108 (G 0 GG 0 ). It is also interesting to compare howthe disulfides in the PDB are distributed amongst the possible conformations. Figure 4shows ascatter plot of the conformations adopted by disulfidesw ith x 3 between 85 and 958 .T his is superimposed on the PES slice for x 3 ¼ 908 .C learly the majority of the disulfidesa re 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 x 3 ¼ 908 the saddle points are also in the lowe nergy region of the PES and significant populations of disulfidesare found with the corresponding sets of dihedral angles.
What is most interestingisthatanappreciable number of disulfidesisalsoseeninthe high energy G 0 GG 0 region (top RH corner). Furtheranalysisofthe structures whichadopt this conformation hasrevealedthat, in almost allcases,the disulfide is fixed in this conformation by thep rotein secondarys tructure.I np articular,m osto ft hese disulfides aref ound to bridge two neighbourings trands in an antiparallel b -sheet.T hiss econdary structurem otif is knownasacross-strand disulfide [4,18,20]. In addition to the strain duetothe unfavourable conformation,the associated b -sheetsare also significantly distortedbyt he presence of these disulfides. They are, therefore, highly strained andp resent very promisingc andidatesf or involvement in redoxp rocesses.F urther analysis of this interesting classofdisulfides hasbeenreportedelsewhere [5].

Use of the 3D-PES to predict disulfide stability
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 asimplethree-dimensional linear interpolation on the calculated PES. The effects of strain in the x 1 and x 1 0 dihedral angles were not takeni nto account in this investigation.
Using our MP2(full) PES, the mean strainenergy of the disulfides in our database was found to be 7.1 kJ mol 2 1 , with as tandardd eviation of 4.8 kJ mol 2 1 .S eventy-nine  percent of disulfides were found to haver elatively low energy ( , 10 kJ mol 2 1 abovethe minimum), with afurther 18% being in the high energy region (between 10 and 15 kJ mol 2 1 ). Only 3% had ar elativee nergy higher than 15 kJ mol 2 1 .
Ah istogram showingt he energyd istributiono ft he disulfides in the PDB can be found in figure5 .E nergy distributions calculated using the Ambertorsional potential and the Amber 3D-PES have beenincluded in addition to the MP2(full) results.T he average disulfide energiesa re also shown along withthe standarddeviations. As discussed in Methods,onlydisulfides with x 3 between 60 and 1308 are included, resulting in aslightlydifferent MP2(full) average comparedt ot hatf or the completed atasets tated above. Usinge nergiesc alculated with theA mber torsional potential, the histogram is seent ob es trongly biased towardsl ower energies; this is an aturalr esult of the predictionbythis potentialthatall the configurations on the PES are equallystable and thusall have arelative energyof zero. Using thee nergiesf romt he Amber3 D-PES, the histogram is much more consistent witht hatf romt he MP2(full)surface, butisstill biased towards lowerenergies. Thisi sc learlys eenb ye xamining the modal energiesf or eachmethod. The histogram peaks at 2.5 kJ mol 2 1 if only the Amber torsional energyi sc onsidered; at 5.0 kJ mol 2 1 when the entireA mberp otential is considered; and at 7.5 kJ mol 2 1 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].

Conclusions
We haves uccessfully constructed aM P2(full)/6-31G(d) PES for the torsiono fD EDS aroundi ts three important dihedral angles. This surface was found to be qualitatively different from that which was predicted usinge ither the Amber torsional energy function or the full Amber force field. In particular,the relative stabilities of the minima on the MP2(full) surface were found to be in good agreement with the G3X benchmarkcalculations, whereas the Amber force field gaven ot only large deviations in the relative energies but also adifferent order for the stabilities of the conformations.T hiso rder is likely to be important in elucidating the mechanisms of reactions that involve a cascadeo fd isulfides.O ne such exampleo ccursi n Staphylococcus aureus Arsenater eductase,i nw hich stepwise formation of the Cys 10 -Cys 82 and Cys 82 -Cys 89 disulfidesf orm part of the reaction cycle to detoxify arsenic [21]. TheC ys 10-Cys 82 disulfide is as hort-lived high energy intermediate trapped in the crystal structure of the Cys89Leu mutant (PDB1lk0). Upon formation, the disulfide likely adopts the SG 0 Tc onformation with arelative energy of 14.9 kJ mol 2 1 (Chain A). Subsequent movements of aflexible region of the backbone (residues8 2-97) twistt he disulfide into an S 0 GG Figure 5. Comparison of relative energies for disulfides in high resolution structures of the PDB as predicted by the (a) Amber torsional potential, (b) full Amber potential including non-bonded terms, and (c) quantum chemical calculations using the MP2(full) level of theory.Note particularly the populations peak in different energy bins. To ensure afair comparison, only disulfides with x 3 between 60 and 1308 were included (see Methods). conformation with ar elativee nergy of 24.0 kJ mol 2 1 , rendering it susceptiblet oa ttackb yt he nearby Cys 89 thiolate (Chain B). Thesubsequently formed Cys 82 -Cys 89 disulfide adopts the G 0 G 0 G 0 conformation with alower energy of 12.2 kJ mol 2 1 ,a nd is sufficiently stable that it must be reduced by thioredoxin in order to regenerate Arsenate reductase for the next reaction cycle (PDB 1lju).
The relative configurational stabilities of the MP2(full) PES were alsof ound to be more consistent with experimentald ata fort he populations of disulfides, which adopt the associated conformations.
Unexpectedly,t he 3D-PES didn ot showaminimum associated witht he G 0 GG 0 conformation for x 3 values of 80 and 908 .Thisisaresult of strongsteric interactions with this particular set of dihedral angles ( x 2 < x 2 0 < 2 608 , x 3 < 908 ). In this conformationt he H a atoms are aligned directlyt owardse acho ther,t hus experiencings trong repulsive forces thatd estabilise the system.Also surprising wasthatasignificantpopulation of disulfides were found to adopt thish ighe nergy conformation.F urther analysis revealedt hati nm ost cases this conformation arose from(andwas requiredfor)anunusual secondary structure motif, the cross-strand disulfide.
The 3D-PES was subsequently used to predict the relative stabilities of all the high resolutiondisulfide bonds reportedi nt he Protein Data Bank. As expected, the vast majority of the disulfides were found to havealow strain energy and are, therefore, likely to be involved solely in structural stabilisation.A pproximately 20%o ft he cystines were of high or very high relative energy and thus havethe potential to be involved in redox processes. Further investigation of these disulfidesiso ngoing.