MM Molecular Dynamics Sampling in Explicit

ARTICLE pubs.acs.org/JPCB

Enhancing QM/MM Molecular Dynamics Sampling in Explicit Environments via an Orthogonal-Space-Random-Walk-Based Strategy Donghong Min,† Mengen Chen,† Lianqing Zheng,† Yonghao Jin,‡ Martin A. Schwartz,‡ Qing-Xiang Amy Sang,‡ and Wei Yang*,†,‡ †

Institute of Molecular Biophysics and ‡Department of Chemistry and Biochemistry, Florida State University, Tallahassee, Florida 32306, United States ABSTRACT: Accurate prediction of molecular conformations in explicit environments, such as aqueous solution and protein interiors, can facilitate our understanding of various molecular recognition processes. Most computational approaches are limited as a result of their compromised choices between the underlying energy model and the sampling length. Taking advantage of a recent second-order generalized ensemble scheme [e.g., the orthogonal space random walk (OSRW) strategy], which can synergistically accelerate the motion of a focused region and its coupled environmental response, we are presenting a QM/MM (combined quantum mechanical/molecular mechanical)-based molecular dynamics sampling technique to explore molecular conformational landscapes in explicit environments. The present QM/MM potential scaling-based OSRW sampling scheme is employed to study the binding of DMSO to the FKBP12 protein, the conformation distribution of a novel mercaptosulfonamide inhibitor in aqueous solution, and its binding poses in zinc-containing matrix metalloproteinase-9 (MMP-9). As demonstrated, the present QM/MM second-order generalized ensemble sampling technique enables feasible usage of the QM/MM model to sample molecular conformations in condensed environments.

’ INTRODUCTION Accurate prediction of molecular conformation in explicit environments, such as aqueous solution and protein interiors, has been an important task in computational chemistry and computational biophysics; such capability is expected to further catalyze our understanding of molecular structure changes and proteinligand recognitions, etc. For the above purpose, two technical aspects are essential: the accuracy of employed energy potentials and the adequacy of sampling length. To acquire sufficient sampling, molecular mechanical (MM) potentials are usually employed. As the matter of fact, applying higher-level potentials, such as combined quantum mechanical/molecular mechanical (QM/MM) potentials,15 has been considered as a natural next step in evaluating proteinligand interactions.69 Thereby, a target molecular region can be described quantum mechanically, and its surroundings can be treated classically, as illustrated in the following energy function form: _

_

Uo ðXt , Xe Þ ¼ ÆΦc jHt ðXt Þ þ Ht=e ðXt , Xe ÞjΦs æQ M=MM þ Ut=e ðXt , Xe ÞMM þ Uee ðXe ÞMM

ð1Þ

where t denotes the target region and e stands for the environment region; ÆΦc|H kt(Xt) þ Hkt/e(Xt, Xe)|ΦsæQM/MM is the QM/MM Hamiltonian that describes all the internal interactions r 2011 American Chemical Society

within the target region and the electrostatic interactions between the target region and the environment region; Ut/e(Xt, Xe)MM is the MM component that describes the remaining interactions between the target region and the environment region; and Uee(Xe)MM is the MM component that describes the internal interactions within the environment region. Despite possibly better accuracy, high computing costs have discouraged QM/MM models from being commonly applied to treat molecular species that are embedded in condensed environments. Facing the above challenge, developing efficient sampling methods for QM/MM-based simulations is a viable choice; then with commonly affordable simulation length, longer-time-scale conformational events could be repetitively captured. To realize accelerated QM/MM MD simulations, especially when no prior knowledge is known on a target system, a direct option is to implement QM/MM potentials into classical first-order generalized ensemble frameworks, such as the multicanonical ensemble method10 and the temperature replica exchange method.11,12 Considering the fact that a target region for sampling is only a localized portion of the whole system, it is desirable to focus the Received: October 1, 2010 Revised: February 28, 2011 Published: March 21, 2011 3924

dx.doi.org/10.1021/jp109454q | J. Phys. Chem. B 2011, 115, 3924–3935

The Journal of Physical Chemistry B

ARTICLE

sampling treatment on an essential conformational space; then the diffusion sampling problem,13,14 which reveals the fact that possible enlargement of diffusion sampling space (where freeenergy surface flattening is realized) causes abolishment of overall sampling efficiency, can be avoided. Along this thought, our group has developed two general QM/MM enhanced-sampling strategies: the hybrid Hamiltonian replica exchange method11 and the simulated-scaling-based method,15 the latter of which was designed to realize a random walk in the scaling parameter (the potential scaling parameter is also called local effective temperature for the fact that when the system potential energy is scaled, the scaling parameter plays the same role as the inverse of the system temperature15,16) space, based on the following modified energy potential: Um ðXt , Xe , λÞ ¼ Uo ðXt , Xe , λÞ þ fm ðλÞ _

_

¼ λÆΦc jHt ðXt Þ þ Ht=e ðXt , Xe ÞjΦs æQM=MM þ λUtedihedral ðXt , Xe ÞMM bond, angle, improper þ ð1 λÞUtt ðXt ÞMM bond, angle, improper, vdw þ Ute ðXt , Xe ÞMM þ Uee ðXe ÞMM þ fm ðλÞ

ð2Þ

Commonly, in a constructed force-scaling potential [Uo(Xt, Xe, λ) in eq 2], (1) the energy terms that determine the conformations of kt(Xt) þ the target region, including the QM/MM term ÆΦc|H H kt/e(Xt, Xe)|ΦsæQM/MM and possible MM dihedral terms bet(Xt, Xe)MM, ween the target and the environment regions Udihedral te are scaled by a dynamic variable λ; (2) to avoid the activation of involved high-frequency vibrations and maintain structural integrity,13 the classical bond, angle, improper, and van der Waals energy terms between the target and the environment (Xt, Xe)MM are not scaled; and (3) regions Ubond,angle,improper,vdw te to counter the activation of internal high-frequency vibrations brought in by the QM/MM energy scaling, the corresponding (Xt)MM are proportionally added MM terms Ubond,angle,improper tt [with a scaling parameter of (1 λ)] as a chaperone component.17 It is noted that the employment of chaperone energy terms is to ensure the acceleration of conformational processes rather than chemical processes.11 The above design of a force-scaling potential follows two principles, which have been discussed in our earlier works: (1) minimizing high-frequency motion activation13,15 and (2) maintaining structural integrity for robust QM self-consistent field calculations.11,15,17,18 Furthermore, to realize a scaling parameter space random walk, a biasing energy potential fm(λ) [λ is constrained in a range between a lower value (than 1) and a larger value, such as “1”] is employed; specifically, the target function form of fm(λ) is Go(λ), which is the negative of the λ-dependent free-energy profile corresponding to the canonical ensemble with Uo(λ) as the energy function. Through the above treatment, when the system visits a lower scaling parameter state, target torsion motions can be effectively activated, and local energy minima can be more easily escaped; then, with the increase of the scaling parameter, newly generated conformations can be brought back to the target scaling parameter state [when λ = 1, Uo(Xt, Xe, λ) in eq 3 is equal to Uo(Xt, Xe) in eq 2]. It should be noted that following our original proposals in the Hamiltonian-scaling-based (also called the energy-scaling or the force-scaling) methods, the coupling between the system dynamics and the move of the scaling

