Molecular insight into pseudolysin inhibition using the MM-PBSA and LIE methods

Journal of

Structural Biology Journal of Structural Biology 153 (2006) 129–144 www.elsevier.com/locate/yjsbi

Molecular insight into pseudolysin inhibition using the MM-PBSA and LIE methods Olayiwola A. Adekoya a, Nils-Peder Willassen b, Ingebrigt Sylte

a,*

a

b

Department of Pharmacology, Institute of Medical Biology, Faculty of Medicine, University of Tromsø, N-9037 Tromsø, Norway Department of Molecular Biotechnology, Institute of Medical Biology, Faculty of Medicine, University of Tromsø, N-9037 Tromsø, Norway Received 26 June 2005; received in revised form 3 November 2005; accepted 4 November 2005 Available online 5 December 2005

Abstract Pseudolysin, the extracellullar elastase of Pseudomonas aeruginosa (EC: 3.4.24.26) plays an important role in the pathogenesis of P. aeruginosa infections. In the present study, molecular dynamics simulations and theoretical affinity predictions were used to gain molecular insight into pseudolysin inhibition. Four low molecular weight inhibitors were docked at their putative binding sites and molecular dynamics (MD) simulations were performed for 5.0 ns, and the free energy of binding was calculated by the linear interaction energy method. The number and the contact surface area of stabilizing hydrophobic, aromatic, and hydrogen bonding interactions appears to reflect the affinity differences between the inhibitors. The proteinaceous inhibitor, Streptomyces metalloproteinase inhibitor (SMPI) was docked in three different binding positions and MD simulations were performed for 3.0 ns. The MD trajectories were used for molecular mechanics-Poisson–Boltzmann surface area analysis of the three binding positions. Computational alanine scanning of the average pseudolysin–SMPI complexes after MD revealed residues at the pseudolysin–SMPI interface giving the main contribution to the free energy of binding. The calculations indicated that SMPI interacts with pseudolysin via the rigid active site loop, but that also contact sites outside this loop contribute significantly to the free energy of association.
2005 Elsevier Inc. All rights reserved. Keywords: Pseudolysin; Molecular dynamics simulations; Free energy of binding; MM-PBSA calculations; SMPI; Protein–protein interactions; LIE calculations; Pseudolysin inhibition; Computational alanine scanning

1. Introduction Theoretical prediction of protein–protein and protein– ligand binding affinities complements experimental analysis by adding molecular insight into the macroscopic properties measured by the experimental analysis. The experimental binding free energy does not refer to a single conformation of the associated complex and the separated molecular partners in solution, but to ensembles of structures representative of the associated and dissociated states. During the last years, free energy models that consider only the initial and final states of the association process have been developed. Compared with the more

*

Corresponding author. Fax: +4777645310. E-mail address: [email protected] (I. Sylte).

1047-8477/$ - see front matter
2005 Elsevier Inc. All rights reserved. doi:10.1016/j.jsb.2005.11.003

rigorous methods like free energy perturbation (FEP) these methods are less computational expensive making them suitable for a variety of molecular systems. In the present study, two such methods, the molecular mechanics-Poisson–Boltzmann surface area (MM-PBSA) and the linear interaction energy (LIE) methods were used to predict the association free energy and the interaction modes of the proteinaceous metalloproteinase inhibitor SMPI and four low molecular weight inhibitors of pseudolysin. Pseudolysin, the extracellullar elastase of Pseudomonas aeruginosa (EC: 3.4.24.26), is a zinc metalloendopeptidase of the thermolysin family, involved in the pathogenesis of P. aeruginosa infections. The three dimensional (3D) structure of pseudolysin is very similar to that of thermolysin (Thayer et al., 1991), the prototype enzyme of the family. The amino terminal domain consists mainly of antiparallel b strands, while the carboxyl terminal domain is mainly

130

O.A. Adekoya et al. / Journal of Structural Biology 153 (2006) 129–144

a-helical. The structure contains a catalytic zinc ion, a calcium ion and two disulfide bonds, (Cys30–Cys58 and Cys270–Cys297). His140, His144, Glu164 and a water molecule ligate the catalytic zinc ion. Site directed mutagenesis studies have identified Glu141 and His223 as essential for the elastase activity. Experimental and clinical studies showed that pseudolysin cleaves casein, elastin, and synthetic peptides (Morihara, 1995), human IgG (Doring et al., 1981; Holder and Wheeler, 1984), collagen types III and IV (Heck et al., 1986), serum a1-proteinase inhibitor (Morihara et al., 1979) and human bronchial mucosal proteinase inhibitor (Johnson et al., 1982). Studies on mice indicated that pseudolysin also has severe hemorrhagic and muscle damaging activities (Komori et al., 2001), and is involved in various lung infections (Mariencheck et al., 2003). Strong evidences suggest that pseudolysin is implicated in chronic ulcers by degradation of human wound fluids and human skin proteins (Schmidtchen et al., 2003) and corneal infection causing corneal liquefaction which could be sight threatening (Hobden, 2002). Experimental attempts were made to treat rabbit cornea infections using antibiotics and P. aeruginosa elastase specific inhibitors as adjuncts to antibiotics (Kessler and Spierer, 1984; Burns et al., 1990), and some inhibitors showed promising effects (Kessler and Blumberg, 1987). A detailed knowledge about the active site geometry and the interaction modes of known inhibitors is the key for designing new inhibitors of clinical use. Streptomyces metalloproteinase inhibitor (SMPI), from Streptomyces nigrescens TK-23, was the first known proteinaceous inhibitor of metalloproteinases. It inhibits the gluzincin metalloproteinases that also include the thermolysin family. The X-ray crystal structure of SMPI is not known, but the NMR structure (Ohno et al., 1998) shows that the protein (102 amino acids) contains two disulphide bridges. The inhibitor function has been connected to the Cys64–Val65 segment of the Arg60–Ala73 loop. This loop also contains a disulphide bridge between Cys64 and Cys69, resulting in a rigid loop structure (Seeram et al., 1997a; Seeram et al., 1997b). The stronger (several thousands fold) binding affinity of SMPI to pseudolysin compared to other available inhibitors of pseudolysin indicates a huge potential in designing inhibitors mimicking the pseudolysin binding region of SMPI. Such inhibitors may have an immersed therapeutic value. Binding studies indicate that SMPI binds stronger to pseudolysin that to thermolysin (a 100-fold higher binding affinity) (Morihara et al., 1979). Thermolysin is the prototype enzyme of the thermolysin family and at least 26 complexes of thermolysin with low molecular weight inhibitors are deposited in the PDB database (Matthews, 1988). However, very little is known about the binding mode of inhibitors to pseudolysin. The X-ray structure of pseudolysin was published in 1991 (Thayer et al., 1991), and in July 2004 the first X-ray structure of pseudolysin (PDB-acquisition: 1u4g) in complex with a small molecule inhibitor was

