key: cord-0061847-j12tjvx1
authors: Krivitskaya, Alexandra V.; Khrenova, Maria G.
title: Boronic Acids as Prospective Inhibitors of Metallo-β-Lactamases: Efficient Chemical Reaction in the Enzymatic Active Site Revealed by Molecular Modeling
date: 2021-04-02
journal: Molecules
DOI: 10.3390/molecules26072026
sha: 7a4969a5dfb660028f0a74c32762ea1ee48da97b
doc_id: 61847
cord_uid: j12tjvx1

Boronic acids are prospective compounds in inhibition of metallo-β-lactamases as they form covalent adducts with the catalytic hydroxide anion in the enzymatic active site upon binding. We compare this chemical reaction in the active site of the New Delhi metallo-β-lactamase (NDM-1) with the hydrolysis of the antibacterial drug imipenem. The nucleophilic attack occurs with the energy barrier of 14 kcal/mol for imipenem and simultaneously upon binding a boronic acid inhibitor. A boron atom of an inhibitor exhibits stronger electrophilic properties than the carbonyl carbon atom of imipenem in a solution that is quantified by atomic Fukui indices. Upon forming the prereaction complex between NDM-1 and inhibitor, the lone electron pair of the nucleophile interacts with the vacant p-orbital of boron that facilitates the chemical reaction. We analyze a set of boronic acid compounds with the benzo[b]thiophene core complexed with the NDM-1 and propose quantitative structure-sroperty relationship (QSPR) equations that can predict IC50 values from the calculated descriptors of electron density. These relations are applied to classify other boronic acids with the same core found in the database of chemical compounds, PubChem, and proposed ourselves. We demonstrate that the IC50 values for all considered benzo[b]thiophene-containing boronic acid inhibitors are 30–70 μM.

Bacterial resistance is known from the very beginning of the penicillin era [1] [2] [3] [4] . One of the mechanisms is attributed to the hydrolytic activity of the bacterial enzymes towards β-lactams. These enzymes are called β-lactamases (BLs), and they hydrolyze the C-N bond of the β-lactam ring of antibiotics forming the inactivated compounds, β-amino acids [5, 6] . The antibiotic hydrolysis is initiated by various nucleophiles depending on the particular type of β-lactamases [6, 7] . The B class BLs are metallo-β-lactamases (MBLs) that comprise one or two zinc cations in the active site [7] . These are divided into three subclasses, B1, B2, and B3 [8, 9] . The most clinically relevant representative in the B1 subclass is the NDM-1 (New Delhi metallo-β-lactamase-1) [10] . It demonstrates substrate promiscuity [11] , including hydrolysis of carbapenems [12] . In addition, it causes particular alertness due to its high spread between different bacterial species by lateral gene transfer favored by globalization and travel [13] .

The active site of NDM-1 contains two zinc cations, Zn1 2+ and Zn2 2+ , bridged by a catalytic hydroxide ion, OH − . The first zinc cation (Zn1 2+ ) is coordinated by His120, His122 and His189, while the second zinc cation (Zn2 2+ ) is coordinated by Asp124, Cys208 and His250 ( Figure 1A ). In the binding site, Zn1 2+ orients and polarizes the carbonyl group of the substrate, facilitating the chemical reaction. The above-mentioned hydroxide anion is a In this work, we focus on two features of boron-containing inhibitors. First, we compare interactions between the nucleophilic OHof the NDM-1 active site and its substrate imipenem with one of the boronic acid inhibitors (cpd5, according to Ref. [26] ). To do this, we calculate the Gibbs energy profiles of the nucleophilic attack step using molecular dynamics (MD) simulations with the combined quantum mechanics/molecular mechanics QM(PBE0-D3/6-31G**)/MM potentials. Additionally, we compare the dynamic behavior of the corresponding reactant complexes with the same compounds in water solution, measuring atomic Fukui electrophilicity indices of the carbonyl carbon (in an imipenem) and boron (in a boronic acid inhibitor) atoms. Next, we consider five boronic acid inhibitors with the benzo[b]thiophene core with known IC50 values (concentration of inhibitor that results in 50% inhibition of the enzyme activity) [26] corresponding to the NDM-1 inhibition. We study covalent adducts, hydroxyboronate anions in the active site of NDM-1 by analysis of equilibrium geometry configurations of the NDM-1 with these compounds and suggest electron density-based criteria responsible for the inhibitor potency. These features were calculated at the bond critical points determined in terms of the quantum theory of atoms in molecules (QTAIM) [38, 39] . It was already shown that such an approach could explain macroscopic features of chemical reactions occurring in the active sites of enzymes [40] [41] [42] [43] [44] . We examined other known compounds with the same core and predicted their inhibition potency based on the correlation found for compounds with the known experimental IC50 values. In this work, we focus on two features of boron-containing inhibitors. First, we compare interactions between the nucleophilic OH − of the NDM-1 active site and its substrate imipenem with one of the boronic acid inhibitors (cpd5, according to Ref. [26] ). To do this, we calculate the Gibbs energy profiles of the nucleophilic attack step using molecular dynamics (MD) simulations with the combined quantum mechanics/molecular mechanics QM(PBE0-D3/6-31G**)/MM potentials. Additionally, we compare the dynamic behavior of the corresponding reactant complexes with the same compounds in water solution, measuring atomic Fukui electrophilicity indices of the carbonyl carbon (in an imipenem) and boron (in a boronic acid inhibitor) atoms. Next, we consider five boronic acid inhibitors with the benzo[b]thiophene core with known IC50 values (concentration of inhibitor that results in 50% inhibition of the enzyme activity) [26] corresponding to the NDM-1 inhibition. We study covalent adducts, hydroxyboronate anions in the active site of NDM-1 by analysis of equilibrium geometry configurations of the NDM-1 with these compounds and suggest electron density-based criteria responsible for the inhibitor potency. These features were calculated at the bond critical points determined in terms of the quantum theory of atoms in molecules (QTAIM) [38, 39] . It was already shown that such an approach could explain macroscopic features of chemical reactions occurring in the active sites of enzymes [40] [41] [42] [43] [44] . We examined other known compounds with the same core and predicted their inhibition potency based on the correlation found for compounds with the known experimental IC50 values.

