Molecular Mechanics for Transition Metal Centers: From Coordination Complexes To Metalloproteins

Robert J. Deeth    Inorganic Computational Chemistry Group, Department of Chemistry, University of Warwick, Coventry CV4 7AL, United Kingdom

I Introduction

Transition metal (TM) systems present a fundamental dilemma for computational chemists. On the one hand, TM centers are often associated with relatively complicated electronic structures which appear to demand some form of quantum mechanical (QM) approach (1). On the other hand, all forms of QM are relatively compute intensive and are impractical for conformational searching, virtual high-throughput screening, or dynamics simulations since all these approaches may require many hundreds of thousands of individual calculations. Consequently, TM computational chemists tend to restrict themselves to smaller “model” systems with limited conformational freedom. This is particularly marked in bioinorganic chemistry where the calculation focuses on the “important” active site region but the bulk of the protein is not treated explicitly (2,3).

In contrast, those interested in purely “organic” systems have long enjoyed the advantages of “cheap,” classical molecular mechanics (MM) and molecular dynamics (MD) to study the entire molecular system including the surrounding solvent. However, conventional MM is not well suited to TM systems since it does not provide a general way of accounting for the important effects arising from the d electrons (4,5). In response, hybrid QM/MM methods have appeared (6). The metal center and its immediate environment is handled by a “high level” QM method, typically based on Density Functional Theory (DFT), with the rest of the system treated by MM. As the many technical difficulties of QM/MM have progressively been solved—most importantly how to couple the quantum region to the classical region—QM/MM has grown in popularity. However, the inclusion of any QM, even on a relatively small piece of the whole system, soon exacts a huge cost in execution time. Just a few minutes per calculation soon equates to years of CPU time. The only options are either to use thousands of computers or to develop a method which is as accurate as QM, but many orders of magnitude more efficient. We have taken the second path by augmenting MM with additional terms designed to provide a physically meaningful description of metal–ligand bonding and thus be able to emulate the behavior of more sophisticated, but expensive, QM methods. However, in order to put our model into perspective, we must first appreciate the nature of “conventional” MM and its shortcomings when applied to TM systems.

II Conventional Molecular Mechanics

Molecular mechanics in its simplest form expresses the total potential energy, Etot, as a sum of terms describing bond stretching, Estr, angle bending, Ebend, torsional twisting, Etor, and nonbonding interactions, Enb(1). The latter can include both van der Waals (vdW) interactions and, by assigning to each atom a partial atomic charge, electrostatics.

tot=∑Estr+∑Ebend+∑Etor+∑Enb

Each term in (1) is represented by a simple mathematical expression as exemplified in (2), where the k are appropriate force constants, θ are bond angles, τ are torsion angles, n is the torsional periodicity parameter, ϕ is the torsion offset, ρ are partial atomic charge, ε is the dielectric constant, A and B are Lennard-Jones vdW parameters, and the summations run over bonded atom pairs (ij), angle triples (ijk) and torsional quadruples (ijkl). The nonbonded terms are summed over the distances, dij, between unique atom pairs excluding bonded pairs and the atoms at either end of an angle triple. For the atoms at the ends of a torsion quadruple, the nonbonded term may be omitted or scaled.

tot=∑i,jkij(rij−r0,ij)2+∑i,j,kkijk(θijk−θ0,ijk)2+∑i,j,k,lkijkl[1+cos(nijklτ−ϕijkl)]+[∑i<jρiρjεdij∑i<j(Aijdij12−Bijdij6)]

These potential energy terms and their attendant empirical parameters define the force field (FF). More complicated FFs which use different and/or more complex functional forms are also possible. For example, the simple harmonic oscillator expression for bond stretching can be replaced by a Morse function, EMorse(3), or additional FF terms may be added such as the stretch-bend cross terms, Estb, (4) used in the Merck molecular force field (MMFF) (7–10) which may be useful for better describing vibrations and conformational energies.

Morse=D{1−ea(r−r0)}2−D

stb=∑i,j,k{kijk(rij−r0,ij)+kkji(rjk−r0,jk)}θ0,ijk

MM is a very successful model but it is clear from expressions such as (2) that the FF may comprise a very large number of parameters and since the quality of the FF will depend crucially on these parameters, developing a truly “universal” FF is an enormous (perhaps impossible) challenge (11). Consequently, FFs tend to be applicable to a specific class of molecular system such as small organic molecules (e.g., MMFF (7–10)), or large biomolecules like proteins and DNA (e.g., AMBER (12) or CHARMM (13)) for which the expressions in (2) are suited. Fortunately, these specific classes encompass an enormous amount of organic chemistry and biology plus the FFs are continually being revised to increase their applicability.

The computational efficiency of a FF approach also enables simulations of dynamical behavior—molecular dynamics (MD). In MD, the classical equations of motion for a system of N atoms are solved to generate a search in phase space, or trajectory, under specified thermodynamic conditions (e.g., constant temperature or constant pressure). The trajectory provides configurational and momentum information for each atom from which thermodynamic properties such as the free energy, or time-dependent properties such as diffusion coefficients, can be calculated.

MD introduces a time dependence but to avoid numerical instability, the time step δt needs to be quite small (∼1 fs) which places strong limitations on the total simulation time. Nevertheless, compared to QM-based schemes such as Carr-Parinello MD, classical MD is currently the only really viable method for computing the dynamical behavior of large systems like biomolecules for any appreciable length of...