Michael Devereux

The importance of multipole moments when describing water and hydrated amino acid cluster geometry

With a view to protein hydration modelling, optimized geometries of pure water clusters, hydrated serine and hydrated tyrosine clusters are compared systematically. Geometries predicted by multipole models according to the theory of Quantum Chemical Topology and by point charge models are contrasted with ab initio geometries obtained at the B3LYP/aug-cc-pVDZ level of theory. The performance of popular point charge models such as AMBER, CHARMM, OPLS, MMFF, TAFF and TIP4P is scrutinized.

Force Field Optimization using Dynamics and Ensemble Averaged Data : Vibrational Spectra and Relaxation in Bound MbCO

Force field parameters are ingredients for realistic atomistic simulations of gas- and condensed-phase systems. Here we discuss the effect of including averaged data from explicit MD simulations in optimizing potential energy functions. It is shown that vibrational frequencies (FeC and CO stretch and FeCO bend) and CO vibrational relaxation times ((v = 1) --> (v = 0) (T(10)) and (v = 2) --> (v = 1) (T(21))) in the active site of CO-bound myoglobin (MbCO) can be well represented with a single set of force field parameters. It is further demonstrated that parameters fitted in a subsystem of MbCO comprising the CO ligand, heme group and proximal histidine, are transferable to investigating the full protein and to providing quantitatively correct results. In particular, it is possible to calculate the CO and FeC stretch and the FeCO bending frequency to within approximately 5%; the relaxation time of the first vibrationally excited state including quantum corrections of T(10) approximately 25 ps is calculated close to the experimental value (17 ps), and the ratio T(10)/T(21) approximately 2 agrees favorably with experimental estimates. In contrast, following the more traditional approach of fitting frequencies from analyzing the Hessian matrix leads to a force field that captures frequencies correctly but not relaxation of vibrations.

Temperature dependence of the heat diffusivity of proteins.

In a combined experimental-theoretical study, we investigated the transport of vibrational energy from the surrounding solvent into the interior of a heme protein, the sperm whale myoglobin double mutant L29W-S108L, and its dependence on temperature from 20 to 70 K. The hindered libration of a CO molecule that is not covalently bound to any part of the protein but is trapped in one of its binding pockets (the Xe4 pocket) was used as the local thermometer. Energy was deposited into the solvent by IR excitation. Experimentally, the energy transfer rate increased from (30 ps)(-1) at 20 K to (8 ps)(-1) at 70 K. This temperature trend is opposite to what is expected, assuming that the mechanism of heat transport is similar to that in glasses. In order to elucidate the mechanism and its temperature dependence, nonequilibrium molecular dynamics (MD) simulations were performed, which, however, predicted an essentially temperature-independent rate of vibrational energy flow. We tentatively conclude that the MD potentials overestimate the coupling between the protein and the CO molecule, which appears to be the rate-limiting step in the real system at low temperatures. Assuming that this coupling is anharmonic in nature, the observed temperature trend can readily be explained.

Application of multipolar charge models and molecular dynamics simulations to study stark shifts in inhomogeneous electric fields.

Atomic multipole moments are used to investigate vibrational frequency shifts of CO and H(2) in uniform and inhomogeneous electric fields using ab initio calculations and Molecular Dynamics (MD) simulations. The importance of using atomic multipole moments that can accurately represent both molecular electrostatics and the vibrational response of the molecule to changes in the local electric field is highlighted. The vibrational response of CO to applied uniform and inhomogeneous electric fields is examined using Density Functional Theory calculations for a range of test fields, and the results are used to assess the performance of different atomic multipole models. In uniform fields, the calculated Stark tuning rates of Deltamu = 0.52 cm(-1)/(MV/cm) (DFT), Deltamu = 0.55 cm(-1)/(MV/cm) (fluctuating three-point charge model), and Deltamu = 0.64 cm(-1)/(MV/cm) (Multipole model up to octupole), compare favorably with the experimentally measured value of 0.67 cm(-1)/(MV/cm). For H(2), which has no permanent dipole moment, CCSD(T) calculations demonstrate the importance of bond-weakening effects in force fields in response to the applied inhomogeneous electric field. Finally, CO in hexagonal ice is considered as a test system to highlight the performance of selected multipolar models in MD simulations. The approach discussed here can be applied to calibrate a range of multipolar charge models for diatomic probes, with applications to interpret Stark spectroscopy measurements in protein active sites.

The quantum topological electrostatic potential as a probe for functional group transferability

The electrostatic potential can be used as an appropriate and convenient indicator of how transferable an atom or functional group is between two molecules. Quantum-chemical topology (QCT) is used to define the electron density of a molecular fragment and the electrostatic potential it generates. The potential generated on a grid by the terminal aldehyde group of the biomolecule retinal is compared with the corresponding aldehyde group in smaller molecules derived from retinal. The terminal amino group in the free amino acid lysine was treated in a similar fashion. Each molecule is geometry-optimized by an ab initio calculation at B3LYP/6-311G+(2d,p)//HF/6-31G(d) level. The amino group in lysine is very little influenced by any part of the molecule further than two C atoms away. However, the aldehyde group in retinal is influenced by molecular fragments six C atoms away. This dramatic disparity is ascribed to the difference in saturation in the carbon chains; retinal contains a conjugated hydrocarbon chain but lysine an aliphatic one.

Toward an ab initio fragment database for bioisosterism: Dependence of QCT properties on level of theory, conformation, and chemical environment