The model systems containing boronic acid inhibitors in the active site of NDM-1 were constructed using crystal structures PDB ID: 6IBV, 6Q2Y [35] . They corresponded to complexes of the NDM-1 with cpd2 and cpd3, respectively. The complexes with other experimentally studied boronic acid inhibitors, cpd1, cpd4 and cpd5 ( Figure 2 , compounds are numbered as in Ref. [26] ) and other boronic acid-containing compounds with the benzo[b]thiophene core, cpd6-cpd15, were constructed by motifs of these X-ray structures. Substituted benzo[b]thiophene-containing boronic acids were found in the database of chemical compounds, PubChem (cpd6-cpd10, Figure 2 ) [45] . Additionally, we proposed ones with different electron-withdrawing groups (cpd11-cpd15, Figure 2 ). 

The model systems containing boronic acid inhibitors in the active site of NDM-1 were constructed using crystal structures PDB ID: 6IBV, 6Q2Y [35] . They corresponded to complexes of the NDM-1 with cpd2 and cpd3, respectively. The complexes with other experimentally studied boronic acid inhibitors, cpd1, cpd4 and cpd5 (Figure 2 , compounds are numbered as in ref [26] ) and other boronic acid-containing compounds with the benzo[b]thiophene core, cpd6-cpd15, were constructed by motifs of these X-ray structures. Substituted benzo[b]thiophene-containing boronic acids were found in the database of chemical compounds, PubChem (cpd6-cpd10, Figure 2 ) [45] . Additionally, we proposed ones with different electron-withdrawing groups (cpd11-cpd15, Figure 2 ).

Boronic acid inhibitors considered in this study. Cpd1-cpd5 are experimentally studied in. [26] , cpd6-cpd10 are found in the PubChem [45] and cpd11-cpd15 are suggested in this study.

The enzyme-substrate complex of NDM-1 and imipenem was obtained from the crystal structure PDB ID: 5YPK [12] with the hydrolyzed imipenem bound to NDM-1. The structure of the imipenem was restored to its non-hydrolyzed state.

In all complexes, hydrogen atoms were added using the Reduce program [46] to reproduce the protonation states of amino acids at neutral pH. The systems were solvated in rectangular water boxes so that the distance from the protein surface to the border of the cell exceeds 12 Å . These systems were neutralized by adding sodium or chloride ions. All atom force fields were utilized; CHARMM36 [47, 48] for the protein, TIP3P [49] for water molecules and CGenFF [50] for the imipenem and inhibitors. Systems were equilibrated for 20 ns; the typical RMSD graph is shown in Figure S1 . It was performed at T = 300 K and p = 1 atm with 1 fs integration time step. The substrate/inhibitor, catalytic OH -, zinc cations, side chains of amino acids forming coordination bonds with Zn 2+ and a sidechain of Asp124 were fixed during this preliminary run. All MD simulations were performed in the NAMD program [51] . Last frames from the MD runs were utilized for the subsequent preparation of the QM/MM models. Water molecules located further than 6 Å from protein or 10 Å from the active site and sodium/chloride ions were removed from the model systems. Equilibrium geometry configurations were obtained at the QM(PBE0-D3/6-31G**)/MM(AMBER) [52] [53] [54] [55] level of theory using the NWChem program [56] . An electronic embedding scheme was utilized with no cutoff on electrostatic interactions. The QM subsystem included the substrate or inhibitor, two Zn 2+ cations and side chains of amino acid residues that form coordination bonds with them (His120, His122, His189, Cys208, His250), catalytic hydroxide anion and side chains of Asp124 and Asn220 that form hydrogen bonds with the substrate or inhibitor. Together, we obtained a set of 15 NDM-1 complexes with hydroxyboronate anion-containing compounds, EI.

Equilibrium geometry configurations of the NDM-1 with cpd5 (E·cpd5) and the ES complex after preliminary classical MD relaxation were utilized as initial structures for the molecular dynamics simulations with the QM/MM potentials. Solvation and neutralization of model systems, QM and MM partitioning and details of MD simulations are discussed above. The NAMD [51] Figure 2 . Boronic acid inhibitors considered in this study. Cpd1-cpd5 are experimentally studied in [26] , cpd6-cpd10 are found in the PubChem [45] and cpd11-cpd15 are suggested in this study.

The enzyme-substrate complex of NDM-1 and imipenem was obtained from the crystal structure PDB ID: 5YPK [12] with the hydrolyzed imipenem bound to NDM-1. The structure of the imipenem was restored to its non-hydrolyzed state.

