Efficient Calculation of Enzyme Reaction Free Energy Profiles Using a Hybrid Differential Relaxation Algorithm

Application to Mycobacterial Zinc Hydrolases

Juan Manuel Romero*,1; Mariano Martin†,1; Claudia Lilián Ramirez*; Victoria Gisel Dumas*,†; Marcelo Adrián Marti†,‡,2    * Instituto de Química Física de los Materiales Medio Ambiente y Energía (INQUIMAE), UBA-CONICET, Facultad de Ciencias Exactas y Naturales, Universidad de Buenos Aires, Ciudad Universitaria, Buenos Aires, Argentina
† Departamento de Química Biológica, Facultad de Ciencias Exactas y Naturales, Universidad de Buenos Aires, Ciudad Universitaria, Buenos Aires, Argentina
‡ Instituto de Química Biológica de la Facultad de Ciencias Exactas y Naturales (IQUIBICEN) CONICET, Ciudad Universitaria, Buenos Aires, Argentina
2 Corresponding author: email address: marcelo@qi.fcen.uba.ar
1 Both authors contributed equally to this work.

1 Introduction

1.1 Free Energy Profiles of Enzymatic Reactions

Understanding the origins of enzyme catalysis is one of the main challenges in current biochemical research, with potential impact in the fields of protein evolution, protein engineering, and drug development, among others. To uncover the underlying physicochemical reasons of a given protein-catalyzed reaction and relate it with the protein's structure and dynamics at a residue-based level, it is of primordial relevance to determine the corresponding reaction free energy profile (FEP). The FEP, also sometimes referred as the potential of mean force, can not only be directly related to experimentally determined properties, such as reaction rates and equilibrium constants, but also constitutes a proof of concept of the reaction dynamics, defined as the motion of atoms, with the concomitant breaking and forming of chemical bonds during the reactive process (Kamerlin & Warshel, 2010).

To be able to describe a (bio)chemical reaction, the system needs to be treated using quantum mechanics (QM) methods, such as density functional theory (DFT). QM methods, however, are computationally expensive, and describing systems larger than a few hundreds of atoms, with that level of accuracy, is very difficult and an active research field (Ferrer et al., 2011). For enzyme-based reactions, those which require consideration of the whole protein-solvent environment—tens of thousands of atoms, hybrid quantum mechanical/molecular mechanical (QM/MM) schemes are the best choice (Warshel & Levitt, 1976). QM/MM methods allow accurate description of the chemical events that take place in the enzyme active site modeled at the QM level, while treating the rest of the system using a less expensive classical force-field based level of theory. Key to these methods is the coupling between the QM and MM regions, which must properly describe the environment electrostatic as well as the steric effect on the reactive subsystem. Although they have been extensively and successfully used during the last decade to study reaction mechanisms, the configurational sampling required to obtain an accurate FEP is still an open challenge faced by QM/MM methods and their application to enzyme (Kamerlin & Warshel, 2010; Van der Kamp & Mulholland, 2013).

1.1.1 MSMD and Jarzynski's Relationship

In plain molecular dynamics (MD), the system under study is usually unable to cross moderate high barriers ( kT), such as those presented by enzyme reactions, thus remaining trapped in the initial (reactant) state, unless driven up the hill. There are several enhanced or biased sampling strategies that nonetheless allow determination of FEP with their associated barriers. The most old and possibly common is the umbrella sampling strategy (Leach, 2001). More recently, newer strategies, such as metadynamics (Laio & Parrinello, 2002), adaptive biasing force (Hénin & Chipot, 2004), free energy perturbation (Pohorille, Jarzynski, & Chipot, 2010), and orthogonal space random walk (Zheng, Chen, & Yang, 2008), have also been developed. One of the most easily implemented, with wide applicability, solid statistical thermodynamic background and multiple computational and experimental validations is the multiple steered molecular dynamics (MSMD) strategy combined with Jarzynski's relationship (JR) (Jarzynski, 1997; Liphardt, Dumont, Smith, Tinoco, & Bustamante, 2002; Park, Khalili-Araghi, Tajkhorshid, & Schulten, 2003). In MSMD (as shown schematically in Fig. 1), the system is driven “multiple” times along the selected reaction coordinate (lambda) under nonequilibrium conditions, by applying an external force onto the reaction coordinate. For each individual trajectory, the work performed by the force is determined (Wi(λ)). Finally, multiple works are exponentially average in JR (Eq. 1) to obtain the corresponding FEP.

λ=−β−1lne−βWiλ

where G(λ) represents the FEP as a function of the reaction coordinate, β = 1/kBT, where kB is Boltzmann constant and T is the system temperature. The brackets represent the average of the function within them. In JR, the exponential average is computed on the work values distribution, and the more narrow the distribution, the more accurate the average. The width of the work distribution is directly proportional to the pulling speed. Thus, in the practice, the main drawback with JR is that in order to obtain a well-converged average, and thus an accurate FEP, either very large number of trajectories and/or very low pulling speeds are needed. These facts result in a high computational cost and sometimes in an insurmountable problem that prevents accurate convergence of the FEP (Pohorille et al., 2010; Xiong, Crespo, Marti, Estrin, & Roitberg, 2006).

1.1.2 Hybrid Differential Relaxation Algorithm

To overcome the above-mentioned difficulties and reduce the overall computational cost of MSMD-JR strategy for FEP determination, we have recently developed a hybrid differential relaxation algorithm (HyDRA) (Ramírez, Zeida, Jara, Roitberg, & Martí, 2014). This scheme allows faster equilibration of the classical environment during the steering process that drives the QM system along the reaction under study. The differential relaxation strategy, inspired in the multiple time step schemes developed earlier, takes advantage of the less expensive calculation of the classical environment and performs multiple pure classical relaxation steps for each QM perturbation step. This allows a better relaxation of the whole system, resulting in closer to equilibrium steering trajectories, with more narrow...