Thermodynamic Analysis of Thrombin Inhibition by Benzamidine and

To accomplish that, we used the finite difference thermodynamic integration (FDTI) algorithm within the Discover program of MSI. We have shown the imp...
0 downloads 0 Views 183KB Size
Download PDF
466

J. Phys. Chem. B 2002, 106, 466-476

Thermodynamic Analysis of Thrombin Inhibition by Benzamidine and p-Methylbenzamidine via Free-Energy Perturbations: Inspection of Intraperturbed-Group Contributions Using the Finite Difference Thermodynamic Integration (FDTI) Algorithm Cristiano Ruch Werneck Guimara˜ es† and Ricardo Bicca de Alencastro* Physical Organic Chemistry Group, Departamento de Quı´mica Orgaˆ nica, Instituto de Quı´mica, UniVersidade Federal do Rio de Janeiro, Cidade UniVersita´ ria, CT, Bloco A, Sala 609, Rio de Janeiro, RJ 21949-900, Brazil ReceiVed: August 2, 2001

Thrombin is a serine protease responsible for blood coagulation. Since thrombin inhibitors appear to be effective in the treatment and prevention of thrombotic and embolic disorders, considerable attention has been focused on the structure and interactions of this enzyme. Here, to evaluate the relative free energies of hydration and binding to thombin between two serine protease inhibitors, p-methylbenzamidine and benzamidine, we employed molecular dynamics simulations in conjunction with free-energy perturbation calculations. To accomplish that, we used the finite difference thermodynamic integration (FDTI) algorithm within the Discover program of MSI. We have shown the importance of the inclusion of intraperturbed-group contributions to the free energy, and demonstrated that the orthogonality problem that occurs in the calculation of these contributions is adequately treated by the FDTI method. We have also demonstrated that problems of singularity and convergence in free-energy calculations can be properly solved by combining the FDTI method with the Gaussian-Legendre quadrature method for numerical integration, associated with the introduction of a physical criterion to determine the breaking point of a bond or angle described by a harmonic potential. The results seem to indicate that the greater affinity of p-methylbenzamidine for thrombin is derived from stronger electrostatic and hydrophobic interactions between this molecule and the enzyme.

1. Introduction Thrombin is a serine protease that plays a central role in blood coagulation through the conversion of fibrinogen in fibrin and platelets activation. The exaggerated clot formation leads to several cardiovascular disorders such as: venous and arterial thrombosis, atrial fibrillation, stroke, and myocardial infarction.1 As it is the most potent stimulator of platelet aggregation, thrombin has become the principal target for the development of new antithrombotic agents.2 Thrombin, related to the serine proteases trypsin and chymotrypsin, consists of an A chain of 36 amino acids and a B chain of 259 amino acids connected by a disulfide bridge.1 Besides the catalytic triad (Asp102-His57-Ser195), the active site has three other important binding pockets that increase the selectivity of the enzyme. The S1 pocket, which contains an aspartate (Asp189), a small hydrophobic site defined by the Tyr60A-Pro60B-Pro60C-Trp60D insertion loop called Ppocket (proximal to the catalytic site), and a larger hydrophobic site called D-pocket (distal from the catalytic site), separated from the P-pocket by the side chain of Leu99.3 Thrombin inhibitors can be divided in two different types: the electrophilic and the noncovalent. The first are characterized by toxicity problems, lack of selectivity and questionable efficaciousness. The second, on the other hand, have the advantage of being more selective as they interact noncovalently with the enzyme.1 A simple reversible inhibitor, benzamidine (Ki ) 300 nM), is one of them.4 When protonated, benzamidine interacts * To whom correspondence should be addressed. E-mail: [email protected]. FAX number: 55-21-25627256. † E-mail: [email protected].

with the active site of thrombin through a salt bridge with the residue Asp189, within the S1 pocket. The free-energy perturbation technique (FEP),5 which is generally attributed to Zwanzig,6 has successfully been applied to an extensive series of works related to the calculation of relative free-energy differences.7 Due to the reliable binding constants available,8 this technique was tested through the calculation of the relative free energies of binding to trypsin by a series of benzamidine derivatives.9 However, the lack of experimental binding constants for the complexes between these benzamidine derivatives and thrombin, except for benzamidine itself, has prevented the test of this methodology in the study of thrombin inhibition by these compounds. As the experimental results of thrombin inhibition by benzamidine (Bz) and pmethylbenzamidine (PMBz) are now available,10 in this work, we have employed molecular dynamics simulations in conjunction with FEP calculations in order to compare the theoretical and experimental relative free energies of binding for Bz and PMBz. Moreover, to verify the influence of hydration in the binding process, we have calculated the relative free energy of hydration for these compounds. 2. Theoretical Aspects Considering two states, 1 and 2, the Helmholtz free-energy variation between them is given by eq 1. Because of the separation of QNVT, the Helmholtz free-energy variation can be expressed as a sum of ideal gas and configurational parts.

Q2NVT

Z2NVT

∆A ) -kT ln 1 ) ∆Aid - kT ln 1 Q NVT Z NVT

10.1021/jp013563l CCC: $22.00 © 2002 American Chemical Society Published on Web 12/13/2001

(1)

Thrombin Inhibition by Benzamidine

J. Phys. Chem. B, Vol. 106, No. 2, 2002 467

Figure 1. Thermodynamic cycles used for the calculation of the relative free energy of hydration (left) and binding to thrombin (right) for p-methylbenzamidine (PMBz) and benzamidine (Bz).

When calculating relative free-energy differences, a thermodynamic cycle is generally applied.11 Here, we have used a thermodynamic cycle, as exhibited in Figure 1. Since the free energy is a thermodynamic state function, the left and right sides of the cycle give the relative free energy of hydration (∆∆Ahyd ) ∆AhydBz - ∆AhydPMBz ) ∆Atr(sol) - ∆Atr(g)) and the relative free energy of binding (∆∆Abind ) ∆AbindBz - ∆AbindPMBz ) ∆Atr(comp) - ∆Atr(sol)), respectively. Because the quantities (∆Ahyd) and (∆Abind) are impracticable to compute, ∆∆Ahyd and ∆∆Abind are obtained through the simulation of a nonphysical process in three different environments. In our case, the mutation of PMBz into Bz in the gas phase, in solution, and inside the solvated complex, gives (∆Atr(g)), (∆Atr(sol)), and (∆Atr(comp)), respectively. This mutation involves the perturbation of a methyl group into a hydrogen atom. Therefore, the differences in ∆A due to kinetic energy, as in eq 1, may be assumed to be identical in calculating, for example, ∆Atr(sol) and ∆Atr(g), or ∆Atr(comp) and ∆Atr(sol) (Figure 1). Given that, it is reasonable to suppose that H ∼ U.5 Thus, the perturbation method equation (PM)5 may be written as n-1

∆A ) -kT

ln 〈e-(U(λ ∑ i)0

i+1)

- U(λi))/kT

〉λi

(2)

In eq 2, U(λi) and U(λi+1) are the potential energy functions for the states λi and λi+1, respectively, and 〈 〉λi refers to an average over the ensemble of configurations generated using U(λi). Equation 2 shows that, to propitiate a good overlap between the successive phase spaces, a parameter λ is introduced to couple the initial and final states. Another equivalent expression for ∆A can be employed in free-energy difference calculations. It is referred to as the thermodynamic integration method (TI),5 where ∆A is obtained by eq 3:

∆A )

∫01dλ〈

〉

∂U(λ) ∂λ

λ

(3)

Boresch and Karplus have dissected the bonded contributions to free-energy simulations in single and dual topologies simulations. They have shown that if the rigid rotor harmonic oscillator approximation is valid, then the bond or the bond angle freeenergy components of single topology simulations can be separated in vibrational, potential-of-mean force (pmf) type and Jacobian factor contributions:12-14 the first arises from changes in the force constant, when going from the initial to the final state; the second results from changes in nonbonded interactions (intramolecular nonbonded or solute-environment nonbonded interactions) if the equilibrium geometry of a molecule is altered between the states; the last comes from changes in the equilibrium bond lengths and bond angles in the absence of nonbonded interactions. In many FEP calculations, it is assumed that the differences in ∆A due to intraperturbed-group interactions or self-terms,12 such as bond stretching and angle bending, may be identical