In all complexes, hydrogen atoms were added using the Reduce program [46] to reproduce the protonation states of amino acids at neutral pH. The systems were solvated in rectangular water boxes so that the distance from the protein surface to the border of the cell exceeds 12 Å. These systems were neutralized by adding sodium or chloride ions. All atom force fields were utilized; CHARMM36 [47, 48] for the protein, TIP3P [49] for water molecules and CGenFF [50] for the imipenem and inhibitors. Systems were equilibrated for 20 ns; the typical RMSD graph is shown in Figure S1 . It was performed at T = 300 K and p = 1 atm with 1 fs integration time step. The substrate/inhibitor, catalytic OH − , zinc cations, side chains of amino acids forming coordination bonds with Zn 2+ and a side-chain of Asp124 were fixed during this preliminary run. All MD simulations were performed in the NAMD program [51] . Last frames from the MD runs were utilized for the subsequent preparation of the QM/MM models. Water molecules located further than 6 Å from protein or 10 Å from the active site and sodium/chloride ions were removed from the model systems. Equilibrium geometry configurations were obtained at the QM(PBE0-D3/6-31G**)/MM(AMBER) [52] [53] [54] [55] level of theory using the NWChem program [56] . An electronic embedding scheme was utilized with no cutoff on electrostatic interactions. The QM subsystem included the substrate or inhibitor, two Zn 2+ cations and side chains of amino acid residues that form coordination bonds with them (His120, His122, His189, Cys208, His250), catalytic hydroxide anion and side chains of Asp124 and Asn220 that form hydrogen bonds with the substrate or inhibitor. Together, we obtained a set of 15 NDM-1 complexes with hydroxyboronate anion-containing compounds, EI.

Equilibrium geometry configurations of the NDM-1 with cpd5 (E·cpd5) and the ES complex after preliminary classical MD relaxation were utilized as initial structures for the molecular dynamics simulations with the QM/MM potentials. Solvation and neutralization of model systems, QM and MM partitioning and details of MD simulations are discussed above. The NAMD [51] program was utilized for MD steps performance and calculation of forces and energies of the MM subsystem. The TeraChem program [57] was used to calculate forces in the QM region, and a NAMD-TeraChem interface [58] was applied. The cutoff distance on the point charges of the MM subsystem contributing to the QM Hamiltonian was 12 Å. Gibbs energy profiles of the OH − nucleophilic attack were calculated using the umbrella sampling approach. For the ES model, the reaction coordinate was the distance between the carbonyl carbon atom of the imipenem and the oxygen atom of OH − , ξ(ES) (Figure 3 ). Harmonic potentials 1 /2·K·(ξ − ξ 0 ) 2 were centered at 1.3, 1.5, 1.7, 1.9, 2.1, 2.4, 2.6 and 2.8 Å with the force constant K = 40 kcal/(mol·Å 2 ); additional runs were performed with potentials centered at 1.7, 1.8, 1.9 Å and K = 80 kcal/(mol·Å 2 ); 1.7, 1.8, 1.9, 2.0, 2.1 Å and K = 120 kcal/(mol·Å 2 ). A complex collective variable was suggested as a reaction coordinate of the nucleophilic addition reaction in the E·cpd5 complex, ξ(E·cpd5); it was set as a sum of two coordination bond distances, d(Zn2 2+ . . . O2) and d(Zn1 2+ . . . O1), and a distance of nucleophilic attack, d(B . . . O w ) ( Figure 3 ). Harmonic potentials were centered at the following ξ 0 values: from 5.3 Å to 10.4 Å with 0.3 Å increment, from 10.5 Å to 10.8 Å with 0.1 Å increment, from 11.0 Å to 11.6 with 0.2 Å increment and 11.9 Å. The force constant was 40 kcal/(mol·Å 2 ) in all simulations. The reaction coordinates were analyzed between 5 Å and 12 Å and divided into 200 bins in the case of E·cpd5 and between 1.3 Å and 3 Å with 50 bins for the complex NDM-1imipenem complex. Each MD trajectory with additional harmonic potential added on the reaction coordinate was not less than 5 ps, and the first 1 ps was excluded from the data analysis. The weighted histogram analysis method (WHAM) was utilized to restore the Gibbs energy profile. The lengths of the QM/MM trajectories with harmonic potentials are considerably shorter than the conventional classical MD trajectories. However, we analyze the distributions that we obtain with each harmonic potential and check whether those changes with the further elongation of MD runs. Figures S2 and S3 depict distributions obtained with different harmonic potentials. MD trajectories of the comparable lengths were calculated in other enzymatic reaction mechanism studies; see, for example, Ref. [59] . Additional 10 ps runs were performed for the ES complex and E·cpd5 with the ξ 0 = 10.1 Å and K = 40 kcal/(mol·Å 2 ) harmonic potential. calculation of forces and energies of the MM subsystem. The TeraChem program [57] was used to calculate forces in the QM region, and a NAMD-TeraChem interface [58] was applied. The cutoff distance on the point charges of the MM subsystem contributing to the QM Hamiltonian was 12 Å . Gibbs energy profiles of the OHnucleophilic attack were calculated using the umbrella sampling approach. For the ES model, the reaction coordinate was the distance between the carbonyl carbon atom of the imipenem and the oxygen atom of OH -, ξ(ES) (Figure 3 ). Harmonic potentials ½·K·(ξ-ξ0) 2 were centered at 1.3, 1.5, 1.7, 1.9, 2.1, 2.4, 2.6 and 2.8 Å with the force constant K = 40 kcal/(mol·Å 2 ); additional runs were performed with potentials centered at 1.7, 1.8, 1.9 Å and K = 80 kcal/(mol·Å 2 ); 1.7, 1.8, 1.9, 2.0, 2.1 Å and K = 120 kcal/(mol·Å 2 ). A complex collective variable was suggested as a reaction coordinate of the nucleophilic addition reaction in the E·cpd5 complex, ξ(E·cpd5); it was set as a sum of two coordination bond distances, d(Zn2 2 + …O2) and d(Zn1 2+ …O1), and a distance of nucleophilic attack, d(B…Ow) (Figure 3 ). Harmonic potentials were centered at the following ξ0 values: from 5.3 Å to 10.4 Å with 0.3 Å increment, from 10.5 Å to 10.8 Å with 0.1 Å increment, from 11.0 Å to 11.6 with 0.2 Å increment and 11.9 Å . The force constant was 40 kcal/(mol·Å 2 ) in all simulations. The reaction coordinates were analyzed between 5 Å and 12 Å and divided into 200 bins in the case of E·cpd5 and between 1.3 Å and 3 Å with 50 bins for the complex NDM-1-imipenem complex. Each MD trajectory with additional harmonic potential added on the reaction coordinate was not less than 5 ps, and the first 1 ps was excluded from the data analysis. The weighted histogram analysis method (WHAM) was utilized to restore the Gibbs energy profile. The lengths of the QM/MM trajectories with harmonic potentials are considerably shorter than the conventional classical MD trajectories. However, we analyze the distributions that we obtain with each harmonic potential and check whether those changes with the further elongation of MD runs. Figures S2 and S3 depict distributions obtained with different harmonic potentials. MD trajectories of the comparable lengths were calculated in other enzymatic reaction mechanism studies; see, for example, ref [59] . Additional 10 ps runs were performed for the ES complex and E·cpd5 with the ξ0 = 10.1 Å and K = 40 kcal/(mol·Å 2 ) harmonic potential. Additional models were constructed for the cpd5 and imipenem (in the corresponding reactant states) solvated in the rectangular water box with the distance to the cell border being larger than 15 Å. The 1 ns preliminary classical MD simulation was performed before the 10 ps QM/MM MD run. Computational protocols were the same as discussed above.