parameter λ can be achieved via either the hybrid Monte Carlo strategy19 or the λ-dynamics strategy.20 The above QM/MM “simulated scaling” strategy falls into the classical generalized ensemble regime, where random walks along a predefined order parameter are realized to activate target events. As discussed in our recent works,21,22 the first-order generalized ensemble treatment generally suffers from the “Hamiltonian lagging” issue, which reveals the fact that necessary environmental relaxations fall behind activated central events; this sampling issue is particularly severe when a local event for the selective sampling treatment is embedded in complex explicit environments. Following the second-order generalized ensemble strategy21,22 [e.g., the orthogonal space random walk (OSRW) method], which can synergistically accelerate the motions of a focused region and their strongly coupled environmental responses, we are reporting an OSRW update of the above QM/ MM potential-scaling-based generalized ensemble scheme to enhance QM/MM molecular dynamics sampling in explicit environments. As demonstrated by the present studies on the DMSO binding to the FKBP12 protein, the conformational distribution of a novel mercaptosulfonamide inhibitor in aqueous solution, and its binding structure in zinc-containing matrix metalloproteinase-9 (MMP-9), the present QM/MM potential-scaling-based OSRW strategy can enable feasible usage of the QM/MM model for the sampling of molecular conformations in condensed environments.

II. THEORETICAL DESIGN The OSRW strategy was designed to simultaneously overcome explicit free-energy barriers along the preset order parameter and “hidden” free-energy barriers, the latter of which usually hinder the collective reorganization of condensed environments that strongly couple with central conformational dynamics. Specifically, on the top of the first-order generalized ensemble treatment [namely, with a biasing energy term fm(λ), as shown in eq 1], an additional biasing energy term, Fm[λ, h(λ)], where h(λ) represents a robust order parameter that describes the direction of slow structural relaxations at each λ state, is included; thereby, required environmental sampling can be synchronously accelerated. In OSRW, (∂Uo(Xt, Xe, λ))/∂λ is applied as such a function form. As discussed in our earlier works,22 this order parameter function is particularly robust for systems in which orthogonal space structural transitions strongly couple with the move of the order parameter, λ. Following the simulated scaling strategy as depicted in eq 2, the QM/MM potentialscaling-based OSRW target energy function can be written as Um ðXt , Xe , λÞ ¼ Uo ðXt , Xe , λÞ þ fm ðλÞ DUo ðXt , Xe , λÞ þ Fm λ, Dλ

ð3Þ

where Uo(Xt, Xe, λ) can be the same as that in eq 2 or can be redesigned, as long as it follows the constraint Uo(Xt, Xe, 1) = Uo(Xt, Xe); the target function of fm(λ) is Go(λ) [Go(λ) is the λdependent free-energy profile in the canonical ensemble with Uo(Xt, Xe, λ) as the potential energy function]; and the target function form of Fm[λ, (∂Uo(Xt, Xe, λ)/∂λ)] is Go0 [λ, (∂Uo(Xt, Xe, λ)/∂λ)] {Go0 [λ, (∂Uo(Xt, Xe, λ)/∂λ)] is the [λ, (∂Uo(Xt, Xe, λ)/∂λ)] space free-energy profile corresponding to the canonical ensemble with Uo(Xt, Xe, λ) Go(λ) as the potential energy function}. The same as that in our original OSRW proposal, to ensure a continuous description of various λ states, the coupling of the scaling parameter 3925

dx.doi.org/10.1021/jp109454q |J. Phys. Chem. B 2011, 115, 3924–3935


ARTICLE

λ is realized via the λ-dynamics method20 (in the extendedLagrangian scheme) (e.g. the scaling parameter λ is treated as a one-dimension particle moving between two preset boundary values, for instance a lower value and 1, and its dynamics can be thermalized via an independent Langevin reservoir). Two technical issues need to be considered in generalized ensemble simulations: (1) how to perform an adaptive recursion simulation to efficiently and robustly obtain a set of biasing energy functions, fm(λ) and Fm[λ, (∂Uo(Xt, Xe, λ)/∂λ)]; (2) how to recover canonical ensemble properties using the samples collected in a production equilibrium simulation, which is based on a fixed set of biasing energy functions, fm(λ) and Fm[λ, (∂Uo(Xt, Xe, λ)/∂λ)]. II.A. OSRW Recursion Simulations To Obtain Biasing Energy Terms and To Perform Initial Structural Refinement and Optimization. Notably, in addition to its major function of obtaining biasing energy terms, the recursion step simulation can also serve as a tool for initial structural refinement and optimization.22 The second generalized ensemble method requires two recursion components: a “recursion kernel” procedure to adaptively obtain Fm[λ, (∂Uo(Xt, Xe, λ)/∂λ)] and a corresponding “recursion slave” relation to derive fm(λ) from Fm[λ, (∂Uo(Xt, Xe, λ)/ ∂λ)]. Specifically, on the “recursion kernel”, we employ a similar recursion strategy as the one in the metadynamics method;23 Fm[λ, (∂Uo(Xt, Xe, λ)/∂λ)] is obtained as the sum of relatively small Gaussian-shaped repulsive potentials: 0 2 1 B λ λðti Þ C B C C h expB B C 2 2w1 B C @ A

its generalized force space can be estimated as Fm[λ0 , (∂Uo(Xt, Xe, λ)/∂λ)|λ ]; correspondingly, the energy derivative probability should be proportional to exp{βFm[λ0 , (∂Uo(Xt, Xe, λ)/∂λ)|λ ]}. Then, the free-energy derivative, as the energy derivative ensemble average at the state λ0 , can be estimated as dGo DUo ðXt , Xe , λÞ ¼ Dλ dλ 0 λ0

0 1 DU ðX , X , λÞ DU ðX , X , λÞ 2 o t e B o t e ðt i Þ C B C Dλ Dλ C expB B C 2w2 2 B C @ A

in which the meaning of each energy term was described in the portion following eq 2. II.B. Data Analysis Based on the OSRW Production Simulation Results. After biasing potential functions Fm[λ, (∂Uo(Xt, Xe, λ)/∂λ)] and fm(λ) are generated from an OSRW recursion simulation, generalized ensemble simulations can be performed on the resulting energy potential as described in eq 3 and with the optimized structure from the OSRW recursion simulation as the starting structure. According to the umbrella sampling reweighting scheme,25 the relative weight of the sample collected at the simulation time tj should be Uo[Xt, Xe, λ(ti)] Uo[Xt, Xe, 1) þ fm [λ(ti)] þ Fm[λ(ti), ((∂Uo(Xt, Xe, λ)/∂λ)(ti))], where as described earlier, Uo(Xt, Xe, 1) is the potential of the target chemical state. Therefore, the potential of mean force along any trial set of reaction coordinates, B, s can be obtained according to the following equation: 8 < δ½ B s B s ðtj Þ Gð B s Þ ¼ RT ln :