when one calculates, for example, ∆Atr(sol) and ∆Atr(g), or ∆Atr(comp) and ∆Atr(sol) (Figure 1). In other words, only contributions from interactions between the perturbed group and the rest of the system are included.12-15 In the case of bond stretching and angle bending terms, vibrational and Jacobian factor contributions can be canceled in the two nonphysical legs of the thermodynamic cycle only if there is no significant coupling between them and pmf-type contributions.12 If such a coupling occurs, due to interactions with the environment (solvent and/or protein), vibrational and Jacobian factor contributions different from the gas phase would be obtained. This would give rise to vibrational and Jacobian factor contributions of the intraperturbed group to the free-energy difference. However, even with no coupling between vibrational and Jacobian factor contributions and pmf-type contributions, in single topology simulations, the last one is projected on to the bonded free-energy component.12,14 Therefore, to obtain accurate results, this type of contribution should be included in FEP calculations. Due to severe convergence problems, as shown by Pearlman and Kollman,15 determining the contribution to the free energy from the harmonic stretching terms can be a problem when one uses the PM formula. Why? Because the harmonic stretching potential energy surface is too steep, and two states differing in equilibrium bond length are too orthogonal. Fortunately, the same is not true for the softer nonbonded, the torsional, and the angle-bending terms. Thus, a large number of windows should be used in eq 2 to propitiate a good overlap between the phase spaces of the states. Still, this procedure would make FEP calculations impracticable. To overcome the orthogonality problem of the harmonic stretching terms in the PM formalism, Pearlman and Kollman proposed a method called bond potential of mean force correction15,16 used in conjunction with the SHAKE algorithm17 to calculate pmf-type contributions to free energy. In the pmf correction method of Pearlman and Kollman, vibrational and Jacobian factor contributions of the bond stretching terms are not calculated, although the latter could be obtained analytically. This is because the Hamiltonian of the system is not dependent on the bond stretching terms (the SHAKE algorithm is applied). In another work, Severance et al.18 have reported a modified derivation of the PM equation, which corrects for the change in equilibrium bond length as a function of the coupling parameter λ. When used to perturb an isolated harmonic oscillator from an initial to a final state, this method converges rapidly to the analytical free-energy difference. The finite difference thermodynamic integration (FDTI) algorithm19 combines the formalisms of PM and TI. FDTI computes numerically the derivatives of the free energy in relation to the coupling parameter at fixed values of λ, followed by numerical integration using a quadrature scheme (eq 4). In other words, differently from PM, the number of λ to be used in the free-energy calculation is not dependent on the degree of overlap between the initial and final phase spaces.

ln〈e-[U(λi + δλ)-U(λi)]/kT〉i ∆λi ∆A ) -kT δλ i)1 n

∑

(4)

In eq 4, n is the number of quadrature points, ∆λi is a parameter that depends on the numerical integration scheme, and δλ is the increment used to compute the numerical derivatives. Since δλ should tend to zero to compute adequately the free-energy derivatives at each λi, FDTI can propitiate a good overlap between orthogonal phase spaces. Therefore, FDTI

468 J. Phys. Chem. B, Vol. 106, No. 2, 2002

Figure 2. The structure of p-amidinophenylpyruvate (APPA).

may afford a good opportunity to calculate correctly the bonded contributions to the free energy. To test this methodology, we first mutated the C-F equilibrium bond length of fluoromethane behaving as an ideal gas, using the FDTI algorithm in the single topology approach. The method was validated by comparing the calculated results with theoretical results derived from classical statistical thermodynamics. After its validation, we applied the FDTI algorithm to the calculation of the relative free energies of hydration and binding to thrombin for Bz and PMBz, including all intraperturbed-group interactions to the free energy. 3. Computational Procedures All calculations were performed in a molecular modeling software from MSI (USA), running on a Silicon Graphics O2 R10000 workstation. The Insight II program (version 97.0) was employed as a graphical interface for the construction and visualization of molecular structures. The all-atom CVFF force field,20 within the DiscoVer program (version 2.9.7), was employed in the energy minimization and the molecular dynamics simulations. The bond stretching was expressed by a simple harmonic function. The cross terms of the force field were not included in the energy expression. As there is no reported X-ray crystal structure of thrombin complexed with PMBz, we used the X-ray crystal structure of the p-amidinophenylpyruvate(APPA).thrombin.hirugen ternary complex with crystal water molecules (structure 1AHT),21 taken from the Brookhaven Protein Data Bank, to build the PMBz‚ thrombin complex model. APPA (Figure 2) is an eletrophilic inhibitor forming a covalent bond with Ser195. The hirugen molecule is bonded to the fibrinogen recognition exosite. In our model, basic residues, such as Arg and Lys are protonated, and acid residues such as Asp and Glu are deprotonated. Due to its normal pKa, the His residues were assumed to be neutral at physiological pH. Since PMBz is a noncovalent inhibitor, we restored the catalytic triad by disrupting the covalent bond between APPA and Ser195, and transferred the proton from His57 to Ser195. Keeping the crystal coordinates fixed, only the hydrogen atoms of this initial configuration were energy minimized. Here, we employed 2000 steps of the steepest descent algorithm followed by 5000 steps of the conjugate-gradient algorithm. At the end of this procedure, the maximum derivative was less than 0.01 kcal mol-1 Å-1. To alleviate bad van der Waals contacts, we then fully optimized the whole system. In this step, we employed 5000 steps of the steepest descent algorithm followed by 10000 steps of the conjugate-gradient algorithm. At the end of this procedure, the maximum derivative was less than 0.01 kcal mol-1 Å-1. Then, the hirugen molecule and all the residues and water molecules of APPA.thrombin.hirugen complex beyond a sphere (R ) 16 Å) centered on the C1 atom of APPA (Figure 2) were removed. This simplification reduced the number of amino acids to 123. Acetyl and N-methylamine blocking groups were used to cap the truncation points. Since we included in our model all the residues with atoms within a radius of 16 Å, our complex model has a radius of approximately 19 Å. Therefore, a cap of water molecules of 19 Å was added to the space of the sphere centered at the C1 atom of APPA (305 water molecules). To

Guimara˜es and Bicca de Alencastro evaluate the influence of the solvent radius in the calculated free energies, the complex was solvated with a smaller cap of water molecules (R ) 16 Å, 155 water molecules). The left side of the thermodynamic cycle, shown in Figure 1, was used to calculate the relative free energy of hydration. To do that, it is only necessary to compute ∆Atr(sol) and ∆Atr(g). To compute ∆Atr(sol), the PMBz structure was centered in a periodic cubic box of 20 Å of length (∼250 water molecules) using periodic boundary conditions. To calculate the relative free energy of binding, which is done by the application of the right part of the thermodynamic cycle, we need to calculate ∆Atr(comp) and ∆Atr(sol). As show by Essex and Jorgensen,22 the use of a spherical cap of water rather than the conventional periodic boundary conditions affects the calculated free energies of hydration in simple systems. Therefore, to cancel the errors introduced by the application of a spherical cap of water in the solvated complex model, we employed an identical approximation for the calculation of ∆Atr(sol), necessary for the calculation of ∆∆Abind. In this case, the aqueous phase was obtained by surrounding PMBz with a water sphere of 16 and 19 Å, centered at the C1 atom of PMBz (same atom position of APPA, see Figure 2), including a total of 558 and 944 water molecules, respectively. The initial preparation of all systems consisted of an energy minimization calculation (2000 steps of the steepest descent algorithm followed by 5000 steps of the conjugate-gradient algorithm), followed by 20 ps of a molecular dynamics (MD) run, to which we applied the NVT ensemble. We kept the solute atoms fixed throughout the whole preparation process. This step was performed in order to relax the water molecules in the aqueous phase and in the solvated complex models to the solute potential. Then, all the coordinates of the systems were fully optimized. Due to the absence of the discarded residues, the methyl carbon atom of the acetyl and N-methylamine blocking groups in the solvated complex model were kept fixed to avoid an anomalous behavior in the protein structure. At this point, we employed 5000 steps of the steepest descent algorithm followed by 10000 steps of the conjugate-gradient algorithm. At the end of this procedure, the maximum derivative was less than 0.01 kcal mol-1 Å-1. Next, we performed the transformation simulations. Four successive mutations inside the complex model transformed APPA into PMBz, generating three new candidates (benzamidine derivatives) for the thrombin inhibition. The results of the relative free energy of hydration and binding of these candidates will be reported elsewhere.23 Each transformation simulation was achieved by a MD run of 1.2 ns, giving a total of 4.8 ns to mutate APPA into PMBz. The 4.8 ns structure was the initial configuration for the transformation of PMBz into Bz inside the solvated complex model. While APPA is neutral at physiological pH, PMBz and Bz have a net charge of +1. Hence, the overall charges of the solvated complex and the aqueous phase systems are the same as those of PMBz and Bz. In every PMBz f Bz transformation, in the gas and aqueous phase, and inside the solvated complex, the FDTI method was employed. This method was also employed to change the C-F equilibrium bond length of fluoromethane in the gas phase. The initial and final states were linearly coupled through the λ parameter. The number of quadrature points, n, used in the calculations was six. The values of ∆λi and the quadrature points (λ1 ) 0.03377, λ2 ) 0.16940, λ3 ) 0.38069, λ4 ) 0.61931, λ5 ) 0.83060, and λ6 ) 0.96623) were automatically calculated by the Gaussian-Legendre quadrature method.24 Good agreement between the values calculated by FDTI and the theoretical