Analysis of all MD trajectories and visualization of the equilibrium geometry configurations were performed in the VMD program [60] .

We performed an analysis of electron density-based and geometry descriptors at the stationary points on potential energy surfaces and along QM/MM trajectories. In the latter case, we selected a set of 100 MD frames equally distributed along each considered MD trajectory. The electron densities were calculated, taking into account partial atomic charges from the MM subsystems that contributed to the one-electron part of the QM Hamiltonian. All electron density-based descriptors were calculated using the Multiwfn program [61] . We calculated electron density-based descriptors at bond-critical points (BCPs) corresponding to the coordination bonds between the zinc cations and hydroxyboronate anion-containing inhibitors (cpd1-cpd15) at minima on potential energy surfaces. These are electron density, ρ(r), and Laplacian of electron density, ∇ 2 ρ(r). In addition, we calculate atomic contribution S(r, Ω i ) [62] [63] [64] evaluated as

where Ω i is the atomic basin of the i-th atom; here, we consider oxygen atoms O1, O2 and O w that are covalently bound to a boron atom (Figure 1 ). More details on the utilization of electron density-based descriptors for biologically relevant systems are described in Refs. [41, 43] . In the case of QM/MM MD trajectories, we calculated atomic Fukui electrophilicity indices, f+, on a carbonyl carbon atom of the imipenem and a boron atom of cpd5. These are evaluated as differences between Hirshfeld charges [65] calculated for the model system with N electrons (as it is set in the calculations) and N + 1 electrons [66] [67] [68] .

3.1. Imipenem and Boronic Acid Inhibitor Cpd5 in Water Solution and in the Active Site of NDM-1

The boronic acid inhibitors form covalent adducts with the catalytic OHion as known from the crystal structures [26] . Here, we study the mechanism of this chemical reaction and compare it with the same nucleophilic addition step during the hydrolysis reaction of the imipenem substrate.