ð4Þ

which is centered at [λ(ti), ((∂Uo(Xt, Xe, λ)/∂λ)(ti))] and thereby discourages the system from often-visited configurations. Here, ti represents the time-step when the ith update is scheduled. With this procedure repeated, the overall biasing potential 0 2 1 B λ λð t i Þ C B C DUo ðXt , Xe , λÞ C ¼ h expB Fm λ, B C 2 Dλ 2w B C 1 ti @ A

∑

0 1 DU ðX , X , λÞ DU ðX , X , λÞ 2 o t e B o t e ðti Þ C B C Dλ Dλ C expB B C 2w2 2 B C @ A

0

0

λ

DUo DUo ðXt , Xe , λÞ exp βFm λ, δðλ λ0 Þ Dλ DUo =Dλ Dλ

¼ Z DUo ðXt , Xe , λÞ exp βFm λ, δðλ λ0 Þ Dλ DUo =Dλ Z

ð6Þ Following the thermodynamic integration formula,24 the freeenergy change between the initial state with λi, which is the lower bound of the scaling parameter range, and any target state with the scaling parameter λ can be estimated as a function of λ: Z λ dGo ðλÞ ð7Þ Go ðλÞ ¼ dλ0 dλ 0 λi λ

Therefore, eqs 6 and 7 can act as the “recursion slave” to obtain Go(λ), which is the target function for fm(λ). It should be noted that on the basis of Uo(Xt, Xe, λ) (eq 2), in the present work, (∂Uo(Xt,Xe,λ)/∂λ) is equal to _

_

ÆΦc jHt ðXt Þ þ Ht=e ðXt , Xe ÞjΦs æQM=MM bond, angle, improper þ Utedihedral ðXt , Xe ÞMM Utt ðXt ÞMM

ð8Þ

∑t j

39

DUo > tj 7> Uo λ tj Uo ð1Þ þ fm λ tj þ Fm λ tj , = 6 6 7 Dλ 7 exp6 4 5> RT > 2

ð5Þ

builds up and eventually flattens the underlying curvature of the target free-energy surface. Then, like in traditional metadynamics simulations, the free-energy profile along the orthogonal space [λ, (∂Uo(Xt, Xe, λ)/∂λ)] can be estimated as Fm[λ, (∂Uo(Xt, Xe, λ)/∂λ)]. Thus, for any state λ0 , the free-energy profile along

;

þ const

ð9Þ

In summary, two steps are needed in the present QM/MM OSRW-based generalized ensemble simulations. First, a recursion simulation is performed to obtain a set of second-order generalized biasing potentials, Fm[λ, (∂Uo/∂λ)] and fm(λ). 3926



ARTICLE

Figure 1. The global minimum binding mode between DMSO and FKBP12. (a) The binding site structure (from PDB code 1d7h); DMSO is shown with the space-filling model; the side chain and Trp59 and the backbone of Ile56 are shown with the stick model; the rest of the protein is shown with the molecular surface model. (b) Schematic illustration of the DMSO and FKBP12 complex.

Second, a production generalized ensemble simulation is carried out on the resulting energy potential (eq 3). In the analysis step, any ensemble property including the free-energy surfaces along some trial reaction coordinate system can be calculated on the basus if the reweighting scheme described by eq 9.

III. COMPUTATIONAL DETAILS The present QM/MM OSRW sampling method was implemented in our customized CHARMM program.26,27 In the current implementation, the QM/MM potential is based on the SCCDFTB (self-consistent charge density functional tight binding)/CHARMM treatment.28 It is noted that achieving the proposed sampling strategy in the regime of higher-level QM models is certainly one of our future goals. To illustrate the presently proposed algorithm, two studies were performed. III.A. Study 1: DMSO Binding to the FKBP12 Protein. FKBP12 is the receptor of the immunosuppressant drug FK506;29 understanding its binding with various small molecule inhibitors, such as dimethyl sulfoxide (DMSO), can be of significance to the future design of more potent ligands.30 As shown in Figure 1, DMSO has no internal degree of freedom. This study is to demonstrate the capability of the present method in sampling relative translations and rotations between a protein and a ligand. The complex structure was obtained based on the PDB code, 1d7h.30 In the present simulation study, the DMSO molecule was treated quantum mechanically; and the rest of the system, including the FKBP12 protein and the water molecules, was described classically. As illustrated in Figure 1, in the crystal structure, the DMSO molecule forms a hydrogen-bonding interaction with the amide group of Ile 56 and a strong packing interaction with the sidechain of Trp59. In the initial structure of the recursion simulation, the DMSO molecule was rearranged so that it lost the native packing contact with Trp59; in this structure, the distance between the DMSO sulfur atom and the CD2 atom of Trp59, as shown by the dashed line in Figure 1b, is 5.3 Å instead of its original value of 3.6 Å. As discussed in the section of Theoretical Design, one set of generalized ensemble conformational sampling requires two consecutive simulations; that is, a recursion simulation and a production simulation. Here, the recursion

simulation was continued by the production simulation when the native binding mode in the crystal structure was visited three times (around 420 ps of the recursion simulation); on the basis of the obtained set of the biasing potentials, a total of 950 ps of the production equilibrium simulation was performed. In this study, the range of the scaling parameter was set between 0.8 and 1. To test the robustness of the recursion simulation in searching the native binding structure, two more independent recursion simulations were performed with the initial DMSO molecule located farther away from the native binding position. Specifically, the initial distance between the DMSO sulfur atom and the CD2 atom of Trp59 was set as 7.0 Å, and the distance between the oxygen of DMSO and the amide hydrogen of Ile56 was set as 3.5 Å. III.B. Study 2: The Conformational Distribution of a Novel Mercaptosulfonamide Inhibitor in Aqueous Solution and Its Binding Pose in the Zinc-Containing Matrix Metalloproteinase-9 (MMP-9). The development of selective inhibitors for matrix metalloproteinase-9 (MMP-9) has attracted a strong interest due to their possible usage as “chemical tools” in understanding relevant biological processes.3133 Recently, a novel family of mercaptosulfonamide inhibitors, as represented by the one in Figure 2, was designed to target MMP enzymes.34 As shown in Figure 2a, this mercaptosulfonamide molecule has several rotating bonds; one would expect that multiple conformations may exist when it is in aqueous solution. In the present study, this system was employed to exam the capability of our QM/MM OSRW strategy in sampling torsional dynamics of flexible molecules in aqueous solution. In the present simulation study, the mercaptosulfonamide molecule was treated quantum mechanically, and the solvent water molecules were described classically. The recursion simulation was continued by the production simulation when the indentified major conformations were repetitively visited (around 300 ps of the recursion simulation). On the basis of the obtained biasing potentials, a total of 350 ps of the production equilibrium simulation was performed. In this study, the range of the scaling parameter was set between 0.9 and 1. MMP-9 is a zinc-containing enzyme.3133 As shown in Figure 2c, upon the binding to MMP-9, the thiol group of the mercaptosulfonamide molecule loses the proton, and the resulting thiolate forms a strong interaction with the zinc center. The 3927