Thrombin Inhibition by Benzamidine

Figure 3. A. The potential energy surfaces at λi, λi - δλ and λi + δλ. B. Graphical representation of the computation of ∆Ai/δλ at each λi, sampling forward (∆A(λi f λi + δλ)/δλ) and backward (∆A(λi f λi - δλ)/δλ). ∆Ai/δλ at λi is the first quantity plus the negative of the second, divided by two.

values (derived from classical statistical thermodynamics) for the mutations of fluoromethane in the gas phase determined the increment δλ to be used to calculate the numerical freeenergy derivatives of all transformations. All transformations were performed using MD simulations in an NVT ensemble (T ) 300 K), and the single topology approach. In the PMBz f Bz transformation, the methyl group was mutated into a hydrogen atom. To keep the number of atoms fixed, the hydrogen atoms of the methyl group were substituted by dummy atoms (Du). The H-Du equilibrium bond length was not shrunk in the simulations. All energy terms (including bonded terms) of the dummy atoms were turned off in the simulations. The equations of motion were integrated every 1.0 fs using the Verlet Leapfrog algorithm.25 All bonds and angles were allowed to move in the simulations. In all simulations, the ensemble average of the Boltzmann factor, e-∆U/kT, was evaluated through a MD run of 100 ps after a 100 ps of equilibration at each λi, giving a total of 1.2 ns of simulation for each transformation. For further analysis, the trajectory was sampled every 1.0 ps. The convergence at each λi was verified by plotting the Boltzmann factor vs time. The FDTI method computes ∆Ai/δλ at each λi sampling both forward (∆A(λi f λi + δλ)/δλ) and backward (∆A(λi f λi - δλ)/δλ). ∆Ai/δλ at λi is the first quantity plus the negative of the second one, divided by two (Figure 3A,B). Similar values for (∆A(λi f λi + δλ)/δλ) and (-∆A(λi f λi - δλ)/δλ) are a second measure of convergence. During the MD runs, the temperature was held at 300 K via the velocity-scaling algorithm26 at the equilibration stage, and via a weak coupling to an external temperature bath with a time constant of 0.1 ps27 at the data collection stage. To save computational time, a spherical residue-based nonbonded interaction cutoff of 10 Å was applied. To “turn off” the interactions smoothly, we employed a quintic spline function

J. Phys. Chem. B, Vol. 106, No. 2, 2002 469 from 8.5 to 10 Å, minimizing, in this manner, discontinuities in the potential energy surface. Whenever an atom moved more than half the length of the buffer region (between 10 and 11 Å), the neighbor list was updated. This ensured that no atoms outside the buffer region were able to move close enough to interact. To prevent the evaporation of water molecules, a halfharmonic restraining potential of 0.5 kcal mol-1 Å-2 was used when the distance between a water oxygen atom and the center of the solvation model exceeded 16 or 19 Å. For the reasons mentioned above, the methyl carbon atoms of the blocking groups were kept fixed in the MD transformation simulations occurring inside the solvated complex model. To check the hysteresis in the calculations, we performed the mutation in the reverse direction, Bz into PMBz, in the aqueous phase (R ) 19 Å) and in the solvated complex model (R ) 19 Å). However, the solvated complex model used in this simulation was constructed from the X-ray crystal structure of the Bz.thrombin‚hirudin(55-65) ternary complex (structure 1DWB),28 taken from the Brookhaven Protein Data Bank. Hirudin(55-65) is an undecapeptide that, like hirugen, binds to the fibrinogen recognition exosite. The construction was performed by selecting the same 123 amino acids used in the forward simulation, and by applying the spherical cap of water molecules of R ) 19 Å centered at the C1 atom of Bz (same atom position of APPA, see Figure 2). The water relaxation, energy minimization, and simulation stages followed the same procedures described above. 4. Results and Discussion Testing the FDTI Method in the Calculation of Jacobian Factor Contributions to the Free Energy. Fluoromethane, behaving as an ideal gas, was selected as a test case. As the internal degrees of freedom of fluoromethane (bond stretching and angle bending) do not correlate at all, the configurational partition function can be written in a decoupled way. Since there are no intra- and intermolecular nonbonded interactions, and the C-F force constant does not change from the initial to the final state, there are no pmf-type and vibrational contributions to the free energy. Therefore, the calculated free energy should contain only the Jacobian factor contributions resulting from changes in the C-F equilibrium bond length. Using eq 5,13 these contributions can be calculated and compared to the calculated values from FDTI.

( )

∆A ) ∆AJ ) -kT ln

req2 req1

2

(5)

Since in the CVFF force field the C-F force constant is very high (KC-F ) 496.0 kcal mol-1 Å-2), the harmonic stretching potential surface is too steep, and thus, two states differing in the equilibrium bond length are extremely orthogonal. Because this case represents an extreme situation, as far as the orthogonality problem is concerned, it is a good test for the FDTI method. All simulations were performed in the vacuum, using the FDTI algorithm in the single topology approach. To test the method in more drastic situations, we increased the final C-F equilibrium bond length (req2) to make more orthogonal the potential energy surfaces. To choose the best overlap between the orthogonal phase spaces, we used a set of different δλ values. Table 1 compares the values calculated by eq 4, setting δλ as 0.005 and 0.0005, to the values obtained from eq 5 in five different final C-F equilibrium bond lengths. As expected, when the C-F equilibrium bond length increases, i.e., the potential

470 J. Phys. Chem. B, Vol. 106, No. 2, 2002

Guimara˜es and Bicca de Alencastro

TABLE 1: Comparison between the Jacobian Factor Contribution to the Free Energy Calculated from Statistical Themodynamics and from the FDTI Method fluoromethane req1(C-F) f req2(C-F) (Å)

∆Acalcda (δλ ) 0.005) (kcal/mol)

∆Acalcd (δλ ) 0.0005) (kcal/mol)

∆Atheorb (kcal/mol)

1.363 f 1.463 1.363 f 1.563 1.363 f 1.663 1.363 f 1.763 1.363 f 1.863