First, we analyze the behavior of cpd5 and imipenem in solution in terms of the electrophilicity of the boron and carbon atoms that form covalent bonds with the OH − during the reaction. We use the Fukui function for the electrophilic attack, f+, calculated for the C (in the imipenem) and B (in the cpd5) atoms along the MD trajectories ( Figure 4) . The f+ distribution is considerably shifted to the larger values in the case of B, indicating that it is a stronger nucleophile in the solution. The f+ value for the carbonyl carbon atom in the imipenem is 0.032 ± 0.015 a.u. and 0.074 ± 0.011 a.u. for B of the cpd5. Next, we study the behavior of these two compounds in the active site of NDM-1. We constructed the ES complex with an imipenem molecule. During the QM/MM MD simulation, the distance of the nucleophilic attack varied in the range of 2.89 ± 0.15 Å . The mean value of the f+(C) increases more than twice compared with water solution being 0.084 a.u (Figure 4 ). This distribution is considerably broader, with the standard deviation being 0.031 a.u. that is twice larger than in solution. It means that not all ES complexes (frames along the MD trajectory) correspond to the reactive species. The fraction with the larger f+ values demonstrates the substrate activation in the active site of the enzyme discussed in Refs [59, 69] . The Gibbs energy profile for the nucleophilic attack step was calculated using the umbrella sampling technique ( Figure 5 ). The minimum corresponding to the ES complex is located at the reaction coordinate value of 2.79 Å . It is slightly shorter than the mean C…Ow distance along the MD trajectory of the ES complex. During the unconstrained MD simulation of the ES complex, we observe both reactive and nonreactive species, whereas only reactive species contribute to the chemical reaction [43] . The C…Ow distance at the transition state is 1.92 Å that is typical for this type of reaction occurring in the zinc-dependent proteins. To compare, it is 1.86 Å in the case of nitrocefin hydrolysis by L1 metallo-β-lactamase [70] and 1.80 Å during the oligopeptide hydrolysis by MMP-2 matrix metalloproteinase [71] . The energy barrier at this step is about 14 kcal/mol, and the first intermediate, Int, is ~4 kcal/mol higher in energy than the ES complex. Next, we study the behavior of these two compounds in the active site of NDM-1. We constructed the ES complex with an imipenem molecule. During the QM/MM MD simulation, the distance of the nucleophilic attack varied in the range of 2.89 ± 0.15 Å. The mean value of the f+(C) increases more than twice compared with water solution being 0.084 a.u (Figure 4 ). This distribution is considerably broader, with the standard deviation being 0.031 a.u. that is twice larger than in solution. It means that not all ES complexes (frames along the MD trajectory) correspond to the reactive species. The fraction with the larger f+ values demonstrates the substrate activation in the active site of the enzyme discussed in Refs. [59, 69] . The Gibbs energy profile for the nucleophilic attack step was calculated using the umbrella sampling technique ( Figure 5 ). The minimum corresponding to the ES complex is located at the reaction coordinate value of 2.79 Å. It is slightly shorter than the mean C . . . O w distance along the MD trajectory of the ES complex. During the unconstrained MD simulation of the ES complex, we observe both reactive and nonreactive species, whereas only reactive species contribute to the chemical reaction [43] . The C . . . O w distance at the transition state is 1.92 Å that is typical for this type of reaction occurring in the zinc-dependent proteins. To compare, it is 1.86 Å in the case of nitrocefin hydrolysis by L1 metallo-β-lactamase [70] and 1.80 Å during the oligopeptide hydrolysis by MMP-2 matrix metalloproteinase [71] . The energy barrier at this step is about 14 kcal/mol, and the first intermediate, Int, is~4 kcal/mol higher in energy than the ES complex. Next, we study the behavior of these two compounds in the active site of NDM-1. We constructed the ES complex with an imipenem molecule. During the QM/MM MD simulation, the distance of the nucleophilic attack varied in the range of 2.89 ± 0.15 Å . The mean value of the f+(C) increases more than twice compared with water solution being 0.084 a.u (Figure 4 ). This distribution is considerably broader, with the standard deviation being 0.031 a.u. that is twice larger than in solution. It means that not all ES complexes (frames along the MD trajectory) correspond to the reactive species. The fraction with the larger f+ values demonstrates the substrate activation in the active site of the enzyme discussed in Refs [59, 69] . The Gibbs energy profile for the nucleophilic attack step was calculated using the umbrella sampling technique ( Figure 5 ). The minimum corresponding to the ES complex is located at the reaction coordinate value of 2.79 Å . It is slightly shorter than the mean C…Ow distance along the MD trajectory of the ES complex. During the unconstrained MD simulation of the ES complex, we observe both reactive and nonreactive species, whereas only reactive species contribute to the chemical reaction [43] . The C…Ow distance at the transition state is 1.92 Å that is typical for this type of reaction occurring in the zinc-dependent proteins. To compare, it is 1.86 Å in the case of nitrocefin hydrolysis by L1 metallo-β-lactamase [70] and 1.80 Å during the oligopeptide hydrolysis by MMP-2 matrix metalloproteinase [71] . The energy barrier at this step is about 14 kcal/mol, and the first intermediate, Int, is ~4 kcal/mol higher in energy than the ES complex. The Gibbs energy profile of the covalent bond formation between the OH − and the cpd5 was studied within the same QM/MM protocol. First, we performed the QM/MM MD run of the product state (E·cpd5) with the covalent bond between the cpd5 and OH − . We gradually increased the reaction coordinate and performed sequential umbrella sampling runs. The first attempt was a simple reaction coordinate, d(B . . . O w ), similarly to the reaction with the imipenem. However, this resulted in the artificial behavior of the model system. Therefore, we proposed a more complex reaction coordinate that accounts not only for the distance of nucleophilic attack but also for two coordination bonds between oxygen atoms, O1 and O2, of the cpd5 and zinc cations (Figure 3) . The reaction coordinate equals 5.8 Å at the minimum, corresponding to the covalent complex E·cpd5 with two coordination bonds, Zn1 2+ . . . O1 and Zn2 2+ . . . O2, and a covalent B-O w bond being formed. The dissociation energy profile has two regions. A smaller reaction coordinate values, the rapid increase of energy results in the cleavage of the covalent bond B-O w and elongation of coordination bonds. Starting from the ξ = 8 Å, the energy increases much slower; the covalent bond is already cleaved, and the dissociation process is mostly related to the complete cleavage of coordination bonds. One should keep in mind that even at the ξ = 12 Å, the NDM-1 and cpd5 complex is not fully dissociated. More precisely, the coordination bonds between the cpd5 and zinc cations are cleaved, but components of this complex are not fully solvated. We expect that at larger ξ values, the Gibbs energy will decrease.

