The Effect of Tensile Stress on the Conformational Free Energy Landscape of Disulfide Bonds

Disulfide bridges are no longer considered to merely stabilize protein structure, but are increasingly recognized to play a functional role in many regulatory biomolecular processes. Recent studies have uncovered that the redox activity of native disulfides depends on their C–C–S–S dihedrals, and . Moreover, the interplay of chemical reactivity and mechanical stress of disulfide switches has been recently elucidated using force–clamp spectroscopy and computer simulation. The and angles have been found to change from conformations that are open to nucleophilic attack to sterically hindered, so–called closed states upon exerting tensile stress. In view of the growing evidence of the importance of C–C–S–S dihedrals in tuning the reactivity of disulfides, here we present a systematic study of the conformational diversity of disulfides as a function of tensile stress. With the help of force-clamp metadynamics simulations, we show that tensile stress brings about a large stabilization of the closed conformers, thereby giving rise to drastic changes in the conformational free energy landscape of disulfides. Statistical analysis shows that native TDi, DO and interchain Ig protein disulfides prefer open conformations, whereas the intrachain disulfide bridges in Ig proteins favor closed conformations. Correlating mechanical stress with the distance between the two –carbons of the disulfide moiety reveals that the strain of intrachain Ig protein disulfides corresponds to a mechanical activation of about 100 pN. Such mechanical activation leads to a severalfold increase of the rate of the elementary redox reaction step. All these findings constitute a step forward towards achieving a full understanding of functional disulfides.


A. Equilibrium and Metadynamics Force Field Simulations
All force field molecular dynamics (MD) simulations were carried out with the software suite GROMACS-4.5.3. 1,2 The OPLS-AA 3,4 all-atom force field was applied for simulating the polypeptide as well as cystine and diethyl disulfide (DEDS) together with using the SPC water model 5 for explicit solvation. The simulations were run in the NVT ensemble with periodic boundary conditions. The temperature of the systems was kept constant at 300 K by coupling to Nosé-Hoover chain thermostats 6 , an integration time step of 2 fs was used, trajectory coordinates were recorded every 200 fs for analysis, and all bond lengths were constrained using LINCS. 7 The electrostatic interactions were incorporated using a reaction field method 8 that has been shown to perform well for the thermodynamics of protein folding 9 and thus for computing conformational probability density maps are required here.
In order to converge the sampling of the conformational χ 2 /χ 2 -space, metadynamics acceleration 10,11 was carried out using the PLUMED plugin 12 together with GROMACS-4.5.5.
During metadynamics, the conformational sampling of the χ 2 and χ 2 torsional angles, being the collective variables for metadynamics, as defined in Fig. 1 of the main text were accelerated for all three systems using Gaussian biasing functions (height: 0.23 kcal/mol, width: 0.1 radians) and deposited every 10 ps for a total simulation length of 300 ns. We assume that by sampling up to a free energy of ≈ 17 kcal/mol via metadynamics almost all possible conformations of the disulfides should be properly sampled as the highest strain energy in a disulfide moiety is found to be of that order or magnitude 13,14 (corresponding to the Cys37-Cys54 bond of 3-oxoacyl-acyl carrier protein reductase from Thermotoga maritima, PDB ID 1O5I). In the case of DEDS and cystine, we assume this value to be even lower as they are free at their termini and hence not that much constrained as disulfide bonds in proteins.
DEDS, cystine and the polypeptide model Ile-Cys-Leu-Ser-Glu-Pro-Asp-Val-His-Cys-Gln (with the disulfide bond between the 2nd and 10th cysteine residue, were generated using PyMOL 15 , minimized (1000 steps using the steepest descent algorithm) and subsequently equilibrated. The species where solvated with SPC water in a periodic box of (3.0 x 3.0 x 3.0) nm 3 for DEDS and cystine whereas a simulation box of (5.0 x 4.0 x 4.0) nm 3 was used for the polypeptide using 886, 880, and 2633 solvent molecules, respectively In addition to the solvent, 10 Na + and 9 Cl − ions were also added to the protein model system in order to keep the system neutral. The termini of the polypeptide were capped with hydrogens as to avoid any Coulomb attraction bias on the conformational free energies when using such a rather short peptide.
Force-clamp MD simulations of these species (that have been equilibrated at zero force) were then performed using the same computational approach as described above by explicitly including constant forces 16  by deuterium masses to allow for larger fictitious orbital mass and thus a larger time step resulting into more efficient AIMD propagation. 23 For each force value, the system was equilibrated initially for 2 ps followed by production runs of 200 ps for analyses.

C. QM Simulations on Disulfide Macrocycle
In previous experiments 29 , chemically engineered macrocycles, so-called "molecular force probes", 30,31 have been used in order to generate internal mechanical stress that strains the disulfide bond that is embedded in the macrocycle. 19 For such systems the tensile force corresponding to the produced strain has been extracted by a validated approach. 29,31 In order to compute the corresponding C α -C α distance at room temperature, an AIMD 23 simulation of an isolated disulfide macrocycle (macrocycle 9 29 , has been performed. Importantly, no external force must be applied to the macrocyle, where strain is generated purely intramolecularly due to the photoswitch, stiff stilbene, being in the cis/trans conformation, see Fig. 1. These AIMD simulations have been conducted as explained in the previous section using a 20Å cubic box, proper hydrogen masses, and thus a timestep of 3 a.u. with µ = 400 a.u. was set to be a cubic box of 20Å. Again, the system was initially equilibrated for 2 ps followed by 120 ps for computing the average C α -C α distance.