-0.0533 ( 0.0094 -0.0900 ( 0.0216 -0.1767 ( 0.0047 -0.1967 ( 0.0236 -0.2733 ( 0.0047

-0.0600 ( 0.001 -0.1167 ( 0.0125 -0.1700 ( 0.0141 -0.2367 ( 0.0047 -0.2700 ( 0.0282

-0.0844 -0.1633 -0.2372 -0.3068 -0.3726

a The FDTI method was employed to compute the potential energy contribution to the Helmholtz free energy. The number of quadrature points n used in the calculations was six. The value of ∆λi is automatically calculated by the Gaussian-Legendre quadrature method for the total number of quadrature points specified. Each calculated value is an average of three MD runs of 1.2 ns for each transformation (100 ps of equilibration followed by 100 ps of data collection at each λi). The reported errors are the standard deviation from the average ∆Acalcd value. b Calculated using eq 4.

energy surfaces become more orthogonal, the agreement between the calculated and the theoretical values deteriorates. However, despite consistently underestimating the Jacobian factor contribution to the free energy, Table 1 shows that, regardless of the δλ value applied, the FDTI method gives a better estimation than other similar methodologies.15 The Problems of Singularities and Convergence in FEP. As mentioned above, to keep the number of particles constant, we transformed each methyl hydrogen atom into a dummy atom. All energy terms (bonded and nonbonded) of the dummy atoms were turned off in the simulations. When the TI formalism is applied in creating or annihilating atoms, its observable, the derivative of the potential energy with respect to the coupling parameter, contains a singularity. As discussed by Beutler et al.,29 under favorable conditions, the singular regions of the observable correspond to high-energy regions of the phase space, which means that these configurations are never sampled in Monte Carlo (MC) and MD simulations. Nonetheless, depending on the path chosen, the ensemble average of the observable of the TI formalism (integrand of eq 3) may diverge as λ goes to zero or one. This occurs when nonbonded interaction sites, bond stretching and angle-bending terms are created, or annihilated. In the first case, when van der Waals (VDW) interaction sites are created or removed (also called VDW end-point problems), the divergence of the integrand can be eliminated by the use of nonlinear coupling ((1 - λ)nU1(r) + λnU2(r), with n g 4) between the initial and final states.29 Yet, even then, the VDW end-point problems remain because of the instability in MD simulations resulting from very small (if VDW interaction sites are created) or very large λ values (if VDW interaction sites are annihilated). This means that when the VDW interactions decrease, the time step of the simulation has to be decreased at these λ values to avoid the sampling of the high-energy regions of the phase space, and consequently, to avoid the sampling of the singular regions of the observable.29 Another way to avoid singularities, when VDW sites are created or annihilated, is to omit the calculation of the integrand of eq 3 for very small or very large λ values, or to replace it by an extrapolation to very small or very large λ values.12,29 These procedures, however, may introduce an unknown error in the calculation. Thus, to solve VDW end-point problems alternative solutions were proposed by Beutler and co-workers29 and by Zacharias and co-workers.30 When VDW and Coulomb terms are scaled together, attractive Coulomb forces may become larger than repulsive VDW forces, sampling in this manner the singular regions of the observable. This may occur when electrostatic and VDW terms of a polar group are created or annihilated. To overcome this problem, electrostatic and VDW terms should be created or annihilated in a decoupled way.31 When breaking or making bonds the integrand of eq 3 diverges logarithmically as λ approaches the state where the

bond is broken (λ ) 0 or 1), regardless of the dependence of the coupling parameter. As stated by Boresch and Karplus,12 an alternative approach to this problem would be the introduction of a criterion to determine the breaking point of a bond described by a harmonic potential. This could be accomplished by integrating eq 3 from a lower limit λ )  to λ ) 1 when a bond is formed, or from λ ) 0 to λ ) 1 -  when it is broken. The values  and 1 -  ( ) 2πkT/Kf, where Kf is the force constant of the final state if a bond is formed, or  ) 2πkT/Ki, where Ki is the force constant of the initial state, if it is broken) have a simple physical meaning as they correspond to the limit value of the force constant for a bond existence (K ) 2πkT).12 In the gas phase, at λ )  or 1 - , the configurational partition function, which has the dimensionality of [V] in three dimensions, equals the unit volume.12 The omitted range of integration for bond forming or breaking accounts for the fact that the free particle has access to the full volume of the system. In the presence of nonbonded (nb) interactions, the intermolecular degrees of freedom between the solute and the environment are correlated to the bond stretching degree of freedom through the potential of mean force. The same is not true for λ values smaller than  or higher than 1 -  because at these λ values the bond is considered broken. Thus, the omitted range of integration for bond forming or breaking is a constant contribution, regardless of the environment, and it is canceled when the two relevant parts of the thermodynamic cycle for the free-energy difference calculation are combined. There remains only the fact that the free particle has access to the full volume of the system, interacting with the environment. Therefore, to describe bond formation in the presence of nonbonded interactions, the end points (λnb ) 0, λbond )  f λnb ) 1, and λbond ) 1) should be applied. However, at very small (λ ) ) or very large values of λ (λ ) 1 - ), besides the van der Waals end-point problems for the non bonded interaction terms, severe convergence problems of the integrand of eq 3 are found.12 The same arguments apply to harmonic angle bending terms. When the PM formalism is applied, singularity and convergence problems in the observable occur when atoms are created or annihilated. The major difference between the PM and the TI formalisms is that in the former the observable is the Boltzmann factor, e-∆U/kT. To circumvent end-point problems in the PM formalism, a simulation at a λ value outside the critical region of λ values would be the reference state necessary to calculate the free-energy difference over the complete critical region.29 However, this procedure can rarely be applied because, to converge, the PM formula needs an overlap of the initial and final phase spaces. Using the FDTI Method in the Calculation of the Relative Free Energies of Hydration and Binding to Thrombin for Benzamidine and p-Methylbenzamidine. After the validation

Thrombin Inhibition by Benzamidine of the FDTI algorithm, we applied the method to the calculation of the relative free energies of hydration and binding to thrombin for benzamidine and p-methylbenzamidine. To have a better overlap between the phase spaces, we used the smaller δλ value (0.0005). In these simulations, we mutated the methyl group into a hydrogen atom. To do so, one must include the bonded contributions from harmonic stretching and angle bending terms (vibrational, Jacobian factor, and pmf-type contributions). Here, as in the fluoromethane case described above, the FDTI method should give reliable results. In the FDTI method19 as implemented in DiscoVer, the derivatives of the free energy in relation to the coupling parameter (∆Ai/δλ) were numerically computed using the PM formalism at fixed values of λ (see eq 4). Next, a GaussianLegendre quadrature method was applied to integrate numerically the generated function from λ ) 0 to λ ) 1. We used six quadrature points to compute numerically the free-energy derivatives (λ1 ) 0.03377, λ2 ) 0.16940, λ3 ) 0.38069, λ4 ) 0.61931, λ5 ) 0.83060, and λ6 ) 0.96623). The ∆λi values (weight of integration) are 0.0857, 0.1804, 0.2340, 0.2340, 0.1804, and 0.0857, respectively. As this more sophisticated integration method avoids the need to calculate the derivatives of the free energy in the vicinity of the end points (very small or large λ values), the introduction of a cutoff for the existence of a bond is not necessary. However, at small or large λ values, convergence and VDW end-point problems may still exist, as for example, at λ6 ) 0.96623, for the mutation of the methyl group of PMBz into the hydrogen atom of Bz. VDW end-point problems may derive from instabilities in MD simulations. As these instabilities may manifest themselves as poor convergence behavior of the Boltzmann factor (e-[U(λi+δλ)-U(λi)]/kT), both the convergence and the VDW end-point problems at each λi were verified by plotting the Boltzmann factor vs time. As the methyl group is not polar, the decoupling procedure of electrostatic and VDW terms, used to avoid VDW end-point problems in freeenergy simulations, was not necessary. Figure 4 shows the Boltzmann factor plots for the PMBz f Bz transformation both in the gas and the aqueous phases, using periodic boundary conditions and spherical caps of water of radius 16 and 19 Å, and also inside the solvated complex, using spherical caps of water of radius of 16 and 19 Å. Poor convergence of the Boltzmann factor can be viewed at the last λ value (λ6 ) 0.96623) for the PMBz f Bz transformation, regardless of the environment. To check hysteresis in the calculations, we simulated the reverse transformation, Bz f PMBz, in the aqueous phase and in the solvated complex model, applying in both cases a spherical cap of water molecules of R ) 19 Å. Poor convergence at the first λ (λ1 ) 0.0337) for the Bz f PMBz transformation was also observed. For the PMBz f Bz transformation in the gas phase, VDW end-point problems can only result from intramolecular VDW interactions. As the conformational rigidity of PMBz avoids the sampling in VDW repulsive energy regions of the phase space, which causes instabilities in MD simulations, the inferior convergence observed at λ6 for the gas-phase simulation (Figure 4A) must be derived from the bonded terms that are annihilated in the transformation. Since the instantaneous Boltzmann factor can be written as an exponential multiplication of the Boltzmann factors of each individual contribution of the total hybrid potential energy function (eq 6), we can dissect each contribution to the poor convergence of the observable.

e-([U(λi+δλ,t)-U(λi,t)]/kT) ) e-([Ubond(λi+δλ,t)-Ubond(λi,t)]/kT) ) e-([Uangle(λi+δλ,t)-Uangle(λi,t)]/kT) ... (6)

J. Phys. Chem. B, Vol. 106, No. 2, 2002 471

Figure 4. Convergence analysis of the Boltzmann factor at each λi for the PMBz f Bz transformation. A. In the gas phase. B. In the aqueous phase using periodic boundary conditions (PBC). C. In the aqueous phase using a spherical cap of water with 16 Å of radius. D. In the aqueous phase using a spherical cap of water with 19 Å of radius. E. Inside the solvated complex using a spherical cap of water with 16 Å of radius. F. Inside the solvated complex using a spherical cap of water with 19 Å of radius.