We failed to obtain a stable prereaction complex of the NDM-1 and the cpd5. However, we analyzed the dynamic behavior of the system in the geometry configuration that is similar to the ES complex. We performed a constrained MD run with the ξ 0 = 10.1 Å and the force constant imposed on the reaction coordinate K = 40 kcal/(mol·Å 2 ). The B . . . O w distance varied in the range of 2.97 ± 0.13 Å. We expected that the f+ on the boron atom is higher than in the solution. However, it turned out that it was considerably smaller and similar to that of a carbon atom of imipenem in the solution being 0.040 ± 0.012 a.u. (Figure 4 ). To analyze this unexpected phenomenon, we calculated the Laplacian of electron density, ∇ 2 ρ(r), along the C . . . O w and B . . . O w at different frames. Figure 6 demonstrates the alignment of ∇ 2 ρ(r) for complexes of the NDM-1 with the cpd5 or the imipenem with the same distances of the nucleophilic attack. The ∇ 2 ρ(r) curves are similar around the oxygen atom of the OH − and are different in the B/C region. In the case of a carbon atom, we observe a minimum at~0.5 Å and a maximum at~0.7 Å from a C nucleus. In the case of a boron atom, we observe a flat region between two atoms with the small deconcentration of electron density (∇ 2 ρ(r) > 0).

We analyzed the ∇ 2 ρ(r) curves at different MD frames. For different NDM-1 with cpd5 complexes, we observe similar behavior of the Laplacian of electron density in the interatomic region: a flat area with the electron density deconcentration, with ∇ 2 ρ(r) being around 0.05 a.u. Such regions are attributed to the presence of the substrate activation effect that was previously demonstrated on the 2D maps of ∇ 2 ρ(r) [42, 43, 69, 72] . From this viewpoint, the boron atom is activated in all presented MD frames. However, its Fukui electrophilicity index is lower when bound to NDM-1 compared with the aqueous solution. Similar abnormal behavior of a boron atom is discussed in Ref. [73] and it is attributed to the strong p-π electronic interactions; that is, the empty p-orbital of a boron center is partly filled by the π-electron of the neighboring atom. In our case, the partner of these interactions is a nucleophilic OH − that partly "donates" its electron lone pair to the empty boron p-orbital. It seems that these interactions, together with the long-range electrostatic interactions between the coordination bond partners, may be the reason for the failure to locate a minimum corresponding to the prereaction complex.

Laplacian of electron density curves for the ES complex at different MD frames are diverse. Depending on the particular frame, a minimum that is 0.5 Å apart from a carbon atom has either a positive (solid green line on Figure 6 ) or negative (green dashed line on Figure 6 ) ∇ 2 ρ(r) value. It means that a substrate is activated only at certain frames. It is in line with the presence of equilibrium between the reactive and nonreactive species [43] . 

Another important issue is the prediction of the inhibitor potency of boronic acid compounds. We calculated equilibrium geometry configurations for a set of covalent complexes of NDM-1 and hydroxyboronate anion-containing compounds with available experimental IC50 values [26] . All complexes demonstrate similar geometry features. The length of the newly formed covalent bond between a boron-based inhibitor and a former nucleophilic hydroxide anion, B-Ow, is 1. 45 Figure 3) . Here, the QTAIM [38, 39] is utilized to identify and characterize interatomic interactions. It was already demonstrated that even in large systems, such as enzymes, one could identify key interatomic interactions that are responsible for the observed macroscopic property [41] . The first step in the search of key interatomic interactions is interatomic distances check. Geometry parameters demonstrate a worse correlation. However, if these exist, this can be the hint for electron density-based descriptor search. We found that electron density-based descriptors at the bond critical points (BCPs) of the coordination bonds between the hydroxyboronate anion and zinc cations correlate with the IC50 values (Table 1, Figure 7 ). These are Zn1 2+ …O1, Zn2 2+ …O2 and Zn1 2+ …Ow coordination bonds with the corresponding BCP1, BCP2 and BCP3. We obtained QSPR equations between the IC50 values and coordination bond length, electron density at BCP, Laplacian of electron density at BCP and atomic contribution of an oxygen atom forming coordination bond to a corresponding BCP ( Table 1) . Sums of corresponding values at BCP1, BCP2 and BCP3 were also examined. The best correlations were obtained for the IC50 dependency on the atomic contributions at BCP and sums of descriptors calculated at three BCPs. These QSPR equations were further utilized for the prediction of the IC50 values of new compounds. 