ARTICLE

Figure 2. The target mercaptosulfonamide inhibitor molecule and its possible binding pose in MMP-9. (a) The structure of the mercaptosulfonamide molecule. (b) The schematic illustration of the mercaptosulfonamide molecule. (c) The schematic illustration of a possible binding pose of the complex of the mercaptosulfonamide molecule and MMP-9.

thiolatezinc interaction has a large covalent component; therefore, the QM/MM potential is ideal to describe the binding of the mercaptosulfonamide molecule in MMP-9. It should be particularly noted that earlier works have demonstrated the feasibility of the SCCDFTB energy function in treating such metalligand interactions.35 Due to the fact that the SCCDFTB is a semiempirical DFT method, careful benchmark studies like the ones in ref 35 are required when other metalligand interactions are to be modeled. When bound to MMP-9, the aromatic rings of the inhibitor molecule are expected to take a conformational pose and fold into the well-defined S10 binding pocket (Figure 2c). Therefore, the analysis of the conformations of the inhibitor molecule will be focused around two dihedral angles (χ1, χ2), the definition of which is illustrated in Figure 2b. In this part of the study, two independent QM/MM OSRW recursion simulations, each of which was started with a distinct docked conformation, were performed to refine the guessed complexes of the target inhibitor and MMP-9. Regarding the QM/MM setup, the mercaptosulfonamide molecule, the zinc center, and the three histidine ligands were treated quantum mechanically. As exemplified by the treatment on His405 in Figure 2c, the cross-bond boundaries between these histidine residues and the rest of the protein were treated via the generalized hybrid orbitals (GHO).3638 On the OSRW setup, the range of the scaling parameter was set between 0.9 and 1. III.C. General Simulation Setups. Each of the protein systems was solvated in a corresponding cubic water (with the TIP3P water model39) box with the distance between the edges of the box and the system initial structure set as 10 Å; the small molecule system was solvated in a 31 Å 31 Å 31 Å box. Counterions (the sodium and chloride ions) were added to neutralize each of the protein systems with the ionic strength set as 0.15 M. The molecular mechanical interactions were based on the CharmmCMAP parameters.40,41 As shown in eq 2, on the target QM region of each study, only bond, angle, and improper energy terms are required; most of these parameters were either in the original CHARMM parameter set40 or obtained from the QuantaCharmm parameter set.42 In the MM description of the zinc center, the bond and angle terms were introduced between zinc and its ligands to keep the structural integrity at the partial QM (scaled by λ) states.

It is noted that these bond terms are scaled off when the λ parameter visits the QM state; therefore, these terms only act as the chaperon energy terms during the sampling, but do not determine the conformations of the target QM state. The longrange electrostatics was treated by the recently integrated SCCDFTB-based QM/MM lattice Ewald method;43 the shortrange interactions were switched starting at 8 Å and totally off at 12 Å. Following the same implementation strategy in our earlier work,18 the present QM/MM OSRW sampling method is coupled with the Ewald treatment and the GHO treatment. The NoseHoover method44 was employed to maintain a constant temperature at 300 K, and the Langevin piston algorithm was used to maintain the constant pressure at 1 atm. The time-step was set as 1 fs. In all the OSRW setups, the height of the Gaussian function h (eq 4) was set as 0.01 kcal/mol; the widths of the Gaussian function, ω1 and ω2, were set as 0.01 and 4 kcal/mol, respectively. In the recursion kernel, Fm(λ,(∂Uo)/∂λ) was updated every 10 time steps; in the recursion slave, fm(λ) was updated every 100 time steps.

IV. RESULTS AND DISCUSSIONS IV.A. Study 1: DMSO Binding to the FKBP12 Protein. The Recursion Simulation. As shown in Figure 3a, it took the recursion

simulation ∼40 ps to travel from the state of λ = 1 to the state of λ = 0.8; and then it took this simulation another 40 ps to complete the first round-trip visit in the scaling parameter space. During this stage of the recursion simulation, DMSO remained at the same binding pose as where it was initially arranged (Figure 3c). As noted in section III.A, in the starting binding mode, the substrate’s native binding interaction with Ile56 was kept (with the distance between the oxygen of DMSO and the amide hydrogen of Ile56 as 2.3 Å), and its native packing with Trp59 was intentionally flipped away (with the distance between the sulfur of DMSO and the CD2 atom of Trp59 as 5.1 Å, instead of the original value of 3.6 Å). With ∼110 ps recursive activation, DMSO lost its initial binding pose and adopted a new binding gesture, in which its native binding interaction with Ile56 was also

3928



ARTICLE

Figure 3. The results of the recursion (docking) simulation of DMSO in FKBP12. (a) The time-dependent scaling parameter changes. (b) The timedependent energy derivative changes. (c) The time-dependent changes of two monitored distances: the one between the DMSO oxygen and the amide hydrogen of Ile56 and the one between the DMSO sulfur and the CD2 atom of Trp59. (d) The obtained set of the biasing potentials.

lost, with the distance between the DMSO oxygen and the Ile56 amide hydrogen increased to ∼6.3 Å (Figure 3c). As shown in Figure 3b, during this structural transition (around 110 ps), the energy derivative (∂Uo(Xt, Xe, λ)/∂λ), the change of which as revealed in eq 8 describes largely the fluctuation of the interaction energy between DMSO and the rest of the environment, increased by about 6 kcal/mol. At about 150 ps, the energy derivative decreased to a lower value range (Figure 3b). Corresponding to this energy derivative decrease, as shown in Figure 3c, DMSO was docked into the native binding mode with both of the monitored distances changed to the same ranges as those in the crystal structure. In the following 50 ps, the interaction between DMSO and the environment became stronger (Figure 3b); however, its binding poses remain in the same global minimum energy well as the native one. The above observations indicate that this 50 ps recursion simulation (150200 ps) is mostly responsible for relaxing the DMSO environment around its global minimum binding mode. It is particularly noted that during the 150200 ps recursion simulation, the environmental responses to the formation of the DMSO native binding mode were enabled through the active changes (lowering) of the energy derivative (∂Uo(Xt, Xe, λ)/∂λ). The result of the first 200 ps recursion simulation confirms the fact that the recursion simulation can serve as a tool for initial structural refinement and optimization. In the following 220 ps