472 J. Phys. Chem. B, Vol. 106, No. 2, 2002 Figure 5 shows that the convergence of the Boltzmann factor at λ6, or at λ1 for Bz f PMBz (not shown), is affected mainly by the poor convergence of the individual Boltzmann factors of the harmonic stretching and angle bending terms that are annihilated in the simulations. This is probably due to very low force constants at λ6; the C-H f H-Du, Cp-C-H f CpH-Du (Cp is the CVFF nomenclature for a sp2 carbon atom) and H-C-H f Du-H-Du transformations change the force constant values from 340.6, 44.4, and 39.5 to 11.5 kcal mol-1 Å-2, and 1.5 and 1.3 kcal mol-1 θ-2 at λ6 ) 0.96632, respectively. Therefore, as the restoring force is very small, the intra and intermolecular nonbonded forces significantly distort the disappearing bonds and angles from their equilibrium values. In other words, at each instant of the MD simulation at λ6, a very different configuration of the methyl group is sampled (Figure 6). At λ6, as the force constants of the harmonic stretching and angle bending terms annihilated in the simulations are very small, their angles and bonds are extremely flexible. Thus, in this state, a weak nonbonded force is enough to keep the disappearing atoms away from the high-energy regions of the phase space. This, consequently, avoids the sampling of the singular regions of the observable, and the poor convergence of the Boltzmann factor at λ6, or at λ1 for Bz f PMBz (not shown), is not derived from VDW end-point problems. As discussed above, the limiting value of the force constant for a bond (or angle) existence is 2πkT (∼3.7 kcal mol-1 Å-2 or kcal mol-1 θ-2). Consequently, the limiting value  for the existence of a bond stretching or an angle bending is 2πkT/Ki if it is annihilated, or 2πkT/Kf if it is created. Therefore, eq 4 should be numerically integrated from 0 to 1 -  in the first case, or from  to 1 in the second one. As the limiting value λ )  or 1 -  is a function of the respective bonded force constant, eq 4 should be numerically integrated from 0 to 0.98900 for the annihilation of C-H, from 0 to 0.91564 for the annihilation of Cp-C-H, and from 0 to 0.90517 for the annihilation of H-C-H. Also, since at λ6 ) 0.96623 the force constants of C-H, Cp-C-H, and H-C-H are 11.5 kcal mol-1 Å-2, and 1.5 and 1.3 kcal mol-1 θ-2, respectively, only the last two terms, the angle terms, are already annihilated at this state. Consequently, we eliminated from our calculations the contributions of the disappearing angle bending terms to the Boltzmann factor at λ6, solving most of the convergence problems, and included all other contributions (bonded and nonbonded terms) of the force field. To accomplish that, we calculated the free-energy derivatives including the contributions of the disappearing angle bending terms at λ1 ) 0.03377, λ2 ) 0.16940, λ3 ) 0.38069, λ4 ) 0.61931, and λ5 ) 0.83060. The generated function was integrated using as weights of integration 0.0857, 0.1804, 0.2340, 0.2340, 0.1804, respectively. At λ6 ) 0.96623, we used as weight of integration ∆λ6 ) 0.0857. In other words, the freeenergy derivative function, including all individual Boltzmann factor contributions of the total hybrid potential energy function, was integrated from 0 to 0.91430. The integration from 0.91340 to 1 excluded the Boltzmann factor contributions of the disappearing angle bending terms. One can see that our integration limit is different from the theoretical value for the annihilation of Cp-C-H (1 -  ) 0.91564). As suggested by Boresch and Karplus,12 the use of the same cutoff for λ may result in a cancellation of errors when the two relevant parts of the thermodynamic cycle for the calculation of the double freeenergy difference are combined. However, when a bond is broken, coupling with nonbonded interactions at values of λ

Guimara˜es and Bicca de Alencastro

Figure 5. Convergence analysis of the total Boltzmann factor and the individual Boltzmann factors of the disappearing harmonic stretching and angle bending terms at λ6 for the PMBz f Bz transformation in different environments. A. In the gas phase. B. In the aqueous phase using periodic boundary conditions (PBC), C. In the aqueous phase using a spherical cap of water with 16 Å of radius. D. In the aqueous phase using a spherical cap of water with 19 Å of radius. E. Inside the solvated complex using a spherical cap of water with 16 Å of radius. F. Inside the solvated complex using a spherical cap of water with 19 Å of radius.

Thrombin Inhibition by Benzamidine

J. Phys. Chem. B, Vol. 106, No. 2, 2002 473 TABLE 2: Free-Energy Changes Calculated by the FDTI Method (kcal/mol)a

Figure 6. Graphical representation of the potential energy surfaces at λ6 and λ6 + δλ for a disappearing harmonic term. In each instant of the MD simulation at λ6, a very different configuration of the methyl group is sampled, giving rise to poor convergence of the individual Boltzmann factors of the disappearing harmonic stretching and angle bending terms.

higher than the cutoff values may give nonzero contributions to the relative free-energy difference. These contributions are probably not important for the annihilation of Cp-C-H because the upper limit of the first part of the integral is very close to the theoretical value. The same may not be true for the annihilation of C-H, as the correct upper limit of the integral is very different. Therefore, to test this assumption, we also decided to exclude from our calculations the disappearing harmonic stretching term contributions to the Boltzmann factor at λ6 because they are the other important source of convergence problems. The same type of corrections was applied to the reverse Bz f PMBz transformation. We included in the free-energy derivatives the contributions of the disappearing angle bending terms and the disappearing harmonic stretching terms from λ2 to λ6. At λ1 ) 0.03377, again, we excluded both the contribution of the disappearing angle bending terms, and these same contributions added to the contributions of the disappearing harmonic stretching terms, and integrated the free-energy derivative function from 0 to 0.0857. The free-energy derivative function, including all contributions to the Boltzmann factor, was integrated from 0.0857 to 1. The use of a spherical cap of water affects the calculated free energies of hydration in simple systems.22 In consequence, we used the conventional periodic boundary conditions (PBC) to calculate ∆∆Ahyd(∆Atr(sol-PBC) - ∆Atr(g)). Since we used a spherical cap of water in the solvated complex model, to cancel the errors, we also used it in the calculation of ∆Atr(sol), necessary for the calculation of ∆∆Abind. To probe the influence of the radius of solvation in the relative free energies of binding, we employed two different radius (R ) 16 and 19 Å) to obtain ∆∆Abind(16 Å) and ∆∆Abind(19 Å). Table 2 gives the calculated ∆Atr(g), ∆Atr(sol-PBC), ∆∆Ahyd, ∆Atr(sol-16 Å), ∆Atr(comp16 Å), and ∆∆Abind(16 Å) for the forward transformation, PMBz f Bz. It also gives ∆Atr(sol-19 Å), ∆Atr(comp-19 Å), and ∆∆Abind(19 Å) for both forward and reverse transformations. Moreover, Table 2 shows these quantities corrected as defined above. The errors presented in the Table are the standard deviation of the Boltzmann factor, obtained during the data collection stage of the simulations. As shown in Table 2, the ∆∆Ahyd values are extremely dependent on the inclusion of the contributions of the disappearing harmonic stretching and angle bending terms at λ6. With the inclusion of all contributions to the Boltzmann factor, PMBz is more solvated in water than Bz by 1.47 kcal/mol. With the exclusion of the contributions of the disappearing angle bending terms, ∆∆Ahyd becomes 3.84 kcal/mol (∆∆Ahyd(-a)). However, the exclusion of both contributions (∆∆Ahyd(-b,-a)) makes Bz

transfer

PMBz f Bz

Bz f PMBz

∆Atr(g) ∆Atr(g)(-a)c ∆Atr(g)(-b,-a)d ∆Atr(sol-PBC) ∆Atr(sol-PBC)(-a) ∆Atr(sol-PBC)(-b,-a) ∆Atr(sol-16 Å) ∆Atr(sol-16 Å)(-a) ∆Atr(sol-16 Å)(-b,-a) ∆Atr(sol-19 Å) ∆Atr(sol-19 Å)(-a) ∆Atr(sol-19 Å)(-b,-a) ∆Atr(comp-16 Å) ∆Atr(comp-16 Å)(-a) ∆Atr(comp-16 Å)(-b,-a) ∆Atr(comp-19 Å) ∆Atr(comp-19 Å)(-a) ∆Atr(comp-19 Å)(-b,-a) ∆∆Ahyd ∆∆Ahyd(-a) ∆∆Ahyd(-b,-a) ∆∆Abind(16 Å) ∆∆Abind(16 Å)(-a) ∆∆Abind(16 Å)(-b,-a) ∆∆Abind(19 Å) ∆∆Abind(19 Å)(-a) ∆∆Abind(19 Å)(-b,-a) 〈∆∆Abind(19 Å)〉f,g 〈∆∆Abind(19 Å)(-a)〉 〈∆∆Abind(19 Å)(-b,-a)〉 ∆∆Abindexp10