deposited in the PDB database. However, neither the interaction mode of SMPI with pseudolysin nor thermolysin is known in detail. Ideally, the X-ray structure of the pseudolysin–SMPI complex would be used to guide a comprehensive functional survey of residues present at the interaction interface. Protein–protein interfaces usually have some residues (
hot spot residues) that give the main contribution to the free energy of binding (Clackson et al., 1998). Identifying these residues is an important starting point for a rational design of small molecule mimics. The aim of the present study was to gain molecular insight into pseudolysin inhibition using molecular dynamics (MD) simulations and theoretical affinity predictions. Three different binding modes of SMPI were explored by 3 ns of MD simulations and MM-PBSA calculations of the free energy of protein–protein association. The three binding modes were also studied by computational alanine scanning of residues at the pseudolysin–SMPI interaction interfaces to determine the hot spot residues responsible for high affinity binding. Further, four small molecule inhibitors were docked at the active site and their interactions were studied by MD simulations for 5 ns and calculations of the free energy of binding using the LIE method. 2. Materials and methods The calculations were performed using four processors on a HP Superdome with 550 MHz CPUs. The X-ray crystal structure of pseudolysin (PDB-acquisition: 1ezm) was used as starting structure. The coordination of the zinc and calcium ions within the pseudolysin structure was described by non-bonded energy terms in the calculations (Terp et al., 2000). Zinc was assigned a formal charge of ˚ and well depth +2.0, van der Waals radius 0.69 A e = 0.014 kcal/mol. Calcium was assigned a formal charge ˚ and well depth of +2.0, van der Waals radius 1.6 A e = 0.1 kcal/mol (Terp et al., 2000). All minimizations and MD simulations were performed using the AMBER 7.0 package (Pearlman et al., 1995; Case et al., 2002) with a TIP3P water model (Jorgensen, 1982). The amber94 force field (Cornell et al., 1995) was used for both ligand bound and free states of the pseudolysin. All molecular systems (bound and free states) were solvated with a cubic box of ˚ qvist, 1996). The water and neutralized by counter ions (A number of counter ions varied in the different systems. The particle mesh Ewald (PME) method was used for the treatment of long range electrostatic interactions (Darden ˚ and et al., 1993). The non-bonded cutoff was set to 9.0 A the SHAKE option (Gunsteren and Berendsen, 1977) was used to constrain the bonds involving hydrogen atoms. The dielectric constant was set to 1 in all calculations. Generally all structures were minimized by the conjugate gradient method for 1000 cycles and equilibrated by 125 ps of MD. The temperature was gradually increased to 300 K during the first 25 ps of the equilibration, and kept constant for the rest of the equilibration. The time step during MD was 1 fs, and the non-bonded pair list was updated

O.A. Adekoya et al. / Journal of Structural Biology 153 (2006) 129–144

every 25 steps. Various parameters (density, temperature, pressure, kinetic energy and potential energy) were monitored during the MD simulations, and were generally stable at the post-equilibration period of the simulations. His223 was treated as protonated since the corresponding histidine in thermolysin (His231) stabilized the transition state during catalytic cleavage by donating a hydrogen bond to the hydrated peptide (Matthews, 1988). 2.1. The LIE method The LIE method uses only the intermolecular interactions of the ligand in the associated and dissociated states ˚ qvist et al., 1994, to predict the free energy of binding (A 2002; Marelius et al., 1998). Two simulations are required, one with the ligand free in solution and one with the ligand bound to the solvated target. The absolute binding free energies of the ligand are calculated as: el DGbind ¼ aDhVvdw l–s i þ bDhVl–s i þ c;

where Æ æ denotes the average non-bonded van der Waals (vdw) and electrostatic (el) interactions between the ligand and its surrounding environment (l–s) obtained from MD or Monte Carlo (MC) trajectories. The D denotes the changes in these averages when transferring the ligand from solution (free state) to the target binding state (bound state). The a and b are weight coefficients for the non-polar and polar binding energy contributions, respectively, while ˚ qvist et al., 1994; Hansson et al., c is a constant term (A 1998). The LIE method was used to calculate the binding free energies using the trajectories between 4.5 and 5.0 ns of MD simulations. 2.2. MM-PBSA method The MM-PBSA method (Kollman et al., 2000) was used to estimate the free energy of association of the pseudolysin–SMPI complexes. For this method, three MD simulations are required for the analysis: (1) pseudolysin, (2) SMPI, and (3) the pseudolysin–SMPI complex. Absolute binding free energies are then estimated by: DG ¼ Gcomplex
Gprotein
Gprotein–ligand . The MM-PBSA method combines explicit solvation models, Poisson–Boltzmann analysis (Gilson and Honig, 1988; Honig and Nicholls, 1995) and non-polar solvation free energy calculations (Sanner et al., 1996) to estimate the free energy of binding. MD simulations of the solute were performed in a periodic box of water with counter ions, using the Particle Mesh Ewald method for representing long-range electrostatics effects (Essmann et al., 1995). Snapshots of the molecular system were sampled during MD, and these were post MD processed after removal of the solvent and counter ions (the catalytic zinc ion and calcium ions were not removed) and used for the MM-PBSA calculations.

131

The average free energy G was calculated by: G ¼ EMM þ GPBSA
TS, where EMM is the average molecular mechanical energy calculated as: EMM ¼ Eint þ Evdw þ Eelec ; Eint ¼ Ebond þ Eangle þ Etors which corresponds to the average internal strain energies in bonds, angles and torsions angles. Evdw is the average van der Waals energy, while Eelec is the average electrostatics energy. GPBSA is the solvation free energy given by: GPBSA ¼ GPB þ GSA . GPB is the electrostatic solvation free energy calculated by solving the Poisson–Boltzmann equation. GSA is the non-polar contribution to the solvation free energy from the solvent-accessible surface area (SASA) (Sitkoff et al., 1994). The electrostatic solvation free energy was calculated with the DelPhi software package (Rocchia et al., 2002), using a dielectric constant of e = 1 for the solute, and e = 80 for the solvent. An 80% boxfill lattice and a ˚ were used in the calculations. grid spacing of 2.0 grids/A The non-polar contribution to the solvation free energy was determined from the equation: GSA ¼ cSA þ b, where SA is the solvent-accessible surface area calculated by the MSMS program. (Sanner et al., 1996) and c and b are ˚ 2 and 0.92 kcal/mol, respectively (Sitkoff 0.00542 kcal/molA ˚. et al., 1994). The probe radius was 1.4 A The solute entropy,
TS was estimated by normal mode calculations (Lamm and Szabo, 1986; Kottalam and Case, 1990) using the Nmode module in the Amber 7.0 package. A distance dependent dielectric function e = 4r was used (r: interatomic distance). A convergence criterion of 0.0001 kcal/mol was used for the energy gradient. 2.3. Calculation of ligand RESP charges The 3D structure of the inhibitors, HONHCOCH2CO– Ala–Gly–NH2 (HRO), benzyloxycarbonyl-L-leucine(ZOH) and benzyloxycarbonyl-L-phenylalanine (POH) (Fig. 1) were constructed using the xLeap program of the AMBER package. The X-ray crystal structure of phosphoramidon (Tlp) (Fig. 1) was extracted from a complex with thermolysin (PDB-acquisition: 1tlp), and used as the starting structure for Tlp. Their electrostatic potentials were calculated using Gaussian 98 (Frisch et al., 2001) and used to derive atomic restrained electrostatic point (RESP) charges (Bayly et al., 1993; Cornell et al., 1993). 2.4. Docking of small molecule inhibitors Docking of small molecule inhibitors was performed with the ICM docking module (Abagyan et al., 1994). The X-ray crystal structure of thermolysin in complex with Tlp, (1tlp) was a guideline for docking of Tlp into pseudolysin. As in the X-ray structure, two oxygen atoms bonded to the phosphor of Tlp coordinated the zinc ion (Fig. 1). The X-ray crystal structure of thermolysin in complex with HONH–benzylmalonyl-L-alanylglycine-p-nitroanilide, (5tln), Cbz–GlyP–(O)–Leu–Leu (ZgP(O)Ll), (6tmn), and Cbz–GlyP–Leu–Leu (ZgPLl), (5tmn), were guiding the

132

O.A. Adekoya et al. / Journal of Structural Biology 153 (2006) 129–144

Fig. 1. The small molecule inhibitors at the active site of pseudolysin. [ ]: Ligand atomic bond; [ ]: Amino acid atomic bond; ]: Hydrogen bond, [ ]: Amino acid residue involved in hydrophobic contact(s); [ ]: Ligand atom involved in hydrophobic contact(s). [ Color coding of inhibitor atoms: oxygen; red, nitrogen; blue; zinc ion; green. The schematic diagrams of pseudolysin–ligand complexes were generated using Ligplot (Wallace et al., 1995).

docking of the hydroxymate containing inhibitors (HRO, ZOH, and POH). As in the X-ray complexes, the two oxygen atoms of the hydroxymate were coordinating the zinc ion. Due to the +2 charge of the zinc ion and several surrounding anionic residues, the zinc atom and its ligated amino acids needed special treatment during MD to prevent distortion of the zinc coordination. Two histidines, a glutamic acid and a water molecule coordinate the zinc ion in the active site of pseudolysin. The X-ray structures of thermolysin–inhibitor complexes show that inhibitors usually replace this water molecule with a carboxylate, hydroxamate or phosphoramidate group. Harmonic ˚ and distance restraints of 2.1– restraints of 50 kcal/A ˚ 2.5 A were therefore used between the two oxygen atoms on the phosphorus of Tlp (Tronrud et al., 1986, 1987), the two oxygen atoms of the hydroxymate (Izquierdo-Martin and Stein, 1992), and the zinc ion (Fig. 1). Further, har-

˚ ) and a distance restraint monic restraints (50 kcal/A ˚ (2.1 A) were applied between zinc and its amino acid ligands: His140 (atom NE2); His144 (atom NE2) and Glu164 (atom OE1). These restraints were applied for the first 100 ps of the simulations and thereafter released. Similar approaches have previously been used by others during MD of zinc containing proteins (Tate et al., 1998; Wasserman and Hodge, 1996). 2.5. Docking of SMPI We have previously studied two binding modes of SMPI to thermolysin (Adekoya et al., 2005). One of those positions, termed
the horizontal arrow head docking (HAHD), was similar to that previously suggested by Tate et al. (1998). This position was based on the X-ray structure of thermolysin in complex with the dipeptide inhibitor Val–Trp (PDB

O.A. Adekoya et al. / Journal of Structural Biology 153 (2006) 129–144

accession code 3TMN) (Holden and Matthews, 1988). The rigid loop of SMPI (Arg60–Ala73) was positioned in the active site of thermolysin with the Cys64–Val65 segment overlapping the dipeptide inhibitor of the complex. The carbonyl oxygen atom of Cys64 (P1) formed tetrahedral coordination with the active site zinc atom, while the methyl group of Val65 (P1 0 ) are located in the S 0 pocket of thermolysin. In the second position, termed
the vertical arrow head docking (VAHD), the carbonyl oxygen of Cys64 still coordinated the zinc ion. However, the number of stabilizing protein–protein surface contacts (hydrophobic contacts and hydrogen bonds) seemed more favorable than in the HAHD complex. The P1 0 residue (Val65) had hydrophobic contacts with the aromatic ring of His142 and with the side chain of Val139. The P2 0 residue (Arg66) was hydrogen bonded with main chain oxygen atoms of Val121 and Val139, while the P4 residue (Lys61) was hydrogen bonded with the side chain oxygen of Asn227. The P5 residue (Arg60) was hydrogen bonded with the main chain oxygen of Gln225. In the present study, SMPI was docked into pseudolysin in similar positions (VAHD and HAHD), but a third position was also considered. In the HAHD position, Val110 (pseudolysin) had hydrophobic interactions with the P1 residue (Cys64) of SMPI. The backbone oxygen of the P1 residue also had hydrogen bonding interaction with Arg198 (pseudolysin). In addition Arg198 interacted by hydrogen bonding with the backbone oxygen of the P2 residue (Thr63). The P2 0 residue (Arg66) interacted by hydrogen bonding with the imidazole ring of His223. The backbone oxygen of the P1 residue (Val65) formed a hydrogen bond with the imidazole ring of His140. Outside the reactive loop, Lys38 (SMPI) was hydrogen bonded with the backbone oxygen atoms of Val110 and Asn112, while Arg44 (SMPI) interacted with the backbone oxygen atom of Ser45. In the VAHD position, the backbone oxygen of the P5 residue (Arg60) interacted by a hydrogen bond with the backbone nitrogen atom of Asp221 (pseudolysin). A terminal nitrogen atom of Arg60 also interacted by hydrogen bonds with the backbone oxygen atom of Gly219. The side chain of the P2 residue (Thr63) interacted by a hydrogen bond with the carboxylic oxygen atom of Asp206. The backbone NH-group of the P3 residue (Val62) interacted by hydrogen bonding with the nitrogen atom in the imidazole ring of His224. Outside the reactive loop, SMPI had several interactions with pseudolysin. The side chain nitrogen atom of Lys31 (SMPI) interacted with the backbone oxygen of Gly107 and was also involved in salt bridges with the carboxylic oxygens of Asp124. The side chain oxygen of Ser27 (SMPI) interacted with the side chain nitrogen atom of Arg108. The backbone oxygen of Asp28 (SMPI) interacted with the backbone NH-group of Arg108. Furthermore, the backbone oxygen atom of Gly36 (SMPI) interacted with the side chain amide nitrogen of Asn112 (pseudolysin). A third position, a variant of the VAHD position, was also explored. Arg66 (SMPI) has been suggested by mutational studies to be important for SMPI–thermolysin

133

complexation (Hiraga et al., 1999). This position was considered to improve the hydrogen bonding interactions of Arg66 observed in the VAHD position. In this position, the terminal nitrogen atom of Arg66 (the P2 0 residue) interacted with the carboxylic acid oxygen atoms of Asp136 and Asp189 (Fig 2). The backbone nitrogen of the P2 residue (Thr63) had a hydrogen bond with the backbone oxygen of His223. The hydroxyl group of Thr63 also had a hydrogen bond with the side chain oxygen of Ser226. The P4 residue (Lys61) was hydrogen bonded with the side chain oxygen of Tyr155. Amino acids outside the reactive loop were also involved in hydrogen bonds with pseudolysin. The backbone oxygen of Asp28 (SMPI) interacted with the NH-group of Arg108. The side chain oxygen atom of Tyr114 (pseudolysin) was hydrogen bonded with Asp35 (SMPI). Lys31 (SMPI) interacted by hydrogen bonds with the backbone oxygen atom of Gly107. 2.6. MD simulations and data analysis MD simulations for a total of 5 ns were performed for the complexes of pseudolysin with the small molecule inhibitors. Coordinates were sampled every 250 steps. The LIE analysis of the inhibitor–protein complexes, and the inhibitors in water (free states) were done by extracting 2000 coordinate sets from the trajectories during 4.5–5.0 ns of MD. The various non-polar (van der Waals) and polar (electrostatic) contributions were calculated by the analyze module in the AMBER suite of programs. These contributions were averaged and inserted into the LIE equation using b = 0.5 and a = 0.16 and c = 0 (Marelius et al., ˚ qvist et al., 2002). 1998; A The HAHD, VAHD and the position 3 complexes of SMPI with pseudolysin were starting complexes for 3.0 ns MD simulations. During the first 100 ps of MD, ˚ ) and a distance restraint harmonic restraints (50 kcal/A ˚ ) were used between the carbonyl oxygen of (2.1–2.5 A Cys64 and the zinc ion. Coordinates were sampled every 250 steps. Three hundred and fifty snapshots of each

Asp189

Arg66 (SMPI)

Asp136

Fig. 2. The position 3 complex showing the hydrogen bonding interactions between Arg66 located at the tip of the rigid loop (Arg60–Ala73) of SMPI (green) and the main chain and side chain oxygen atoms (red) of Asp136 and Asp189 of pseudolysin. The nitrogen atoms are colored blue.

134

O.A. Adekoya et al. / Journal of Structural Biology 153 (2006) 129–144

MD simulations (between 2.3 and 3.0 ns of simulation) were used for the MM-PBSA calculations and for the normal mode analysis. The CARNAL module of the AMBER program package was used to determine possible hydrogen bonds and ˚ between the salt bridges within an atomic distance of 3 A atoms of the inhibitors and pseudolysin in all the complexes. The analysis was performed for the snapshots sampled between 2.3 and 3.0 ns (2800 coordinate sets) for the protein–SMPI complex and for the snapshots between 4 and 5 ns for the small molecule-pseudolysin complexes (4000 coordinate sets). 2.7. Computational alanine scanning Computational alanine scanning is useful for identifying important amino acid residues at protein–protein interfaces that contribute significantly to the free energy of binding. The interactions of such residues can be used as starting points for designing small molecule mimics and are important for the design experimental studies. Binding free energy changes upon alanine mutation were calculated using the ROBETTA server at: http:// robetta.bakerlab.org/alascansubmit.jsp (Kortemme and Baker, 2002; Kortemme et al., 2004). The energy function for calculating the effects of alanine mutation is dominated by Lennard–Jones interactions, solvation and hydrogen bonding interactions. Average structures of the various pseudolysin–SMPI complexes (2800 coordinate sets of each) between 2.3 and 3.0 ns were used for the computational alanine scanning. The free energy function, DG, for calculating effects of alanine mutation on binding free energy of a protein–protein complex is represented by: DG ¼ W attr ELJattr þ W rep ELJrep þ W HBðsc–bbÞ EHBðsc–bbÞ þ W HBðsc–scÞ EHBðsc–scÞ þ W sol Gsol þ W U=W EU=W þ

20 X

ðaaÞ

naa Eref aa ;

aa
1

where ELJattr = Attractive part of a Lennard–Jones potential; ELJrep = A linear distance-dependent repulsive term; EHB(sc–bb) = Orientation-dependent side chain–backbone hydrogen bond potential; EHB(sc–sc)= Orientation-dependent side chain–side chain hydrogen bond potential; Gsol = Implicit solvation model; W = Relative weights of the different energy terms; EU/W (aa) = Amino-acid type (aa) dependent backbone torsion angle propensity; Eref aa = Amino-acid type dependent reference energy, which approximates the interactions made in the unfolded state ensemble; naa = the number of amino acids of a certain type. The effects of alanine replacement were computed both for the protein complex and for the corresponding uncomplexed partners to yield the change in binding energy DDGbind:

WT DDGbind ¼ DGMUT bind
DGbind MUT ¼ ðDGMUT complex
DGpartner

A

WT
ðDGWT complex
DGpartner


DGMUT partner B Þ A


DGWT partner B Þ;

where DGcomplex, DGpartner A, and DGpartner B are the free energy functions of the complex and the unbound partners, and WT and MUT describe wild type and mutant proteins. 3. Results and discussion Pseudolysin is implicated as a virulence factor of P. aeruginosa, an opportunistic pathogen able to cause both local and disseminated infection. Adequate inhibition of pseudolysin is therefore necessary, especially with the emergence of resistance to antibiotics. The basic mechanisms behind the biological effects of endogenous ligands and therapeutic effects of drugs involve an initial process of molecular recognition followed by complex formation between the ligand and the target. Knowledge about the molecular features underlying specific molecular recognition and binding is important for rational design of ligands with potential therapeutic effects. A detailed knowledge of the inhibitor binding site of pseudolysin and the identification of the amino acid residues giving the main contributing to the binding free energies are very important for rational design of new pseudolysin inhibitors. Theoretical prediction of binding free energies is helpful for gaining such information. The FEP approach is still the most important technique for free energy calculations by MD or MC simulations. However, the perturbations (transformations) involved in FEP analysis cannot be too drastic, which limits the structural diversity of the ligands being compared. Extensive conformational sampling is also required to obtain convergent results, and structural flexibility may be a major source for uncertainties. The LIE and MM-PBSA methods consider only the initial and final states of the association process, and are therefore not that dependent of extensive conformational sampling as the FEP approach. But to obtain reliable thermally averaged energies for being used in the analysis, stable complexes must be obtained. Therefore, to assess the stability of the simulations, the RMSD was calculated relative to the average structure calculated from the coordinate sets sampled during the production phase of the simulations (data not shown). The RMSD values indicated that stable complexes were obtained, and that the windows used for the LIE (4.5–5 ns, 2000 sets of coordinates) and MM-PBSA (2.3–3 ns, 350 sets of coordinates) were representative. The RMSD was between 0.4 and ˚ during the period of the simulation used for the 0.85 A analysis. 3.1. The MM-PBSA analysis A more complex relationship between entropic and enthalpic contributions to the free energy of binding is

O.A. Adekoya et al. / Journal of Structural Biology 153 (2006) 129–144

generally anticipated for protein–protein association than for protein–small molecule association. Using the LIE approach requires convergent values for the interactions energies. The quantities of protein–protein interactions can be in the order of several thousand kcal/mol, indicating that extremely long simulations are required to get stable averages. That would result in an inefficient approach with ˚ qvist corresponding high uncertainties in the energies (A et al., 2002). The LIE approach also requires empirical parameters that depend on the nature of the ligand-protein interaction. Protein–protein association very often involves more than one site of association between the molecules that contribute to the free energy of binding. These association sites may be of different chemical nature, and require different parameterization. The MM-PBSA method applies no empirical parameters in the calculations of the free energy. Based on all this, we desired to use the MM-PBSA method to validate the binding free energies of the pseudolysin–SMPI complexes. Binding free energy terms (kcal/mol) of the HAHD, VAHD and position 3 pseudolysin–SMPI complexes including rotational, translational and vibrational entropies of SMPI are shown in Tables 1 and 2. The MM-PBSA calculations suggested that SMPI binds strongest to pseudolysin in the VAHD orientation and weakest in the HAHD orientation. This is in agreement with our previous calculations for thermolysin (Adekoya et al., 2005). However, a complex corresponding to the position 3 complex was not considered in the thermolysin–SMPI study. A combination of the electrostatic contributions calculated by the molecular mechanical energy term and the electrostatic contribution to the solvation free energy calculated by DelPhi, seemed to contribute most to the binding free energy. The electrostatic contribution to the solvation free energy seemed to be the main factor for the stronger binding energy of SMPI in the VAHD complex than in the HAHD and the position 3 complexes (Table 2). The free energies calculated with the MM-PBSA method were overestimated compared with the experimental inhibition constants (Ki)

135

for pseudolysin–SMPI complex, which is 2.54 · 10
12 M (Oda et al., 1979; Seeram et al., 1997a) (Table 2). Several explanations are possible for the overestimation: (1) Incorrect parameterization and geometry of the zinc coordination. (2) The charge of the catalytic zinc. A formal charge of +2 was used for the catalytic zinc atom, which also has been used by others during MD of zinc containing proteins (Donini and Kollman, 2000; Hou et al., 2001; Hou et al., 2002). This might also contribute to errors, since it is obvious that the atomic charge of zinc will depend on the orbitals of the coordinating groups. (3) Inaccurate harmonic approximations of solute entropy in nmode calculations. Nmode analysis give only qualitative estimates of solute entropy (Cheatham et al., 1998). 3.2. Molecular interactions in the pseudolysin–SMPI complexes The binding free energy terms indicated that SMPI binds stronger to pseudolysin in the VAHD than in the HAHD and position 3 complexes, which is in agreement with our previous calculations of the SMPI–thermolysin interactions (Adekoya et al., 2005). The analysis of the complexes revealed that hydrophobic interactions at the interface between the pseudolysin and SMPI are important contributors to the higher affinity in the VAHD orientation. The position 3 complex did not have such interactions. In the VAHD complex, SMPI had a hydrophobic interface involving the P1 0 residue (Val65), the P3 0 residue (Phe67), and the P3 residue (Val62). Pseudolysin interacted with these residues via the aromatic amino acids Phe129 and Phe67 (with the P3 0 residue), and the hydrophobic amino acids Val137 and Leu197 (with the P1 0 and P3 residues of SMPI). Similar hydrophobic interactions of Val137 and Leu197 were also seen in the X-ray structure of pseudolysin in complex with N-(1-carboxy-3-phenylpropyl)phenylalanyl-alpha-asparagine (1u4g). Further, aromatic interactions were seen between His223 (pseudolysin) and Tyr70 (P6 0 residue of SMPI). The imidazole ring of

Table 1 Energy terms (kcal/mol) of pseudolysin, SMPI and the pseudolysin–SMPI complexes Energy term

Position 3 complex

HAHD complex

VAHD complex

Pseudolysin

SMPI

Eelec Evdw Eint EMM GSA GPB GPBSA GPB þ Eelec EMM þ GPBSA TStrans TSrot TSvib TS total


12960.38
1735.60 6379.30
8316.68 89.91
5126.86
5036.96
18087.25
13353.64


12935.61
1103.10 7615.52
6423.19 92.09
5052.83
4960.74
17988.44
11383.93


12724.18
1725.41 6394.75
8054.84 91.70
5448.90
5357.20
18173.08
13412.04


10799.03
1385.05 4766.49
7417.58 63.92
3376.45
3312.53
14175.48
10730.12


1974.35
302.67 1588.76
688.26 31.82
1980.64
1948.82
3954.99
2637.08 16.05 16.01 1037.92 1069.98

(4.2) (1.7) (2.8) (4.4) (0.0) (5.0) (5.0) (2.6) (3.4)

(3.6) (1.8) (3.3) (4.0) (0.0) (4.7) (4.7) (3.3) (3.9)

(4.1) (1.5) (2.9) (4.7) (0.0) (6.9) (6.9) (4.8) (5.1)

(3.5) (1.4) (2.4) (4.3) (0.0) (4.4) (4.4) (2.6) (2.9)

(1.9) (0.7) (1.4) (2.4) (0.0) (1.8) (1.8) (0.8) (1.4) (0.0) (0.0) (0.0) (0.0)

A description of the energy terms is given in Section 2. Standard error of the mean is given in parentheses. The values are calculated from the 350 snapshots extracted at regular time-intervals during 2.3–3.0 ns of MD.

136

O.A. Adekoya et al. / Journal of Structural Biology 153 (2006) 129–144

Table 2 Free energy terms (kcal/mol) for the pseudolysin–SMPI complexes Energy term

DEelec DEvdw DEint DEMM DGSA DGPB DGPBSA DðGPB þ Eelec Þ DðEMM þ GPBSA Þ
TStotal DGcalc

Energies of pseudolysin–SMPI interactions (kcal/mol)a HAHD complex

Position 3 complex

VAHD complex


162.24 584.62 1260.27 1682.66
3.65 304.27 300.62 142.03 1983.27
1069.98 914.29


187.01
47.88 24.05
210.84
5.84 230.23 224.40 43.22 13.56
1069.98
1056.42

49.19 (5.0)
37.69 (2.2) 39.5 (3.9) 51.01 (6.5)
4.04 (0.0)
91.81 (5.7)
95.85 (5.7)
42.61 (4.1)
44.84 (5.2)
1069.98 (0.0)
1114.82

(5.2) (2.3) (4.2) (5.9) (0.0) (5.4) (5.4) (3.8) (4.9) (0.0)

(4.9) (2.3) (3.9) (5.6) (0.0) (4.9) (4.9) (2.9) (4.0) (0.0)

A description of the energy terms is given in Section 2. The experimentally measured inhibition constant (Ki) for SMPI is: 2.54 · 10
12 M (Oda et al., 1979; Seeram et al., 1997a). a The energies have been calculated as: Energy of complex
energy of pseudolysin
energy of SMPI.

His223 was also close to the disulphide bridge between Cys64 and Cys69 in SMPI. The disulphide bridge has been found to be important for the interactions of SMPI with thermolysin (Hiraga et al., 1999). Table 3 shows the hydrogen bonding interactions obtained by analyzing the trajectories between 2.3 and 3.0 ns of MD simulations (2800 coordinate sets for each pseudolysin–SMPI complex). The pseudolysin–SMPI hydrogen bonds of the positions 3 complex had the highest occupancy during MD, especially those involving Arg66 of SMPI and Asp136 and Asp189 of pseudolysin. These hydrogen bonds had occupancies ranging from 91 to 99.4% (Table 3). Asp136 and Asp189 have previously not been suggested to be involved in the binding interface, and there are no available experimental results verifying their importance. Computational alanine scanning of the positions 3 complex also confirmed that they were important for the free energy of SMPI binding in this position (Table 4). However, the VAHD complex seemed more realistic due to the existence of both polar and non-polar interactions at its interfaces and the stronger SMPI–pseudolysin interactions (Table 2). In spite of that, Asp136 and Asp189 should be the subject of further experimental studies as a validation of the suggested complexes. Relatively fewer hydrogen bond interactions were observed in the HAHD complex (Table 3). The X-ray structure of pseudolysin in complex with N-(1-carboxy-3-phenylpropyl)-phenylalanyl-alpha-asparagine (1u4g) indicates that Glu164, His223, Asn112, Ala113, and Arg198, are involved in hydrogen bonds with the ligand, while Leu197, Val137, Phe129, and Ile186 are involved in hydrophobic interactions. Interestingly, the computational alanine scanning suggested that the residues involved in hydrogen bonds (except for Ala113) all are predicted to be important for the interaction of pseudolysin with SMPI (Tables 3 and 4). It has been proposed that Tyr155 and Asp221 are important for substrate binding (Morihara, 1995). However, in the X-ray structure of pseudolysin in complex with N-(1-carboxy-3-phenylpro-

pyl)-phenylalanyl-alpha-asparagine, (1u4g), these amino acids (Tyr155 and Asp221) are not directly involved in ligand binding, but Asp221 is connected to the side chain of His223 by a hydrogen bond. This hydrogen bond was seen in all the average pseudolysin–SMPI complexes after MD. However, in the VAHD complex, Asp221 interacted with the P6 0 residue of SMPI (Table 3). 3.3. Computational alanine scanning Protein–protein interactions are essential to many processes within living cells. The present pseudolysin–SMPI complexes provide a system to understand the active site of pseudolysin and are important for the identification of the hot spot binding residues at the protein–protein interface that could be mimicked by inhibitor design. Systematic replacement of amino acid residues by alanine at the pseudolysin–SMPI interfaces may reveal hot spot residues (Kortemme and Baker, 2002; Kortemme et al., 2004). There are altogether almost 50 residues at the pseudolysin–SMPI interfaces contributing to the complex building. For each of the pseudolysin–SMPI complexes there are 23–30 at the pseudolysin side, and about 20–22 at the SMPI side. Residues with a DDG greater than 1 kcal/mol were considered hot spot residues for the complex formation. The alanine scanning indicates that only 7–12 residues (on the pseudolysin side) and 4–7 residues on SMPI side are the main contributors to free energy of binding (Tables 4 and 5). Most of the amino acids residues having hydrogen bonding interactions between 2.3 and 3 ns were also suggested to be important for the free energy of binding by the alanine scanning (Tables 3 and 4). We have previously suggested that residues outside the reactive loop of SMPI interacts with thermolysin (Adekoya et al., 2005). In the present study, residues outside the rigid loop were also found important for binding to pseudolysin (Fig. 3). Thr11, Thr22, and Lys31 in the VAHD complex, Ser13 and Arg44 in the HAHD and Lys31 in the position 3 complex were all identified as hot spot residues for binding

O.A. Adekoya et al. / Journal of Structural Biology 153 (2006) 129–144

137

Table 3 Observed hydrogen bonds between SMPI and pseudolysin during 2.3–3 ns of MD simulation SMPI

Position 3 complex

HAHD complex

Residue

Atom

Residue

Atom

Dist

Occu

Asp28 Asp28

OD1 OD2

Arg108 Arg108

HH11 HH11

2.8 2.8

26.4 22.6

Glu30 Glu30 Glu30 Glu30

O OE1 OE1 OE2

Thr46 Thr46 Thr46

H HG1 HG1

2.7 2.7 2.7

21.6 32.1 39.7

Lys31 Lys31 Lys31

HZ2 HZ3 HZ1

Val110 Tyr106 Tyr106

O O O

2.8 2.8 2.8

29.5 34.3 28.1

Asp35 Asp35

OD1 OD1

Asn112 Tyr114

HD22 HH

2.8 2.7

94.1 98.5

Gly36

O

Asn112

HD21

2.8

82.4

Val37

O

Residue

Atom

VAHD complex Dist

Occu

Tyr114

HH

2.6

98.6

Ser44

O

2.8

81.8

Residue

Atom

Dist

Occu

Ser109

HG

2.6

82.1

Tyr114

HH

2.7

51.1

Val110 Val110

O O

2.8 2.8

30.5 29.5

Arg44

HH11

Arg60 (P5) Arg60 Arg60

O HH12 HH22

Val62 (P3)

H

His224

NE2

2.9

29.2

Thr63 (P2)

HG1

Asp206

OD2

2.7

59.9

Tyr114 Ala113

O O

2.8 2.8

79.3 82.2

Asp221

OD1

2.6

56.9

0

Arg66 (P2 ) Arg66 Arg66 Arg66

HH11 HH12 HH22 HH21

His223 Gly219 Gly219

HD1 O O

2.9 2.8 2.8

32.5 51.3 71.8

Asp189 Asp136

OD2 OD2

2.8 2.7

93.9 99.4

Asp189

OD2

2.8

91.9

0

O

Trp115

H

2.9

73.0

0

Cys69 (P5 )

O

Asn112

HD22

2.9

80.0

Tyr70 (P6 0 )

HH

Phe67 (P3 )

0

Gln71 (P7 ) Gln 71

HE21 OE1

Asn112

OD1

2.9

37.6 Asn112

HD21

2.9

45.1

˚ ). Occu: The The number of sampled coordinate sets was 2800 for each of the psedulysin–SMPI complex. Dist: The average hydrogen bonding distances (A ˚ . Only the hydrogen bonds with occupancy percent of coordinate sets between 2.3 and 3.0 ns with the specific hydrogen bonding distance shorter than 3 A above 20% are included in the table. The binding site nomenclature for SMPI is similar to that used by Schechter and Berger (1967).

(Table 5). Although the VAHD complex seems most realistic, the SMPI residues outside the rigid loop predicted as important for pseudolysin binding in all three complexes should be subjected to mutational analysis. Such studies would verify their role in pseudolysin inhibition, and the quality of the present predictions. Within the reactive loop of SMPI, Thr63, Val65, Arg66, Phe67 (the VAHD complex), Arg60, Thr63, Arg66, Phe67, Tyr70 (the position 3 complex) and Phe67, Gln71, and Thr95 (the HAHD complex) were predicted as hot spot residues. No mutational studies of the SMPI–pseudolysin complex have been published, but mutational analysis of the SMPI–thermolysin complex showed that Arg60, Arg61, and Arg66 within the rigid active site loop are important for thermolysin inhibition. The double mutant, R60/61A, had 5 times weaker inhibition than the wild type, while the single mutant R66A showed 100 times reduced inhibition, whereas triple mutant R60/61/66A showed 200 reduced inhibition compared to the wild type (Hiraga et al., 1999).

At the pseudolysin side of the protein–protein interface, Asn112, Tyr114, His140, Leu197, Asp206, His224 (the VAHD complex), Asn112, Tyr114, Ser134, Asp136, His140, Glu164, Asp189, Arg198, His 223, His224 (the position 3 complex) and Asn112, Tyr114, Arg198 (the HAHD complex) were predicted as important for SMPI binding (Table 4). These results are in agreement with the inhibitor binding site seen in the X-ray structure of pseudolysin in complex with N-(1-carboxy-3-phenylpropyl)phenylalanyl-alpha-asparagine, (1u4g). 3.4. The SMPI interactions with pseudolysin and thermolysin Experimental binding studies have indicated that pseudolysin binds the inhibitors 10 to 100-fold stronger than thermolysin (Seeram et al., 1997a). Thermolysin is irreversibly inhibited by ClCH2CO–HOLeu–OCH3 (HOLeu: N-hydroleucine), while pseudolysin is not, whereas

138

O.A. Adekoya et al. / Journal of Structural Biology 153 (2006) 129–144

Table 4 Hot spot and neutral residues of pseudolysin at the protein–protein interactions surface predicted by computational alanine scanning of the average pseudolysin–SMPI complexes between 2.3 and 3.0 ns of MD simulation Pseudolysin mutant S44A S45A T46A D47A N67A Y106A R108A S109A V110A E111A N112A Y114A W115A L121A T127A M128A F129A L132A S134A D136A V137A H140A E141A Y155A E164A I186A D189A I190A L197A R198A D206A R208A Y217A D221A H223A H224A

VAHD complex DDGbind

Position 3 complex DDGbind 0.00 1.16

0.05 0.94
0.07 0.51 1.47
0.17 0.80 2.68 0.49 0.00 0.84 0.11

0.62 1.35


0.04

1.02 0.46 2.04 0.57 0.11
0.31 0.87 2.57

0.15
0.10 0.89
0.08 2.43 1.87 0.32 0.02 0.06 0.51 0.98 0.45 1.62 1.14 0.90 2.72 0.13 6.98 0.16 3.27 0.33 1.81 3.91

HAHD complex DDGbind
0.04
0.03 1.13 0.00

0.09
0.17 0.54 0.01 3.37 2.58 0.22 0.05 0.03 0.23

0.24
4.64 1.64

0.70 3.58

0.07
0.16 1.87 2.19

0.53

DDGbind: Predicted changes in binding energy upon alanine mutation (kcal/mol). A positive value means that the replacement by alanine is predicted to destabilize the complex. A negative value predicts a stabilizing effect of the alanine mutation. Residues for which the alanine mutation has a DDGbind of more than 1 kcal/mol are considered hot spot residues. Residues with DDGbind < 1 kcal/mol are considered as neutral amino acid residues. Residues in bold type are predicted hot spot residues in at least two of the average pseudolysin–SMPI complexes.

ClCH2CO–HOLeu–Ala–Gly–NH2 can irreversibly inhibit both enzymes (Rasnick and Powers, 1978). This may indicate that the binding pocket of pseudolysin is bigger than the binding pocket of thermolysin, since an increase in the length of the ligand by two amino acids (Ala and Gly) results in inhibition of both pseudolysin and thermolysin. Most of the amino acids involved in active site interactions and ligand binding are highly conserved between thermolysin and pseudolysin. Previous studies of thermolysin have indicated the following amino acids as important for substrate/inhibitor binding or catalyses (Matthews,

Table 5 Hot spot and neutral residues of SMPI at the protein–protein interaction surface predicted by computational alanine scanning of the average pseudolysin–SMPI complexes between 2.3 and 3.0 ns of MD simulation SMPI mutant

VAHD complex DDGbind

T 11A S13A T22A I24A D28A M29A E30A K31A T34A D35A V37A K38A S42A R44A R60A K61A V62A T63A C64A V65A R66A F67V C69A Y70A Q71A T93A T95A

1.16
0.08 1.42 0.66 0.56 0.08
0.02 3.20
0.16 0.05

0.10 0.29 0.03
0.02 1.59
0.04 1.65 1.88 1.00
0.08 0.31

Position 3 complex DDGbind

HAHD complex DDGbind

0.05 1.13 0.04 0.52 0.02 0.00 0.80 2.45
0.03 0.36 0.23
0.02 0.20 2.10 0.07 1.25
0.09 0.58 4.29 3.84
0.04 1.24 0.50

0.13 0.36
0.05 2.44 0.00
0.14 0.58 0.18
0.08 0.11
4.40 1.49
0.05 5.51 0.96 3.44

Residues in bold type are predicted as hot spot residues in at least two of the binding positions. The meaning of DDGbind (kcal/mol) is explained in Table 4.

1988) (the corresponding amino acids in pseudolysin are given in brackets): Asn112 (Asn112), Ala113 (Ala113), Trp115 (Trp115), Arg203 (Arg198), Tyr157 (Tyr155), Phe130 (Phe129), Leu133 (Leu132), Val139 (Val137), Glu143 (Glu141), Ile188 (Ile186), Val192 (Ile190), Leu202 (Leu197) and His231 (His223) (Fig. 4). The present study indicates that the corresponding amino acids in pseudolysin also are important for binding of both SMPI (Tables 3–5) and the other small molecule inhibitors (Table 6). The experimentally determined inhibition constants (Ki) of SMPI is 2.54 · 10
12 M for pseudolysin, and 1.14 · 10
10 M for thermolysin. In a previous study, we used the MM-PBSA method to study the interactions of SMPI with thermolysin, and an orientation of SMPI as in the VAHD complex also gave the most favorable free energy of binding to thermolysin. The free energy of SMPI binding was
984.49 kcal/mol for thermolysin (Adekoya et al., 2005), and
1114.82 kcal/mol for pseudolysin (Table 2). These results confirm the experimental observation that SMPI binds stronger to pseudolysin than to thermolysin. The calculations indicated that some amino acids not conserved between pseudolysin and thermolysin are important for SMPI binding to pseudolysin. These amino acids are Tyr114, Asp206 and His224, (Tables 3 and 4, Fig. 4).

O.A. Adekoya et al. / Journal of Structural Biology 153 (2006) 129–144

139

Fig. 3. The average pseudolysin–SMPI complexes between 2.3 and 3.0 ns of MD simulation. A closer view of the auxiliary binding interactions outside the rigid loop of SMPI is shown to the right. (A) The VAHD complex. (B) The HAHD complex. (C) Position 3 complex. Color coding: blue cpk; zinc ion; purple cpk; calcium ion; red; pseudolysin backbone, green; SMPI backbone.

The corresponding amino acids of thermolysin are Phe114, Tyr206, and Ile224, respectively. These differences partly explain the stronger affinities of SMPI to pseudolysin than to thermolysin. For both thermolysin and pseudolysin, the main contribution to the free energy of SMPI binding is a combination of the electrostatic contribution from the molecular mechanical energy term and the electrostatic contribution to the solvation free energy calculated by the DelPhi program. The VAHD complexes also showed that pseudolysin had more contact sites with SMPI outside the Cys64– Val65 segment of the rigid loop than had thermolysin. However, the corresponding region in both enzymes interacted outside the Cys64–Val65 segment of SMPI. This region consisted of Arg108, Ser109, Val110, Glu111, Asn112, Ala113, and Tyr114 in pseudolysin and Gln108, Gly109, Tyr110, Asn111, Asn112, Ala113, and Phe114 in thermolysin (Fig. 4). Within these regions, the contacts sites outside the Cys64–Val65 segment of SMPI were Gln108 and Tyr110 for thermolysin and Ser109, Val110, and Tyr114 for pseudolysin. The occupancy of the observed hydrogen bonds during MD was also higher for the pseudolysin–SMPI complex (Table 3) than for the thermolysin–SMPI (Adekoya et al., 2005).

3.5. Pseudolysin–small molecule inhibitor interactions The LIE method is a good alternative to the more time consuming MM-PBSA analysis for ranking the binding constants of small molecules. The method has been applied to several different systems of protein–small molecule complexes and the results are reported in good agreement with experimental data (Ljungberg et al., 2001; Wang et al., 1999). To obtain reliable thermally averaged energies for the LIE calculations, the MD simulation must be sufficiently long enough. In the present study, the vdw and electrostatic binding contributions between 4.5 and 5.0 ns of MD were used for the LIE calculations. The estimated rank order of binding affinity was: Tlp > POH > ZOH > HRO with the calculated binding free energy in the range of
67.8 to
58.0 kcal/mol (Table 7). The ranking according to the experimental inhibition constants was: Tlp > ZOH > POH > HRO (Table 7). The number, and the contact surface area of stabilizing hydrophobic, aromatic and hydrogen bonding interactions appears to reflect the affinity differences between the inhibitors (Table 6). The average complexes obtained between 4.0 and 5.0 ns of MD (4000 coordinate sets of each) indicated that Tlp had the highest number of hydrogen bond interactions as well as stabilizing hydrophobic interactions (Table 6 and Fig. 1).

140

O.A. Adekoya et al. / Journal of Structural Biology 153 (2006) 129–144

Fig. 4. Alignments of the sequences of thermolysin (bathe) and pseudolysin (pseu).

These hydrophobic and hydrogen bonding interactions by Tlp appear to be the most distinguishing factors for the stronger pseudolysin affinity than the other inhibitors. The four inhibitors appeared to have about the same number of destabilizing hydrophobic–hydrophilic interaction, except for POH having slightly less of such interaction than the others. POH also had the highest number of aromatic interactions, followed by Tlp. Such interactions were not observed for HRO and ZOH. The experimental binding affinity of ZOH is higher than that of POH, which is contradictory to the predicted affinities (Table 7) The higher number of aromatic interaction and less hydrophobic– hydrophilic (destabilizing) interactions of POH with pseudolysin might explain why POH was predicted to bind stronger to pseudolysin than ZOH. ZOH had more hydrophobic interactions with pseudolysin than HRO. The reason for that is that HRO lacks the hydrophobic side chain of ZOH. Some of the amino acids in pseudolysin contributed both to stabilizing and destabilizing interactions with the inhibitors, while others were detected as stabilizing or destabilizing (Table 6). The number of amino acids contacts found to be entirely stabilizing without any destabilizing (hydrophobic–hydrophilic) contribution was much higher for Tlp and POH than for the other ligands. The residues involved in stabilizing hydrophobic interactions to HRO and ZOH were also involved in destabilizing interactions (Table 6). For Tlp, the interactions with Leu132 and Leu197 were all defined as entirely stabilizing,

while Leu132 and Arg198 formed corresponding stabilizing interactions with POH. The pseudolysin–inhibitor interactions (Table 6 and Fig. 1) were in agreement with the X-ray structure of the pseudolysin–inhibitor complex (1u4g). Also in this complex, aromatic–aromatic interactions were observed between the inhibitor and Tyr114 and His223, while aromatic–hydrophobic interactions were seen between the inhibitor and Phe129, Ile186, Leu197, Leu137, and Glu141. These amino acids were involved in similar interactions in the averaged pseudolysin-Tlp complex after the MD simulation (Fig. 1). Further, in the X-ray complex, Arg198, Glu164, His223, and Asn112 are involved in hydrogen bonding interactions with the inhibitor. In the present pseudolysin–Tlp complex Arg198, His223, Glu141, and Trp115 formed hydrogen bonds with Tlp. Glu164 was also involved in hydrogen bonding interactions with POH (Table 6). The amino acids found to be important for binding of the small molecular inhibitors were also the main amino acids for SMPI binding (Table 4). 3.6. Towards rational design of pseudolysin inhibitors As for thermolysin, the VAHD complex had the most favourable free energy of association. The identification of hot spot residues for the SMPI–pseudolysin contacts in the VAHD complex could be used as a guide for a structure based design of pseudolysin inhibitors. Table 3 shows

O.A. Adekoya et al. / Journal of Structural Biology 153 (2006) 129–144

141

Table 6 The most important residual contacts with the inhibitors in average structures between 4.0 and 5.0 ns of MD Interactions

Ligand HRO Pseudolysin residue

Zn Aromatic–aromatic

ZOH Dist.

Surf.

1.9

33.8

POH

Pseudolysin residue

Dist.

Surf

1.6

41.1

Tlp

Pseudolysin residue

Dist.

Surf.

Phe129 His140 His224

1.6 3.7 3.8 3.7

38.0 31.2 45.9 22.4

Pseudolysin residue

Dist.

Surf.

Phe129

2.9 4.1

25.1 18.4

Hydrophobic–hydrophobic

Ala113 Leu137 Leu197

4.2 3.3 3.5

30.6 40.8 55.1

Val137 Ile186 Leu197 Arg198

3.8 3.3 3.8 3.7

39.5 30.3 65.5 59.4

Leu132 Leu197 Arg198

3.5 3.5 3.3

20.9 56.1 39.4

a Ala113 Tyr114 b Leu132 a Val137 a His140 a Glu141 b His223 a Leu197 b Leu197

4.0 3.7 4.0 3.2 3.6 4.3 3.5 4.4 3.6

18.8 26.0 14.6 42.5 32.7 11.9 32.4 6.1 30.1

Hydrophobic–hydrophilic (destabilizing)

Glu141 Ala113 Arg198 Glu164 Leu137 Leu197 Asp206

2.1 4.2 2.9 2.6 3.3 3.5 4.4

44.1 30.6 60.1 5.2 40.8 55.1 13.1

Val137 Glu141 Glu164 Asp168 Ile186 Leu197 Arg198

3.8 3.0 2.5 4.1 3.3 3.8 3.7

39.5 28.6 9.9 4.3 30.3 65.5 59.4

Val137 His140 Glu164 Ile190 Leu197

3.8 3.8 2.6 3.6 3.5

20.9 45.9 19.2 17.6 56.1

Tyr114 Trp115 a Ala113 a Val137 a His140 a Glu141 b His223

3.7 3.2 4.0 3.2 3.6 4.3 3.5

26.0 30.6 18.8 42.5 32.7 11.9 32.4

Hydrogen bonds

Glu141

2.1

44.1

Glu141

3.0

28.6

Glu164

2.6

19.2

Trp115 His223 a Arg198

3.2 3.1 2.7

30.6 21.0 39.4

˚ ). Surf.: Surface Ligand–protein contacts were analysed with the ligand–protein contacts (LPC) software (Sobolev et al., 1999). Dist.: Atomic distance (A ˚ 2). contact area (A a b

Indicates residue interacting with Leu of Tlp. Indicates residue interacting with Trp of Tlp.

Table 7 Average vdw and electrostatic interactions (kcal/mol) from the bound and unbound (free) MD simulations, and the calculated free energy of binding (DGb) Ligand

ÆVvdwæbound

ÆVvdwæfree

ÆVelæbound

ÆVelæfree

DGb

Ki (lM)

HRO ZOH POH Tlp


16.8
19.9
21.1
37.8


0.02
0.16
0.08
0.11


110.9
122.6
126.1
128.9


0.3 (1.3)
0.17 (0.7)
0.2 (1.2)
4.5 (2.3)


58.0
64.4
66.3
67.8

340 11 21 0.2

(4.4) (3.7) (4.1) (3.3)

(0.3) (0.3) (0.3) (0.3)

(29.5) (8.1) (7.1) (12.1)

The coordinate sets sampled between 4.5 and 5.0 ns of the MD simulations were used for the LIE calculations. Standard deviations are given in parentheses. The experimentally observed inhibition constants, (Ki) are given (Nishino and Powers, 1980; Kessler et al., 1982).

˚ ) between the residues of pseudolysin and the distances (A SMPI. Identification of the molecular interaction fields (hydrophobic, hydrogen bond donor and acceptor, etc) of the residues in the binding pocket and the effects of flexibility of these residues on their interaction fields upon inhibitor binding are also important in the de novo ligand design. Table 5 indicates that the hot spot residues in the reactive loop of SMPI for pseudolysin binding (VAHD complex) were Thr63 (P2 residue), Val65 (P1 0 residue), Arg66 (P2 0 residue) and Phe67 (P3 0 residue). All these residues interact very close to the active site residues of pseudolysin. Their interactions with active site residues could be re-constructed by de novo ligand design.

3.7. Treatment of catalytic zinc Zinc in the binding site of metalloproteinases performs essential biological functions and contributes considerably to the binding affinity of both small molecule inhibitors and protein inhibitors. Docking and MD of zinc-containing enzymes is a challenge due to the multiple coordination geometries of zinc and the lack of appropriate force field parameters. Zinc is most often found to be four-coordinated with a tetrahedral geometry. However, five and six coordinated zinc with an octahedral, trigonal or square-base geometry are also seen in zinc metalloproteinases (Alberts et al., 1998). The multiple possibilities

142

O.A. Adekoya et al. / Journal of Structural Biology 153 (2006) 129–144

for coordination geometry indicate that parameters developed for one system are not directly transferable to another system. It is therefore notoriously difficult to evaluate metal–ligand interactions (Ferrara et al., 2004). A study assessing the scoring function for protein–ligand interactions indicates that neither a more accurate treatment of solvation nor a more sophisticated charge model for zinc improved the quality of the result (Ferrara et al., 2004). The correlation values for the metalloproteinases were very low for all scoring functions using +2 charge for zinc. A charge transfer between the zinc ion and its ligand atoms has therefore been suggested. In previous free energy prediction studies of zinc metalloproteinases both a bonded (Hou et al., 2001) and non-bonded (Donini and Kollman, 2000; Hou et al., 2002) approach for zinc coordination have given results in agreement with experimental studies. In bonded approaches, the coordination between the zinc ion and ligated groups are described by usual terms including bond stretching, angle bending and torsional contributions. In such an approach the coordination number is well preserved, but the zinc atom is almost frozen during the MD. In a non-bonded approach, non-bonded van der Waals and electrostatic interactions are used for the coordination between zinc and the ligated groups. Such an approach may give a much more flexible molecular system. In the present study, we used a non-bonded approach for the catalytic zinc atom. However, to ensure that the coordination was retained during MD, restraints between the zinc ion and coordinating atoms were used. A similar approach has also been used by others (Wasserman and Hodge, 1996; Tate et al., 1998). Similar zinc parameters were assigned both for four- and five coordinated zinc, which is regarded as a simplification and might contribute to energetic and conformational errors. The energies were overestimated compared with the values calculated from experimentally detected inhibition constants. As discussed above, an incorrect parameterization and geometry of the zinc coordination may be one of the reasons for that. A formal charge of +2 was used for the catalytic zinc atom, which also has been used by others during MD of zinc containing proteins (Donini and Kollman, 2000). This might also contribute to errors, since it is obvious that the atomic charge of zinc will depend on the orbitals of the coordinating groups. In spite of the possible limitations in zinc parameterization, the relative ranking of the ligand free energies was in good correlation with experimental inhibition constants. Acknowledgment This work has received support from The Research Council of Norway (program for Supercomputing) through a grant for computer time.

References Abagyan, R., Maxim, T., Kuznetsov, D., 1994. ICM—A new method for protein modelling and design: applications to docking and structure prediction from the distorted native conformation. J. Comp. Chem. 15, 488–506. Adekoya, O.A., Willassen, N.P., Sylte, I., 2005. The protein–protein interactions between SMPI and thermolysin studied by molecular dynamics and MM/PBSA calculations. J. Biomol. Struct. Dyn. 22, 521–532. Alberts, I.L., Nadassy, K., Wodak, S.J., 1998. Analysis of zinc binding sites in protein crystal structures. Protein Sci. 7, 1700–1716. ˚ qvist, J., 1996. Calculation of absolute binding free energies for charged A ligands and effects of long-range electrostatic interactions. J. Comput. Chem. 17, 1587–1597. ˚ qvist, J., Luzhkov, V.B., Brandsdal, B.O., 2002. Ligand binding affinities A from MD simulations. Acc. Chem. Res. 35, 358–365. ˚ qvist, J., Medina, C., Samuelsson, J.E., 1994. A new method for A predicting binding affinity in computer-aided drug design. Protein Eng. 7, 385–391. Bayly, C.I., Cieplak, P., Cornell, W.D., Kollman, P.A., 1993. A wellbehaved electrostatic potential based method using charge restraints for deriving atomic charges—the RESP model. J. Phys. Chem. 97, 10269–10280. Burns, F.R., Paterson, C.A., Gray, R.D., Wells, J.T., 1990. Inhibition of Pseudomonas aeruginosa elastase and Pseudomonas keratitis using a thiol-based peptide. Antimicrob. Agents Chemother. 34, 2065–2069. Case, D.A., Pearlman, D.A., Caldwell, J.W., Cheatham III, T.E., Wang, J., Ross, W.S., Simmerling, C., Darden, T., Merz, K.M., Stanton, R.V., et al., 2002. AMBER 7. AMBER 7. University of California, San Francisco. Cheatham III, T.E., Srinivasan, J., Case, D.A., Kollman, P.A., 1998. Molecular dynamics and continuum solvent studies of the stability of polyG–polyC and polyA–polyT DNA duplexes in solution. J. Biomol. Struct. Dyn. 16, 265–280. Clackson, T., Ultsch, M.H., Wells, J.A., deVos, A.M., 1998. Structural and functional analysis of the 1:1 growth hormone:receptor complex reveals the molecular basis for receptor affinity. J. Mol. Biol. 277, 1111–1128. Cornell, W.D., Cieplak, P., Bayly, C.I., Gould, I.R., Merz, K.M., Ferguson, D.M., Spellmeyer, D.C., Fox, T., Caldwell, J.W., Kollman, P.A., 1995. A 2nd generation force-field for the simulation of proteins, nucleic-acids, and organic-molecules. J. Am. Chem. Soc. 117, 5179–5197. Cornell, W.D., Cieplak, P., Bayly, C.I., Kollman, P.A., 1993. Application of RESP charges to calculate conformational energies, hydrogen-bond energies, and free-energies of solvation. J. Am. Chem. Soc. 115, 9620– 9631. Darden, T., York, D., Pedersen, L., 1993. Particle mesh Ewald: An N Æ log (N) method for Ewald sums in large systems. J. Chem. Phys. 98, 10089–10092. Doring, G., Obernesser, H.J., Botzenhart, K., 1981. Extracellular toxins of Pseudomonas aeruginosa. II. Effect of two proteases on human immunoglobulins IgG, IgA and secretory IgA (authors transl.). Zentralbl. Bakteriol. A 249, 89–98. Donini, O.A.T., Kollman, P.A., 2000. Calculation and prediction of binding free energies for the matrix metalloproteinases. J. Med. Chem. 43, 4180–4188. Essmann, U., Lalith, Perera, Berkowitz, M.L., Darden, T., Hsing, Lee, Pedersen, L.G., 1995. A smooth particle mesh Ewald method. J. Chem. Phys. 103, 8577–8593. Ferrara, P., Gohlke, H., Price, D.J., Klebe, G., Brooks 3rd, C.L., 2004. Assessing scoring functions for protein–ligand interactions. J. Med. Chem. 47, 3032–3047. Frisch, M.J., Trucks, G.W., Schlegel, H.B., Scuseria, G.E., Robb, M.A., Cheeseman, J.R., Zakrzewski, V.G., Montgomery, J.A.J., Stratmann, R.E., Burant, J.C., et al., 2001. Gaussian 98 (Revision A.11). Gaussian, Pittsburgh, PA.

O.A. Adekoya et al. / Journal of Structural Biology 153 (2006) 129–144 Gilson, M.K., Honig, B., 1988. Calculation of the total electrostatic energy of a macromolecular system: solvation energies, binding energies, and conformational analysis. Proteins 4, 7–18. Gunsteren, W.F.v., Berendsen, H.T.C., 1977. Algorithms for macromolecular dynamics and contraints dynamics. Mol. Phys. 34, 1311. ˚ qvist, J., 1998. Ligand binding affinity Hansson, T., Marelius, J., A prediction by linear interaction energy methods. J. Comput.-Aid. Mol. Des. 12, 27–35. Heck, L.W., Morihara, K., McRae, W.B., Miller, E.J., 1986. Specific cleavage of human type III and IV collagens by Pseudomonas aeruginosa elastase. Infect. Immun. 51, 115–118. Hiraga, K., Seeram, S.S., Tate, S., Tanaka, N., Kainosho, M., Oda, K., 1999. Mutational analysis of the reactive site loop of Streptomyces metalloproteinase inhibitor, SMPI. J. Biochem. (Tokyo) 125, 202–209. Hobden, J., 2002. Pseudomonas aeruginosa proteases and corneal virulence. DNA Cell Biol. 21, 391–396. Holden, H.M., Matthews, B.W., 1988. The binding of L-valyl–L-tryptophan to crystalline thermolysin illustrates the mode of interaction of a product of peptide hydrolysis. J. Biol. Chem. 263, 3256–3260. Holder, I.A., Wheeler, R., 1984. Experimental studies of the pathogenesis of infections owing to Pseudomonas aeruginosa: elastase, an IgG protease. Can. J. Microbiol. 30, 1118–1124. Honig, B., Nicholls, A., 1995. Classical electrostatics in biology and chemistry. Science 268, 1144–1149. Hou, T., Guo, S., Xu, X., 2002. Predictions of binding of a diverse set of ligands to gelatinase-A by a combination of molecular dynamics and continuum solvent models. J. Phys. Chem. B 106, 5527–5535. Hou, T.J., Zhang, W., Xu, X.J., 2001. Binding affinities for a series of selective inhibitors of gelatinase-A using molecular dynamics with a linear interaction energy approach. J. Phys. Chem. B 105, 5304– 5315. Izquierdo-Martin, M., Stein, R.L., 1992. Mechanistic studies on the inhibition of thermolysin by a peptide hydroxamic acid. J. Am. Chem. Soc. 114, 325–331. Johnson, D.A., Carter-Hamm, B., Dralle, W.M., 1982. Inactivation of human bronchial mucosal proteinase inhibitor by Pseudomonas aeruginosa elastase. Am. Rev. Respir. Dis. 126, 1070–1073. Jorgensen, W.L., 1982. Revised TIPS for simulations of liquid water and aqueous solutions. J. Chem. Phys. 77, 4156. Kessler, E., Blumberg, S., 1987. Specific inhibitors of Pseudomonas aeruginosa elastase as potential drugs for the treatment of Pseudomonas keratitis. Antibiot. Chemother. 39, 102–112. Kessler, E., Israel, M., Landshman, N., Chechick, A., Blumberg, S., 1982. In vitro inhibition of Pseudomonas aeruginosa elastase by metalchelating peptide derivatives. Infect. Immun. 38, 716–723. Kessler, E., Spierer, A., 1984. Inhibition by phosphoramidon of Pseudomonas aeruginosa elastase injected intracorneally in rabbit eyes. Curr. Eye Res. 3, 1075–1078. Kollman, P.A., Massova, I., Reyes, C., Kuhn, B., Huo, S., Chong, L., Lee, M., Lee, T., Duan, Y., Wang, W., Donini, O., Cieplak, P., Srinivasan, J., Case, D.A., Cheatham, T.E., 2000. Calculating structures and free energies of complex molecules: combining molecular mechanics and continuum models. Acc. Chem. Res. 33, 889–897. Komori, Y., Nonogaki, T., Nikai, T., 2001. Hemorrhagic activity and muscle damaging effect of Pseudomonas aeruginosa metalloproteinase (elastase). Toxicon 39, 1327–1332. Kortemme, T., Baker, D., 2002. A simple physical model for binding energy hot spots in protein–protein complexes. PNAS 99, 14116–14121. Kortemme, T., Kim, D.E., Baker, D., 2004. Computational alanine scanning of protein–protein interfaces. Sci. STKE 2004, pl2. Kottalam, J., Case, D.A., 1990. Langevin modes of macromolecules— applications to Crambin and DNA hexamers. Biopolymers 29, 1409– 1421. Lamm, G., Szabo, A., 1986. Langevin modes of macromolecules. Biophys. J. 49, A57. Ljungberg, K.B., Marelius, J., Musil, D., Svensson, P., Norden, B., ˚ qvist, J., 2001. Computational modelling of inhibitor binding to A human thrombin. Eur. J. Pharm. Sci. 12, 441–446.

143

˚ qvist, J., 1998. Calculation of ligand binding Marelius, J., Hansson, T., A free energies from molecular dynamics simulations. Int. J. Quant. Chem. 69, 77–88. Mariencheck, W.I., Alcorn, J.F., Palmer, S.M., Wright, J.R., 2003. Pseudomonas aeruginosa elastase degrades surfactant proteins A and D. Am. J. Respir. Cell Mol. Biol. 28, 528–537. Matthews, B.W., 1988. Structural basis of the action of thermolysin and related zinc peptidases. Acc. Chem. Res. 21, 333–340. Morihara, K., 1995. Pseudolysin and other pathogen endopeptidases of thermolysin family. Methods Enzymol. 248, 242–253. Morihara, K., Tsuzuki, H., Oda, K., 1979. Protease and elastase of Pseudomonas aeruginosa: inactivation of human plasma alpha 1proteinase inhibitor. Infect. Immun. 24, 188–193. Nishino, N., Powers, J.C., 1980. Pseudomonas aeruginosa elastase. Development of a new substrate, inhibitors, and an affinity ligand. J. Biol. Chem. 255, 3482–3486. Oda, K., Koyama, T., Murao, S., 1979. Purification and properties of a proteinaceous metallo-proteinase inhibitor from Streptomyces nigrescens TK-23. Biochim. Biophys. Acta 571, 147–156. Ohno, A., Tate, S., Seeram, S.S., Hiraga, K., Swindells, M.B., Oda, K., Kainosho, M., 1998. NMR structure of the Streptomyces metalloproteinase inhibitor, SMPI, isolated from Streptomyces nigrescens TK-23: another example of an ancestral beta gamma-crystallin precursor structure. J. Mol. Biol. 282, 421–433. Pearlman, D.A., Case, D.A., Caldwell, J.W., Ross, W.S., Cheatham III, T.E., DeBolt, S., Ferguson, D.M., Seibel, G.L., Kollman, P., 1995. AMBER, a package of computer programs for applying molecular mechanics, normal mode analysis, molecular dynamics and free energy calculations to simulate the structural and energetic properties of molecules. Comput. Phys. Commun. 91, 1–41. Rasnick, D., Powers, J.C., 1978. Active site directed irreversible inhibition of thermolysin. Biochemistry 17, 4363–4369. Rocchia, W., Sridharan, S., Nicholls, A., Alexov, E., Chiabrera, A., Honig, B., 2002. Rapid grid-based construction of the molecular surface and the use of induced surface charge to calculate reaction field energies: applications to the molecular systems and geometric objects. J. Comp. Chem. 23, 128–137. Sanner, M., Olson, A.J., Spehner, J.C., 1996. Reduced surface: An efficient way to compute molecular surfaces. Biopolymers 38, 305– 320. Schechter, I., Berger, A., 1967. On the size of the active site in proteases. I. Papain. Biochem. Biophys. Res. Commun. 27, 157–162. Schmidtchen, A., Holst, E., Tapper, H., Bjorck, L., 2003. Elastaseproducing Pseudomonas aeruginosa degrade plasma proteins and extracellular products of human skin and fibroblasts, and inhibit fibroblast growth. Microb. Pathogen. 34, 47–55. Seeram, S.S., Hiraga, K., Oda, K., 1997a. Resynthesis of reactive site peptide bond and temporary inhibition of Streptomyces metalloproteinase inhibitor. J. Biochem. (Tokyo) 122, 788–794. Seeram, S.S., Hiraga, K., Saji, A., Tashiro, M., Oda, K., 1997b. Identification of reactive site of a proteinaceous metalloproteinase inhibitor from Streptomyces nigrescens TK-23. J. Biochem. (Tokyo) 121, 1088–1095. Sitkoff, D., Kim, A., Sharp, Honig, B., 1994. Accurate calculation of hydration free energies using macroscopic solvent models. J. Phys. Chem. 98, 1978–1988. Sobolev, V., Sorokine, A., Prilusky, J., Abola, E.E., Edelman, M., 1999. Automated analysis of interatomic contacts in proteins. Bioinformatics 15, 327–332. Tate, S., Ohno, A., Seeram, S.S., Hiraga, K., Oda, K., Kainosho, M., 1998. Elucidation of the mode of interaction of thermolysin with a proteinaceous metalloproteinase inhibitor, SMPI, based on a model complex structure and a structural dynamics analysis. J. Mol. Biol. 282, 435–446. Terp, G.E., Christensen, I.T., Jorgensen, F.S., 2000. Structural differences of matrix metalloproteinases. Homology modeling and energy minimization of enzyme–substrate complexes. J. Biomol. Struct. Dyn. 17, 933–946.

144

O.A. Adekoya et al. / Journal of Structural Biology 153 (2006) 129–144

Thayer, M.M., Flaherty, K.M., McKay, D.B., 1991. Three-dimensional ˚ resolustructure of the elastase of Pseudomonas aeruginosa at 1.5-A tion. J. Biol. Chem. 266, 2864–2871. Tronrud, D.E., Holden, H.M., Matthews, B.W., 1987. Structures of two thermolysin–inhibitor complexes that differ by a single hydrogen bond. Science 235, 571–574. Tronrud, D.E., Monzingo, A.F., Matthews, B.W., 1986. Crystallographic structural analysis of phosphoramidates as inhibitors and transitionstate analogs of thermolysin. Eur. J. Biochem. 157, 261–268.

Wallace, A.C., Laskowski, R.A., Thornton, J.M., 1995. LIGPLOT: a program to generate schematic diagrams of protein–ligand interactions. Protein Eng. 8, 127–134. Wang, W., Wang, J., Kollmann, P.A., 1999. What determine the van der Waals coefficient b in the LIE (Linear Interaction Energy) method to estimate binding free energies using molecular dynamics simulations? Proteins: Struct. Funct. Gen. 34, 395–402. Wasserman, Z.R., Hodge, C.N., 1996. Fitting an inhibitor into the active site of thermolysin: a molecular dynamics case study. Proteins 24, 227–237.