Challenges in Modelling Metalloenzymes
513
([Fe 2 S 2 (His) 2 (Cys) 2 ]) as well as its hydrogen-bonded partners (with dielectric constant ε 1), rest of the protein (ε 4) and the solvent (ε 80). DFT gas-phase
calculations for the cluster model allowed establishing geometries and energies of
the active site cluster with various numbers of electrons and protons bound in order
to account for every combination of protonation and metal oxidation state. In this
case, 12 different variants were considered resulting from combination of 4 possible protonation states of two His ligands and three different electronic states of the
Fe 2 S 2 cluster (fully oxidized, with Fe coordinated by His reduced, with Fe coordinated by Cys reduced). The protein and solvent environment were modelled with a
dielectric continuum model with embedded charges representing the protein utilizing
macroscopic electrostatics with atomic detail (MEAD) model [51–54]. The calculations were iterative and consisted of the following main steps. First, ESP atomic
charges on the cluster were fitted with the CHELPG protocol. Next, the Poisson
or Poisson–Boltzmann equation was solved yielding the reaction field. The protein
field was computed in the same manner using partial charges from the surrounding
protein. Subsequently, the protein and reaction fields were added to the cluster Hamiltonian and the DFT calculations were repeated until convergence was reached. The
numerical solution of the Poisson–Boltzmann equation (conducted with the program
MEAD [52]) was used to calculate protonation and interaction energies and Monte
Carlo simulation for protonation state sampling [55]. Finally, the microscopic protonation equilibrium constants were calculated by methods described by Li et al. [56,
57] using a thermodynamic cycle that represents deprotonation in gas/protein/solvent
phase and solvation process. The application of such approach allowed Ullmann et al.
to predict apparent pK a values of the oxidized Rieske centre as 6.9 and 8.8, i.e. values fairly close to those obtained from the experiment (7.5 and 9.2, respectively).
However, it turned out that these pK a values cannot be associated with deprotonations of individual His residues as both of them titrated over a wide range of pH. In
accordance with experimental data, the computational results demonstrated that it is
the Fe ion coordinated by His that receives electron upon the site reduction and that
such redox change results in shift of apparent pK a to 11.3 and 12.8, i.e. above 10.0,
as was also found experimentally. Computational results showed that after reduction
of the FeS centre, the His161 exhibits a higher proton affinity than His141. As a
result, the His161 residue will be able to transfer a proton from ubiquinol faster than
His141. Finally, the authors considered the possibility of coupling between pK a of
the second shell residues and the redox state of the Rieske centre. However, it turned
out that for most residues, pK a are not strongly affected by the redox state of the
cluster, with the only exception of Tyr157, whose pK 1/2 shifted by 0.6 pH unit upon
cluster reduction.
7 Molecular Dynamics Simulations
Molecular dynamics simulations are widely recognized as a very powerful tool that
allows one to study various processes. For example, it is frequently used to: (a) “re-
513
([Fe 2 S 2 (His) 2 (Cys) 2 ]) as well as its hydrogen-bonded partners (with dielectric constant ε 1), rest of the protein (ε 4) and the solvent (ε 80). DFT gas-phase
calculations for the cluster model allowed establishing geometries and energies of
the active site cluster with various numbers of electrons and protons bound in order
to account for every combination of protonation and metal oxidation state. In this
case, 12 different variants were considered resulting from combination of 4 possible protonation states of two His ligands and three different electronic states of the
Fe 2 S 2 cluster (fully oxidized, with Fe coordinated by His reduced, with Fe coordinated by Cys reduced). The protein and solvent environment were modelled with a
dielectric continuum model with embedded charges representing the protein utilizing
macroscopic electrostatics with atomic detail (MEAD) model [51–54]. The calculations were iterative and consisted of the following main steps. First, ESP atomic
charges on the cluster were fitted with the CHELPG protocol. Next, the Poisson
or Poisson–Boltzmann equation was solved yielding the reaction field. The protein
field was computed in the same manner using partial charges from the surrounding
protein. Subsequently, the protein and reaction fields were added to the cluster Hamiltonian and the DFT calculations were repeated until convergence was reached. The
numerical solution of the Poisson–Boltzmann equation (conducted with the program
MEAD [52]) was used to calculate protonation and interaction energies and Monte
Carlo simulation for protonation state sampling [55]. Finally, the microscopic protonation equilibrium constants were calculated by methods described by Li et al. [56,
57] using a thermodynamic cycle that represents deprotonation in gas/protein/solvent
phase and solvation process. The application of such approach allowed Ullmann et al.
to predict apparent pK a values of the oxidized Rieske centre as 6.9 and 8.8, i.e. values fairly close to those obtained from the experiment (7.5 and 9.2, respectively).
However, it turned out that these pK a values cannot be associated with deprotonations of individual His residues as both of them titrated over a wide range of pH. In
accordance with experimental data, the computational results demonstrated that it is
the Fe ion coordinated by His that receives electron upon the site reduction and that
such redox change results in shift of apparent pK a to 11.3 and 12.8, i.e. above 10.0,
as was also found experimentally. Computational results showed that after reduction
of the FeS centre, the His161 exhibits a higher proton affinity than His141. As a
result, the His161 residue will be able to transfer a proton from ubiquinol faster than
His141. Finally, the authors considered the possibility of coupling between pK a of
the second shell residues and the redox state of the Rieske centre. However, it turned
out that for most residues, pK a are not strongly affected by the redox state of the
cluster, with the only exception of Tyr157, whose pK 1/2 shifted by 0.6 pH unit upon
cluster reduction.
7 Molecular Dynamics Simulations
Molecular dynamics simulations are widely recognized as a very powerful tool that
allows one to study various processes. For example, it is frequently used to: (a) “re-