The rupture mechanism of rubredoxin is more complex than previously thought†

The surprisingly low rupture force and remarkable mechanical anisotropy of rubredoxin have been known for several years. Exploiting the first combination of steered molecular dynamics and the quantum chemical Judgement of Energy DIstribution (JEDI) analysis, the common belief that hydrogen bonds between neighboring amino acid backbones and the sulfur atoms of the central FeS4 unit in rubredoxin determine the low mechanical resistance of the protein is invalidated. The distribution of strain energy in the central part of rubredoxin is elucidated in real-time with unprecedented detail, giving important insights into the mechanical unfolding pathway of rubredoxin. While structural anisotropy as well as the contribution of angle bendings in the FeS4 unit have a significant influence on the mechanical properties of rubredoxin, these factors are insufficient to explain the experimentally observed low rupture force. Instead, the rupture mechanism of rubredoxin is far more complex than previously thought and requires more than just a hydrogen bond network.


Benchmark of Density Functionals Against Experimental Data
Various density functionals were benchmarked against experimental data. The tested density functionals comprise B3LYP, 1-3 B97, 4 B97M-rV, 5 BHHLYP, 1-3 , BLYP, 1,2 BP86VWN, 1,6 M06-L, 7 PBE, 8 TPSS,9 ωB97M-V, 10 and ωB97X-V. 11 Although multi-reference ab initio methods are of course to be preferred over Density Functional Theory (DFT) in the case of metal complexes, the large number of calculations needed to determine the rupture forces and the Hessian matrices for the JEDI analysis necessitate the use of a cost-efficient density functional. The def2-TZVP 12 basis set was used in all calculations for all atoms to provide a reliable description of the Fe−S bonds and the hydrogen bonds.
The performance of the density functionals was assessed by comparison to the experimental data provided by Day et al., 13 who reported the bond lengths between Fe(III)/Fe(II) and the surrounding sulfur atoms in rubredoxin. Day et al. find that the bond lengths in the Fe(III) complex are shorter than in the Fe(II) complex, which agrees with chemical intuition. Also, those Fe−S bond lengths in which the sulfur atom is involved in two hydrogen bonds with neighboring amino acids are observed to be longer than if only one hydrogen bond is formed. Given that the protein environment has a significant influence on these structural parameters and only a small model system can be treated quantum chemically, it is unlikely that a single density functional can reproduce all experimental bond lengths. Nevertheless, it is desirable that the general trend is reproduced. Another crucial point is that the same density functional shall be used for all investigated systems so that the results for Fe(III) and Fe(II) can be compared to each other.
To model the protein environment more realistically, a total of six formamide molecule were added to the model systems [Fe(III)(SCH 3 ) 4 ] − and [Fe(II)(SCH 3 ) 4 ] 2− , which form hydrogen bonds to the sulfur atoms. As a result, two sulfur atoms are involved in two hydrogen bonds (mimicking Cys5 and Cys38) and the other two sulfur atoms are involved in only one hydrogen bond (mimicking Cys8 and Cys41).
The results of the benchmark are given in Figures S1 and S2. As expected, the performance of the density functionals is very diverse and no single functional reproduces all experimental parameters correctly. Particularly in the case of Cys8 and Cys41, in which one hydrogen bond is involved, all tested functionals overestimate the experimentally observed Fe−S bond length. Nevertheless, most functionals reproduce the general trend that those Fe−S bonds in which the sulfur atom is involved in two hydrogen bonds are longer than in those cases in which only one hydrogen bond is formed. We found that the BP86VWN functional offers an attractive compromise between agreement with the experiment and computational cost, so we chose this functional for calculations reported in the paper. Moreover, BP86VWN was used successfully before to reproduce experimental Fe−S bond lengths in model complexes. 14 It is noted in passing that the application of a solvent model for water (C-PCM, 15,16 = 78) generally leads to a deterioration of the agreement with the experimental data ( Figures S3 and S4), so we ran all our calculations in the gas phase.     . Formamide molecules were added to the sulfur atoms in the order S1, S2, S3, S4, S1, and S2. Lines were included to guide the eye. with no hydrogen bonds are given in Tables S1 and S2. The numbering scheme is the same as in the main paper.  The mechanical anisotropy of the [Fe(III)(SCH 3 ) 4 ] − pseudo-tetrahedron and the importance of the bendings for the rupture process can be observed in the progression of the harmonic strain energies stored in the different internal coordinates with increasing stretching force ( Figure S11A). The bond angle bendings Fe−S1−C1 and Fe−S2−C2 store most strain energy throughout the entire stretching coordinate, signifying that they are softer than the central bending S1−Fe−S2. Hence, in the static strain analysis, the bond angles in the FeS 4 unit store more strain energy than in the dynamic calculations, signifying that different points of view on the rupture process of rubredoxin are provided by the two approaches.
The strain in the Fe−S bonds strongly increases with increasing force. This effect is more pronounced for the Fe−S1 bond as a result of the observed anisotropy and preconditions this bond for rupture. Not surprisingly, the C−S bonds are stronger and do not store as much strain energy as the Fe−S bonds.

CHELPG atomic charges for MD simulations
The atomic charges of the rubredoxin active site, comprising four deprotonated cysteine residues and the Fe(III) ion were re-parametrized by means of quantum-chemical calculations. The procedure included a geometry optimization of the minimal model system (see main text), representing the side chains of the respective amino acids, at the PBE 8 /6-31G* 19 level of theory as implemented in the ORCA program package. 20 Afterwards, CHELPG charges 21 were calculated and assigned to the corresponding CHARMM atom types for subsequent MD simulations. The used parameters are summarized in Table S5.