Another important issue is the prediction of the inhibitor potency of boronic acid compounds. We calculated equilibrium geometry configurations for a set of covalent complexes of NDM-1 and hydroxyboronate anion-containing compounds with available experimental IC50 values [26] . All complexes demonstrate similar geometry features. The length of the newly formed covalent bond between a boron-based inhibitor and a former nucleophilic hydroxide anion, B-O w , is 1. 45 Figure 3 ). Here, the QTAIM [38, 39] is utilized to identify and characterize interatomic interactions. It was already demonstrated that even in large systems, such as enzymes, one could identify key interatomic interactions that are responsible for the observed macroscopic property [41] . The first step in the search of key interatomic interactions is interatomic distances check. Geometry parameters demonstrate a worse correlation. However, if these exist, this can be the hint for electron density-based descriptor search. We found that electron density-based descriptors at the bond critical points (BCPs) of the coordination bonds between the hydroxyboronate anion and zinc cations correlate with the IC50 values (Table 1, Figure 7 ). These are Zn1 2+ . . . O1, Zn2 2+ . . . O2 and Zn1 2+ . . . O w coordination bonds with the corresponding BCP1, BCP2 and BCP3. We obtained QSPR equations between the IC50 values and coordination bond length, electron density at BCP, Laplacian of electron density at BCP and atomic contribution of an oxygen atom forming coordination bond to a corresponding BCP ( Table 1) . Sums of corresponding values at BCP1, BCP2 and BCP3 were also examined. The best correlations were obtained for the IC50 dependency on the atomic contributions at BCP and sums of descriptors calculated at three BCPs. These QSPR equations were further utilized for the prediction of the IC50 values of new compounds. Figure 7 . The IC50 = f(S Σ ) dependency; red dots correspond to cpd1-cpd5 inhibitors; cpd6-cpd15 are marked with lavender stars. Prospective compounds, cpd14 and cpd15, are named. Coordination bonds between the hydroxyboronate anion and zinc cations are shown with violet dashed lines, and corresponding bond critical points (BCPs) are highlighted yellow.

NDM-1 complexes were modeled with compounds cpd6-cpd10 presented in the PubChem [45] and having the same benzo[b]thiophene core. However, their predicted IC50 values were found to be higher than 60 µM. Therefore, we returned to the experimentally studied set. The cpd5 has the lowest value of IC50 among those studied experimentally, so it was chosen as a lead compound to design new inhibitors. Namely, we proposed that electron-withdrawing groups are responsible for the IC50 value. We suggested and examined another set of compounds (cpd11-cpd15), replacing the carboxyl group in cpd5 with other electron-withdrawing groups or adding them at different positions. All of them demonstrated IC50 values in the same range. Among them, the most promising are cpd14 and cpd15 (Figure 2) , which have predicted IC50 values 40 µM and 47 µM, respectively (Figure 7) . It seems that a further decrease of the IC50 value can be achieved if changing the core of boronic acid compounds.

Application of molecular dynamics simulations with the QM/MM potentials complemented with on-the-fly calculations of electron density-based descriptors allowed us to study and compare nucleophilic attacks of the catalytic species in the NDM-1 active site. We compared the reaction mechanism with a substrate imipenem and a boronic acid inhibitor. The nucleophilic attack occurs with the energy barrier of 14 kcal/mol in the case of the imipenem and simultaneously upon binding of an inhibitor. Activation of the carbonyl carbon atom of the imipenem is observed in the enzyme-substrate. The boronic acid compound is already activated, even being surrounded by water molecules. These properties are quantified by the atomic Fukui electrophilicity indices. Binding to the active site of NDM-1 comprises both formations of the prereaction interactions between the boron atom of an inhibitor and a catalytic OH − and coordination bonds between OH groups of an inhibitor and zinc cations. The latter interactions are predominantly of an electrostatic nature and are long-range. Upon the formation of the prereaction complex between NDM-1 and inhibitor, the lone electron pair of the nucleophile interacts with the vacant p-orbital of boron that facilitates the chemical reaction. This explains the barrierless process of covalent bond formation between an OH − and a boronic acid if a latter comes to the active site of NDM-1.

We analyzed a set of boronic acid compounds with the benzo[b]thiophene core complexed with the NDM-1 and proposed QSPR equations to predict IC50 values from the calculated descriptors of electron density. Among them, the prospective ones are quantities obtained as sums of electron densities, Laplacian of electron densities or atomic contributions calculated at BCPs of three coordination bonds between oxygen atoms of hydroxyboronate anion and zinc cations. We examined a set of compounds with the same benzo[b]thiophene from the PubChem database and proposed others ourselves. The predicted IC50 values were found to be larger than 30 µM, similar to those studied experimentally. We suppose that a further decrease of this value can be achieved if changing the core of boronic acid compounds. 

An Enzyme from Bacteria able to Destroy Penicillin

Three Decades of β-Lactamase Inhibitors

Antimicrobial Properties on Non-Antibiotic Drugs in the Era of Increased Bacterial Resistance

Bacterial Resistance: Antibiotics of Last Generation used in Clinical Practice and the Arise of Natural Products as New Therapeutic Alternatives

Mechanisms of Antibiotic Resistance

Revised Ambler classification of β-lactamases

Metallo-β-Lactamases: Influence of the Active Site Structure on the Mechanisms of Antibiotic Resistance and Inhibition

Where Do We Stand? Curr. Drug Targets

Metallo-β-lactamase structure and function

Ten Years with New Delhi Metallo-β-lactamase-1 (NDM-1): From Structural Insights to Inhibitor Design

NDM-1, the ultimate promiscuous enzyme: Substrate recognition and catalytic mechanism

The mechanism of NDM-1-catalyzed carbapenem hydrolysis is distinct from that of penicillin or cephalosporin hydrolysis

New Delhi metallo-beta-lactamase (NDM-1): Towards a new pandemia?

bla NDM-1 -positive Klebsiella pneumoniae from Environment

Molecular Basis of NDM-1, a New Antibiotic Resistance Determinant

Shared evolutionary origin of vertebrate neural crest and cranial placodes

Specific NDM-1 Inhibitor of Isoliquiritin Enhances the Activity of Meropenem against NDM-1-positive Enterobacteriaceae in vitro

Discovery of Cyclic Boronic Acid QPX7728, an Ultrabroad-Spectrum Inhibitor of Serine and Metallo-β-lactamases

