Thermodynamic Cycle-Perturbation Study of the Binding of Trifluoroacetyl Dipeptide Anilide Inhibitors with Porcine Pancreatic Elastase B O G D A N LESYNG* and EDGAR F. MEYER, JR.' T r x a s A X M University, Biochemistry and Biophysics Department, College Station, Texas 77843
SYNOPSIS
The variety of results of crystallographic studies of the serine proteases complexed with isocoumarin inhibitors presents a challenging problem to modeling methods and molecular energetics. Therefore, the thermodynamic cycle-perturbation technique has been used to study a model system of elastase and two peptidic inhibitors. Using the program AMBER, the technique correctly predicts changes of the binding constants for the trifluoroacetyl dipeptide inhibitors in comparison with available experimental ( kinetic and crystallographic) data. However, the absolute values obtained are shown to be sensitive to the specific electrostatic interaction poiential parameters used in the simulations. The reader and user are cautioned that thermodynamic cyle-perturbation results may be too optimistic by underestimating the accuracy of free energy values. This is especially a matter of concern for those cases where a direct comparison with experimental values is not possible, viz., ( 1)the simulation of binding of novel compounds, ( 2 ) structurally uncertain binding sites, or ( 3 ) structurally different binding modes. With our best 4-31G* ESP (electrostatic potential) charges we were able to reproduce experimentally determined free energy differences ( A A A ) with an accuracy of about 1.5 kcal/mol. Dynamically induced structural changes in the binding site of elastase, and particularly changes in hydrogen-bond patterns of the binding site, are also reported.
INTRODUCTION Modeling and drug design efforts are increasingly attracting the attention of synthetic chemists as a means for creating new or improved compounds. While generally successful, such efforts are not without difficulties derived both from the methods themselves and then compounded by the choice of the system being analyzed. A well-studied but still enigmatic system is the serine proteinases and especially elastase. Recent reviews discuss the chemistry, biochemistry, and structure of elastase and related serine proteinases. In addition to the beneficial physiological properties of elastase ( digestion, tissue remodeling and turnover, defense c 1990 .John Wiley & Sons, Inc. CCC OOOfi-3525/90/7-80773-08 $04.00 Biopolymers, Vol. 30, 773-780 (1990) * On leave from University of Warsaw, Department of Biophysics, Poland. T o whom correspondence should be addressed. +
against bacterial infection), pathological conditions (pancreatitis, pulmonary emphysema, etc.) are the direct result of out-of-control elastases; no drugs are currently available to treat these diseases. It therefore is the objective of both experimental and computational / simulation investigations that the mode of action and detailed specificity of binding be better understood. Crystallographic structure analyses of human leucocyte elastase ( HLE ) and porcine pancreatic elastase ( P P E ) ' demonstrate that while amino acid sequences agree to only about 40%, the active site regions of both enzymes can be superimposed to a rms fit of 0.3 A.4 Many more P P E complexes have been studied crystallographically than those of HLE. With this growing data base, P P E is a good subject for theoretical investigations, particularly when we want to compare experimental and simulation results, and among others if we are interested in the determination of binding constants using molecular dynamics ( M D ) simulation techniques. One should point out that MD techniques 773
774
LESYNG AND MEYER
are the object of intense development, particularly when using these methods for a new class of systems /problems and with the increasing speed and availability of supercomputers or dedicated workstations. It is important, therefore, that computational results be directly compared with experimental results, as is done here. In a continuation of our previous study in which we simulated the Asp102 + AsnlO2 mutation in the active site of PPE," we have used the thermodynamic-cycle perturbation technique see (e.g., Refs. 7 and 17), to simulate a mutation in a P P E inhibitor complex. In particular, the purpose of this study was
perturbation of the inhibitor-this should provide information on the rigidity and robustness of the binding site of elastases, subject t o MD treatment. Because regular MD simulations for elastase with and without the amide form of the dipeptide inhibitor had already been performed in this laborathe present study is a natural extension of our previous investigations, and could be performed with moderately low additional computational costs on the CRAY X - M P computer.
+
SYSTEMS A N D METHODS 1. to simulate binding free energy changes for
two representatives of a known class of trifluoroacetyl dipeptide anilides ( TFA inhibitors), previously used for mapping of reception subsites of P P E and HLE'; 2 . to check the reliability of the simulation technique when applied to a drug design project; 3 . and specifically to predict how an ester form of a trifluoroacetyl-dipeptide-dimetylanilide inhibitor would bind to elastase (crystallographic data not available); as well as 4. to answer the question whether there is any noticeable structural change in the binding site of the enzyme in the presence of a strong
Two systems have been investigated. T h e first one consists of P P E with the amide form of the TFA inhibitor (see Figure 1) and a cup of 402 water molecules that covers the binding site and simulates the full solvation of the active site region. Crystallographic data" have been used to generate initial coordinates of elastase and the inhibitor in the binding site. T h e water cup has a radius of 18.5 A and was centered a t His-57 N d . T h e oxygen atom positions of the most external (3.5-A thickness) shell of the cup have been fixed in order to prevent evaporation of the water molecules in the course of the MD simulation; these molecules were allowed to rotate freely.
0
;i /
Mutation ,____._.__ , --.------
F
I
F-
region, N-H replaced wlth 0
CHCH,
F
: c i II 0
H
II 0
7'
i"
CH3
Figure 1. Amide form of the TFA dipeptide anilide inhibitor and its mutation site.
P P E THERMODYNAMIC CYCLE PERTURBATION
Independently, a simulation has been performed for the inhibitor within a water droplet. T h e spherical water cluster with the same radius and a “frozen” shell with the same thickness contains a total of 759 water molecules. In the present study a n united atom mode113,14for carbon atoms was employed. Bound hydrogen atoms are omitted in the interaction potential function, as they contribute only to the masses of the atoms and to their electrostatic and van der Waals interactions. Within this model the enzyme inhibitor cup of water system contains 3958 atoms, and the inhibitor plus the water droplet comprises 2313 atoms. The free energy difference has been determined using a thermodynamic-cycle perturbation approach, slow growth method7sL5;the backbone ) N H group of the inhibitor ( II) has been changed into --0--,which resulted in generation of (I2) (Figure 1) . Only the potential energy term was made a function of a coupling parameter, A. The ratio of binding constants has been determined by means of AA3 and AA4 free energy calculations (compare Figure 2; note that AAl - AA2 = AA:{ - AA4; the latter difference can be more readily calculated; real “horizontal” processes are diffusion dependent, thus too slow for current computational methods). The AMBER 3.0 library with its standard parametrization6 was employed; a residual based cutoff distance was chosen. The original 2.5-A resolution P P E structure” was used for minimization; the norm of the energy gra-
+
+
~~
I
AA.,
-
1
AA3
AA
=
AA3
-
AA4
AA4
Figure 2. Thermodynamic cycle used to compute the relative free energy binding for the amide ( I l ) and ester (I,) forms of the TFA inhibitor to PPE.
775
dient was less than 0.2 kcal/mol/A. Prior to the MD perturbation simulation the model systems were thermalized and equilibrated using the following procedure: Atoms were assigned initial velocities from a Maxwellian distribution at a temperature 285 K, i.e., below the reference thermal bath temperature (300 K ) . Pulses of kinetic energy were provided to the system, simulating the dynamics in the microcanonical ensemble 15 times for 0.2-ps intervals ( i n total 4.0 p s ) . At this moment, additional 2.0-ps check simulations for each of the systems were performed, showing that the temperature of the systems was stabilized a t about 295 K. In the next step the systems were coupled to a thermal bath (300 K ) with a temperature relaxation time of 0.05 ps. Equilibration a t a constant temperature required 40 ps. The potential energy of the systems stabilized after ca. 30 ps. The equations of motion were integrated with a time step of 0.001 ps (1.0 f s ) . The masses of the hydrogen atoms were increased to 10; this lowers librational motions of these atoms and makes the numerical integration procedure more stable when used with the above time step. Before starting regular simulation the final equilibration coordinates were inspected on the PS-300 graphics display, and in the case of the enzyme inhibitor complex, were compared with the available crystallographic data.” The molecule conserved its shape and no pronounced backbone deformations were noticed. In the case of the inhibitor-enzyme complex, a 40-ps thermodynamic perturbation procedure was started after a regular 100-ps MD simulation (not reported in this study, described in Ref. 10; it is being analyzed and was undertaken to seek a detailed analysis of the binding and active site dynamics of elastase). In the case of the inhibitor in water, the perturbation was performed using the same thermalization and equilibration procedure described above. After the perturbation procedures were completed, which resulted in the “mutation” of the amide group of the inhibitor into a n ester, the direction of the perturbations was reversed using coordinates and velocities from the restart files. The backward simulation also required 40 ps. The electrostatic interactions have been calculated based on the charges that reproduce the molecular electrostatic potential (ESP charges). In order to have a set of electrostatic interaction parameters complementary to the standard AMBER parametrization, the SCF-LCAO method with its ab initio formulation has been app1ied.l7 The STO-3G E S P charges for the whole inhibitor have been re-
776
LESYNG AND MEYER
H
0.313
Changes
i
,c,~ .,,,
of the ESP 4-31G. charges
o.6,6
during "mutatlon"
0.248
\
i
0.526
F: 0 . 5 1 3
,$ 0 . 6 1 6
Changes
f: 0 . 6 2 3
/I
of the ESP STO-3G charges
I durlng "mutatlon"
Figure 3. Main changes of the ESP charges in the mutation site of TFA. Absolute values of the changes a t the other atoms are smaller than 0.002.
ported in our previous study." In addition, in the present study, the 4-31G* E S P as well as the STO3G ESP charges have been calculated for a molecular fragment, which modeled the mutation site displayed in the window in Figure 1. The molecular fragment was substituted with =O, -H, -H, and =CH2, starting from the carbonyl group of the inhibitor and going clockwise. T h e calculations have been performed for the amide and the ether forms. Then, the charges of the inhibitor mutation site have been corrected for the difference between the 4-31G* and the STO-3G charges of the molecular fragment under consideration. This procedure resulted in the model of the inhibitor, which reproduces the 4-31G * ESP changes in the mutation site of the inhibitor when going from the amide to the ether form. The changes of the 4-31G* ESP charges as well as changes of the original STO-3G ESP charges in the mutation site, when perturbing I I into 12,are given in Figure 3. During the equilibration, all the bonds and data collection periods were constrained using the SHAKE a l g ~ r i t h m ,with ' ~ a very small tolerance set to 0.0001 A.
RESULTS AND DISCUSSION The free energy differences for the elastase-inhibitor complex, as well as for the inhibitor in the water
droplet (when mutating this inhibitor from X1 to A,), are displayed in Figure 4. The X is the perturbational parameter; X = 0 corresponds to the initial amide form of TFA and X = 1 corresponds to the final ester form. The results show that the enzyme-inhibitor complex is destabilized during mutation of the inhibitor. This destabilization is mostly due to the loss of a n hydrogen bond between the ligand and Val216 in the binding site (compare Figure 1, cf. Table
4
1 00
.
, 0 2
.
, 0 4
.
, 0 6
.
, 0 6
.
, 1 0
1
4
1 2
Lambda
Figure 4. Free energy changes for transformation of the amide form of the TFA inhibitor into its ester form. The upper plot describes the complex of PPE with TFA, forward ( F ) and backward ( B ) transformations; the plot below describes TFA in the water droplet, forward ( F ) and backward ( B ) transformations.
P P E THERMODYNAMIC CYCLE PERTURBATION
2 of Ref. 6 for a summary of active site residues in PPE) . The free energy difference between the initial and the final states is equal to 1.23 kcal/mol. The mutation also results in better hydration of the inhibitor, with the free energy difference equal to -3.20 kcal/mol. In total, the free energy difference is equal to 4.43 kcal/mol (compare Figure 3 ) . This indicates weaker receptor binding of the ester form than the amide form. In both cases the perturbation procedures is reversible. The free energy difference reproduces its initial value with a deviation smaller than 0.2 kcal/mol. Breaking of the hydrogen bond of the )N-H group of the inhibitor and the )C=O group of Val216 is not the only structural change that occurs in the enzyme inhibitor complex. T o our surprise, the mutation also results in changing of interactions of the trifluoroacetyl ( T F A ) acetyl group, which is Zocated at the other end of the inhibitor. In the amide form the carbonyl group of the TFA substituent had formed a n hydrogen bond with H-Nt of Gln-192 (see Figure 5 a ) . This bond was broken and a new bond with the hydroxyl group of Ser-195 was formed (Figure 5 b ) . It is especially noteworthy that perturbation in the reverse direction reproduced the initial hydrogen bond pattern (Figure 5c). The hydrogen-bond lengths are collected in Table I. T h e experimental data show6 that the ratio of the binding constants is equal to about lo2, which implies a free energy difference of about 2.8 kcal/ mol. Although the theoretical method correctly predicts the direction of the AAA magnitude, it overestimates its absolute value. This is partially due to an uncertainty of the interaction parameters in the potential energy function, and to a lesser extent to possible limits of precision of experimental data. TFA dipeptide anilides are strong inhibitors. Their binding constants Kl are of the order of 10-7-10~8 M . and spectrofluorometric techniques with Dixon plot analyses carry with them a n uncertainty that can reach one order of magnitude in the ratio of the binding constants for such strong inhibitors (although smaller errors may be expected in the experimental work referenced here). We attribute the fact that the ester form of the enzyme interacts more strongly with water molecules than the amide form partially to a higher absolute value of the 4-31G* E S P charge a t -0-, with respect to the )N-H group. Another reason is that the perturbation center is located in the immediate neighborhood of a highly hydrophobic group. This is a different situation compared to what has been observed in the case of thermolysin complexed with a pair of phosphonamide and phospho-
777
nate ester inhibitors (where the perturbational center was bound to a hydrophylic group). In the case discussed here, water molecules, when interacting as hydrogen-bond acceptors ( amide case, I ]) , have more limited access to the amide group in comparison t o the situation when they provide a proton to the hydrogen bond (ester case with the oxygen atom characterized by the isotropic ESP charge, 1 2 ) .In this way, the ester form is better hydrated. With the potential force field used in this study, this effect is probably overestimated because of a lack of a lone electron pair for H bonding a t PO-. As the ester oxygen atom contains the lone electron pair, its directional properties should be similar to the directional properties of the amide group when forming hydrogen bonds with water. In this case the interaction energy of the water molecule with the ester domain would depend more strongly on the specific direction in which the molecule approaches the domain. The optimal (linear) direction would be influenced by the presence of the propylphenyl group of TFA (Figures 1 and 5 ) ; the hydration energy would be somewhat smaller than the value obtained in this simulation. This problem is therefore more general. It relates back to the parametrization of the interaction potential function, and suggests that in future simulations all atoms and lone pair orbitals should be included in the simulations. Regardless of the influence of the lone electron pair interactions on macroscopic thermodynamical functions, there is also a question to what extent the basis set selection in SCF computations influences the E S P charges and the electrostatic interactions, and of course the subsequent dynamics of macromolecules and thermodynamical properties of these systems. In order to check this effect we repeated the thermodynamics-perturbation calculations for the inhibitor in water using the STO-3G charges (Figure 3 ) . The free energy difference obtained in this simulation is equal to - 1.15 kcal/mol. It retains the direction of the changes obtained with the 4-31G * basis set, but its absolute value is about three times smaller. In addition, in both cases, the conformation of the inhibitor in water is different. Whereas the 4-31G* case predicts an elongated form of the lysine residue, the STO-3G leads to a more compact form with a weak intramolecular hydrogen bond of the -N ( H3)+ group and the Val carbonyl group, with a distance between the nitrogen and oxygen atoms oscillating around 3.05 A. This is qualitatively a different situation in comparison to what has been observed performing dynamics studies for a small-model N-methylalanylacetamide molecule,
778
LESYNG AND MEYER
I/
VAL
0 9
9
VAL
9
VAL
HIS
SER
‘r/
‘/
-99
VAL
-9
VAL
I
9
VAL
v5 SER
Figure 5 . The TFA inhibitor in the PPE binding site. ( a ) Initial configuration of the TFA amide form before perturbation was started; perturbation parameter Z = 0. This model is drawn as background (thin lines) in figures ( b ) and ( c ) for reference purposes. ( b ) TFA after transformation into its ester form; I = 1. ( c ) TFA after retransformation into the amide form, i.e., a “before-and-after” ( F - B ) description of the simulation of the complex.
PPE THERMODYNAMIC CYCLE PERTURBATION
Table I Hydrogen-Bond Distances in the Binding Site of PPE"
779
(A), After Minimization and During Perturbation Perturbation
Inhibitor Group
Enzyme Group
Minimization Data
Xh (anilide) 0 (LAYS) H - N (LYs) 0 (trifluoroacetyl)
0 (Val 216) H-N (Val 216) 0 (Ser 214) Gln 192)
1.92 1.87 1.83 1.81
X = O(F)
X = l
X = O(B)
2.11 1.76 1.76 1.85
3.34' 1.83 2.33
2.35 1.91 1.79 1.84
d
a The distances are given between acceptor atoms and the proton of donors. The X is the perturbation parameter. F and B denote the forward and backward directions of the perturbation, respectively. X represents the H-N group in the amide form and the -0-atom in the ester form. ' H - N group was mutated into -0-and the proton of the donor thus does not exist. This is therefore the distance between the heavy atoms. The hydrogen bond with Gln 192 was broken, and the new bond with H - 0 Ser 195 was formed. The distance 0-H-0 is equal to 2.21 A
where no important influence of the electrostatic term on dynamical properties of this molecule has been noticed. Most likely for larger, more flexible systems, differences in electrostatic charges and/or different initial conditions (separate thermalization and equilibration procedures for both cases have been employed here) can drive the systems into different regions of configurational space. This is not the place to discuss this case in detail; it is generally acknowledged that the STO-3G basis set does not reproduce electrostatic interactions in a satisfactory way (see, e.g., Ref. 1 9 ) ; however, we wish to stress that this type of computation is sensitive to the parametrization of the electrostatic term, and probably some results reported in the literature may be too optimistic in underestimating of the free energy error. In addition, in the case when the inhibitor/drug is very flexible, its conformation should be probably restricted to a given conformational region, which is interesting from a structural and biological point of view, preventing occupation of other regions that can be reached in a random way during thermalization /equilibration procedure ( s ) ; i.e., this method may not be robust enough t o treat alternate binding modes of several linked functional groups.
CONCLUSIONS Starting with a known structure of the complex P P E : amide form of TFA, the transformation of the inhibitor into its ester form has been simulated computationally using a thermodynamic-cycle perturbation method. The same procedure has been used for the inhibitor in a water droplet.
In both cases, the thermodynamic cycle was reversible, with a small ( < 0.2 kcal/mole) error of closure; this demonstrates the conformational changes of the inhibitor and the binding properties of its ester form, and perhaps changes to limits of precision of the experimental data. The AAA free energy difference is overestimated in comparison to the experimental kinetic data reported by Renaud et al.' most likely due to an overestimation of the hydration free energy for the ester form and perhaps to limits of precision of the experimental data. The presence of the lone electron pair representation in the standard ( monopole-monopole term plus Lennard-Jones type energy contribution ) potential energy function should improve these and other simulation data. We conclude that future drug design projects will require more precise potential energy functions. The free energy data seem t o be sensitive to the electrostatic charges in the potential energy function. We checked this using STO-3G ESP charges in addition to the main simulations that employed 4-31G* E S P charges. The hydration free energy difference changed from -3.20 to -1.15 kcal/mol. It is especially worth noting that the inhibitor adopted a different conformation when using the STO-3G charges. This suggests that in some cases the perturbation procedure should be performed while restricting the system to a given conformational region; otherwise, the system is not able to scan all local conformational minima in a relatively short span of time, introducing a n additional influence on the free energy difference that makes valid comparisons increasingly more difficult as the structures diverge.
780
LESYNG AND MEYER
We did not notice any important global structural changes in the binding site of elastase after perturbation; this supports our hypothesis about the rigidity of elastase, which has been confirmed by a comparison of experimental results.4 In addition, we should like t o call attention t o the fact that present MD methods are unable to simulate the influence of p H upon binding constants, which sometimes is crucial for comparison of theoretical and experimental data, particularly near those p H values a t which polar groups are partially ionized/protonated; this requires a generalization of present MD methods, including in the MD formalism a description of proton dissociation / association processes. When some of the above characteristics are neglected, some thermodynamic-cycle perturbation results reported so far in the literature may be too optimistic in underestimating the error of free energy values.
This research has been supported by the U.S. National Science Foundation (DMB-8517286). Computational resources have been provided by the Pittsburgh Supercomputer Center (CRAY X - M P ) , Dr. J. Dinkel, Associate Provost for Computing at Texas A&M University, and by the Department of Biochemistry and Biophysics at the Texas Agricultural Experimentation Station. We acknowledge Professor J. A. McCammon and Dr. S. Swanson for their insights and useful comments. BL wishes to acknowledge travel support from Ministry of Science and Higher Education (Poland) within the project CPBP 01.06.
REFERENCES 1. Bieth, J . G. (1986) in Regulation of Matrix Accumulation, Mechan, R. D., Ed., Academic Press, New York, pp. 217-320.
2. Bode, W. & Huber, R. ( 1984) in Molecular and Cellular Basis of Digestion, Desnuelle, P., Sjostrom, H. & Noren, O., Eds., Elsevier, New York, pp. 213-234. 3. Powers, J. C. & Harper, J. W. (1986) in Proteinase Inhibitors, Barret, A. J. & Salvensen, G., Eds., Elsevier, Amsterdam. 4. Bode, W., Meyer, E. F., Jr. & Powers, J. C. (1989) Biochemistry 28, 1951-1963. 5. Bode, W., Wei, A.-Z., Huber, R., Meyer, E., Travis, J. & Neumann, S. (1986) E M B O J. 10, 2453-2458. 6. Meyer, E., Cole, G., Radhakrishnan, R. & Epp, 0. (1988) Acta Cryst. B 44, 26-38. 7. McCammon, J. A. & Harvey, S. C. ( 1987) Dynamics of Proteins and Nucleic Acids, Cambridge University Press, Cambridge, England. 8. Lybrand, T., McCammon, J. A. & Wipff, G. (1986) Proc. Natl. Acad. Sci. U S A 83,833-835. 9. Renaud, A,, Lestienne, P., Hughes, D. L., Bieth, J. G. & Dimicoli, J.-L. (1983) J . Biol. Chem. 258,83128316. 10. Carlson, G. M. (1987) Ph.D. dissertation, Texas A&M University. 11. Lesyng, B. & Meyer, E. F. (1987) J. Cornput.-Aided Mol. Design 1, 211-217. 12. Hughes, D. L., Sieker, L. C., Bieth, J . & Dimicoli, J.-L. (1982) J. Mol. Biol. 162,645-658. 13. Weiner, S. J., Kollman, P. A,, Nguen, D. T. & Case, D. A. (1986) J. Comput. Chem. 7,230-252. 14. Singh, U. C., Weiner, P. K., Caldwell, J. & Kollman, P. (1987) A M B E R 3.0: Assisted Model Building with Energy Refinement, Department of Pharmaceutical Chemistry, University of California, San Francisco. 15. Bash, P. A., Singh, U. C., Brown, F. K., Langridge, R. & Kollman, P. (1987) Science 235,574-575. 16. Singh, U. C. & Kollman, P. A. (1984) J . Comput. Chem. 5, 129-145. 17. Ryckaert, J. P., Ciccotti, G. & Berendsen, H. J. C. (1977) J. Comput. Phys. 23,327-341. 18. Pettitt, B. M. & Karplus, M. (1985) J. Am. Chem. SOC.107,1166-1173. 19. Koehler, J . E. H., Saenger, W. & Lesyng, B. (1987) J. Comput. Chem. 8 , 1090-1098. Received February 8, 1990 Accepted March 19, I990
SYNOPSIS
The variety of results of crystallographic studies of the serine proteases complexed with isocoumarin inhibitors presents a challenging problem to modeling methods and molecular energetics. Therefore, the thermodynamic cycle-perturbation technique has been used to study a model system of elastase and two peptidic inhibitors. Using the program AMBER, the technique correctly predicts changes of the binding constants for the trifluoroacetyl dipeptide inhibitors in comparison with available experimental ( kinetic and crystallographic) data. However, the absolute values obtained are shown to be sensitive to the specific electrostatic interaction poiential parameters used in the simulations. The reader and user are cautioned that thermodynamic cyle-perturbation results may be too optimistic by underestimating the accuracy of free energy values. This is especially a matter of concern for those cases where a direct comparison with experimental values is not possible, viz., ( 1)the simulation of binding of novel compounds, ( 2 ) structurally uncertain binding sites, or ( 3 ) structurally different binding modes. With our best 4-31G* ESP (electrostatic potential) charges we were able to reproduce experimentally determined free energy differences ( A A A ) with an accuracy of about 1.5 kcal/mol. Dynamically induced structural changes in the binding site of elastase, and particularly changes in hydrogen-bond patterns of the binding site, are also reported.
INTRODUCTION Modeling and drug design efforts are increasingly attracting the attention of synthetic chemists as a means for creating new or improved compounds. While generally successful, such efforts are not without difficulties derived both from the methods themselves and then compounded by the choice of the system being analyzed. A well-studied but still enigmatic system is the serine proteinases and especially elastase. Recent reviews discuss the chemistry, biochemistry, and structure of elastase and related serine proteinases. In addition to the beneficial physiological properties of elastase ( digestion, tissue remodeling and turnover, defense c 1990 .John Wiley & Sons, Inc. CCC OOOfi-3525/90/7-80773-08 $04.00 Biopolymers, Vol. 30, 773-780 (1990) * On leave from University of Warsaw, Department of Biophysics, Poland. T o whom correspondence should be addressed. +
against bacterial infection), pathological conditions (pancreatitis, pulmonary emphysema, etc.) are the direct result of out-of-control elastases; no drugs are currently available to treat these diseases. It therefore is the objective of both experimental and computational / simulation investigations that the mode of action and detailed specificity of binding be better understood. Crystallographic structure analyses of human leucocyte elastase ( HLE ) and porcine pancreatic elastase ( P P E ) ' demonstrate that while amino acid sequences agree to only about 40%, the active site regions of both enzymes can be superimposed to a rms fit of 0.3 A.4 Many more P P E complexes have been studied crystallographically than those of HLE. With this growing data base, P P E is a good subject for theoretical investigations, particularly when we want to compare experimental and simulation results, and among others if we are interested in the determination of binding constants using molecular dynamics ( M D ) simulation techniques. One should point out that MD techniques 773
774
LESYNG AND MEYER
are the object of intense development, particularly when using these methods for a new class of systems /problems and with the increasing speed and availability of supercomputers or dedicated workstations. It is important, therefore, that computational results be directly compared with experimental results, as is done here. In a continuation of our previous study in which we simulated the Asp102 + AsnlO2 mutation in the active site of PPE," we have used the thermodynamic-cycle perturbation technique see (e.g., Refs. 7 and 17), to simulate a mutation in a P P E inhibitor complex. In particular, the purpose of this study was
perturbation of the inhibitor-this should provide information on the rigidity and robustness of the binding site of elastases, subject t o MD treatment. Because regular MD simulations for elastase with and without the amide form of the dipeptide inhibitor had already been performed in this laborathe present study is a natural extension of our previous investigations, and could be performed with moderately low additional computational costs on the CRAY X - M P computer.
+
SYSTEMS A N D METHODS 1. to simulate binding free energy changes for
two representatives of a known class of trifluoroacetyl dipeptide anilides ( TFA inhibitors), previously used for mapping of reception subsites of P P E and HLE'; 2 . to check the reliability of the simulation technique when applied to a drug design project; 3 . and specifically to predict how an ester form of a trifluoroacetyl-dipeptide-dimetylanilide inhibitor would bind to elastase (crystallographic data not available); as well as 4. to answer the question whether there is any noticeable structural change in the binding site of the enzyme in the presence of a strong
Two systems have been investigated. T h e first one consists of P P E with the amide form of the TFA inhibitor (see Figure 1) and a cup of 402 water molecules that covers the binding site and simulates the full solvation of the active site region. Crystallographic data" have been used to generate initial coordinates of elastase and the inhibitor in the binding site. T h e water cup has a radius of 18.5 A and was centered a t His-57 N d . T h e oxygen atom positions of the most external (3.5-A thickness) shell of the cup have been fixed in order to prevent evaporation of the water molecules in the course of the MD simulation; these molecules were allowed to rotate freely.
0
;i /
Mutation ,____._.__ , --.------
F
I
F-
region, N-H replaced wlth 0
CHCH,
F
: c i II 0
H
II 0
7'
i"
CH3
Figure 1. Amide form of the TFA dipeptide anilide inhibitor and its mutation site.
P P E THERMODYNAMIC CYCLE PERTURBATION
Independently, a simulation has been performed for the inhibitor within a water droplet. T h e spherical water cluster with the same radius and a “frozen” shell with the same thickness contains a total of 759 water molecules. In the present study a n united atom mode113,14for carbon atoms was employed. Bound hydrogen atoms are omitted in the interaction potential function, as they contribute only to the masses of the atoms and to their electrostatic and van der Waals interactions. Within this model the enzyme inhibitor cup of water system contains 3958 atoms, and the inhibitor plus the water droplet comprises 2313 atoms. The free energy difference has been determined using a thermodynamic-cycle perturbation approach, slow growth method7sL5;the backbone ) N H group of the inhibitor ( II) has been changed into --0--,which resulted in generation of (I2) (Figure 1) . Only the potential energy term was made a function of a coupling parameter, A. The ratio of binding constants has been determined by means of AA3 and AA4 free energy calculations (compare Figure 2; note that AAl - AA2 = AA:{ - AA4; the latter difference can be more readily calculated; real “horizontal” processes are diffusion dependent, thus too slow for current computational methods). The AMBER 3.0 library with its standard parametrization6 was employed; a residual based cutoff distance was chosen. The original 2.5-A resolution P P E structure” was used for minimization; the norm of the energy gra-
+
+
~~
I
AA.,
-
1
AA3
AA
=
AA3
-
AA4
AA4
Figure 2. Thermodynamic cycle used to compute the relative free energy binding for the amide ( I l ) and ester (I,) forms of the TFA inhibitor to PPE.
775
dient was less than 0.2 kcal/mol/A. Prior to the MD perturbation simulation the model systems were thermalized and equilibrated using the following procedure: Atoms were assigned initial velocities from a Maxwellian distribution at a temperature 285 K, i.e., below the reference thermal bath temperature (300 K ) . Pulses of kinetic energy were provided to the system, simulating the dynamics in the microcanonical ensemble 15 times for 0.2-ps intervals ( i n total 4.0 p s ) . At this moment, additional 2.0-ps check simulations for each of the systems were performed, showing that the temperature of the systems was stabilized a t about 295 K. In the next step the systems were coupled to a thermal bath (300 K ) with a temperature relaxation time of 0.05 ps. Equilibration a t a constant temperature required 40 ps. The potential energy of the systems stabilized after ca. 30 ps. The equations of motion were integrated with a time step of 0.001 ps (1.0 f s ) . The masses of the hydrogen atoms were increased to 10; this lowers librational motions of these atoms and makes the numerical integration procedure more stable when used with the above time step. Before starting regular simulation the final equilibration coordinates were inspected on the PS-300 graphics display, and in the case of the enzyme inhibitor complex, were compared with the available crystallographic data.” The molecule conserved its shape and no pronounced backbone deformations were noticed. In the case of the inhibitor-enzyme complex, a 40-ps thermodynamic perturbation procedure was started after a regular 100-ps MD simulation (not reported in this study, described in Ref. 10; it is being analyzed and was undertaken to seek a detailed analysis of the binding and active site dynamics of elastase). In the case of the inhibitor in water, the perturbation was performed using the same thermalization and equilibration procedure described above. After the perturbation procedures were completed, which resulted in the “mutation” of the amide group of the inhibitor into a n ester, the direction of the perturbations was reversed using coordinates and velocities from the restart files. The backward simulation also required 40 ps. The electrostatic interactions have been calculated based on the charges that reproduce the molecular electrostatic potential (ESP charges). In order to have a set of electrostatic interaction parameters complementary to the standard AMBER parametrization, the SCF-LCAO method with its ab initio formulation has been app1ied.l7 The STO-3G E S P charges for the whole inhibitor have been re-
776
LESYNG AND MEYER
H
0.313
Changes
i
,c,~ .,,,
of the ESP 4-31G. charges
o.6,6
during "mutatlon"
0.248
\
i
0.526
F: 0 . 5 1 3
,$ 0 . 6 1 6
Changes
f: 0 . 6 2 3
/I
of the ESP STO-3G charges
I durlng "mutatlon"
Figure 3. Main changes of the ESP charges in the mutation site of TFA. Absolute values of the changes a t the other atoms are smaller than 0.002.
ported in our previous study." In addition, in the present study, the 4-31G* E S P as well as the STO3G ESP charges have been calculated for a molecular fragment, which modeled the mutation site displayed in the window in Figure 1. The molecular fragment was substituted with =O, -H, -H, and =CH2, starting from the carbonyl group of the inhibitor and going clockwise. T h e calculations have been performed for the amide and the ether forms. Then, the charges of the inhibitor mutation site have been corrected for the difference between the 4-31G* and the STO-3G charges of the molecular fragment under consideration. This procedure resulted in the model of the inhibitor, which reproduces the 4-31G * ESP changes in the mutation site of the inhibitor when going from the amide to the ether form. The changes of the 4-31G* ESP charges as well as changes of the original STO-3G ESP charges in the mutation site, when perturbing I I into 12,are given in Figure 3. During the equilibration, all the bonds and data collection periods were constrained using the SHAKE a l g ~ r i t h m ,with ' ~ a very small tolerance set to 0.0001 A.
RESULTS AND DISCUSSION The free energy differences for the elastase-inhibitor complex, as well as for the inhibitor in the water
droplet (when mutating this inhibitor from X1 to A,), are displayed in Figure 4. The X is the perturbational parameter; X = 0 corresponds to the initial amide form of TFA and X = 1 corresponds to the final ester form. The results show that the enzyme-inhibitor complex is destabilized during mutation of the inhibitor. This destabilization is mostly due to the loss of a n hydrogen bond between the ligand and Val216 in the binding site (compare Figure 1, cf. Table
4
1 00
.
, 0 2
.
, 0 4
.
, 0 6
.
, 0 6
.
, 1 0
1
4
1 2
Lambda
Figure 4. Free energy changes for transformation of the amide form of the TFA inhibitor into its ester form. The upper plot describes the complex of PPE with TFA, forward ( F ) and backward ( B ) transformations; the plot below describes TFA in the water droplet, forward ( F ) and backward ( B ) transformations.
P P E THERMODYNAMIC CYCLE PERTURBATION
2 of Ref. 6 for a summary of active site residues in PPE) . The free energy difference between the initial and the final states is equal to 1.23 kcal/mol. The mutation also results in better hydration of the inhibitor, with the free energy difference equal to -3.20 kcal/mol. In total, the free energy difference is equal to 4.43 kcal/mol (compare Figure 3 ) . This indicates weaker receptor binding of the ester form than the amide form. In both cases the perturbation procedures is reversible. The free energy difference reproduces its initial value with a deviation smaller than 0.2 kcal/mol. Breaking of the hydrogen bond of the )N-H group of the inhibitor and the )C=O group of Val216 is not the only structural change that occurs in the enzyme inhibitor complex. T o our surprise, the mutation also results in changing of interactions of the trifluoroacetyl ( T F A ) acetyl group, which is Zocated at the other end of the inhibitor. In the amide form the carbonyl group of the TFA substituent had formed a n hydrogen bond with H-Nt of Gln-192 (see Figure 5 a ) . This bond was broken and a new bond with the hydroxyl group of Ser-195 was formed (Figure 5 b ) . It is especially noteworthy that perturbation in the reverse direction reproduced the initial hydrogen bond pattern (Figure 5c). The hydrogen-bond lengths are collected in Table I. T h e experimental data show6 that the ratio of the binding constants is equal to about lo2, which implies a free energy difference of about 2.8 kcal/ mol. Although the theoretical method correctly predicts the direction of the AAA magnitude, it overestimates its absolute value. This is partially due to an uncertainty of the interaction parameters in the potential energy function, and to a lesser extent to possible limits of precision of experimental data. TFA dipeptide anilides are strong inhibitors. Their binding constants Kl are of the order of 10-7-10~8 M . and spectrofluorometric techniques with Dixon plot analyses carry with them a n uncertainty that can reach one order of magnitude in the ratio of the binding constants for such strong inhibitors (although smaller errors may be expected in the experimental work referenced here). We attribute the fact that the ester form of the enzyme interacts more strongly with water molecules than the amide form partially to a higher absolute value of the 4-31G* E S P charge a t -0-, with respect to the )N-H group. Another reason is that the perturbation center is located in the immediate neighborhood of a highly hydrophobic group. This is a different situation compared to what has been observed in the case of thermolysin complexed with a pair of phosphonamide and phospho-
777
nate ester inhibitors (where the perturbational center was bound to a hydrophylic group). In the case discussed here, water molecules, when interacting as hydrogen-bond acceptors ( amide case, I ]) , have more limited access to the amide group in comparison t o the situation when they provide a proton to the hydrogen bond (ester case with the oxygen atom characterized by the isotropic ESP charge, 1 2 ) .In this way, the ester form is better hydrated. With the potential force field used in this study, this effect is probably overestimated because of a lack of a lone electron pair for H bonding a t PO-. As the ester oxygen atom contains the lone electron pair, its directional properties should be similar to the directional properties of the amide group when forming hydrogen bonds with water. In this case the interaction energy of the water molecule with the ester domain would depend more strongly on the specific direction in which the molecule approaches the domain. The optimal (linear) direction would be influenced by the presence of the propylphenyl group of TFA (Figures 1 and 5 ) ; the hydration energy would be somewhat smaller than the value obtained in this simulation. This problem is therefore more general. It relates back to the parametrization of the interaction potential function, and suggests that in future simulations all atoms and lone pair orbitals should be included in the simulations. Regardless of the influence of the lone electron pair interactions on macroscopic thermodynamical functions, there is also a question to what extent the basis set selection in SCF computations influences the E S P charges and the electrostatic interactions, and of course the subsequent dynamics of macromolecules and thermodynamical properties of these systems. In order to check this effect we repeated the thermodynamics-perturbation calculations for the inhibitor in water using the STO-3G charges (Figure 3 ) . The free energy difference obtained in this simulation is equal to - 1.15 kcal/mol. It retains the direction of the changes obtained with the 4-31G * basis set, but its absolute value is about three times smaller. In addition, in both cases, the conformation of the inhibitor in water is different. Whereas the 4-31G* case predicts an elongated form of the lysine residue, the STO-3G leads to a more compact form with a weak intramolecular hydrogen bond of the -N ( H3)+ group and the Val carbonyl group, with a distance between the nitrogen and oxygen atoms oscillating around 3.05 A. This is qualitatively a different situation in comparison to what has been observed performing dynamics studies for a small-model N-methylalanylacetamide molecule,
778
LESYNG AND MEYER
I/
VAL
0 9
9
VAL
9
VAL
HIS
SER
‘r/
‘/
-99
VAL
-9
VAL
I
9
VAL
v5 SER
Figure 5 . The TFA inhibitor in the PPE binding site. ( a ) Initial configuration of the TFA amide form before perturbation was started; perturbation parameter Z = 0. This model is drawn as background (thin lines) in figures ( b ) and ( c ) for reference purposes. ( b ) TFA after transformation into its ester form; I = 1. ( c ) TFA after retransformation into the amide form, i.e., a “before-and-after” ( F - B ) description of the simulation of the complex.
PPE THERMODYNAMIC CYCLE PERTURBATION
Table I Hydrogen-Bond Distances in the Binding Site of PPE"
779
(A), After Minimization and During Perturbation Perturbation
Inhibitor Group
Enzyme Group
Minimization Data
Xh (anilide) 0 (LAYS) H - N (LYs) 0 (trifluoroacetyl)
0 (Val 216) H-N (Val 216) 0 (Ser 214) Gln 192)
1.92 1.87 1.83 1.81
X = O(F)
X = l
X = O(B)
2.11 1.76 1.76 1.85
3.34' 1.83 2.33
2.35 1.91 1.79 1.84
d
a The distances are given between acceptor atoms and the proton of donors. The X is the perturbation parameter. F and B denote the forward and backward directions of the perturbation, respectively. X represents the H-N group in the amide form and the -0-atom in the ester form. ' H - N group was mutated into -0-and the proton of the donor thus does not exist. This is therefore the distance between the heavy atoms. The hydrogen bond with Gln 192 was broken, and the new bond with H - 0 Ser 195 was formed. The distance 0-H-0 is equal to 2.21 A
where no important influence of the electrostatic term on dynamical properties of this molecule has been noticed. Most likely for larger, more flexible systems, differences in electrostatic charges and/or different initial conditions (separate thermalization and equilibration procedures for both cases have been employed here) can drive the systems into different regions of configurational space. This is not the place to discuss this case in detail; it is generally acknowledged that the STO-3G basis set does not reproduce electrostatic interactions in a satisfactory way (see, e.g., Ref. 1 9 ) ; however, we wish to stress that this type of computation is sensitive to the parametrization of the electrostatic term, and probably some results reported in the literature may be too optimistic in underestimating of the free energy error. In addition, in the case when the inhibitor/drug is very flexible, its conformation should be probably restricted to a given conformational region, which is interesting from a structural and biological point of view, preventing occupation of other regions that can be reached in a random way during thermalization /equilibration procedure ( s ) ; i.e., this method may not be robust enough t o treat alternate binding modes of several linked functional groups.
CONCLUSIONS Starting with a known structure of the complex P P E : amide form of TFA, the transformation of the inhibitor into its ester form has been simulated computationally using a thermodynamic-cycle perturbation method. The same procedure has been used for the inhibitor in a water droplet.
In both cases, the thermodynamic cycle was reversible, with a small ( < 0.2 kcal/mole) error of closure; this demonstrates the conformational changes of the inhibitor and the binding properties of its ester form, and perhaps changes to limits of precision of the experimental data. The AAA free energy difference is overestimated in comparison to the experimental kinetic data reported by Renaud et al.' most likely due to an overestimation of the hydration free energy for the ester form and perhaps to limits of precision of the experimental data. The presence of the lone electron pair representation in the standard ( monopole-monopole term plus Lennard-Jones type energy contribution ) potential energy function should improve these and other simulation data. We conclude that future drug design projects will require more precise potential energy functions. The free energy data seem t o be sensitive to the electrostatic charges in the potential energy function. We checked this using STO-3G ESP charges in addition to the main simulations that employed 4-31G* E S P charges. The hydration free energy difference changed from -3.20 to -1.15 kcal/mol. It is especially worth noting that the inhibitor adopted a different conformation when using the STO-3G charges. This suggests that in some cases the perturbation procedure should be performed while restricting the system to a given conformational region; otherwise, the system is not able to scan all local conformational minima in a relatively short span of time, introducing a n additional influence on the free energy difference that makes valid comparisons increasingly more difficult as the structures diverge.
780
LESYNG AND MEYER
We did not notice any important global structural changes in the binding site of elastase after perturbation; this supports our hypothesis about the rigidity of elastase, which has been confirmed by a comparison of experimental results.4 In addition, we should like t o call attention t o the fact that present MD methods are unable to simulate the influence of p H upon binding constants, which sometimes is crucial for comparison of theoretical and experimental data, particularly near those p H values a t which polar groups are partially ionized/protonated; this requires a generalization of present MD methods, including in the MD formalism a description of proton dissociation / association processes. When some of the above characteristics are neglected, some thermodynamic-cycle perturbation results reported so far in the literature may be too optimistic in underestimating the error of free energy values.
This research has been supported by the U.S. National Science Foundation (DMB-8517286). Computational resources have been provided by the Pittsburgh Supercomputer Center (CRAY X - M P ) , Dr. J. Dinkel, Associate Provost for Computing at Texas A&M University, and by the Department of Biochemistry and Biophysics at the Texas Agricultural Experimentation Station. We acknowledge Professor J. A. McCammon and Dr. S. Swanson for their insights and useful comments. BL wishes to acknowledge travel support from Ministry of Science and Higher Education (Poland) within the project CPBP 01.06.
REFERENCES 1. Bieth, J . G. (1986) in Regulation of Matrix Accumulation, Mechan, R. D., Ed., Academic Press, New York, pp. 217-320.
2. Bode, W. & Huber, R. ( 1984) in Molecular and Cellular Basis of Digestion, Desnuelle, P., Sjostrom, H. & Noren, O., Eds., Elsevier, New York, pp. 213-234. 3. Powers, J. C. & Harper, J. W. (1986) in Proteinase Inhibitors, Barret, A. J. & Salvensen, G., Eds., Elsevier, Amsterdam. 4. Bode, W., Meyer, E. F., Jr. & Powers, J. C. (1989) Biochemistry 28, 1951-1963. 5. Bode, W., Wei, A.-Z., Huber, R., Meyer, E., Travis, J. & Neumann, S. (1986) E M B O J. 10, 2453-2458. 6. Meyer, E., Cole, G., Radhakrishnan, R. & Epp, 0. (1988) Acta Cryst. B 44, 26-38. 7. McCammon, J. A. & Harvey, S. C. ( 1987) Dynamics of Proteins and Nucleic Acids, Cambridge University Press, Cambridge, England. 8. Lybrand, T., McCammon, J. A. & Wipff, G. (1986) Proc. Natl. Acad. Sci. U S A 83,833-835. 9. Renaud, A,, Lestienne, P., Hughes, D. L., Bieth, J. G. & Dimicoli, J.-L. (1983) J . Biol. Chem. 258,83128316. 10. Carlson, G. M. (1987) Ph.D. dissertation, Texas A&M University. 11. Lesyng, B. & Meyer, E. F. (1987) J. Cornput.-Aided Mol. Design 1, 211-217. 12. Hughes, D. L., Sieker, L. C., Bieth, J . & Dimicoli, J.-L. (1982) J. Mol. Biol. 162,645-658. 13. Weiner, S. J., Kollman, P. A,, Nguen, D. T. & Case, D. A. (1986) J. Comput. Chem. 7,230-252. 14. Singh, U. C., Weiner, P. K., Caldwell, J. & Kollman, P. (1987) A M B E R 3.0: Assisted Model Building with Energy Refinement, Department of Pharmaceutical Chemistry, University of California, San Francisco. 15. Bash, P. A., Singh, U. C., Brown, F. K., Langridge, R. & Kollman, P. (1987) Science 235,574-575. 16. Singh, U. C. & Kollman, P. A. (1984) J . Comput. Chem. 5, 129-145. 17. Ryckaert, J. P., Ciccotti, G. & Berendsen, H. J. C. (1977) J. Comput. Phys. 23,327-341. 18. Pettitt, B. M. & Karplus, M. (1985) J. Am. Chem. SOC.107,1166-1173. 19. Koehler, J . E. H., Saenger, W. & Lesyng, B. (1987) J. Comput. Chem. 8 , 1090-1098. Received February 8, 1990 Accepted March 19, I990
