Article pubs.acs.org/JPCB
Hybrid Potential Simulation of the Acylation of Enterococcus faecium L,D‑Transpeptidase by Carbapenems Nicholus Bhattacharjee,† Martin J. Field,*,† Jean-Pierre Simorre,‡ Michel Arthur,§ and Catherine M. Bougault*,‡ †
DYNAMO/DYNAMOP, UMR 5075, Université Grenoble 1, CNRS, CEA, Institut de Biologie Structurale, 71 Avenue des Martyrs, CS 10090, 38044 Grenoble Cedex 9, France ‡ RMN, UMR 5075, Université Grenoble 1, CNRS, CEA, Institut de Biologie Structurale, 71 Avenue des Martyrs, CS 10090, 38044 Grenoble Cedex 9, France § Centre de Recherche des Cordeliers, Equipe 12, UMR S 872, Université Pierre et Marie Curie-Paris 6, INSERM, Université Paris Descartes, Sorbonne Paris Cité, 15 rue de l’Ecole de Médecine, 75006 Paris, France S Supporting Information *
ABSTRACT: The L,D-transpeptidases, Ldts, catalyze peptidoglycan cross-linking in β-lactam-resistant mutant strains of several bacteria, including Enterococcus faecium and Mycobacterium tuberculosis. Although unrelated to the essential D,D-transpeptidases, which are inactivated by the β-lactam antibiotics, they are nevertheless inhibited by the carbapenem antibiotics, making them potentially useful targets in the treatment of some important diseases. In this work, we have investigated the acylation mechanism of the Ldt from E. faecium by the carbapenem, ertapenem, using computational techniques. We have employed molecular dynamics simulations in conjunction with QC/MM hybrid potential calculations to map out possible reaction paths. We have focused on determining the following: (i) the protonation state of the nucleophilic cysteine of the enzyme when it attacks; (ii) whether nucleophilic attack and β-lactam ring-opening are concerted or stepwise, the latter occurring via an oxyanion intermediate; and (iii) the identities of the proton acceptors at the beginning and end of the reaction. Overall, we note that there is considerable plasticity in the mechanisms, owing to the significant flexibility of the enzyme, but find that the preferred pathways are ones in which nucleophilic attack of cysteine thiolate is concerted with β-lactam ring-opening.
■
INTRODUCTION β-Lactams were discovered in the 1920s and first used as therapeutic agents in the 1940s. Since then, the β-lactam family of drugs has played a major role in antibiotherapies.1 These molecules target ubiquitous enzymes, the penicillin-binding proteins (PBPs), which are essential in the extracytoplasmic assembly of disaccharide-pentapeptide precursors into peptidoglycan. The latter is a giant biopolymer that surrounds the bacterial cytoplasmic membrane and provides the bacterial cell with the mechanical resistance required for sustaining the intracellular osmotic pressure. By preventing the transpeptidation reaction between peptide stems of adjacent peptidoglycan motifs,2 β-lactams block bacterial growth and eventually lead to bacterial cell death and lysis. The in vivo efficacy of these drugs is nevertheless limited by β-lactamases in a large variety of bacterial species, which are associated with an increasing number of resistance pathways.3,4 There has been renewed interest in β-lactams in recent years with the demonstration of the efficacy of Meropenem, a member of the carbapenem subfamily, in the inhibition of extensively drug-resistant strains of Mycobacterium tuberculosis.5 The bactericidal in vitro and in vivo activity of Meropenem © XXXX American Chemical Society
against both exponentially growing and nonreplicating bacilli proved to be further increased when the drug was administered in combination with clavulanic acid, a recognized inhibitor of the mycobacterial β-lactamase BlaC. Nevertheless, it was apparent that PBPs were not the sole target for carbapenems in these bacilli, as only 20% of the mycobacterial peptidoglycan showed the classical PBP-catalyzed stem peptide cross-linking pattern.6,7 The additional, new inhibitory target for carbapenems were the L,D-transpeptidases, Ldts. These enzymes, which were initially isolated and characterized from an in vitro selected ampicillin-resistant mutant of Enterococcus faecium, Ldtfm,8,9 form 3 → 3 cross-links in the peptidoglycan, between amino acids at the third position of two different stem peptides, instead of the commonly found 4 → 3 cross-links. Five Ldt paralogues have been identified in M. tuberculosis, LdtMt1 to LdtMt5. Four of these, with the exception of LdtMt3, are active in in vitro peptidoglycan cross-linking assays, whereas all but LdtMt5 are inhibited by carbapenems.10 Furthermore, LdtMt2 is Received: March 18, 2016 Revised: May 10, 2016
A
DOI: 10.1021/acs.jpcb.6b02836 J. Phys. Chem. B XXXX, XXX, XXX−XXX
Article
The Journal of Physical Chemistry B essential for virulence in an acute infection mouse model.11 Its efficient inhibition by carbapenems may thus suffice in the case of mycobacterial infection. In addition, LdtMt1 may play a role in the transition to the persistent state of the bacilli,6 and its inhibition at early stages of the infection may prove crucial, while bacilli lacking a functional copy of LdtMt5 show aberrant growth and enhanced susceptibility to crystal violet and osmotic shock.12 If the role of the other Ldts remains more elusive,10 their efficient inhibition may also be of importance. The in vitro inhibition of L,D-transpeptidases from E. faecium13−15 and M. tuberculosis11,16,17 by carbapenems or penems and, to a lesser extent, cephalosporins, has been demonstrated by mass spectrometry. The reaction proceeds through nucleophilic attack of the antibiotic’s carbonyl group by the enzyme’s catalytic cysteine, yielding a stable covalent acylation product.13,15 The resulting acyl-enzyme is stable for hours, and hydrolysis, leading to the initial enzyme and the opened β-lactam ring, only occurs with a slow rate constant (typically less than 10−3 min−1 for carbapenems). Using a combination of fluorescence and UV−vis spectroscopy approaches, kinetic studies have indicated a two-step acylation mechanism, with the formation of an intermediate with decreased intrinsic tryptophane fluorescence properties in LdtMt116 and Ldtfm.14,18 These studies also allowed the determination of the individual kinetic constants for a series of carbapenems and cephems, identifying antibiotic features essential for efficient inactivation of the enzyme and stability of the resulting acyl-enzyme. A more detailed understanding of the inhibition of Ldts by carbapenem may nevertheless contribute to the design of new drugs that are active against M. tuberculosis, and anticipate the emergence of carbapenemresistant strains of this bacterium. In this respect, simulation has proved a useful complement to experiment for the elucidation of enzymatic reaction mechanisms. Much work has been done in particular on enzymes that undergo acylation reactions with β-lactams and/or peptidoglycan fragments, such as the β-lactamases19−21 and penicillinbinding proteins,22,23 respectively. These enzymes are nevertheless structurally unrelated to Ldts, and contain a catalytic serine in the active site instead of a cysteine residue. Alternatively, the acylation of cysteine proteases, such as papain24 and cathepsin K,25 has been extensively investigated experimentally and computationally during the hydrolysis of their peptide substrates. These enzymes can be considered as a possible model for Ldts, despite the drastically different organization of their catalytic site and overall architecture. Apart from the previously described kinetic results emphasizing the formation of a carpabenem/Ldt intermediate, experimental data on the acylation of Ldts by β-lactams are scarce. Information is limited to X-ray or NMR structures of Ldts in two alternative states: (i) the apoenzyme in the absence of β-lactam antibiotic (PDB codes 4JMN, 3VYO/4GSQ/ 4QRA, 4Z7A, and 3ZG4) or (ii) the acyl-enzyme covalently bound to a carbapenem (PDB codes 4JMX, 3VYP/4GSU/ (4QR7, 4QTF), and 3ZGP) for LdtMt1,26 LdtMt2,27,28 LdtMt5,12 and Ldtfm,29 respectively. To the best of our knowledge, computational studies so far only include a comprehensive series of studies by Silva and co-workers that have targeted M. tuberculosis LdtMt2. These authors have employed a number of approaches, including docking, free-energy calculations, and quantum chemical/molecular mechanical (QC/MM) hybrid potentials, to investigate the formation of 3 → 3 cross-links between two peptidoglycan motifs, 30 and the protein
inactivation mechanism with two carbapenems, imipenem and Meropenem.31,32 They found a reaction scheme that consisted of three steps, namely: (i) the energetically favored generation of a cysteine thiolate/histidine imidazolium pair;30 (ii) nucleophilic attack by the thiolate on the carbapenem, together with a concerted ring-opening of the β-lactam ring; and (iii) protonation of the lactam amide nitrogen by the same histidine that accepted the proton in the first step.32 Although the energy barriers that these workers found for the last two steps are consistent with experimental data on LdtMt1,16 the first step contradicts the relative pKa’s of the catalytic cysteine (around 4) and histidine (around 11) that were experimentally and computationally determined for Ldtfm29 and LdtMt2,33 respectively. Furthermore, in their study, the antibiotic in the acylenzyme product lies in the substrate acyl-acceptor cavity (pocket 2) instead of the substrate acyl-donor cavity (pocket 1) that was determined in the X-ray structures of the LdtMt2 (PDB codes 3VYP,28 4GSU,27 4QR7 for Meropenem, and 4QTF for imipenem). See ref 34 for a description of the role of the different cavities, and Figure S3 of ref 32 for the position of the antibiotic in the product. Given the aforementioned uncertainities, our aim in this work has been to investigate the acylation mechanism of the E. faecium L,D-transpeptidase, Ldtfm, by a carbapenem, ertapenem, using computational techniques. We have employed molecular dynamics simulations in conjunction with QC/MM hybrid potential calculations to map out possible reaction paths. Like previous authors in the field of β-lactam hydrolysis by enzymes with catalytic serine or cysteine, we have been interested in identifying the nucleophilic species (thiol or thiolate), in determining whether the nucleophilic attack by cysteine and ring-opening of the β-lactam are concerted or stepwise, occurring via an oxyanion intermediate, and in identifying the proton acceptors in the first and last steps of the reaction.
■
METHODS Initial Preparation of the Acyl-Enzyme Models. The enzyme model was built using the NMR structure of the catalytic domain (residues 338−446) of Ldtfm acylated by ertapenem (PDB code 3ZGP).29 The structure of the model of lowest energy (model 1) among the 20 refined deposited structures was used in the current study. The protonation states of the titratable residues were set in accordance with the NMR results that were obtained at the physiological pH of 6.4. These were confirmed by performing Poisson−Boltzmann electrostatic calculations combined with Monte Carlo titration. The calculations were done using the continuum electrostatics module35 of the pDynamo molecular modeling library,36 coupled to the MEAD Poisson−Boltzmann equation solver,37 and to the Monte Carlo sampling program, GMCT.38 A starting solvated acyl-enzyme structure was generated using the psfgen module of VMD.39 The initial NMR structure was solvated in an orthorhombic water box of approximate dimensions 65.6 × 51.8 × 55.4 Å3 to ensure at least a 10 Å layer of water on each side of the protein, using the solvate module of VMD. In addition, Na+ cations were included to make the total system neutral. The complete model consisted of 129 residues, 5126 water molecules, and 12 Na+ cations (17 428 atoms in total). MM parameters for all residues, except the ertapenem-modified cysteine residue, were taken from the CHARMM27 force field,40 and the TIP3P model for water.41 The parameters and atomic charges for the ertapenem-modified cysteine were assigned by analogy to preexisting compounds in B
DOI: 10.1021/acs.jpcb.6b02836 J. Phys. Chem. B XXXX, XXX, XXX−XXX
Article
The Journal of Physical Chemistry B the force field, and its topology file is provided in Table S1 of the Supporting Information. A first acyl-enzyme protein model was generated through a 5 ns MD simulation of the previous starting solvated acyl-enzyme structure, using the NAMD software,42 with harmonic positional restraints on all non-hydrogen protein atoms. This step was performed to facilitate the diffusion of water molecules toward the ertapenem-binding site, as some of them may play a key role in the acylation mechanism. The obtained structure was further optimized through our QC/MM protocol (see following section) to yield Model1. A second structure, Model2, was obtained by carrying out a 51 ns MD simulation with no restraints on the atoms with the aim of exploring for alternative acyl-enzyme conformations. Both MD simulations were performed at a temperature of 300 K using a constant temperature algorithm. QC/MM Model Setup. All QC/MM calculations were done using Python scripts prepared within the framework of the pDynamo library,36 coupled to the ORCA QC program.43 Unless otherwise mentioned, the QC region for all calculations comprised the carbapenem nucleus of ertapenem (Ert), the catalytic cysteine Cys442, the δ-protonated His421 residue, the backbone atoms of Asp422, the N and HN atoms of Ile443, and the C, O, CA, HA1, and HA2 atoms of Gly441. The side-chain atoms of ertapenem beyond the sulfur atom of the carbapenem nucleus (SAX in Table S1) were placed in the MM region, rather than the QC region, because they are far from the reaction center (greater than 7 Å from the sulfur of Cys442 for both Model1 and Model2) and because they were found to be very flexible in the NMR structures. This is also consistent with overlays of X-ray structures of LdtMt2 acylated by Meropenem. This basic model consisted of 68 atoms and is shown in Figure 1. Some of the calculations on Model1 and Model2 were done with an additional water molecule in the QC region, raising the total number of QC atoms to 71.
QC/MM boundary were treated using the link atom approach.36 During the QC/MM calculations, different levels of restraint were applied to atoms of the MM region, depending on their proximity to the QC region. No restraints were applied to atoms that were within 8 Å of the QC region. At distances between 8 and 16 Å, the MM atoms were harmonically restrained with a force constant that increased linearly with distance from 0 to 12 kcal mol−1 Å−2. Beyond 16 Å from the QC region, all MM atoms were restrained with a maximum force constant of 12 kcal mol−1 Å−2. The numbers of unrestrained (restrained) atoms in Model1 and Model2 were 1098 (13 011) and 1012 (13 137), respectively. Reaction-Path Search. The only experimental structures available with antibiotics bound to the enzyme are for the acylenzyme. There are some modeled structures of Ldts docked with carbapenems, but some of them31,32 show localizations of the carbapenems that contradict experimental data, and others29 may not be completely reliable due to the arbitrary restraint of 5 Å between the β-lactam carbonyl and cysteine sulfur used during the CNS minimization protocol. For this reason, we studied the reaction in the reverse direction by starting with the known acyl-enzyme structure. Possible reaction mechanisms were studied with a mixture of potential energy surface (PES) scans and reaction-path location using the nudged elastic band (NEB) method. 43,47,48 Exploratory calculations were done with the PES scan technique in which normally one geometric variable, usually a distance, was chosen as the reaction coordinate variable. These variables were systematically incremented or decremented by regular amounts (for example, ±0.1 Å for a distance), and each new structure was fully reoptimized, except for restraints that kept the reaction coordinate variables around their current reference values. Paths coming from these PES scans typically consisted of between 5 and 30 structures, depending on the reaction step, and were stopped when a terminal, minimumenergy structure was identified. The latter was then fully optimized, without any restraints, to serve as a possible intermediate for subsequent steps in the mechanism. It is worth mentioning that many such scans, using different starting structures and reaction coordinate variables, were often necessary to generate reasonable reaction steps. Once the PES scans had given rise to feasible mechanisms with well-defined intermediates, the reaction paths between reactants, products, and intermediates were refined using the NEB method as implemented in the pDynamo library.47,48 The NEB method requires a trajectory of structures as an initial guess that approximates the path between two terminal minimum-energy structures. For simple reaction steps, structures generated by linear interpolation between end-points were adequate, but it was normally much more efficient to employ starting trajectories containing structures taken directly from the PES scans. Typical NEB pathways consisted of between 11 and 21 structures. As noted previously,48 the highest energy point interpolated from a NEB pathway is often a very good approximation to the transition-state structure for the reaction. After full NEB optimization, single point QC/MM calculations were performed with the higher-level QC/MM method on each path structure to obtain the final energy profile.
Figure 1. QC region employed in the QC/MM calculations. Some calculations also included an additional water molecule.
The atoms in the QC region were treated with a DFT method using the B3LYP functional,44−46 whereas the atoms in the MM region employed the same force field as in the previous section. For all QC/MM optimizations a medium size basis set, 6-31G(d), was used, whereas final energies were evaluated by performing single point calculations on the optimized structures using a larger TZVP basis set, which is of triple-ζ quality and has both polarization and diffuse functions. All QC/MM calculations were carried out with full electrostatic embedding, which incorporates MM charges in the QC wave function, and severed bonds between QC and MM atoms at the C
DOI: 10.1021/acs.jpcb.6b02836 J. Phys. Chem. B XXXX, XXX, XXX−XXX
Article
The Journal of Physical Chemistry B
■
RESULTS Acyl-Enzyme Models. Two different structures of the acylenzyme, Model1 and Model2, were employed during this study. The final solvated QC/MM optimized Model1 structure can be overlaid with the starting NMR structure (PDB code 3ZGP) with an RMSD of 0.96 and 1.49 Å on backbone and all protein atoms, respectively. These values are slightly higher than the corresponding 0.30 and 0.50 Å values for the RMSD between the 20 NMR models, but are well within the RMSD values between the Ldtfm and LdtMt2 acyl-enzymes (1.95, 1.97, 2.00, and 2.24 Å with the B chains of 4GSU, 3VYP, 4QR7, and 4QTF, respectively). The RMSD values between the LdtMt2 acyl-enzymes themselves are in the range 0.17−0.28 Å. Differences between the initial NMR structure and Model1 mainly arise from the so-called flap region between the acyldonor and acyl-acceptor cavities, residues 388−409 of Ldtfm, which are dynamic in solution (see below and Figure S2), as well as in the crystal.49 The overall orientation of the carbapenem nucleus in the acyl-donor cavity of Ldtfm is comparable to that in the initial NMR structure. It is stabilized in this position essentially through two hydrogen bonds: the pyrroline amino proton (HAV) interacts with the oxygen of the side-chain of Tyr403 of the flap at 1.99 Å, and the carboxylate oxygen (OAI) is hydrogen bonded (1.73 Å) to the hydroxyl hydrogen of the same tyrosine side-chain. In this conformation, the hydroxyl proton (HAK) of the β-lactam hydroxyethyl substituent is only distant by 1.89 Å from the Lys394 carbonyl oxygen (see Table S1 in the Supporting Information for the nomenclature of the ertapenem atoms). These stabilizing interactions are consistent with experimental results, as the first hydrogen bond is also present in the initial NMR structure (with a distance of 2.11 Å in 3ZGP), and both interactions are conserved in the 3VYP LdtMt2 acyl-enzyme structure, where Tyr318 plays the same role as Tyr403 in Ldtfm. Conversely, the interaction between the βlactam carbonyl oxygen (OAH) and the Gly441-Cys442 amide proton is significantly increased in Model1 (4.40 Å) in comparison to the initial Ldtfm NMR structure (2.06 Å) or its equivalent in the LdtMt2 structures (1.63−2.47 Å). It is replaced by the close proximity of this oxygen to the pyrroline amino proton (1.89 Å). The second solvated QC/MM optimized Model2 structure obtained after the 51 ns unconstrained MD simulation differs more from the NMR structure than the previous model, with RMSDs with respect to 3ZGP of 1.95 and 2.54 Å for backbone and all protein atoms, respectively. It nevertheless shows similar RMSDs to the LdtMt2 acyl-enzyme structures as the previous model with RMSD values varying from 1.82 to 1.92 Å for all protein atoms. Furthermore, the RMSD value on backbone atoms between the two QC/MM optimized models is 1.92 Å. In the Model2 structure, the orientation of the pyrroline ring is no longer stabilized by its interaction with Tyr403 from the flap. Instead, the carboxylate oxygen OAI establishes a hydrogen bond with a water molecule and the side-chain of Lys372 after an approximately 180° rotation. Although the acylated form of Ldtfm is found to be dynamically quite stable during the 51 ns unconstrained simulation used to generate Model2, significant fluctuations (see Figure S2) appear during the molecular dynamics on a nanosecond time scale. Most of these fluctuations are localized in the loop and turn regions (residues 371−378 and 430−441) and the flap (residues 394− 404 including Tyr403) of the protein that define access of the
antibiotic to the catalytic acyl-donor site. We note that the hydrogen bond between the pyrroline carboxylate and Tyr403 breaks and forms several times during the course of the simulation, indicating that this structural change is reversible and transient. A plot of the corresponding Tyr403-ertapenem distance during the molecular dynamics simulation is shown in Figure S3. This conformational flexibility on the nanosecond time scale is consistent with the presence of two distinct Meropenem conformations in published X-ray structures of acylated LdtMt2 (PDB codes 3VYP and 4GSU). One of these is similar to the Ldtfm NMR structure, whereas the other has an alternative orientation in which the carboxylate group of the 2pyrroline ring is solvated by a water molecule. The latter is, however, different from that found in Model2. Also shown, for comparison, in Figure S2 are the fluctuations calculated from a 51 ns simulation of the apo-form of Ldtfm. Interestingly, this is shown to be even more flexible than the acyl-enzyme form of Ldtfm, which is in accordance with 15N-relaxation measurements by NMR. Figure S4 plots the RMSDs from the 51 ns simulations of the acyl- and apo-forms of Ldtfm, but considering only the nonhydrogen atoms of the catalytic center (for a full definition see the caption to the figure). As expected, the RMSDs are smaller than the full-protein RMSDs of Figure S2, and indicate that the catalytic center is quite stable throughout both simulations. By contrast, though, the acyl-enzyme RMSDs are larger than those of the apoenzyme, due to the greater mobility of the ertapenem atoms. For completeness, we note that the all-atom (backboneatom) catalytic center RMSDs between the starting NMR structure (PDB code 3ZGP) and Model1 and Model2 are 0.96 (0.75) and 1.71 (0.51) Å, respectively. The corresponding values between Model1 and Model2 are 1.50 (0.57) Å. Model1 and Model2 (see Figure 2) show marked differences within the carbapenem β-lactam core. In Model1, the β-lactam carbonyl group points toward the 2-pyrroline ring and the nitrogen of the latter ring is in an S-configuration, with the proton pointing in the same direction as the ring methyl substituent. In Model2, the β-lactam carbonyl group points away from the 2-pyrroline ring (after a 180° crankshaft rotation around the sulfur−carbonyl carbon bond) to make an interaction with the Gly419 amide proton, and the nitrogen of the 2-pyrroline ring is in the opposite R-configuration as in the initial NMR structure. DFT calculations by us on model compounds show that the R- and S-nitrogen configurations of the 2-pyrroline ring are approximately isoenergetic and that the energy barrier to their interconversion is low,