The goal of this work is to assess the scope and suitability of atomic and bond properties for use in a bioisostere fragment database. This database will contain fragment descriptors that can be used to represent portions of larger molecules and similarity in properties between fragments, which will then be used to find bioisosteric replacements in future work. Seventeen common organic fragments relevant to drug design featured as “linker groups” that were capped by two terminal groups. Each terminal group could be one of the set of 12 possible sets: 10 aromatic heterocycles, a phenyl ring, or an ethyl. This enabled a systematic investigation of the chemical environment, enriched with conformational flexibility within the linker group, for a total of 307 different atoms. Five different levels of theory were investigated. This work paves the way to the construction of a quantum mechanical bioisosteric fragment database, for which transferability of stored fragment properties is of fundamental importance.

Anharmonic coupling in molecular dynamics simulations of ligand vibrational relaxation in bound carbonmonoxy myoglobin

Vibrational relaxation of CO bound to myoglobin (MbCO) following photoexcitation is investigated using nonequilibrium molecular dynamics (MD) simulations. It is found that harmonic potential energy functions for bond vibrations are not suited to simultaneously and accurately describe vibrational de-excitation and the vibrational spectroscopy of the bound ligand. Only when anharmonic (e.g. Morse) potentials are introduced for both the C-O and the adjacent Fe-C(CO) bonds to allow anharmonic coupling, rapid (tens of ps) relaxation of the vibrationally excited CO is possible. To capture both relaxation and vibrational spectroscopy, the parameters of the potential energy functions are fitted by an interactive, nonlinear least-squares procedure using averages over multiple MD trajectories. The sensitivity of cooling rate to the difference in vibrational frequency between coupled modes is demonstrated. Potential cooling mechanisms are suggested, based on the sensitivity of the CO relaxation rate to changes in the force field parameters of local degrees of freedom. Accounting for quantum correction leads to relaxation rates around 20 ps, in good agreement with experiment. Finally, the importance of electronic effects is explored by fitting a 2D potential energy surface to ab initio data to describe the strengthening and weakening of the CO bond as a function of Fe-C(CO) bond length, and vice versa.

Structural assignment of spectra by characterization of conformational substates in bound MbCO

Residue motions of the distal heme pocket and bound CO ligand of carbonmonoxy Myoglobin are studied using a combination of molecular dynamics simulations and quantum chemical methods. Using mixed quantum mechanics/molecular mechanics calculations together with sampling from molecular dynamics simulations (QM/MM(MD)), the experimentally observed spectroscopic A(0) and A(1) substates of the bound CO ligand are assigned to the open and closed conformation of His(64) and the His(epsilon)(64) tautomer, respectively. Several previously proposed origins of the A(3) substate, including rotamers of the doubly protonated His(64)H(+) side chain, His(64)H(+) inside the distal pocket, and cooperative motions with Arg(45), are investigated with QM/MM(MD). However, the signatures of the calculated infrared spectra do not agree with the experimentally observed ones. For additional insight on this, extensive molecular dynamics simulations are used together with improved electrostatics for the bound ligand. A CO fluctuating charge model is developed to describe the ab initio dipole and quadrupole moments of the bound ligand. CO absorption spectra are then obtained directly from the dynamics simulations. Finally, the electrostatics of the heme pocket is examined in detail in an attempt to determine the structural origins of the observed spectroscopic A-states from MD simulations. However, contrary to related simulations for unbound CO in myoglobin, the shifts and splittings for carbonmonoxy Myoglobin are generally small and difficult to relate to structural change. This suggests that coupling of the CO motion to other degrees of freedom, such as the Fe-CO stretching and bending, is important to correctly describe the dynamics of bound CO in myoglobin.

A supervised fitting approach to force field parametrization with application to the SIBFA polarizable force field

A supervised, semiautomated approach to force field parameter fitting is described and applied to the SIBFA polarizable force field. The I‐NoLLS interactive, nonlinear least squares fitting program is used as an engine for parameter refinement while keeping parameter values within a physical range. Interactive fitting is shown to avoid many of the stability problems that frequently afflict highly correlated, nonlinear fitting problems occurring in force field parametrizations. The method is used to obtain parameters for the H2O, formamide, and imidazole molecular fragments and their complexes with the Mg2+ cation. Reference data obtained from ab initio calculations using an auc‐cc‐pVTZ basis set exploit advances in modern computer hardware to provide a more accurate parametrization of SIBFA than has previously been available.

Toward force fields for atomistic simulations of iridium‐containing complexes

The structural and energetic characterization of metal complexes is important in catalysis and photochemical applications. Unraveling their modes‐of‐action can be greatly assisted by computation, which typically is restricted to computationally demanding methods including electronic structure calculations with density functional theory. Here, we present an empirical force field based on valence bond theory applicable to a range of octahedral Ir(III) complexes with different coordinating ligands, including iridium complexes with a chiral P,N ligand. Using an approach applicable to metal‐containing complexes in general, it is shown that with one common parametrization 85% of the 116 diastereomers—all within 21 kcal/mol of the lowest energy conformation of each series—can be correctly ranked. For neutral complexes, all diastereomers are ranked correctly. This helps to identify the most relevant diastereomers which, if necessary, can be further investigated by more demanding computational methods. Furthermore, if one specific complex is considered, the root mean square deviation between reference data from electronic structure calculations and the force field is ≈1  kcal/mol  . This, together with the possibility to carry out explicit simulations in solution paves the way for an atomistic understanding of iridium‐containing complexes in catalysis.