1. Introduction
Drug resistance in bacteria remains a major concern for public health [
1,
2]. Antibiotics have been the tool of choice in the fight against bacterial infections following their inception with penicillin. Due to widespread overuse and poor waste management, resistance mechanisms have spread throughout the global bacterial community, rendering many drugs clinically ineffective [
3,
4]. Despite increasing resistance, beta-lactams remain the most commonly prescribed class of antibiotics [
5]. Beta-lactams are small molecules that contain the 4-membered beta-lactam ring motif, and include penicillin, cephalosporin, monobactam, and carbapenem-like molecules [
6]. These drugs act by binding to transpeptidase (penicillin binding protein), preventing them from linking together the peptidoglycans which hold a bacteria’s cell wall together, ultimately leading to bacteria cell death [
7].
One of the leading causes of resistance to beta-lactams in bacteria is the production of an enzyme called beta-lactamase [
8]. Beta-lactamase proteins are broadly categorized into four major classes, namely, A, B, C, and D, based on structural similarity [
9]. Classes A, C, and D share the same serine-mediated catalytic mechanism for inactivating beta-lactams, while class B is a family of metalloenzymes that use a unique zinc-mediated hydrolysis mechanism. The classes are further divided into families with a high sequence identity. One of the earliest class A beta-lactamase enzymes to be encountered was TEM-1, which is considered the wild-type variant [
10]. In recent decades many new TEM beta-lactamases have been found worldwide, partly due to new antibiotics that induce adaptations in beta-lactamase enzymes to protect their parent bacteria.
Alternatively, beta-lactamases can be classified based on the class of drugs they provide resistance to [
11]. Narrow Spectrum Beta-Lactamase (NSBL) primarily exhibits penicilinase and limited cephalosporinase activity, which is considered the wild-type substrate profile (Bush-Jacoby classification 2a, 2b). Additional classes of beta-lactamase include Extended Spectrum Beta-Lactamase (ESBL), which confers resistance against third-generation cephalosporins and beyond (2be, 2de, 2e, 2ber) [
12]; inhibitor-resistant beta-lactamase (IRBL) [
13], able to confer resistance in the presence of traditional beta-lactamase inhibitors (2br, 2ber); and more recently, carbapenemase beta-lactamase (CBL) [
14], which today is the main source of concern regarding beta-lactamase-mediated antibiotic resistance (2f, 3a, 3b).
Outside of its clinical significance, beta-lactamase provides an interesting model for studying the role of dynamics on ligand binding. Despite hydrolysis of many substrates, the catalytic mechanism is relatively conserved between class A, C, and D enzymes. The structure of wild-type TEM-1 beta-lactamase is shown in
Figure 1, with the active site highlighted in magenta. The core residues responsible for hydrolysis are Ser70 (the primary catalytic residue), Lys73, Ser130, Asn132, Glu166, Lys234, and Ala237 [
15]. Ser70 is activated by Lys73 and covalently bonds to the beta-lactam ring, while Glu166 activates a catalytic water to deacylate this bond. The result is that the beta-lactam ring is broken, leaving the drug inactive. Other catalytic residues help transport protons for hydrolysis and maintain the structure of the catalytically active site [
16]. As this mechanism is well conserved, mutations outside the catalytic pocket control substrate specificity.
Mutations in TEM beta-lactamase at residues such as Arg164 and Gly238, known as primary mutations, are often associated with novel resistance through the reduction of substrate specificity [
15]. Many secondary mutations can fine-tune resistance spectra, and mutations far from the active site (such as residue Met182) support primary mutations by increasing the thermodynamic stability of the enzyme. Many mutations have no or unknown effects on the function of beta-lactamase [
17], such as Gln39Lys, which distinguishes wild-type TEM-1 and TEM-2 [
18,
19].
In general, the effects of a mutation may induce changes in structure or motions far away from the vicinity of the substitution [
20]. Mutations can impact binding efficacy by changing how different substrates are accommodated into the binding pocket by inhibiting or enhancing functionality related to changes in active site dynamics. Additionally, dynamic allostery can modulate binding affinities after the binding of a secondary ligand or cofactor at another binding site perturbs the dynamics in the active site [
21].
In this work, the dynamics of four class A TEM beta-lactamases are investigated by Molecular Dynamics (MD) simulation. The dataset contains simulations of TEM-1 and TEM-2 (representing NSBL) and TEM-10 and TEM-52 (representing ESBL) both in apo form and bound to four antibiotic substrates: ampicillin, amoxicillin, cefotaxime, and ceftazidime. The functional dynamics critical for substrate recognition are investigated using a model of dynamic allostery developed via machine learning. Dynamic allostery is a form of functional dynamics that deals with how a secondary binding event that is distant from the catalytic pocket alters vibrations [
22]. Mechanistic explanations involving correlated motions are identified through a multivariate comparative analysis of molecular dynamics trajectories using Supervised Projective Learning for Orthogonal Completeness (SPLOC).
In the rest of this paper,
Section 2 describes the MD simulations, the model for predicting dynamic allostery, and how functional dynamics were extracted using the SPLOC machine learning method.
Section 3 describes the results of the MD simulations by analyzing dynamic allostery and functional dynamics.
Section 4 explains the significance of the results and how the results of this work are consistent with other previously reported results in the literature. Finally,
Section 5 concludes by highlighting the successes of the methods used to gain insight into the dynamic mechanisms of substrate binding, and proposes a strategy for novel inhibitors that utilize long-range dynamic allostery.
2. Materials and Methods
2.1. System Preparation
Eight X-ray crystal structures (PDB codes: 1ERM, 1ERO, 1ERQ, 1HTZ, 1JWP, 1LHY, 1XPB, and 3JYI [
23,
24,
25,
26,
27,
28]) of TEM beta-lactamases were collected and downloaded from the Protein Data Bank [
29]. The all-atom structure of each enzyme was extracted by removing any crystallization waters, co-factors, or ligands. In pymol [
30], each crystal structure was mutated to represent TEM-1, and then further mutated so that each of the eight crystal structures had a structure representing each of the four TEM mutants: TEM-1 (wild type), TEM-2 (Gln39Lys), TEM-10 (Arg164Ser/Glu240Lys), and TEM-52 (Glu104Lys/Met182Thr/Gly238Ser). In total this resulted in 32 starting apo structures.
For the ligands, the PDB structures (PDB codes: 3KP3, 6I1E, 4PM5, and 5TWE [
31,
32,
33,
34]) provided an example of each of the ligands in complex with beta-lactamase or another enzyme. The coordinates of the all-atom models for each ligand were extracted and used as a base model. When necessary, bonds were reformed in Avogadro [
35] to represent the intact ligand prior to enzyme hydrolysis. Special care ensured that each molecule was in its bioactive protonation state [
36], as shown in
Figure 2.
Four ligands were docked into the binding pocket of a beta-lactamase structure; ampicillin and amoxicillin were docked with an apo TEM-1 molecule (1XPB) and cefotaxime and ceftazidime were docked with an apo TEM-52 molecule (1HTZ), following the functional substrate profiles of the enzymes. Rigid receptor docking was performed using AutoDock Vina [
37] and DockThor [
38,
39], and the resulting top-scoring poses were examined for physical relevance. In particular, the beta lactam ring was positioned in the active site such that interactions with Ser70 were favorable. In order to obtain the sixteen starting structures for MD simulation, the four previously-obtained docked structures were used as templates to align each of the four beta-lactamase mutants with the four ligand types by superposition; this methodology was chosen because all of the crystal structures had high structural similarity (RMSD <
Å). Each ligand was then observed in order to determine how it naturally reacted to the presence of the beta-lactamase active site and mutations.
2.2. Molecular Dynamics Simulation
Simulations were performed following the protocol described in prior work [
40]. The amber99sb-ildn protein force field [
41] was used to compute molecular interactions at the all-atom level. Each ligand was parameterized using Antechamber [
42] and ACPYPE [
43] was used to convert the ligand parameters to the GROMACS file format, which were then added to the itp and topology files. All simulations were performed using GROMACS [
44]. For each system, the starting structure was placed in a cubic simulation box with 1 nm distance between the edge of the protein and the side of the box, then solvated in TIP3P water solvent. Sodium ions were added in order to neutralize the net electrostatic charge of the system. Energy minimization was applied until the net force on any atom did not exceed 1000 kJ/(mol nm). The system was then coupled to a Berendenson thermostat and equilibrated to 300 K for 1 ns, then to a Parrinello–Rahman barostat and equilibrated to 1 bar for another 1 ns, followed by a 500 ns production run. The first 100 ns in the production run for each mutant was used for further equilibration of the protein.
There are 32 apo systems, consisting of the TEM-1, TEM-2, TEM-10 and TEM-52 mutants, each starting with eight different crystal structures, yielding a combined total of 16 s of dynamics. There are 16 holo systems where ampicillin, amoxicillin, cefotaxime, and ceftazidime are bound to mutants TEM-1, TEM-2, TEM-10, and TEM-52 using crystal structures 1XPB and 1HTZ for broad and extended spectrum antibiotics, respectively. For holo systems, a combined total of 8 s of dynamics was simulated. In all cases, the trajectories were centered in the simulation box, and no nonphysical artifacts appeared when visualizing the motion.
Using Java Essential Dynamics Inspector (JEDi) [
45], each trajectory was aligned to a common reference frame (the first frame of the 1ERM TEM-1 simulation) and the carbon alpha coordinates were extracted into a plain text file that could be read into other analysis programs. Further validation for convergence in the MD trajectories was performed by computing the RMSD and essential dynamics for the simulations. RMSD was calculated for each trajectory according to the equation
where
represents the alpha carbon coordinates of the reference conformation. The RMSD as a function of time shows how the global motion of an enzyme fluctuates during a simulation. The essential dynamics of the apo and holo trajectories were found by constructing and diagonalizing the pooled covariance matrix of all MD simulations. The distributions of the first two principal components were then analyzed using a PDF estimator in MATLAB [
46,
47].
2.3. Dynamic Allostery Model
Dynamic allostery is quantified by observing changes in correlated motions due to coupled binding events. Normal modes with the lowest frequencies represent the largest-scale correlated motions of a protein, which contribute most to the transmission of dynamic allostery. Under the quasi-harmonic approximation, the unperturbed Hessian matrix
is proportional to the inverse of the covariance matrix
[
48]. From an MD trajectory of an apo mutant, the unperturbed Hessian matrix is provided by
where
is the eigenvalue for the
k-th eigenvector,
, of the covariance matrix, and
p is the number of eigenvectors (modes) equal to the state space dimension.
The modes with the lowest eigenvalues of the covariance matrix dominate the sum in Equation (
2), yet they characterize noise with a small amplitude motion. In previous work [
49] a sound approach to regularizing the process of taking the inverse of the covariance matrix was provided. However, the method employed here removes the previously-introduced hyper-parameter and shifts the spectrum of the covariance matrix to a much lesser degree than described earlier [
49]. In order to decorrelate noise in the covariance matrix, all eigenvalues that fall below a minimum value,
, are replaced by this minimum value using the floor operation
prior to using Equation (
2). This minimum value corresponds to the square of the estimated uncertainty in the atomic coordinates. Here, the minimum uncertainty in the atomic coordinates is estimated at
Å.
The effect of ligand binding is modeled by adding pairwise harmonic restraints (springs) in a local spherical region that encapsulates a binding site. Each spring with spring constant
modifies the unperturbed Hessian matrix. The model used here deals only with carbon alpha atoms. The region is defined as connecting any carbon alpha atom within 10 Å from a tagged carbon alpha atom that is a member of the binding site. Under quasi-harmonic approximation, the Hessian matrix is symmetric and mathematically equivalent to a matrix of second-order partial derivatives of an unperturbed effective harmonic potential energy,
. Likewise, the Hessian matrix of the perturbation,
, is obtained from the potential energy of the perturbation,
. The elements of the Hessian matrices for the unperturbed system and the perturbation are provided by
Notice that Equation (
2) allows the unperturbed Hessian matrix elements in Equation (
3) to be calculated without having to specify the unperturbed
. The spring placed between atoms
i and
j represents a harmonic potential provided by
where the position vectors are obtained from the reference structure used for structural alignments. All pairwise harmonic potentials that represent the perturbations are added to obtain
. Then, using Equation (
3), the perturbed Hessian matrix is provided by
. In this work,
kcal/(mol Å
2); however, as previously shown [
49], the results are independent of this value in the linear response limit, where
is small. These perturbations to the Hessian matrix change the normal mode frequencies, which are obtained using first order perturbation theory. Following the theory of dynamic allostery set out by Cooper and Dryden [
22], the change in the binding free energy (
) of an endogenous ligand at the active site induced by the binding of an effector molecule can be calculated from the resulting changes in normal mode frequencies.
The
k-th normal mode has three frequencies of interest: the apo frequency,
, the frequency with a single endogenous ligand bound at the active site,
, and the frequency with both an endogenous ligand and an effector molecule bound,
. The change in binding free energy can be estimated as
where two sets of external springs are used to model the two binding events. A sum over all modes unrelated to rigid body translations and rotations is performed after removing the first six lowest frequency modes. In practice, a sum over several modes is sufficient to identify putative dynamic allostery sites;
quickly converges, as high-frequency modes yield little contributions. It should be noted that because this is a coarse-grained model that uses only carbon alpha atoms, and because the perturbation is small, the
are generally much smaller than what would be expected physically; however, it is the shape of the variations as a propensity which is being considered here. This method allows us to perform a single MD simulation on an apo structure and then perform a perturbation scan that calculates changes in active site affinity as different regions of the protein experience the simulated effector molecule binding. When local regions of the effector and endogenous ligand overlap, fewer springs are used, avoiding any double counting of springs. In this work, there is no interesting signal caused by this local interference, and this special case is ignored in all subsequent discussions that focus on distal effects.
2.4. Comparative Multivariate Analysis
Through a comparative and exploratory analysis, the functional mechanisms within beta-lactamase were extracted by supervised projection pursuit using Supervised Projective Learning for Orthogonal Completeness (SPLOC) [
50,
51,
52]. The state space consisted of the
coordinates of the carbon alpha atoms representing 263 residues, yielding
degrees of freedom. As only substitution mutations were considered, all proteins had an identical state space. All conformations were structurally aligned using JEDi [
45]. The input to SPLOC was a collection of data packets, each quantified by the mean position of the carbon alpha atoms and the covariance matrix describing atomic fluctuations about the mean, respectively denoted as
and
, where
r is a data packet index. Multiple data packets (called replicas) were constructed by subsampling conformations in a trajectory. The output of SPLOC provided
p modes defining a complete orthornormal basis set of vectors, where the
k-th mode is denoted as
. It is worth mentioning that the PCA and SPLOC modes are unrelated, and are used in easily recognizable contexts in this work. In
Section 4.1, a quantitative comparison of these two modes is shown and discussed.
A feature space was constructed from the mean and standard deviation of the data projected into each mode, such that
=
·
and
=
·
·
. The ordered pair
defines a point for the
r-th data packet within a mode feature space plane (MFSP) for mode
k. As there are two emergent properties (mean and standard deviation) per mode, the feature space dimension is
. Considering all data packets for functional and nonfunctional systems, statistical measures were defined per MFSP in order to arrive at a single objective function, called mode efficacy,
. The SPLOC process is summarized qualitatively in the next section, with complete details provided in [
50].
For the k-th mode, there are three measures: (1) selection power, ; (2) consensus power, ; and (3) quality of clustering, . Selection is a measure of signal-to-noise, consensus is a measure of statistical consistency across the data (including replicas), and cluster quality is a scale-invariant geometrical measure of how data packets cluster. All these measures are invariant upon the exchange of class labels, which are called functional and nonfunctional in this work. A decision tree involving was used to classify the k-th mode as discriminant, undetermined, or indifferent, respectively called d-, u-, or i-mode. If any one of the three measures fell below a minimum threshold, the mode was classified as u-mode. Cross-validation was conducted within SPLOC through the consensus measure as the basis set was optimized. Thus, if a mode was not identified as a u-mode, only two cases were possible; (d-mode, i-mode) could be identified when functional and nonfunctional data packets were linearly separable and statistically indistinguishable within an MFSP. Each mode was assigned a rectifying adaptive nonlinear unit (RANU) in order to bifurcate mode identity towards either d-mode or i-mode. As the parameters of a RANU are data-driven based on the properties of its mode projection, a heterogeneous network of p RANUs can determines the mode efficacy, .
The total efficacy of a proposed orthonormal basis set is provided by . As the total efficacy of the network of RANU-type activation functions is linearly separable, the subspaces of the modes can be optimized independently. As such, a generalized Jacobi method was implemented, in which pairs of modes were rotated in different pair combinations through a random search process. During this process, mode efficacy can only increase monotonically, and the algorithm continues to rotate the basis vectors until the total efficacy of the network saturates. As each mode provides a means of determining emergent properties (as means and standard deviation) along the projection, and no information is lost due to completeness, this process provides interpretability of functional dynamics.
2.5. Comparative Studies and Data Packet Generation
In this study, two main categories were considered in the course of elucidating the functional dynamics of beta-lactamase by quantifying similarities and differences in dynamics using SPLOC. Each category represents functional versus nonfunctional classes. The first category provides a protein-centric comparison, where differences and similarities in protein dynamics depend on NSBL versus ESBL. The second category provides a ligand-centric comparison, where differences and similarities in protein dynamics depend on broad versus extended spectrum antibiotics. While the protein-centric perspective takes into account all MD trajectories, the ligand-centric perspective is only relevant for proteins when ligands are present. The data packet statistics used in subsequent analyses are summarized in
Table 1.
A bootstrap approach was used in which all trajectories that shared a common characteristic were combined. From a given trajectory pool of trajectories consisting of M frames each, a subsample of m frames were randomly drawn and concatenated together. Through this method, replica data packets for functional and nonfunctional classes and for both categories were constructed. All data packets had samples and an observations per variable (OPV) of . When a functional or nonfunctional class consisted of different subclasses (e.g., TEM-1 and TEM-2 are both NSBL), replica data packets were constructed exclusively with MD frames from each subclass. This subsampling process allowed for data packets that provided unbiased representation of all data while at the same time characterizing any fluctuations that might exist.
2.6. Extracting Functional Motions
A common measure obtained from MD simulations is the root mean square fluctuations (RMSF). It is common to project out small amplitude motions using a small number of PCA-modes (with largest variances) compared to p. In such cases, the actual RMSF and the PCA-projected RMSF will generally be very similar, which is the desired result because it is of interest to capture most of the dynamical motion using a relatively small number of PCA-modes. Using the same projection concept, the RMSF measure can be generalized using the d-modes and i-modes obtained from SPLOC. Specifically, it is possible to observe distinctly different conserved motions between two systems by using (d-modes, i-modes) to project out unwanted motions. When comparing functional versus nonfunctional systems, the observed differences in motion may not be related to function, and conserved motions can be important for function.
By knowing which parts of motion differ between systems and which parts are similar, it is possible to glean insights into which motions support a function and why. A hypothesis involving functional dynamics can then be refined by inductive reasoning as additional comparisons are made. After being projected to be within the (discriminant, indifferent) subspace, dynamical fluctuations are called (discriminant RMSF, indifferent RMSF). The discriminant RMSF (dRMSF) and indifferent RMSF (iRMSF) are calculated using the standard procedure for calculating RMSF, except that the first step is to respectively construct a projection operator from the d-modes and i-modes. The second step then projects out all motions outside the subspace of interest. Where there is a data matrix
having the dimensions
and where
n is the number of frames, the projected data matrices for the discriminant subspace,
, and for the indifferent subspace,
, are provided by
4. Discussion
TEM beta-lactamase is dynamically constrained on a global scale while exhibiting variations in loop motions sensitive to mutations and ligand binding. In agreement with other recent studies [
25,
57], the present study further reveals the importance of loop dynamics in controlling beta-lactamase function. In particular, the H9-H10 loop facilitates structural stability, and the H10-H11 loop stabilizes the ligand in the active site for catalysis. Moreover, different long-range effects between the various mutants provide a means for adaptability in beta-lactamase, promoting further control of substrate specificity.
4.1. Comparison of PCA and SPLOC Modes
Dimension reduction using PCA to describe essential dynamics has long been used as a standard method for extracting functional motions from molecular dynamics simulation data [
58,
59]. Despite its popularity, PCA captures functional dynamics that correspond to the largest amplitude motions described by collective modes with the greatest variance. However, this approach is actually based on an assumption that is not universally true; function can involve motions in localized regions that are not extended. The usual practice of throwing out the low-variance PCA modes thus risks missing the functionally important motions altogether. In the analysis provided herein, PCA does not successfully cluster known functional data classes, as shown in
Figure 3. In this case, a supervised discriminant analysis method is more appropriate. Due to the data-driven efficacy score of SPLOC, motions can be labeled as discriminant or indifferent regardless of the amplitude of motion. It is thus interesting to determine how much information contained in orthogonal subspaces of SPLOC can be reconstructed with PCA modes ordered from the most to the least variance.
In
Figure 9, the SPLOC modes from the results in
Section 3.2.1 are compared with the global PCA modes obtained in
Section 3.1 in order to determine how the information extracted by SPLOC is distributed across the PCA basis. Each basis vector in a given subspace was reconstructed using the top
n PCA modes and averaged over all modes in the SPLOC subspace for both SPLOC subspaces (comprising only d-modes, i-modes, and u-modes) and the full SPLOC basis set (
Figure 9a–d, respectively). The percentage of the reconstructed mode is the cumulative overlap (CO) between the SPLOC basis vector and the PCA subspace, provided by
Figure 9a shows that, on average, the d-modes were not even 20 percent reconstructed using the top 400 PCA modes. This means that information about the differences in functional dynamics between two beta lactamase mutants is lost due to dimension reduction. The i-mode reconstruction shown in
Figure 9b increases linearly in CO until 500-600 PCA modes are used. The u-mode reconstruction shown in
Figure 9c initially increases rapidly, especially compared to TEM-1 and TEM-2. Finally, if we use the complete set of SPLOC modes, a perfect linear relationship emerges with a slope of
, as seen in
Figure 9d. Therefore, it makes sense to compare the reconstruction rate with the reference line shown in
Figure 9d.
The i-modes describe the functional similarity between mutant pairs that is most conserved. This motion is expected to be associated with the dominant enzyme dynamics through large-scale motions, which are mainly determined by the structure of the protein. Due to the enzymes’ high degree of structural similarity, the CO of i-modes expressed in terms of PCA modes increases faster than the reference line until the dimension of the i-mode subspace is reached. The motions described by u-modes are not statistically different or the same between a pair of mutants, and thus capture outlier motions that may have been expressed differently between simulations of the same enzyme. Being sensitive to outliers, PCA often captures this high variation in the top modes, which explains why the u-modes initially increase rapidly in CO until the dimension of the subspace for u-modes is reached. Finally,
Figure 9a clearly shows the statistically significant trend that d-modes contain very different information from the most important PCA modes. This analysis suggests that the SPLOC analysis provides a means of capturing different aspects of protein motions that are important for function that PCA cannot provide.
4.2. The H9-H10 Linker Loop Stabilizes Beta-Lactamase against Structural Changes
The MD/SPLOC results suggest that the H9-H10 loop plays a key stabilizing role, which is in agreement with previous work [
57] postulating that loops transmit allosteric signals through the H9-H10 loop, resulting in the observed high fluctuations. This loop (shown in
Figure 10) connects two helices that are almost perpendicular to each other. The H9 helix at the beginning of the loop shields the back of the H2 helix and the catalytic serine 70 from solvent. This loop sterically counters a resultant push back of the H2 loop when it binds to an antibiotic. From the simulation data, Leu194 braces Leu51 on the nearby
-
connector loop. Mutations such as Leu51Pro, Leu51Phe, or Leu51Ile (e.g., in TEM-60) have been shown to increase the stability of TEM-1 beta-lactamase [
17,
60], which supports this hypothesis. Additionally, residue Lys192 has been observed to form an intermittent interaction with Gly196, and disruption of this favorable interaction should provide the means for greater flexibility.
The results in
Section 3.2.1 show that both specificity-modulating and functionally-silent mutations can induce a dynamic response in the H9-H10 loop. Notably, all three mutation sets were able to induce at least a small change in the top of the omega loop, which directly connects to the H9 loop. The smallest change in the omega loop dynamics was observed in TEM-2. In
Section 3.2.2, it was shown that this loop remains flexible when NSBL enzymes bind to a drug, whereas it becomes more rigid when ESBL enzymes bind to a drug. This is in response to the flexibility induced in the omega loop and active site by ESBL-conferring mutations. The active site takes the more flexible form when it accommodates a ligand, as it takes up slack from the H9-H10 loop. Furthermore, this loop is an intermediate connector between the omega loop and the H10-H11 loop, which forms the top and bottom of the catalytic pocket. These three loops are all key locations where the ligand can anchor itself to the enzyme and where flexibility is altered by mutation or effector binding.
4.3. The Role of the H10-H11 Loop in Stabilizing Ligands
Extended spectrum ligands such as cefotaxime and ceftazidime are typically larger than their broad spectrum counterparts, and are therefore harder to accommodate in the active site. Analysis comparing all three mutation sets shows a change in motion in the H10-H11 loop for TEM-10 and TEM-52 (
Figure 5c,d). A comparison of apo and holo enzymes shows that ESBL enzymes lose flexibility here, while NSBL enzymes exhibit similar dynamics before and after the ligand enters the pocket. Comparing the dynamics in the holo simulations (
Figure 8a) shows that although the loop motions change upon binding, they remain relatively mobile. As shown in
Figure 7, the carboxyl group of cefotaxime and ceftazidime points toward the H10-H11 loop, as does a carbonyl group in cefotaxime and pyridine group in ceftazidime. These functional groups form transient contacts with Val216, which stabilizes the ligand, beyond the primary anchors of the ligand at the
-
turn and the omega loop. The conformational freedom of the H10-H11 loop enables the ligand to move towards the active site and positions Val216 near the carboxyl group of the ligand. This underscores the importance for ligand stability of having Val216 in the correct conformation, which explains why motions drastically change upon binding to an extended-spectrum ligand.
Contact with the H10-H11 loop does not appear necessary for ligand hydrolysis, rather promoting the capture, chaperoning, and stabilization of the protein ligand complex. For example, cefotaxime was able to resist being pulled into the solvent, although this H10-H11 loop contact was only intermittent during the simulation. In the TEM-1 holo simulation, the cefotaxime was able to reposition itself in the pocket to interact with the H10-H11 loop while remaining anchored to the omega loop. Ampicillin and amoxicillin can reside in the binding pocket, despite not being able to reach this loop. This secondary anchor forces the extended spectrum ligands to stretch across the active site, putting them in a more favorable position for catalysis. Although not required for function, these results show that this loop plays an important role in expanding substrate specificity for TEM beta-lactamase.
4.4. Dynamic Allostery as a Mechanism of Inhibition
The dynamic allostery perturbation scans of the four TEM beta-lactamase mutants suggest that substrate binding is coupled to distal loop regions. In combination with previous results from the literature, a mechanism involving the known allosteric binding pocket for reducing antibiotic affinity is examined here in detail.
According to the scan, binding an effector to H11 is likely to enhance the binding affinity of antibiotics at the active site. In prior experimental work, destabilization of this region was shown to inhibit activity [
53,
61]; here, a stabilizing perturbation was found to enhance activity. The mechanism for H10-H11 ligand stabilization above relies on Val216 taking the appropriate conformation to reach the ligand. A plausible explanation of this allosteric response would be that removing flexibility in this region induces the H10-H11 loop into a favorable conformation for ligand binding. By contrast, forcing H11 and H12 away from each other is likely to move the loop away from its favorable conformation, thereby reducing the enzyme’s activity.
Binding an effector near residues 227-229 and nearby to residues 254-256 may inhibit binding of ligands at the active site. This mechanism is directly related to the primary allosteric region of beta-lactamase, as it houses the pi-stacking Pro226-Trp252-Pro229 triad, which is thought to stabilize this region. Strain in this region due to effector binding increases the likelihood of disrupting this interaction and further disrupting the structure of the active site, including the H10-H11 loop via long-range dynamic signal propagation. SPLOC identified large conformational fluctuations in this loop, enabling the enzyme to adapt the shape of the binding pocket to various beta-lactam molecules.
Another allosteric _target area was found at residues 104-106 and 111-115. This allosteric signal is strong for all mutants; however, in TEM-52 this reflects the location of the Glu104Lys mutation. Tyr105 on this loop has been previously identified as an important residue via mutation studies [
62] and NMR experiments [
63]. Mutation study showed that this position prefers an aromatic residue to cap the active site and prevent steric interactions with the ligand in the pocket. Moreover, NMR showed that Tyr105 dynamics are correlated with both the motions of other local active site residues and with distant residues. These residues include Asn132, a catalytic residue found on the loop just below Tyr105, and Lys234, found on the
-
loop on the opposite side of the active site. The same study found a correlation between the motions of Tyr105 and Val216, suggesting a link between the allosteric mechanism described above and the H3-H4 loop. This work further suggests that this allosteric signal can be controlled to inhibit the enzyme by binding events at either H3 or H4, which form the beginning and end of this loop.
It appears that the H9-H10 loop resists structural changes to beta-lactamase, although it was not identified as allosteric in the perturbation scan. While loop dynamics may not directly couple to active site residues, they are influenced by distant residues. Therefore, the H9-H10 loop may be a viable inhibitor _target, as its flexibility is well conserved across mutation sets and ligand binding differences in TEM beta-lactamase. As an alternative approach to dynamic allostery, SPLOC results can be used to gain insights into how to optimize dynamic properties that improve enzyme binding. Binding a ligand to or near this loop may force the H9-H10 loop into a conformation with more rigidity, preventing the loop from maintaining stability through its dynamic mechanism. The region between H10 and the
-
turn, shown in
Figure 11, forms a pocket roughly
Å long by
Å wide.