-10.51 ( 0.10 -8.56 ( 0.08 -2.81 ( 0.07 -9.04 ( 0.35 -4.72 ( 0.27 -2.83 ( 0.24 -9.83 ( 0.36 -5.45 ( 0.28 -3.69 ( 0.25 -10.66 ( 0.37 -5.88 ( 0.29 -4.14 ( 0.25 -11.14 ( 0.37 -6.70 ( 0.29 -4.84 ( 0.26 -9.56 ( 0.37 -4.69 ( 0.30 -2.62 ( 0.26 1.47 ( 0.45 3.84 ( 0.35 -0.02 ( 0.34 -1.31 ( 0.73 -1.25 ( 0.57 -1.15 ( 0.51 1.10 ( 0.74 1.19 ( 0.59 1.52 ( 0.51 0.56 ( 0.54 0.65 ( 0.54 0.92 ( 0.60 0.68

n.c.b n.c. n.c. n.c. n.c. n.c. n.c. n.c. n.c. 10.06 ( 0.33 5.58 ( 0.26 3.70 ( 0.24 n.c. n.c. n.c. 14.49 ( 0.35 (10.04 ( 0.34)e 9.88 ( 0.30 (5.47 ( 0.28) 7.71 ( 0.27 (3.38 ( 0.25) n.c. n.c. n.c. n.c. n.c. n.c. 4.43 ( 0.68 (-0.02 ( 0.67) 4.30 ( 0.56 (-0.11 ( 0.54) 4.01 ( 0.51 (-0.32 ( 0.49)

-0.68

a The FDTI method was employed to compute the potential energy contribution to the Helmholtz free energy. The number of quadrature points n used in the calculations was six. The value of ∆λi was automatically calculated by the Gaussian-Legendre quadrature method for the total number of quadrature points specified. Each calculation was performed using MD runs of 1.2 ns for each transformation (100 ps of equilibration followed by 100 ps of data collection at each λi). The reported errors are the standard deviation of the Boltzmann factor obtained during the data collection stage of the simulation. The calculated ∆Atr(g), ∆Atr(sol-PBC), ∆∆Ahyd, ∆Atr(sol-16 Å), ∆Atr(comp16 Å), and ∆∆Abind(16 Å) values were computed for the forward PMBz f Bz transformation, and ∆Atr(sol-19 Å), ∆Atr(comp-19 Å), and ∆∆Abind(19 Å) for both forward and reverse transformations. (g), free energy computed in the gas phase; (PBC), in the aqueous phase using periodic boundary conditions; (sol 16) and (sol 19) in the aqueous phase using a spherical cap of water with radius of 16 and 19 Å; (comp 16) and (comp 19) inside the solvated complex, using a spherical cap of water with radius of 16 and 19 Å. b n.c.: not calculated. c Excluding contributions of the disappearing angle bending terms at λ6. d Excluding contributions the disappearing angle bending and harmonic stretching terms at λ6. e The values in parentheses were computed using the 1.7 ns equilibrated structure as the initial structure (see text for details). f Averaging the calculated ∆∆Abind for the forward transformation and the value in parentheses for the reverse transformation. g The reported errors were obtained computing (forward value reverse value)/2.

more solvated than PMBz by 0.02 kcal/mol. The analysis of the trajectory in the gas phase shows that the amidino and tolyl groups of PMBz spin around the molecular symmetry axis C2 in opposite directions. This behavior is due to the N-Cp-CpCp torsional forces and to the absence of intermolecular nonbonded forces. At λ6, this movement leads to configurations with similar C-H bond distances that are much larger than the equilibrium bond distance. The potential energy surface of the C-H harmonic stretching is steeper than that of the harmonic angle bending, Cp-C-H and H-C-H. Furthermore, configurations close to the equilibrium bond angle, with the Cp-C-H angles around 90° (equilibrium bond angle, 110°) and the

474 J. Phys. Chem. B, Vol. 106, No. 2, 2002 H-C-H angles around 120° (equilibrium bond angle, 106.4°) are preferentially sampled. Hence, in contrast to our other simulations, in the gas phase, the sampling of these configurations makes the contribution of the disappearing harmonic stretching terms to the Boltzmann factor more significant than the contributions of the disappearing angle bending terms (Figure 5). As a consequence, when the contributions of the disappearing angle bending terms are excluded at λ6, ∆Atr(g) changes from -10.51 to -8.56 kcal/mol, but ∆Atr(sol-PBC) changes from -9.04 to -4.72 kcal/mol, and ∆∆Ahyd changes from 1.47 to 3.84 kcal/mol. If the same behavior is observed for values of λ between the cutoff (0.91340) and the theoretical limit to the existence of Cp-C-H (0.91564), the difference between ∆∆Ahyd and ∆∆Ahyd(-a) will be derived from a nonzero contribution of the disappearing Cp-C-H terms to ∆∆Ahyd. In this case, if the free-energy derivative (including the contributions of the disappearing angle bending terms) is calculated at 0.91564 and integrated from 0 to this value, ∆Atr(g) will change from -10.51 to a value slightly higher than -8.56 kcal/mol, but ∆Atr(sol-PBC) will change from -9.04 to a value far higher than -4.72 kcal/mol. This will reduce ∆∆Ahyd(-a) from 3.84 kcal/mol to a more reasonable level because there are no differences between PMBz and Bz that justify this value to the relative free energy of hydration for these compounds. It should be noted that the noise in the individual Boltzmann factor of the disappearing angle bending terms at λ6, especially in ∆Atr(sol-PBC) (Figure 5B), also contributes to the difference between ∆∆Ahyd and ∆∆Ahyd(-a). Figure 5A,B show that the contributions of the harmonic stretching terms disappearing in the gas-phase transformation are more important to the Boltzmann factor at λ6 than the same type of contribution in the aqueous phase. Moreover, in ∆Atr(g), Figure 5A shows very little noise in the individual Boltzmann factor of the disappearing C-H bonds at λ6. In the same factor of ∆Atr(sol-PBC), Figure 5B shows a small noise especially when compared to the noise in the individual Boltzmann factor of the disappearing angle bending terms. Therefore, these findings suggest that the difference between ∆∆Ahyd(-a) and ∆∆Ahyd(-b,-a) result mainly from nonzero contributions of the disappearing C-H terms to ∆∆Ahyd(-a) at values of λ higher than the cutoff. Even if ∆∆Ahyd includes all nonzero contributions of the disappearing bonded terms, it also includes a significant noncanceled noise at λ6, derived from the disappearing angle bending terms for the transformation in the aqueous phase. ∆∆Ahyd(-a), on the other hand, eliminates this noise but neglects the nonzero contribution of the disappearing Cp-C-H terms between 0.91340 and 0.91564 to ∆∆Ahyd. Therefore, the relative free energy of hydration for PMBz and Bz should lie somewhere between 1.47 and 3.84 kcal/mol. Jorgensen and Nguyen have calculated the hydration of substituted benzenes using Monte Carlo and the FEP technique.32 They showed that the penalty for creating a larger cavity in transforming an aromatic hydrogen into a methyl group is matched by the enhanced VDW interactions between the methyl group and water. If the same applies to our case, the difference in hydration between PMBz and Bz should be due to other reasons. The classical dipole moments calculated with respect to coordinate axes centered on each molecule’s center of mass for Bz and PMBz are 6.94 and 12.18 D, respectively. This suggests that, as a result of stronger electrostatic interactions with the water molecules, PMBz is more solvated than Bz. Regardless of the radius of the spherical cap of water for ∆∆Abind(16 Å) and for ∆∆Abind(19 Å), the results are less

Guimara˜es and Bicca de Alencastro dependent on the inclusion of the contributions of the disappearing angle bendings at values of λ higher than the cutoff of 0.91340. In this case, we may say that the nonzero contributions of the disappearing Cp-C-H to ∆∆Abind are probably not important because the upper limit of the first part of the integral is very close to the theoretical value for the annihilation of these terms. Thus, the difference between ∆∆Abind and ∆∆Abind(-a) is due to noise in the individual Boltzmann factor of the disappearing angle bending terms, which tends to be canceled when the two relevant parts of the thermodynamic cycle are combined (Figure 5C-F). The small discrepancies between ∆∆Abind(-a) and ∆∆Abind(-b,-a) may be due to nonzero contributions (coupling with nonbonded interactions) of the C-H bonds to ∆∆Abind(-a) at values of λ higher than the cutoff of 0.91340, or to noise in the individual Boltzmann factor of the disappearing C-H bonds at λ6. However, it is difficult to determinate with precision the importance of each term to ∆∆Abind(-a). Another important question is the influence of the radius of the spherical cap of water on ∆∆Abind(-a), which eliminates most of the convergence problems and does not exclude bonded contributions to the free energy. The analysis of the PMBz f Bz transformation in the aqueous phase, ∆Atr(sol)(-a), shows very similar results, regardless of the radius employed. These values are, nevertheless, more negative than the transformation using periodic boundary conditions, ∆Atr(sol-PBC)(-a). This is to be expected since the use of a spherical cap of water affects the calculated free energies of hydration in simple systems.22 On the other hand, in the case of the PMBz f Bz transformation in the solvated complex model, ∆Atr(comp)(-a), using a different radius is considerably different. This is reflected in the calculated ∆∆Abind(-a), and is probably due to a worse description of the solvent when a smaller radius is used. The use of a spherical cap of water of 19 Å gives the best calculated free energy of binding. Therefore, we decided to apply this model to check hysteresis. The results show very little hysteresis for the transformation in the aqueous phase, regardless of the inclusion or not, of the disappearing bonded contributions at λ6. This suggests similar noises for the forward and reverse transformations. Moreover, hysteresis is even smaller if the major sources of noise in the calculations are eliminated. For the transformation inside the solvated complex model, hysteresis is very high, even after the elimination of the noise contributions to the Boltzmann factor (Table 2). The derivation of the PM or TI formula assumes that the system is in equilibrium. Besides, the ensemble averages of the Boltzmann factor are reliable only when all relevant configurations of the system are incorporated in the ensemble. This last requirement should not be a problem since the sampling time was very long (100 ps for each λ value). As discussed by van Gunsteren and Mark,11 it is difficult to choose an initial configuration that is a typical state of an equilibrium ensemble. Therefore, an equilibration period is generally applied. This period should be as long as the longest intrinsic relaxation time of the system, which for biomolecular systems is usually very long. Consequently, only a metastable equilibrium is normally reached. As we discussed above, four successive mutations inside the complex model were used to transform APPA (Figure 2) into PMBz. Since the MD time for each transformation was 1.2 ns, a total of 4.8 ns was simulated for the APPA f PMBz transformation. The structure so obtained was the initial configuration for the PMBz f Bz transformation inside the solvated complex model. Figure 7A shows that this system reached at least a metastable equilibrium condition. On the other hand, the initial configuration for the Bz f PMBz transformation was the energy-

Thrombin Inhibition by Benzamidine

Figure 7. A. Time evolution of the potential energy during the forward transformation (PMBz f Bz). B. Time evolution of the potential energy during the reverse transformation (Bz f PMBz), using the energyminimized solvated complex model as the initial structure (see text for details). C. Time evolution of the potential energy during the reverse transformation (Bz f PMBz), using the 1.7 ns equilibrated structure as the initial structure (see text for details).

minimized solvated complex model. As shown in Figure 7B, the equilibration of the structure is very slow, even after 1.2 ns. This result suggests that in order to get reliable calculated free-energy values, longer simulations than those normally applied in the literature should be used. To conclude, the lack of equilibration in the biosystem during the reverse transformation causes the large hysteresis observed. To reduce hysteresis, we selected the last configuration of the Bz f PMBz transformation inside the complex (1.2 ns of MD time); we, then, converted PMBz back again into Bz, and performed an MD run of 500 ps. The last configuration of this MD run (after 1.7 ns of simulation) was used as the new initial configuration for another Bz f PMBz transformation. Figure 7C shows that the potential energy profile is convergent, indicating that the system has also reached a metastable equilibrium. It is interesting to note that, notwithstanding the inclusion of the disappearing bonded contributions at λ6, this procedure leads to reduced hysteresis for the transformation inside the complex (Table 2, values in parentheses). When the two major sources of noise were eliminated, we observed a slightly larger hysteresis in the transformation inside the solvated complex than in the one in the aqueous phase (see Table 2). Comparison of Figure 7A with Figure 7C shows that the potential energy is lower for the forward transformation than for the reverse transformation, even though the system has reached equilibrium in both paths. This suggests that the system is in two different metastable states. We arrived at the same conclusion by the calculation of the RMSd of the CR atoms of the averaged structures at the corresponding λ values for the

J. Phys. Chem. B, Vol. 106, No. 2, 2002 475 forward and reverse simulations; i.e., the calculated RMSd of the averaged structures at λ1 and λ6, λ2 and λ5, λ3 and λ4, λ4 and λ3, λ5 and λ2, and λ6 and λ1 for the forward and the reverse simulations were 4.96, 4.98, 4.91, 4.98, 4.94, and 4.95 Å, respectively. Doniach and Eastman33 have discussed this type of behavior and concluded that as the simulation proceeds, the number of states visited diminishes progressively; thus, the system spends more time in each state before moving to the next one. Since we used two different crystal structures to perform the forward and reverse simulations, and since our initial structures were previously relaxed for different lengths of time, the system assumed different metastable states. Good results, however, were obtained for ∆∆Abind(19 Å) in both cases (Table 2). In the simulations of crambin, Caves and co-workers34 found that up to 5 ns individual trajectories sample only a fraction of the conformational distribution generated by 10 independent simulations, each one spanning 120 ps. In our case, the average of the calculated ∆∆Abind(19 Å) for the forward and reverse simulations includes contributions of the different regions of the phase space sampled in each individual simulation. This means that we are including contributions of the metastable states visited in the forward and reverse simulations. Consequently, as seen in Table 2, the averaged values show an excellent agreement with the experimental value. In addition, as expected, the best agreement is obtained with (〈∆∆Abind(19 Å)(-a)〉), which eliminates most of the convergence problems and does not exclude bonded contributions to the free energy. Analysis of the Interaction Ligand-Enzyme. Experimental results on the trypsin inhibition8 show that Bz binds more strongly to trypsin than PMBz by 0.27 kcal/mol. On the other hand, the experimental10 and the calculated results (Table 2) show that Bz has a lower affinity for thrombin than PMBz. Using free-energy perturbation simulations, Essex et al.9 studied the trypsin inhibition by some benzamidine derivatives and attributed their different affinities to electrostatic effects. These authors suggested that molecules with a higher dipole moment bind more strongly to water than to trypsin, so lowering the affinity constant. Here, we calculated that PMBz is more solvated than Bz by 1.47 to 3.84 kcal/mol. Therefore, PMBz should have a lower affinity constant to thrombin than Bz. However, the binding of PMBz to thrombin is stronger than the binding of Bz by 0.68 kcal/mol (experimental), indicating that PMBz has a net interaction of 2.15 to 4.52 kcal/mol with the protein (0.68 kcal/mol + 1.47 to 3.84 kcal/mol of desolvation penalty). Figure 8A,B shows that the amidino groups of PMBz and Bz are hydrogen bonded to the carbonyl oxygen of Gly219, and form a water-mediated salt bridge with Asp189. Moreover, the carboxylate group of Asp189 interacts with the carbonyl oxygen of Gly219 through a water-mediated hydrogen bond. The methyl group of PMBz is closer to the backbones of Glu192 and Ser195 and to the side chains of His57 and Val213, than the hydrogen atom of Bz is. This suggests that the interaction of PMBz to thrombin has a strong electrostatic character, but also has an important hydrophobic contribution, which is absent in Bz. These hydrophobic interactions should account for the inverted affinity order. 5. Conclusions In this work, we have employed molecular dynamics simulations in conjunction with the finite difference thermodynamic integration (FDTI) method to compute the relative free energy of hydration and binding to thrombin for Bz and PMBz. First, we discussed how the orthogonality problem could be treated

476 J. Phys. Chem. B, Vol. 106, No. 2, 2002

Guimara˜es and Bicca de Alencastro Prof. Elaine Maia for the use of her computational facilities at the University of Brası´lia (Brası´lia, DF, Brazil). References and Notes

Figure 8. A. The last configuration of the PMBz f Bz transformation at λ1 (the hybrid group is ∼97% of a methyl group). B. The last configuration of the PMBz f Bz transformation at λ6 (the hybrid group is ∼97% of a hydrogen atom).

adequately by the FDTI method, showing that this algorithm worked well in the calculation of Jacobian factor contributions to the free energy in a limit case: the transformation of the equilibrium bond distance with a too high force constant, as in the case of the C-F bond. Then, we applied the FDTI method to the calculation of the relative free energies of hydration and binding to thrombin for benzamidine and p-methylbenzamidine. We have demonstrated that problems of singularity and convergence in free-energy calculations can be properly solved by combining the FDTI method with the Gaussian-Legendre quadrature method for numerical integration, associated with the introduction of a physical criterion to determine the breaking point of a bond or angle described by a harmonic potential. This may be checked by the excellent agreement between the experimental and calculated relative free energies of binding. The results show that, despite being more solvated than Bz by 1.47 to 3.84 kcal/mol, PMBz binds stronger to thrombin by 0.68 kcal/mol. We propose that the greater binding to thrombin is a result of stronger electrostatic and hydrophobic interactions between this molecule and the enzyme. Acknowledgment. This research has received partial financial support from the Brazilian agencies CNPq and FAPERJ (Grants E26/151494/99 and E26/151898/2000). C.R.W.G. would like to acknowledge FAPERJ/Brazil for a scholarship. We are indebted to Prof. Edyr Rogana and Prof. Marcos dos MaresGuia for providing p-methylbenzamidine. The authors also thank

(1) Kimbal, S. D. Curr. Pharm. Des. 1995, 1, 441. (2) Tucker, T. J.; Brady, S. F.; Lumma, W. C.; Lewis, S. D.; Gardell, S. J.; Naylor-Olsen, A. M.; Yan, Y.; Sisko, J. T.; Stauffer, K. J.; Lucas, B. J.; Lynch, J. J.; Cook, J. J.; Stranieri, M. T.; Holahan, M. A.; Lyle, E. A.; Baskin, E. P.; Chen, I.-W.; Dancheck, K. B.; Krueger, J. A.; Cooper, C. M.; Vacca, J. P. J. Med. Chem. 1998, 41, 3210. (3) Jones-Hertzog, D. K.; Jorgensen, W. L. J. Med. Chem. 1997, 40, 1539. (4) Babine, R. E.; Bender, S. L. Chem. ReV. 1997, 97, 7 (5), 1359. (5) Kollman, P. A. Chem. ReV. 1993, 93, 2395. (6) Zwanzig, R. J. Chem. Phys. 1954, 22, 1420. (7) Singh, U. C.; Brown, F. K.; Bash, P. A.; Kollman, P. A. J. Am. Chem. Soc. 1987, 109 (6), 1607. Jorgensen, W. L.; Buckner, J. K.; Boudon, S.; Tirado-Rives, J. J. Chem. Phys. 1988, 89 (6), 3742. Jorgensen, W. L.; Blake, J. F.; Buckner, J. K. Chem. Phys. 1989, 129, 193. Essex, J. W.; Reynolds, C. A.; Richards, W. G. J. Chem. Soc., Chem. Commun. 1989, 1152. Dang, L. X.; Merz, K. M., Jr.; Kollman, P. A. J. Am. Chem. Soc. 1989, 111, 8505. Hirono, S.; Kollman, P. A. J. Mol. Biol. 1990, 212, 197. Jorgensen, W. L.; Briggs, J. M.; Contreras, M. L. J. Phys. Chem. 1990, 94, 1683. Caldwell, J. W.; Agard, D. A.; Kollman, P. A. Proteins 1991, 10, 140. Dunn, W. J., III; Nagy, P. I. J. Comput. Chem. 1992, 13 (4), 468. Sun, Y.; Spellmeyer, D.; Pearlman, D. A.; Kollman, P. A. J. Am. Chem. Soc. 1992, 114, 6798. Myamoto, S.; Kollman, P. A. Proteins 1993, 16, 226. Miyamoto, S.; Kollman, P. A. Proc. Natl. Acad. Sci. U.S.A. 1993, 90, 8402. Pang, Y.-P.; Kollman, P. A. Perspect. Drug DiscoVery Des. 1995, 3, 106. Jorgensen, W. L.; Tirado-Rives, J. Perspect. Drug DiscoVery Des. 1995, 3, 123. Rao, B. G.; Kim, E. E.; Murcko, M. A. J. Comput.-Aided Mol. Des. 1996, 10, 23. Rastelli, G.; Costantino, L.; Vianello, P.; Barlocco, D. Tetrahedron 1998, 54, 9415. Pitera, J.; Kollman, P. A. J. Am. Chem. Soc. 1998, 120, 7557. McCarrick, M. A.; Kollman, P. A. J. Comput.-Aided Mol. Des. 1999, 13, 109. Best, S. A.; Merz, K. M., Jr.; Reynolds, C. H. J. Phys Chem. B 1999, 103, 714. Eriksson, M. A. L.; Pitera, J.; Kollman, P. A. J. Med. Chem. 1999, 42, 868. (8) Mares-Guia, M.; Nelson, D. L.; Rogana, E. J. Am. Chem. Soc. 1977, 99, 2331. (9) Wong, C. F.; McCammon, J. A. J. Am. Chem. Soc. 1986, 108, 3830. Essex, J. W.; Severance, D. L.; Tirado-Rives, J.; Jorgensen, W. L. J. Phys. Chem. B 1997, 101, 9663. Radmer, R. J.; Kollman, P. A. J. Comput.Aided Mol. Des. 1998, 12, 215. (10) Guimara˜es, C. R. W.; Fraga, C. A. M.; Barreiro, E. J.; de Alencastro, R. B. Presented at the 40th Sanibel Symposium, St. Augustine, Florida, 2000. (11) van Gunsteren, W.; Mark, A. E. Eur. J. Biochem. 1992, 204, 947. (12) Boresch, S.; Karplus, M. J. Phys. Chem. A 1999, 103, 103. (13) Boresch, S.; Karplus, M. J. Chem. Phys. 1996, 105, 5145. (14) Boresch, S.; Karplus, M. J. Phys. Chem. A 1999, 103, 119. (15) Pearlman, D. A.; Kollman, P. A. J. Chem. Phys. 1991, 94, 4532 (16) Gough, C. A.; Pearlman, D. A.; Kollman, P. A. J. Chem. Phys. 1993, 99, 9103. (17) Ryckaert, J. P.; Cicotti, G.; Berendsen, H. J. C. J. Comput. Phys. 1977, 23, 327. (18) Severance, D. L.; Essex, J. W.; Jorgensen, W. L. J. Comput. Chem. 1995, 16, 311. (19) Mezei, M. J. Chem. Phys. 1987, 86, 7084. (20) Dauber-Osguthorpe, P.; Roberts, V. A.; Osguthorpe, D. J.; Wolff, J.; Genest, M.; Hagler, A. T. Proteins: Struct., Funct., Genet. 1988, 4, 31. (21) Chen, Z.; Li, Y.; Mulichak, A. M.; Lewis, S. D.; Shafer; J. A. Arch. Biochem. Biophys. 1995, 332 (1), 198. (22) Essex, J. W.; Jorgensen, W. L. J. Comput. Chem. 1995, 16, 951. (23) Guimara˜es, C. R. W.; Alencastro, R. B. Int. J. Quantum Chem. 2001, 6, 713. (24) Press, W. H.; Flanney, B. P.; Teukolsky, S. A.; Vetterling, W. T. Numerical Recipes. The art of Scientific Computing; Cambridge University Press: CV Cambridge, UK, 1986. (25) Hockney, R. W. Methods Comput. Phys. 1970, 9, 136. (26) Woodcock, L. V. Chem. Phys. Lett. 1971, 10, 257. (27) Berendsen, H. C. J.; Postma, J. P. M.; van Gunsteren, W. F.; DiNola, A.; Haak, J. R. J. Chem. Phys. 1984, 81, 3684. (28) Banner, D. W.; Hadva´ry, P. J. Biol. Chem. 1991, 266, 20085. (29) Beutler, T. C.; Mark, A. E.; van Schaik, R. C.; Gerber, P. R.; van Gunsteren, W. F. Chem. Phys. Lett. 1994, 222, 529. (30) Zacharias, M.; Straatsma, T. P.; McCammon, J. A. J. Chem. Phys. 1994 100 (12), 9025. (31) Rao, B. G.; Singh, U. C. J. Am. Chem. Soc. 1989, 111, 3125. (32) Jorgensen, W. L.; Nguyen, T. B. J. Comput. Chem. 1993, 14, 195. (33) Doniach, S.; Eastman, P. Curr. Opin. Struct. Biol. 1999, 9, 157. (34) Caves, L. S. D.; Evanseck, J. D.; Karplus, M. Protein Sci. 1998, 7, 649.