of the recursion simulation, the DMSO molecule left and revisited the crystal structure binding pose repetitively. A total of 420 ps of the recursion simulation led to a set of the biasing potentials as shown in Figure 3d, where dfm(λ)/dλ is shown as a white line and Fm[λ, (∂Uo/∂λ)] is shown as a twodimension contour plot. As revealed in Figure 3d, the center of (∂Uo(Xt, Xe, λ)/∂λ) [Æ(∂Uo(Xt, Xe, λ)/∂λ)æλ ] at the state of λ = 0.8 increases by about 5 kcal/mol in comparison with that of the state of λ = 1. Thereby, when the system randomly walks from the target state (λ = 1) to the other end state (λ = 0.8), ∼5 kcal/mol interaction weakening can be realized to accelerate the relative motion between DMSO and FBKP12. In addition, ∼7 kcal/mol activation associated with the biasing term Fm[λ, (∂Uo/ ∂λ)] (the contour portion of Figure 3) was introduced in the orthogonal space to simultaneously accelerate the environmental responses to the enhanced DMSO motions at each λ state. The Production Simulation. On the basis of the set of biasing potential functions that were obtained in the recursion step, as shown in Figure 4a, the equilibrium production simulation could lead DMSO to (1) locate its global minimum binding pose (starting at 130 ps) and (2) repetitively leave and relocate this native binding mode. It is noted that the initial binding pose of DMSO in this equilibrium simulation was similar to the one arranged for the recursion simulation, the structure of which was collected at around 240 ps of the recursion simulation. Specifically, 0

3929



ARTICLE

Figure 4. The results of the production simulation of DMSO in FKBP12. (a) The time-dependent changes of two monitored distances: the one between the DMSO oxygen and the amide hydrogen of Ile56 and the one between the DMSO sulfur and the CD2 atom of Trp59. (b) The free-energy surface via the projection along two monitored distances.

Figure 5. The results of two longer recursion simulations of DMSO in FKBP12. Here, the initial position of DMSO was set farther from the binding site. The red and blue lines show the time-dependent changes of two monitored distances. The color scheme is the same as the one in Figure 3c.

in ∼900 ps production simulation, DMSO visited the native binding mode three times (at around 140, 300, and 850 ps). As discussed earlier, the native binding mode has two native interactions, which are represented by two characteristic distances. On the basis of eq 9, the samples collected from the production simulation were employed to generate the ligand-binding, free-energy landscape (Figure 4b), which is projected along two monitored distances. As discussed earlier, these two distances represent two essential native interactions between DMSO and FBKP12: the hydrogen bonding interaction with the amide group of Ile56 and the packing interaction with the side chain of Trp59, specifically through the sulfur atom of DMSO. As shown in Figure 4b, the binding structure of the global freeenergy minimum (with the distance between the oxygen of DMSO and the amide hydrogen of Ile56 as 2.3 Å and the distance between the sulfur of DMSO and the CD2 atom of Trp59 as 3.6 Å) is the same as the one in the crystal structure. In addition, there are three local free-energy minima: the first one preserves the interaction with the amide hydrogen of Ile56 and loses the packing interaction with Trp59 (with two distances of 2.3 and 4.7 Å respectively); the second one loses the interaction with the amide hydrogen of Ile56 and preserves the packing interaction with Trp59 (with two distances of 5.0 and 5.0 Å, respectively); the

third one further leaves the binding site, with the distance between the oxygen of DMSO and the amide hydrogen of Ile56 increased to around 7.0 Å, although still preserving the packing interaction with Trp59 because in this binding pose, Trp59 leaves the binding site simultaneously with the escaping of DMSO. On the edges of these free-energy wells (the green regions), the height is ∼2.53.0 kcal/ mol, which is consistent with the measured binding affinity of DMSO to FKBP12 (around 2.5 kcal/mol). It should be noted that although the native binding mode has been repetitively visited, like any equilibrium simulation, whether the simulation length has fulfilled the ergodicity requirement needs continuous enlargement of simulation time to exam. To test the robustness of the recursion simulation in searching the native binding structure, the initial DMSO molecule was located farther away from the native binding position. As shown in Figure 5, two independent simulations can robustly locate and relocate the native binding structure, even after the ligand molecule moves very far from the binding site. IV.B. Study 2: The Conformation Distribution of a Novel Mercaptosulfonamide Inhibitor in Aqueous Solution and Its Binding Poses in Zinc-Containing Matrix Metalloproteinase9 (MMP-9). The Recursion Simulation of the Mercaptosulfonamide Inhibitor in Aqueous Solution. As shown in Figure 6a, it took 3930



ARTICLE

Figure 6. The results of the recursion simulation of the mercaptosulfonamide molecule in aqueous solution. (a) The time-dependent scaling parameter changes. (b) The time-dependent energy derivative changes. (c) The time-dependent changes of two monitored dihedral angles: χ1 (red) and χ2 (blue); the definition of χ1 and χ2 is in Figure 2b. (d) The obtained set of biasing potentials.

the recursion simulation ∼25 ps to travel from the state of λ = 1 to the other end state (λ = 0.9). At 25 ps, the energy derivative (∂Uo(Xt, Xe, λ)/∂λ) went through a local peak and thereafter relaxed to a lower value region (Figure 6b). During this recursion period, the molecule left the initial free-energy well (χ1 = 40; χ2 = 30) and underwent a conformational transition to another free-energy well (χ1 = 150; χ2 = 40) (Figure 6c). Then, the inhibitor molecule dwelled in the same free-energy well of (χ1 = 150; χ2 = 40) for 70 ps, in the middle of which (at around 70 ps) it briefly moved into and back from the region of (χ1 = 150; χ2 = 150). This dwelling behavior indicates the fact that the conformational region around (χ1 = 150; χ2 = 40) is a deep, free-energy minimum. At around 100 ps, the molecule moved to the region of (χ1 = 150; χ2 = 150). After another 20 ps, the molecular conformation quickly passed by the region of (χ1 = 150; χ2 = 40) and jumped back to the very beginning conformational well (χ1 = 40; χ2 = 30). In the following recursion simulation time, the system repetitively visited the above three major conformational regions. As shown in Figure 6b, during 300 ps of the recursion simulation, the energy derivative fluctuated, with its range being continuously increased (Figure 6b). Eventually, the recursion simulation led to a set of the biasing potentials (Figure 6d), where dfm(λ)/dλ is shown as a white line and Fm[λ, (∂Uo/∂λ)] is shown as a two-dimension contour plot. As revealed in Figure 6d, the center of (∂Uo(Xt, Xe, λ)/∂λ) [Æ(∂Uo(Xt, Xe, λ))/∂λæλ ] at the state of λ = 0.9 increases by ∼28 kcal/mol in comparison with 0