The assemblage of covalent and metal binding dual functional scaffold for cross-class metallo-β-lactamases inhibition

SAR Studies Leading to the Identification of a Novel Series of Metallo-β-lactamase Inhibitors for the Treatment of Carbapenem-Resistant Enterobacteriaceae Infections That Display Efficacy in an

Triazole Derivatives as New Metallo-β-Lactamase Inhibitors

Halogen-Substituted Triazolethioacetamides as a Potent Skeleton for the Development of Metallo-β-Lactamase Inhibitors

Synthesis and antibacterial evaluation against resistant Gram-negative bacteria of monobactams bearing various substituents on oxime residue

Synergistic effect of doripenem in combination with cefoxitin and tetracycline in inhibiting NDM-1 producing bacteria

Syntheses and studies of new forms of N-sulfonyloxy β-lactams as potential antibacterial agents and β-lactamase inhibitors

Computational and biological profile of boronic acids for the detection of bacterial serine-and metallo-β-lactamases

New Delhi metallo-β-lactamase-1: Structure, inhibitors and detection of producers

β-Lactamases and β-Lactamase Inhibitors in the 21st Century

Will morphing boron-based inhibitors beat the β-lactamases?

Recent developments in the medicinal chemistry of single boron atom-containing compounds

Boron in drug design: Recent advances in the development of new therapeutic agents

The Fifth Element in Drug Design: Boron in Medicinal Chemistry

PS-341): A Novel, First-in-Class Proteasome Inhibitor for the Treatment of Multiple Myeloma and Other Cancers

Cyclic Boronates Inhibit All Classes of β-Lactamases

X-ray Crystallography Deciphers the Activity of Broad-Spectrum Boronic Acid β-Lactamase Inhibitors

Covalent Inhibition of New Delhi Metallo-β-Lactamase-1 (NDM-1) by Cefaclor

Bicyclic Boronates as Potent Inhibitors of AmpC, the Class C β-Lactamase from Escherichia coli

Atoms in Molecules-A Quantum Theory

A quantum theory of molecular structure and its applications

Revealing electronic features governing hydrolysis of cephalosporins in the active site of the L1 metallo-β-lactamase

The QM/MM-QTAIM approach reveals the nature of the different reactivity of cephalosporins in the active site of L1 metallo-β-lactamase

Discrimination of enzyme-substrate complexes by reactivity using the electron density analysis: Peptide bond hydrolysis by the matrix metalloproteinase-2

Dynamical properties of enzyme-substrate complexes disclose substrate specificity of the SARS-CoV-2 main protease as characterized by the electron density descriptors

The N···H hydrogen bond strength in the transition state at the limiting step determines the reactivity of cephalosporins in the active site of L1 metallo-β-lactamase

PubChem in 2021: New data content and improved web interfaces

Asparagine and glutamine: Using hydrogen atom contacts in the choice of side-chain amide orientation

Optimization of the Additive CHARMM All-Atom Protein Force Field Targeting Improved Sampling of the Backbone φ, ψ and Side-Chain χ 1 and χ 2 Dihedral Angles

Impact of 2 -hydroxyl sampling on the conformational properties of RNA: Update of the CHARMM all-atom additive force field for RNA

Comparison of simple potential functions for simulating liquid water

CHARMM general force field (CGenFF): A force field for drug-like molecules compatible with the CHARMM all-atom additive biological force fields

Scalable molecular dynamics with NAMD

Toward reliable density functional methods without adjustable parameters: The PBE0 model

A consistent and accurate ab initio parametrization of density functional dispersion correction (DFT-D) for the 94 elements H-Pu

Self-Consistent Molecular Orbital Methods. XII. Further Extensions of Gaussian-Type Basis Sets for Use in Molecular Orbital Studies of Organic Molecules

The influence of polarization functions on molecular orbital hydrogenation energies

NWChem: A comprehensive and scalable open-source solution for large scale molecular simulations

TeraChem: A graphical processing unit-accelerated electronic structure package for large-scale ab initio molecular dynamics

NAMD goes quantum: An integrative suite for hybrid simulations

Light-Induced Change of Arginine Conformation Modulates the Rate of Adenosine Triphosphate to Cyclic Adenosine Monophosphate Conversion in the Optogenetic System Containing Photoactivated Adenylyl Cyclase

VMD: Visual molecular dynamics

Multiwfn: A multifunctional wavefunction analyzer

The N· · · I halogen bond in substituted pyridines as viewed by the source function and delocalization indices

Chemical information from the source function

A Green's function for the density

Bonded-atom fragments for describing molecular charge densities

Information Conservation Principle Determines Electrophilicity, Nucleophilicity, and Regioselectivity

Stockholders Charge Partitioning Technique. A Reliable Electron Population Analysis Scheme to Predict Intramolecular Reactivity Sequence

Condensed Fukui Functions Derived from Stockholder Charges: Assessment of Their Performance as Local Reactivity Descriptors

Proof of concept for poor inhibitor binding and efficient formation of covalent adducts of KRAS G12C and ARS compounds

Modeling the Transient Kinetics of the L1 Metallo-β-Lactamase

Mechanism of proteolysis in matrix metalloproteinase-2 revealed by QM/MM modeling

Nucleophilic addition to activated double bonds: Predictions of reactivity from the Laplacian of the charge density

Density functional studies on photophysical properties and chemical reactivities of the triarylboranes: Effect of the constraint of planarity