that of the state of λ = 1. Thereby, when the system randomly walks from the target state (λ = 1) to the other end state (λ = 0.9), ∼28 kcal/mol essential energy (∂Uo(Xt, Xe, λ)/∂λ) weakening can be realized to accelerate the conformational changes of the inhibitor molecule. In addition, ∼11 kcal/mol activation associated with the biasing term Fm[λ, (∂Uo/∂λ)] (the contour portion of Figure 6d) was introduced in the orthogonal space to simultaneously accelerate the environmental responses to various dihedral angle rotations of the inhibitor molecule at each λ state. The Production Simulation of the Mercaptosulfonamide Inhibitor in Aqueous Solution. On the basis of the set of the biasing potential functions that were obtained in the recursion step, as shown in Figure 7a, the equilibrium production simulation could lead the inhibitor molecule to (1) locate its global conformational minima (χ1 = 150; χ2 = 40) and (χ1 = 150; χ2 = 150), (2) repetitively leave and relocate three major minima that were identified in the recursion simulation step, and (3) additionally visit two higher free-energy minima (χ1 = 15; χ2 = 50) and (χ1 = 150; χ2 = 150). It is noted that the initial conformation in this equilibrium simulation was the same as the one arranged for the recursion simulation. On the basis of eq 9, the samples collected from the production simulation were employed to generate the conformational free-energy landscape (Figure 7b), which is through the projection along two monitored dihedral angles. As discussed earlier, these two dihedral angles determine the binding pose of the 3931



ARTICLE

Figure 7. The results of the production simulation of DMSO in FKBP12. (a) The time-dependent changes of two monitored dihedral angles: χ1 (red) and χ2 (blue); the definition of χ1 and χ2 is in Figure 2b. (b) The free-energy surface via the projection along two monitored dihedrals; the free-energy surface is shown in the periodic fashion.

Figure 8. The results of the two recursion (docking) simulations of the mercaptosulfonamide molecule in MMP-9. (a) The time-dependent changes of two monitored dihedral angles, χ1 (red) and χ2 (blue), in the first docking simulation. (b) The time-dependent changes of two monitored dihedral angles, χ1 (red) and χ2 (blue), in the second docking simulation.

target inhibitor inside the S0 binding pocket of MMP-9; therefore, the related conformation distribution is our target of interest. As shown in Figure 7b, there are two equally populated global free-energy minima, (χ1 = 150; χ2 = 40) and (χ1 = 150; χ2 = 150), each of which uniquely has its surface area of the aromatic ring portion minimized. As shown in Figure 7b, the initial input conformation is also located in a local free-energy region, (χ1 = 40; χ2 = 30), the free energy of which is ∼3 kcal/mol higher than those of two global free-energy minima. One of the two newly identified free-energy wells, (χ1 = 15; χ2 = 50), has its minimum higher than the global minima by ∼7 kcal/ mol, and the other newly identified free-energy well, (χ1 = 15; χ2 = 50), has its free-energy minimum higher than the global minima by ∼9 kcal/mol. It should be noted that the free-energy barriers that separate various free-energy minima are relatively large: they range from 15 to 22 kcal/mol. It is consistent with the magnitude of the energy boosting employed in the biasing energy potentials. The related details were discussed in the Recursion Simulation subsection. The Docking Simulation of the Mercaptosulfonamide Inhibitor in Zinc-Containing Matrix Metalloproteinase-9 (MMP-9). As shown in Figure 2, when a mercaptosulfonamide inhibitor molecule binds to the MMP-9 active site, one would expect that

its mercapto moiety and its sulfonamide functional group form strong specific interactions in the MMP-9 catalytic site; however, its aromatic rings are flexible, and how they are stacked in the S10 pocket is a key issue to address in the effort of predicting its binding pose to MMP-9. In this docking study, we performed two recursion-only simulations, which, as discussed earlier, can act as corresponding docking procedures to enable structural refinement and optimization. Two of the starting conformations of the mercaptosulfonamide inhibitor in the two docking (recursion) simulations are one of the global conformational minima with the aromatic rings’ conformation of (χ1 = 150; χ2 = 40) and one of the local minima with the aromatic ring conformation of (χ1 = 40; χ2 = 30), the latter of which was the starting conformation in the aqueous solution simulations. As shown in Figure 8a, during the first docking simulation, the aromatic rings stay around the initial conformation, even after 350 ps OSRW-recursion sampling acceleration, and during the second docking simulation (Figure 8b), two aromatic rings moved to the conformational region of (χ1 = 150; χ2 = 150) in around 50 ps and, after another 20 ps, moved to and stayed in the conformational region of (χ1 = 150; χ2 = 40), which is the same as the one located in the first docking simulation. It should be especially noted that the S10 pocket is a very condensed 3932



ARTICLE

Figure 9. The docked binding pose of the mercaptosulfonamide molecule in MMP-9. (a) The structural changes of the binding site before and after the docking simulations. (b) The front view of the binding gesture. (c) The top view of the binding gesture. All the atoms of the inhibitor, the zinc center, the side chains of the zinc-bound histines and Glu402, and the backbones of Ala188 and Leu189 are shown via the stick model. The protein is shown in the cartoon model. As shown, the S10 pocket is formed between a loop and helix structural moieties.

environment; the efficient realization of the conformational transitions as enabled in the second docking (OSRW recursion) simulation is otherwise very challenging. In addition to the conformational arrangement of the aromatic rings in the S10 pocket, the other portion of the inhibitor molecule and the MMP-9 protein binding site also went through the corresponding structural adaptations in 200 ps of the QM/MM OSRW recursion simulations. As shown in Figure 9a, the MMP-9 binding site changed from the initial unbound crystal structure to the tighter one. For instance, the aromatic rings bind deeper into the S10 pocket, and the hydrophobic five-member ring of the inhibitor forms a tighter and complementary contact with the hydrophobic portion of MMP-9. As shown in Figure 8b and c, in the catalytic portion of the MMP-9 binding site, a relatively rigid binding mode was formed: the thiolate group of the ligand forms a strong covalent-like interaction with the zinc(II) center that is also bound by three other histine residues; the amide group of the ligand forms a hydrogen bond with the catalytic residue Glu402; and one of the oxygen atoms in the sulfonamide moiety forms two hydrogen bonds with the amide groups of Ala188 and Leu189, and the other oxygen atom is exposed to the bulk. Obviously, such good characterization of the binding pose benefits from the usage of the QM/MM energy model.

V. CONCLUDING REMARKS Accurate prediction of molecular conformations in explicit environments, such as aqueous solution and protein interiors, can facilitate our understanding of various molecular recognition

processes. Traditional computational approaches are usually limited as a result of their compromised choices between the underlying energy model and the sampling length. Taking advantage of a recent second-order generalized ensemble scheme (e.g., the OSRW strategy), which can synergistically accelerate the motion of a focused region and its coupled environmental dynamics, we present a QM/MM (combined quantum mechanical/molecular mechanical)-based molecular dynamics sampling technique to explore molecular conformational landscapes in explicit environments. As discussed in our early work, the OSRW sampling is a general approach that can be applied to any design of Uo(λ). For instance, when Uo(λ) = λUo, the scaling parameter plays the role as an effective temperature.15,16 In the present work, Uo(Xt, Xe, λ) in eqs 2 and 3 can be the same as the one we employed in the present work or can be redesigned, as long as it follows the constraint of Uo(Xt, Xe, 1) = Uo(Xt, Xe). When an OSRW simulation is set up, the range of the potential scaling parameter λ needs to be preset. On the basis of the nature of the target event, if environmental collective motions pose the major sampling challenge, the lower bound of the potential scaling parameter should set higher so that the orthogonal space sampling can be dominant. Certainly, how to optimally set the potential scaling parameter λ range needs to be further explored. To illustrate the present QM/MM potential-scaling-based OSRW sampling scheme, the model studies on the binding of DMSO to the FKBP12 protein, the conformation distribution of a novel mercaptosulfonamide inhibitor in aqueous solution and its binding poses in zinc-containing matrix metalloproteinase-9 3933


The Journal of Physical Chemistry B (MMP-9) were performed. As demonstrated, the OSRW recursion simulation can serve as a tool for initial structural refinements and optimizations, such as docking and locating various possible solution conformations. The OSRW production simulation allows quantitatively mapping of various related free-energy surfaces. The present QM/MM second-order generalized ensemble sampling technique enables feasible usage of the QM/ MM model to sample molecular conformations in condensed environments.

’ AUTHOR INFORMATION Corresponding Author

*E-mail: [email protected].

’ ACKNOWLEDGMENT Funding support (MCB 0919983 to W.Y.) from the National Science Foundation is acknowledged. We thank the Florida State University High Performance Computing (HPC) Center and the Institute of Molecular Biophysics Computing facility (Dr. Michael Zawrotny) for computing support. ’ REFERENCES (1) Warshel, A.; Levitt, M. Theoretical studies of enzymic reactions Dielectric, electrostatic and steric stabilization of carbounium-ion in reaction of lysozyme. J. Mol. Biol. 1976, 103, 227–249. (2) Field, M. J.; Bash, P. A.; Karplus, M. A combined quantum mechanical and molecular mechanical potential for molecular-dynamics simulations. J. Comput. Chem. 1990, 11, 700–733. (3) Gao, J. L. Hybrid quantum and molecular mechanical simulations: An alternative avenue to solvent effects in organic chemistry. Acc. Chem. Res. 1996, 29, 298–305. (4) Bakowies, D.; Thiel, W. Hybrid models for combined quantum mechanical and molecular mechanical approaches. J. Phys. Chem. 1996, 100, 10580–10594. (5) Monard, G.; Merz, K. M. Combined quantum mechanical/ molecular mechanical methodologies applied to biomolecular systems. Acc. Chem. Res. 1999, 32, 904–911. (6) Hensen, C.; Hermann, J. C.; Nam, K. H.; Ma, S. H.; Gao, J. L.; Holtje, H. D. A combined QM/MM approach to proteinligand interactions: Polarization effects of the HIV-1 protease on selected high affinity inhibitors. J. Med. Chem. 2004, 47, 6673–6680. (7) Cho, A. E.; Guallar, V.; Berne, B. J.; Friesner, R. Importance of accurate charges in molecular docking: Quantum mechanical/molecular mechanical (QM/MM) approach. J. Comput. Chem. 2005, 26, 915–931. (8) Grater, F.; Schwarzl, S. M.; Dejaegere, A.; Smith, J. C. Proteinligand binding free energies calculated with quantum mechanics/ molecular mechanics. J. Phys. Chem. B 2005, 109, 10474–10483. (9) Khandelwal, A.; Lukacova, V.; Comez, D.; Kroll, D. M.; Raha, S.; Balaz, S. A combination of docking, QM/MM methods, and MD simulation for binding affinity estimation of metalloprotein ligands. J. Med. Chem. 2005, 48, 5437–5447. (10) Jono., R.; Watanabe, Y.; Shimizu, K.; Terada, T. Multicanonical ab initio QM/MM molecular dynamics simulation of a peptide in an aqueous environment. J. Comput. Chem. 2010, 31, 1168–1175. (11) Li, H. Z.; Yang, W. Sampling enhancement for the quantum mechanical potential based molecular dynamics simulations: A general algorithm and its extension for free energy calculation on rugged energy surface. J. Chem. Phys. 2007, 126, 114104. (12) Seabra, G. D.; Walker, R. C.; Roitberg, A. E. Are current semiempirical methods better than force fields? A study from the thermodynamics perspective. J. Phys. Chem. A 2009, 113, 11938–11948.

ARTICLE

(13) Li, H.; Min, D.; Liu, Y.; Yang, W. Essential energy space random walk via energy space metadynamics method to accelerate molecular dynamics simulations. J. Chem. Phys. 2007, 127, 094101. (14) Min, D. H.; Yang, W. A divide-and-conquer strategy to improve diffusion sampling in generalized ensemble simulations. J. Chem. Phys. 2008, 128, 094106. (15) Li, H. Z.; Fajer, M.; Yang, W. Simulated scaling method for localized enhanced sampling and simultaneous “alchemical” free energy simulations: A general method for molecular mechanical, quantum mechanical, and quantum mechanical/molecular mechanical simulations. J. Chem. Phys. 2007, 126, 024106. (16) de Koning, M.; Antonelli, A.; Yip, S. Optimized free-energy evaluation using a single reversible-scaling simulation. Phys. Rev. Lett. 1999, 83, 3973–3977. (17) Yang, W.; Bitetti-Putzer, R.; Karplus, M. Chaperoned alchemical free energy simulations: A general method for QM, MM, and QM/MM potentials. J. Chem. Phys. 2004, 120, 9450–9453. (18) Min, D. H.; Zheng, L. Q.; Harris, W.; Chen, M. G.; Lv, C.; Yang, W. Practically efficient QM/MM alchemical free energy simulations: The orthogonal space random walk strategy. J. Chem. Theor. Comput. 2010, 6, 2253–2266. (19) Duane, S.; Kennedy, A. D.; Pendleton, B. J.; Roweth, D. Hybrid Monte-Carlo. Phys. Lett. B 1987, 195, 216–222. (20) Kong, X. J.; Brooks, C. L. λ-Dynamics: A new approach to free energy calculations. J. Chem. Phys. 1996, 105, 2414–2423. (21) Zheng, L. Q.; Chen, M. G.; Yang, W. Random walk in orthogonal space to achieve efficient free-energy simulation of complex systems. Proc. Natl. Acad. Sci. U.S.A. 2008, 105, 20227–20232. (22) Zheng, L. Q.; Chen, M. G.; Yang, W. Simultaneous escaping of explicit and hidden free energy barriers: Application of the orthogonal space random walk strategy in generalized ensemble based conformational sampling. J. Chem. Phys. 2009, 130, 234105. (23) Laio, A.; Parrinello, M. Escaping free-energy minima. Proc. Natl. Acad. Sci. U.S.A. 2002, 99, 12562–12566. (24) Kirkwood, J. G. Statistical mechanics of fluid mixtures. J. Chem. Phys. 1935, 3, 300–313. (25) Torrie, G. M.; Valleau, J. P. Non-physical sampling distributions in Monte-Carlo free-energy estimationUmbrella Sampling. J. Comput. Phys. 1977, 23, 187–199. (26) Brooks, B. R.; Bruccoleri, R. E.; Olafson, B. D.; States, D. J.; Swaminathan, S.; Karplus, M. Charmm A program for macromolecular energy, minimization, and dynamics calculations. J. Comput. Chem. 1982, 4, 187–217. (27) Brooks, B. R.; Brooks, C. L.; Mackerell, A. D.; Nilsson, L.; Petrella, R. J.; Roux, B.; Won, Y.; Archontis, G.; Bartels, C.; Boresch, S.; Calfisch, A.; Caves, L.; Cui, Q.; Dinner, A. R.; Feig, M.; Fischer, S.; Gao, J.; Hodoscek, M.; Im, W.; Kuczera, K.; Lazaridis, T.; Ma, J.; Ovchinnikov, V.; Paci, E.; Pastor, R. W.; Post, C. B.; Pu, J. Z.; Schaefer, M.; Tidor, B.; Venable, R. M.; Woodcock, H. L.; Wu, X.; Yang, W.; York, D. M.; Karplus, M. CHARMM: The biomolecular simulation program. J. Comput. Chem. 2009, 30, 1545–1614. (28) Cui, Q.; Elstner, M.; Kaxiras, E.; Frauenheim, T.; Karplus, M. A QM/MM implementation of the self-consistent charge density functional tight binding (SCC-DFTB) method. J. Phys. Chem. B 2001, 105, 569–585. (29) Brown, E. J.; Albers, M. W.; Shin, T. B.; Ichikawa, K.; Keith, C. T.; Lane, W. S.; Schreiber, S. L. A mammalian protein targeted by G1arresting rapamycin-receptor complex. Nature 1994, 369, 756–758. (30) Burkhard, P.; Taylor, P.; Walkinshaw, M. D. X-ray structures of small ligandFKBP complexes provide an estimate for hydrophobic interaction energies. J. Mol. Biol. 2002, 295, 953–962. (31) Sang, Q. X.; Jin, Y.; Newcomer, R. G.; Monroe, S. C.; Fang, X.; Hurst, D. R.; Lee, S.; Cao, Q.; Schwartz, M. A. Matrix metalloproteinase inhibitors as prospective agents for the prevention and treatment of cardiovascular and neoplastic diseases. Curr. Top. Med. Chem. 2006, 6, 289–316. (32) Muroski, M. E.; Roycik, M. D.; Newcomer, R. G.; Van den Steen, P. E.; Opdenakker, G.; Monroe, H. R.; Sahab, Z. J.; Sang, Q. X. 3934



ARTICLE

Matrix metalloproteinase-9/gelatinase B is a putative therapeutic target of chronic obstructive pulmonary disease and multiple sclerosis. Curr. Pharm. Biotechnol. 2008, 9, 34–46. (33) Hu, J.; Van den Steen, P. E.; Sang, Q. X.; Opdenakker, G. Matrix metalloproteinase inhibitors as therapy for inflammatory and vascular diseases. Nat. Rev. Drug Discovery 2007, 6, 480–98. (34) Schwartz, M. A.; Jin, Y.; Sang, Q. X. Substituted heterocyclic mercaptosulfonamide metalloprotease inhibitors. PCT Int. Patent US2009/055742, 2010, Pub. No.: WO/2010/028051. (35) Elstner, M.; Cui, Q.; Munih, P.; Kaxiras, E.; Frauenheim, T.; Karplus, M. Modeling zinc in biomolecules with the self consistent charge-density functional tight binding (SCC-DFTB) method: Applications to structural and energetic analysis. J. Comput. Chem. 2003, 24, 565–581. (36) Gao, J. L.; Amara, P.; Alhambra, C.; Field, M. J. A generalized hybrid orbital (GHO) method for the treatment of boundary atoms in combined QM/MM calculations. J. Phys. Chem. A 1998, 102, 4714–4721. (37) Amara, P.; Field, M. J.; Alhambra, C.; Gao, J. L. The generalized hybrid orbital method for combined quantum mechanical/molecular mechanical calculations: Formulation and tests of the analytical derivatives. Theor. Chem. Acc. 2000, 104, 336–343. (38) Pu, J. Z.; Gao, J. L.; Truhlar, D. G. Combinding self-consistentcharge density-functional tight-binding (SCC-DFTB) with molecular mechanics by the generalized hybrid orbital (GHO) method. J. Phys. Chem. A 2004, 108, 5454–5463. (39) Jorgensen, W. L.; Chandrasekhar, J.; Madura, J. D.; Impey, R. W.; Klein, M. L. Comparison of simple potential functions for simulating liquid water. J. Chem. Phys. 1983, 79, 926–935. (40) MacKerell, A. D.; Bashford, D.; Bellott, M.; Dunbrack, R. L.; Evanseck, J. D.; Field, M. J.; Fischer, S.; Gao, J.; Guo, H.; Ha, S.; JosephMcCarthy, D.; Kuchnir, L.; Kuczera, K.; Lau, F. T. K.; Mattos, C. Michnick, S.; Ngo, T.; Nguyen, D. T.; Prodhom, B.; Reiher, W. E.; Roux, B.; Schlenkrich, M.; Smith, J. C.; Stote, R.; Straub, J.; Watanabe, M.; Wiorkiewicz-Kuczera, J.; Yin, D.; Karplus, M. All-atom empirical potential for molecular modeling and dynamics studies of proteins. J. Phys. Chem. B 1998, 102, 3586–3616. (41) Mackerell, A. D.; Feig, M.; Brooks, C. L. Extending the treatment of backbone energetics in protein force fields: Limitations of gas-phase quantum mechanics in reproducing protein conformational distributions in molecular dynamics simulations. J. Comput. Chem. 2004, 25, 1400–1415. (42) Momany, F. A.; Rone, R. Validation of the general-purpose Quanta/CHARMM force-field. J. Comput. Chem. 1992, 13, 888–900. (43) Nam, K.; Gao, J. L.; York, D. M. An efficient linear-scaling Ewald method for long-range electrostatic interactions in combined QM/MM calculations. J. Chem. Theory Comput. 2005, 1, 2–13. (44) Nose, S. A unified formulation of the constant temperature molecular-dynamics methods. J. Chem. Phys. 1984, 81, 511–519.

3935


MM Molecular Dynamics Sampling in Explicit

Recommend Documents