Introduction

Pseudomonas aeruginosa is a ubiquitous Gram-negative bacterium from the family Pseudomonadaceae (Pang et al., 2019). This pathogen is commonly found in hospitals and other healthcare settings, where it can cause infections in people who are immunocompromised or have chronic conditions (Kerr and Snelling, 2009). Pseudomonal infections are associated with high morbidity and mortality in many groups, including patients with cystic fibrosis, pneumonia, and chronic obstructive pulmonary disease (COPD) (Curran et al., 2018; Jurado-Martin et al., 2021; Malhotra et al., 2019; Reynolds and Kollef, 2021). β-lactam antibiotics, characterized by the presence of a β-lactam ring in their chemical structure, are often the first-line of treatment for bacterial infections as they tend to have fewer side effects and are less toxic than other antibiotics (Mora-Ochomogo and Lohans, 2021). Mechanistically, β-lactams work by inhibiting the synthesis of the bacterial cell wall, which is necessary for the survival and growth of bacteria (Lima et al., 2020). β-lactam antibiotics are usually effective against a wide range of bacteria, but can lose their effectiveness if the bacteria develop resistance to them. P. aeruginosa is known for its ability to develop resistance to multiple classes of antimicrobial drugs (Horcajada et al., 2019; Pang et al., 2019; Spagnolo et al., 2021). Therefore, the greatest challenge to eradicating P. aeruginosa infections are multidrug-resistant (MDR) and extensively drug-resistant (XDR) isolates. The World Health Organization has identified P. aeruginosa as a top-priority pathogen for research and development of new antibiotics due to its ability to cause serious infections and its increasing resistance to currently available treatment options (Tacconelli et al., 2018).

The production of antibiotic-inactivating enzymes is one of the major mechanisms of intrinsic resistance in bacteria (Munita and Arias, 2016). P. aeruginosa can express a class C β-lactamase, named Pseudomonas-derived cephalosporinase (PDC), which is an antibiotic-inactivating enzyme (Colque et al., 2022). PDC-3 is a serine β-lactamase that can inactivate a broad range of β-lactam antibiotics, including penicillins, cephalosporins, monobactams, and carbapenems, by breaking the amide bond of the β-lactam ring through a catalytic site serine (Figure 1A and B) (Pang et al., 2019; Tripathi and Nair, 2013, 2016). Cephalosporins are known to be highly susceptible to PDC-3 inactivation (Barnes et al., 2018). The active site of PDC-3 is located at the intersection of the enzyme’s α-helical and α/β domains (Figure 1D). The active site can be further divided into two distinct regions: the R1 site and the R2 site. These regions are defined by the specific binding interactions they facilitate with the R1 and R2 side chains of the cephalosporins, respectively (Figure 1A). The R1 site is surrounded by the Ω-loop, while the R2 site is encased by the R2-loop, which is comprises of the α helix H-10. The Ω-loop and the R2-loop are located at opposite ends of the active site, with the catalytic serine residue positioned in the middle (Jacoby, 2009). However, the Ω-loop and R2-loop, are particularly prone to amino acid substitutions, insertions, and deletions that expand the active site and accommodate larger R1 and R2 groups of the cephalosporins (Nordmann and Mammeri, 2007). These modifications have been observed to enhance the capacity of the enzyme to hydrolyze a wider range of β-lactam antibiotics (Barnes et al., 2018). The evolution of PDC-3 β-lactamase and its amino acid variants, which often result in enhanced catalytic activity and expanded-spectrum of cephalosporin hydrolysis, has garnered considerable interest in scientific research (Ruedas-Lopez et al., 2022). As previously reported, several PDC-3 Ω-loop variants, including V211A, V211G, G214A, G214R, E219A, E219G, E219K, Y221A and Y221H were found in highly drug-resistant P. aeruginosa clinical isolates (Barnes et al., 2018). All residues in this study are numbered based on the Structural alignment-based numbering of class C β-lactamase scheme (SANC) (Mack et al., 2020).

Figure 1.

Molecular Dynamics (MD) simulations provide valuable insights into the time-evolved dynamic behavior of biomolecules such as proteins (Hollingsworth and Dror, 2018). Among various enhanced sampling methods, Adaptive Bandit Molecular Dynamics (AB-MD) stands out as a powerful, reinforcement learning (RL)-based adaptive sampling strategy (Pérez et al., 2020). Compared to a single long equilibrium MD simulation (which can become trapped in a metastable state) or bias-based enhanced sampling techniques like accelerated MD (aMD) and Gaussian accelerated MD (GaMD) that add boost potentials to smooth the energy landscape and facilitate barrier crossing, AB-MD employs a reinforcement learning-inspired multi-armed bandit to adaptively guide multiple short, unbiased trajectories toward regions of conformational space (Bhattarai and Miao, 2018). In a multi-armed bandit, each ‘arm’ yields a reward, and an agent must balance exploration (trying less-sampled options) and exploitation (focusing on the best-known options). AB-MD applies this idea to MD by treating different regions of conformational space (or states) as the bandit ‘arms’. At each iteration, the algorithm decides which state to launch a new MD simulation from, aiming to maximize sampling efficiency while still obtaining an accurate, unbiased representation of the system’s equilibrium behavior.

In this study, the AB-MD approach was employed to expore the conformational landscape of PDC-3 and its variants at the atomistic level. Additionally, constant pH MD simulations were performed to examine the protonation-state behavior of key residues in the catalytic site. By analyzing the conformational ensembles and kinetics of PDC-3 through these simulations, the project aimed to uncover underlying mechanisms governing the function of PDC-3 and its variants. Understanding the molecular mechanisms underlying PDC-3 function and the development of resistance is of paramount importance in combating P. aeruginosa infections. This knowledge can aid in the development of more effective treatments to combat these bacteria.

Results and Discussion

Amino acid substitutions change the lexibility of Ω-loops and R2-loops

Both root-mean-square deviation (RMSD) and root-mean-square fluctuations (RMSF) analysis provide insights into the dynamic behavior and structural differences of biomolecules (Martinez, 2015; Prabantu et al., 2022). Firstly, the structural stability of the overall conformation of the wild-type PDC-3 and its variants was investigated using pairwise RMSD analysis. The results of the RMSD analysis reveal that the wild-type PDC-3 displays a relatively low degree of structural variability, indicating that the protein maintains a relatively consistent overall conformation over time (Figure 2A and figure Supplement 1). Similarly, the V211G, G214A, and Y221H variants also exhibit low RMSD values, which suggests that these structures are less flexible. In contrast, the V211A and E219G variants exhibit the highest RMSD values among the set of structures, indicating a high level of structural variability. This implies that these substitutions lead to increased conformational fluctuations in the protein over time. The G214R, E219A, E219K and Y221A variants exhibit RMSD values that are intermediate between the wild-type and the most flexible variants, indicating that these amino acid substitutions have a moderate effect on the structural stability of the protein’s conformation, i.e., not as significant as the V211A and E219G substitutions.

Figure 2.

To identify regions that contribute the most to the conformational changes in the wild-type PDC-3 and its variants, the RMSF values of Cα atoms were calculated. High RMSF values indicate a high degree of flexibility or mobility for the corresponding atoms, while low RMSF values indicate a significant degree of rigidity (Boot et al., 2011). The wild-type PDC-3 and the G214A, G214R, E214G, and Y221A variants exhibit high flexibility in their Ω-loops, as evidenced by the relatively large per-residue RMSF values observed (approximately 4 Å). In contrast, the V211A, V211G, E219K, and Y221H variants display more constrained conformations in the Ω-loops. Notably, the V211A and V211G variants demonstrate the highest stability in the Ω-loop, with average RMSF values around 1.5 Å, whereas the E219K and Y221H variants exhibit intermediate flexibility, with RMSF values averaging between 2 and 2.5 Å. In terms of the R2-loop, wild-type PDC-3 displays a relatively low degree of flexibility while all variants exhibit an increase in structural flexibility. This suggests that these substitutions have a significant impact on the stability of the R2-loops, potentially affecting enzyme function. A detailed observation of the individual variants reveals that the V211A variant exhibits a particularly high degree of flexibility, as evidenced by the comparatively higher RMSF values, followed by the E219G variant. On the other hand, the Y221A and Y221H variants exhibit a relatively lower degree of flexibility, as inferred by the lower RMSF values observed (Figure 2B, C and figure Supplement 2). Therefore, the flexibility of these proteins is mainly caused by the fluctuations in the Ω-loops and R2-loop.

The importance of Ω-loops and R2-loops in class C β-lactamases has been previously confirmed (Philippon et al., 2022). These loops play a crucial role in the binding and activity of the class C β-lactamases (Philippon et al., 2022). Specifically, residues V211 and Y221 within the Ω-loops have been identified to engage in hydrophobic interactions with the R1 side chains of cephalosporins. The substitution of V211A has been reported to be associated with acquired resistance to cefepime or cefpirome (Rodriguez-Martinez et al., 2010). Additionally, the characteristic aminothiazole ring found in most third-generation cephalosporins interacts with Y221 in an edge-to-face manner, which represents typical quadrupole-quadrupole interactions. However, Y221 can sometimes create steric clashes that prevent ligands from entering the binding sites (Barnes et al., 2018; Powers and Shoichet, 2002). Deletion of Y221 has been observed to broaden substrate specificity and confer resistance to ceftazidime-avibactam (Lahiri et al., 2015). Moreover, the expanded Ω-loop of P99, another member of class C β-lactamases, exhibits conformational flexibility that may facilitate the hydrolysis of oxyimino β-lactams by making the acyl intermediate more accessible to attack by water (Crichlow et al., 1999). In terms of the R2-loops, it has been observed that the N289 (N287 in PDC-3) forms hydrogen-bonding interactions with the C4 carboxylate hydrogen bonds directly in the AmpC/13 (moxalactam) complex (Crichlow et al., 2001). Furthermore, the T289, A292, and L293 residues within the R2-loops of class C β-lactamases have been found to exhibit hydrophobic contacts with the dimethyl group in the R2 chains of cephalosporins (Powers and Shoichet, 2002). Additional research suggests that the removal of the R2 group in cephalosporins occurs, while the R1 group remains intact (Chaudhry et al., 2019; Perez-Inestrosa et al., 2005). This observation indicates that the high flexibility of the R2-loops could be a crucial factor in stabilizing substrates during both the acylation and deacylation steps simultaneously. Overall, the flexibility or mobility of the Ω-loops and R2-loops allows the PDC-3 active site cavity to adopt different size and shapes, thus affecting the binding of different β-lactams and allowing for extended-spectrum activity of some class C β lactamases.

E219K and Y221A mutations facilitate the proton transfer

The substitution of E219 to K219 or Y221 to A221 leads to conformational changes in the active site. Notably, crucial residues S64 and K67 undergo clear conformational alterations. In the E219K substitution, nearly half of the structures display gauche(-) (60°) conformations for S64, which is uncommon in the wild-type PDC-3 and other variants. The Y221A variant also shows limited gauche(-) (60°) conformations for S64. Regarding K67, gauche(+) (-60°) conformations dominate in the wild-type PDC-3 and other variants, but almost 50% of the structures exhibit trans (180°) conformations with the E219K variant (Figure 3A). The conformational changes of S64 and K67 can impact the formation of hydrogen bonds (Table 1 and Figure 3 figure Supplement 1). Specifically, the hydrogen bond formation between K67(NZ) and S64(OG) is observed in most structures of wild-type PDC-3 and variants, but it is broken in most structures with the E219K and Y221A mutants (Figure 3B). The fractions of formation for this hydrogen bond in E219K and Y221A mutants are only 30.36% and 40.05%, respectively. In addition, the formation of K67(NZ)-N152(OD1) and K67(NZ)-G220(O) interactions are also greatly decreased in the E219K and Y221A variants. This could potentially have implications for the catalytic activity of the enzyme, as S64 is a catalytic residue involved in the acylation step of the β-lactamase. K67 is believed to act as a general base in the acylation step of the β-lactamase catalytic mechanism, abstracting a proton from the hydroxyl group of S64, which in turn facilitates the nucleophilic attack of the β-lactam ring (Tripathi and Nair, 2013, 2016). Therefore, K67 is thought to toggle between protonated and deprotonated states to facilitate proton transfer in the catalytic cycle (Figure 1B). However, when K67 is involved in persistent and energetically favored hydrogen bonding interactions with S64, N152, and G220, these stable interactions can lock it in the protonated state. This conformational arrest makes it less available to undergo the protonation state toggling required for catalysis.

Figure 3.

Table 1.

The protonation state (pK_a) of K67 in the enzyme is therefore critical: a lowered pK_a could allow K67 to exist as a neutral ‘general base’ at physiological pH, ready to accept a proton in catalysis (Chen et al., 2009). Indeed, analogies from related enzymes suggest that catalytic lysines often have depressed pK_a values (e.g. K47 in PBP5 and K73 in TEM-1) to enable catalytic function (Golemi-Kotra et al., 2004; Zhang et al., 2007; Shi et al., 2008; Meroueh et al., 2005b). However, directly measuring or computing the pK_a of a buried lysine in a large enzyme is challenging. Constant pH molecular dynamics simulations (CpHMD) provide a powerful in silico approach to estimate pK_a by allowing protonation states to fluctuate according to a chosen pH (Kim et al., 2015). Here, we employ CpHMD to compute the pK_a of K67 in wild-type PDC-3 and compare it with the E219K and Y221A variants. Titration curves generated from pH-replicated simulations are analyzed to extract K67 pK_a values. Our results indicate that, in wild-type PDC-3, K67 exhibits a pK_a in the range of approximately 8.50-8.79 (Figure 3C and figure Supplement 2). By contrast, the E219K mutation dramatically reduces the pK_a of K67 to approximately 6.32-6.71, causing K67 to be predominantly deprotonated (neutral) at physiological pH. This deprotonated form is conducive to K67 functioning as a ‘general base’ readily accepting a proton and thereby facilitating the nucleophilic attack of the S64 hydroxyl group on the β-lactam ring (Tripathi and Nair, 2013, 2016). The Y221A mutation also shifts pK_a of K67 down (7.60-8.06), though to a lesser degree than E219K, consistent with the hydrogen-bond occupancies observed in our simulations. As previously noted, the E219K and Y221A mutations weakens the tridentate hydrogen-bond networks (K67(NZ)-S64(OG), K67(NZ)-N152(OD1), K67(NZ)-G220(O)), thereby enabling K67 to more flexibly adjust both its conformation and its protonation state-a change that promotes efficient proton transfer (Table 1). Experimentally, these mutations also confer increased sensitivity to cephalosporin antibiotics, which aligns with the conformational and protonation-state shifts observed in the simulations (Barnes et al., 2018). Collectively, these findings reveal that the E219K and Y221A substitutions disrupt the tridentate hydrogen-bond network, which lowers the pK_a of K67 and enhances its ability to act as a general base. This elevated proton-transfer efficiency, in turn, improves the enzyme’s catalytic performance.

Substitutions enlarge the active-site pocket to accommodate bulkier R1 and R2 groups of β-lactams

From Table 1, we can observe that the wild-type PDC-3 forms a hydrogen bond K67(NZ)-G220(O) in 63.20% of simulations. However, in the variants G214R, E219K, E219G, and Y221A, the percentages of forming this hydrogen bond are greatly decreased, with fractions of formation ranging from 24.11% in G214R to 31.70% in E219G (Figure 4A, B and figure Supplement 1). In the Y150(N)-A292(O) interaction, the wild-type PDC-3 forms this hydrogen bond in 72.94% of simulations. The V211A, V211G, E219A, E219G, E219K, Y221A, and Y221H variants have significantly lower percentages of formation, ranging from 15.67% to 29.36%. Therefore, substitution of the V211, E219, and Y221 residues will affect this hydrogen bond considerably. In addition, approximately 55.09% of the structures form the hydrogen bond N287(ND2)-N314(OD1) in wild-type PDC-3. However, in the variant structures, most do not form this hydrogen bond, only a few do with a range of 18.72% to 36.42%. Residues G220 is positioned in the Ω-loop, while A292, N287, and N314 are located in the R2-loop. Therefore, the disruption of the K67(NZ)-G220(O) interaction in PDC-3 may lead to the expansion of the R1 site, allowing for substrates with larger R1 side chains such as ceftolozane and ceftazidime to be accommodated. Similarly, the loss of the Y150(N)-A292(O) and N287(ND2)-N314(OD1) interactions could potentially enlarge the R2 site of PDC-3, providing additional space for R2 side chains of cephalosporins.

Figure 4.

The confirmation of these observations can be supported by analyzing the pocket volumes and solvent accessible surface areas (Table 2 and Figure 4 figure Supplement 2). In the simulations with the Y221A, E219K, and E219G variants, it was found that all three interactions were disrupted in the majority of structures. On the other hand, the G214A variant displayed relatively lower fractions of all three interactions compared to other variants. Consequently, nearly all variants exhibit significantly larger active-site volumes than the wild-type (WT) (all p < 0.0001). The most pronounced expansions occur in Y221A (1381.32 ± 302.32 ų), E219K (1278.62 ± 293.01 ų), and E219G (1211.74 ± 250.02 ų). By contrast, G214A (1075.46 ± 259.80 ų; p = 0.0373) only slightly exceeds WT (1071.81 ± 157.55 ų). The solvent-accessible surface area (SASA) shows a similar pattern, with each variant differing significantly from WT (p < 0.0001). Notably, Y221A (156.89 ± 2.95 Ų) and G214R (156.65 ± 2.84 Ų) exhibit some of the largest increases in SASA, further suggesting an expanded active-site region that may accommodate bulkier R1 and R2 side chains in β-lactam substrates. As a result, we infer that the R1 site of the G214R, E219G, E219K, and Y221A variants should undergo a noticeable increase in size. Additionally, the R2 site is expected to experience an increase volume in all variants except G214A and G214R.

Table 2.

A previous study by Barnes et al., has shown that G214R, E219G, E219K, and Y221A variants exhibit greater resistance to ceftolozane and ceftazidime (Barnes et al., 2018). In contrast, their resistance to cefotaxime and cefepime, which have smaller R1 side chains, was not significantly affected (Barnes et al., 2018). Additionally, metadynamics studies have also indicated that there will be a lower free energy barrier for the acylation reaction when this interaction between A220 (the equivalent position of G220 in PDC-3) and K67 is broken in the complex of CBL and aztreonam (Tripathi and Nair, 2013). However, only E219G, E219K, and Y221A demonstrate an obvious increase in resistance to ceftolozane while other variants only show minor changes. This difference in resistance might be attributed to the relatively smaller R1 site. The other variants had more structures capable of forming interactions between K67(NZ) and G220(O), which could potentially reduce the impact on their resistance to ceftolozane.

Markov State Modeling Characterizes Key Hydrogen Bond Transitions

The utilization of Markov state models (MSMs) enabled the analysis of long-term conformational alterations of wild-type PDC-3 and its variants by filtering out local fluctuations related to thermal motion and focusing on underlying conformational transformations (Bowman et al., 2009; Husic and Pande, 2018; Scherer et al., 2015; Trendelkamp-Schroer and Noe, 2013). Previous analyses demonstrated that, in addition to the catalytic site, the most significant structural changes occur in the Ω- and R2-loops; consequently, hydrogen bonds and salt bridges in those loops and in the catalytic site were identified for MSM construction. Distances for all relevant interactions were computed in (i) the active motifs (S⁶⁴XXK⁶⁷, Y¹⁵⁰SN¹⁵², K³¹⁵TG³¹⁷), (ii) the Ω-loop (residues G183–S226), and (iii) the R2-loop (residues L280–Q310). To establish the correlation between structural dynamics and active-site pockets, correlation coefficients were computed between the distances of these interactions and the volume of the active-site pocket. A correlation coefficient exceeding 0.3 or falling below -0.3 indicates a positive or negative relationship, respectively. Only the distances of salt bridges and hydrogen bonds that exhibited a positive or negative relationship with the volume of the active site pockets were selected as features to construct the MSMs. This resulted in the selection of 8 salt bridges and 24 hydrogen bonds (Table 3).

Table 3.

From the metastable state results, we observe that E219K adopts a highly stable conformation in which all the tridentate hydrogen-bonding interactions (K67(NZ)-S64(OG), K67(NZ)-N152(OD1) and K67(NZ)-G220(O) mentioned above are broken (Figure 5). Notably, this ‘fully broken’ configuration appears only in E219K and Y221A variants (states 1, 6, and 7 in E219K, and state 3 in Y221A). The mean first passage time (MFPT) data indicate that once the E219K variant forms one of these bonds-broken states, it remains there for thousands of nanoseconds (8,262.0 ± 2,573.0 ns to 12,769.0 ± 3,327.0 ns) before transitioning to a bonds-formed state (state 3). Likewise, the reverse process also occurs on a microsecond timescale, demonstrating that both directions are kinetically stabilized in E219K. As a result, E219K displays two distinct energy minima in its free-energy landscape, whereas other variants typically show only one (Figure 5A). Prolonged residence in a bonds-broken conformation implies that K67 is more likely to remain deprotonated (consistent with its reduced pK_a in constant-pH MD simulations), enhancing catalytic function. By contrast, in Y221A, the equivalent bonds-broken state (state 3) shifts more readily into other metastable states, including the most stable state 7 (bonds-formed), in only 780.2 ± 46.8 ns. Although the reverse transition from bonds-formed to bonds-broken in Y221A requires a somewhat longer 1,737.6 ± 139.7 ns, this timescale remains far shorter than E219K’s multi-microsecond range. Consequently, Y221A can dynamically switch between ‘formed’ and ‘broken’ conformations with much greater ease. This difference in conformational kinetics helps explain differences in how each mutant enhances hydrolysis rates—E219K achieve it through stable, long-lived ‘active’ conformations, while Y221A rely on faster switching and conformational plasticity.

Figure 5.

The hydrogen bond between Y150(N)-A292(O) is commonly broken in variants, but certain metastable states exist where it re-forms, e.g., V211A (state 7), G214A (states 5, 7, 8), G214R (state 7), and E219K (state 6). This hydrogen bond is likely to play a crucial role in maintaining the R2-loop’s proximity to the active site and stabilizing other hydrogen bonds. Although the fraction of time that Y150(N)–A292(O) is formed in V211A is comparable to other variants, the R2-loop in this mutant is more flexible because the bond can form and break with relative ease (Table 1). Specifically, transitioning in V211A from a bonds-broken state (state 6) to a bonds-formed state (state 7) takes 876.8 ± 78.0 ns, whereas the reverse transition requires 1,874.8 ± 177.6 ns. These timescales show that V211A’s R2-loop toggles between ‘broken’ and ‘formed’ states on the order of nanoseconds to low microseconds, signifying its flexibility. By contrast, G214A exhibits an asymmetry in its transitions: moving from bonds-broken (state 3) to bonds-formed (states 5, 7, 8) takes 1,608.2–2,089.0 ns, but returning to bonds-broken requires 13,175.0–13,657.1 ns. This disparity implies that once the bond is formed, G214A’s loop conformation becomes kinetically stabilized in a ‘closed’ state. G214R shows a similar pattern, requiring 1,701.6 ± 551.9 ns to transition from bonds-broken state (state 4) to bonds-formed state (state 7), but 24,040.6 ± 7,186.8 ns to revert, again indicating that forming the hydrogen bond is far easier than breaking it. Meanwhile, E219K encounters high energy barriers in both directions, based on its long mean first passage times (MFPT), suggesting stable conformations whether the bond is formed or broken.

We also observe that V211A, V211G, G214R, E219G, and E219K variants show metastable states with N287(ND2)-N314(OD1) completely broken. In the case of V211A, state 7 is dominated by conformations that have formed these hydrogen bonds. However, state 3 and state 4 show that the hydrogen bonds are lost. The transition from bonds-formed states (state 2) to bonds-broken states (state 3) is facile and takes only 133.3 ± 12.8 ns, while the transition from bonds-broken states (state 3) to bonds-formed states (state 5) takes 238.4 ± 27.2 ns. In V211G, state 3 and state 4 states exhibit the loss of these hydrogen bonds. The formation of them takes only 48.3 ± 4.1 ns, whereas breaking them takes 613.7 ± 56.6 ns. Despite the Ω-loop being highly stable in V211A and V211G, these states can easily adopt different conformations for the R2-loop, thereby accommodating β-lactams with larger R2 sites. In E219G, state 6 and state 7 show these hydrogen bonds formed and broken, respectively. The MFPT from bonds-formed states (state 6) to bonds-broken states (state 7) is 220.9 ± 23.5 ns, while the MFPT from bonds-broken states (state 7) to bonds-formed states (state 6) is 333.1 ± 19.7 ns. Therefore, high flexibility also can be observed in the R2-loop of the E219G variant. In the case of G214R, the transformation from bonds-broken states (state 3 and 4) to bonds-formed states (state 6), which corresponds to the formation of the hydrogen bond, takes 411.3 ± 21.5 ns and 53.8 ± 5.8 ns, respectively. Conversely, the transition from bonds-formed states (state 6) to bonds-broken states (state 3 and 4) takes much longer, requiring 1393.1 ± 108.4 ns and 1557.8 ± 223.1 ns, respectively. This imbalance means G214R can ‘lock in’ the formed conformation but can still reorganize if given enough time. In the E219K variant, the hydrogen bond is broken in state 1 and state 5, and the transition to state 6, where the hydrogen bond is formed, takes 1621.4 ± 415.5 ns and 4360.4 ± 972.0 ns, respectively. Conversely, the transition from bonds-formed states (state 6) to bonds-broken states (state 1 and state5) takes 3408.7 ± 724.0 ns and 11162.4 ± 3016.3 ns. Therefore, V211A, V211G, and E219G can rapidly form or disrupt the N287(ND2)–N314(OD1) hydrogen bond, whereas G214R and E219K exhibit more substantial kinetic barriers to these transitions.

Conclusions

This investigation of the effects of substitutions in the PDC-3 β-lactamase has provided valuable information on the protein’s dynamics. The study indicates that substitutions can have a significant impact on the stability and flexibility of the Ω-loop and R2-loop, both of which are critical for proper β-lactamase function. Specifically, the G214A, G214R, E214G, and Y221A variants, as well as the wild-type PDC-3, exhibit high flexibility in the Ω-loop, while the V211A and E219G variants show the highest flexibility in the R2-loop. Moreover, a three-pronged hydrogen-bond network around K67—specifically K67(NZ)–S64(OG), K67(NZ)–N152(OD1), and K67(NZ)–G220(O)—governs the proton-transfer step essential for β-lactam hydrolysis. When these three bonds remain intact, K67 tends to stay protonated, limiting its availability to accept protons from S64. Mutations such as E219K and Y221A disrupt this tridentate network, reducing K67’s pK_a and facilitating efficient proton transfer in hydrolysis. Additionally, K67(NZ)–G220(O), the Y150(N)–A292(O) and N287(ND2)–N314(OD1) interactions further modulate R1/R2-loop conformations. Breaking these hydrogen bonds typically shifts the active site to an ‘open’ configuration, accommodating larger cephalosporin substrates. Conversely, their re-formation favors a ‘closed’ state that may stabilize smaller or structurally simpler β-lactams. Overall, the findings of this study provide significant insight into the dynamics of the PDC-3 β-lactamase, revealing the critical roles played by the Ω-loop and R2-loop in its function. These insights gained from this study will aid in the design of more potent antibiotics and β inhibitors for treating bacterial infections.

Methodology

Initial structure preparation

Ten independent all-atom MD simulations of wild-type PDC-3 and its variants were performed. First, the simulations of their conformations were initiated from the X-ray crystallographic structure (PDB ID: 4HEF) at 1.86 Å, after modification of T79A (Lahiri et al., 2013). PDC-3 wild-type and nine variants (V211A, V211G, G214A, G214R, E219A, E219G, E219K, Y221A, and Y221H) were constructed in silico using the ICM mutagenesis program (Abagyan et al., 1994). To ensure the accurate protonation states of the protein, PROPKA 3.0 was employed to assign the protonation states of N-terminus, C-terminus, cationic residues, and anionic residues based on a neutral pH local environment and the protonation states and side chain orientations were also checked by visual inspection (Olsson et al., 2011). In addition, all acidic residues were negatively charged, while alkaline Lys and Arg residues remained positively charged. His was protonated based upon the suggestion by PROPKA 3.0 analysis and also checked by visual inspection.

AdaptiveBandit simulations

Adaptive Bandit MD (AB-MD) simulation is a reinforcement learning based enhanced sampling method that offers a more efficient exploration of the protein’s conformational space while maintaining unbiased, thermodynamically accurate ensembles (Pérez et al., 2020). The advantage of using AB-MD is that it does not alter the underlying potential energy surface, it retains physically realistic dynamics and eliminates the need for reweighting of biased trajectories. As a result, AB-MD can attain a similar or greater depth of conformational sampling with significantly less total simulation time (i.e. lower computational cost) than either extended conventional MD or other enhanced sampling approaches.

The WT and variant structures served as the starting point for a later molecular dynamics (MD) simulation. Multi-microsecond MD simulations of wild-type PDC-3 and its variants were conducted using the Amberff14SB force-field (Maier et al., 2015a). All simulations were run using the ACEMD engine (Doerr et al., 2016; Harvey et al., 2009). Each structure was solvated in a pre-equilibrated periodic cubic box of water molecules represented by the three-point charge TIP3P model, whose boundary is at least 10 Å from any atoms so that the protein does not interact with its periodic images. Periodic boundary conditions in all directions were utilized to reduce finite system size effects. The potassium ions were added to make each system electrically neutral. Long-range electrostatic interactions were computed using the particle mesh Ewald summation method (Cerutti et al., 2009). Subsequently, each system was energy minimized for 5000 steps by conjugate gradient to remove any local atomic clashes and then equilibrated for 5 ns at 1 atmospheric pressure using Berendsen barostat (Feenstra et al., 1999).

Initial velocities within each simulation were sampled from Boltzmann distribution at temperature of 300 K. Isothermic-isobaric NVT ensemble using a Langevin thermostat with a damping of 0.1 ps-1 and hydrogen mass repartitioning scheme to achieve time steps of 4 fs. Multiple short MSM-based adaptively sampled simulations were run using the ACEMD molecular dynamics engine (Doerr et al., 2016; Harvey et al., 2009). The standard adaptive sampling algorithm performs several rounds of short parallel simulations. To avoid any redundant sampling, the algorithm generates an Markov state model (MSM) and uses the stationary distribution of each state to obtain an estimate of their free energy. It then selects any sampled conformation from a low free energy stable state and respawns a new round of simulations. In this context, the MetricSelfDistance function was set to consider the number of native Cα contacts formed for all residues, which were then used to build the MSMs. The exploration value was 0.01 and goal-scoring function was set to 0.3. For each round, 4 simulations of 300 ns were run in parallel until the cumulative time exceeded 30 µs. The trajectory frames were saved every 0.1 ns. 100 trajectories for each system were collected with each trajectory counting 3000 frames.

Constant pH molecular dynamics

To investigate the protonation behavior of K67 in wild-type PDC-3 and its E219K and Y221A variants, K67, Y150, E/K219, and K315 were selected as titratable sites in CpHMD simulations because they lie in proximity to the site of interest (K67) and together form its immediate electrostatic network (Kim et al., 2015; Lahiri et al., 2013). The simulations performed in the Amber suite with the ff99SB force field and an implicit Generalized Born (GB) solvent model (Maier et al., 2015b; Mongan et al., 2004). The protein–solvent complex was energy-minimized for a total of 5,000 steps—10 steps of steepest-descent followed by 4,990 steps of conjugate-gradient—with harmonic positional restraints (10 kcal/mol·Å²) on the backbone atoms to relax side-chain clashes (Brooks et al., 1983; Schlegel, 1982). The system was then heated from 10 K to 300 K over 1 ns, followed by another 1 ns of equilibration at 300 K under Langevin dynamics (Schneider and Stoll, 1978), with a 2 fs time step and SHAKE constraints on hydrogen-containing bonds. During heating, a weaker restraint (2 kcal/mol·Å²) was applied to the backbone atoms to maintain the overall fold while allowing side-chain relaxation, and protonation states were kept fixed in this phase. Equilibrium simulations under constant pH conditions were subsequently conducted at pH 5, 6, 7, 8, 9, 10, and 11 by periodically (every 10 steps) attempting protonation-state changes via a Monte Carlo protocol (Kim et al., 2015). Each pH condition consisted of 50 ns of equilibration followed by 200 ns of production. To confirm convergence, the deprotonation fraction of each residue was monitored over time and found to reach a stable plateau within the final portion of the production simulations. Consequently, the last 50 ns of production for each pH value were used to calculate the deprotonation fractions, ensuring that the analyzed region reflected a converged state. Each system was simulated in triplicate, ultimately providing consistent K67 pK_a estimates for the wild-type, E219K, and Y221A variants. All analyses were done with AmberTools cphstatsand in-house Python scripts (Case et al., 2023).

Markov state models

The PyEMMA software (version 2.5.9) was employed to construct the Markov state models (Husic and Pande, 2018; Scherer et al., 2015). The software determines the kinetically relevant metastable states and their inter-conversion rate from all trajectories of the all-atom molecular dynamics of the wild-type PDC-3 and its variants. Firstly, to evaluate the MSM construction, the conformations defining each frame of the MD trajectories were converted into an intuitive basis. In this step, the features which can represent the slow dynamical modes of these systems were selected. Then, the conformational space was projected to a two-dimensional space using time-lagged independent component analysis (TICA) (Perez-Heandez and Noe, 2016). Using the k-means clustering technique, all conformations from MD simulations were grouped into microstates based on the TICA embedding (Peng et al., 2018). The conformations in the same cluster are geometrically similar and interconvert quickly. After that, the transition matrix between the microstates was built using Bayesian estimation at the appropriate lag time (Trendelkamp-Schroer and Noe, 2013). The lag time was selected where the implied time scales converged, and the transitions between the microstates became the Markovian process. Each indicated time scale represents the average transition time between two groups of states. The microstates were then clustered into a few metastable states using Perron cluster cluster analysis (PCCA) based on their kinetic similarities (Bowman et al., 2009). Additionally, the Chapman-Kolmogorov (CK) test was performed to validate the constructed model further (Barendregt et al., 2019). The CK test measures the reliability of the Markov state models by comparing the predicted residence probability of each microstate obtained from MSMs with those directly computed from MD simulations at longer timescales. Furthermore, the free energies for each metastable state (S_i) were computed from its stationary MSM probability π using the relation:

where π_i denotes the MSM stationary weight of the jth metastable state, k_B is the Boltzmann constant, and T is the temperature. Subsequently, the mean first passage times (MFPT) out of and into the macrostate Si were computed using the Bayesian MSM (Polizzi et al., 2016).

Structural analysis

MDTraj is a robust software package that facilitates the analysis of molecular dynamics (MD) simulations by enabling the manipulation of MD trajectory data from a variety of files (McGibbon et al., 2015). The package provides Python-based tools that allow for the efficient computation of structural and dynamic properties of biomolecules. In this study, MDTraj was utilized to compute several important metrics that are critical in the analysis of MD simulations. Specifically, we used MDTraj to calculate the root-mean-square-fluctuations (RMSF) and root-mean-square-deviation (RMSD) of the protein structure, as well as the x1 torsion angle, hydrogen bonds, salt bridges and the solvent-accessible surface area (SASA). RMSF quantifies the average positional fluctuation of the Cα atom of each residue during the MD simulation relative to its position in the equilibrated reference structures. RMSD quantifies the average displacement of the protein’s Cα atoms from their positions in an equilibrated reference structure over time. In addition, MDLovoFit was used to show the graphical representation of RMSD results (Martínez, 2015). Furthermore, pairwise RMSD analyses were performed using the pytraj package, which allowed us to assess the structural similarities and differences among the conformations sampled during the simulation (Roe and Cheatham, 2013). The x1 torsion angles were calculated for each residue except for alanine and glycine. Hydrogen bonds were defined as a distance of less than 3.5 Å between a hydrogen bond donor and acceptor, with a hydrogen-donor-acceptor angle greater than 30°. Salt bridges were defined as a distance of less than 4.0 Å between a positively charged amino acid side chain (lysine or arginine) and a negatively charged side chain (aspartate or glutamate). The volume of the active site pocket was calculated using ParkVFinder (Guerra et al., 2020). Visualization of the structures of protein was performed using PyMOL (Schrödinger and DeLano, 2020).

Data Availability

All files required to run the simulations (topology, coordinates, input), processed trajectories (xtc), corresponding coordinates (pdb), metastable PDB files for each system described in this manuscript can be downloaded from the DOI https://doi.org/10.57760/sciencedb.15876.

Acknowledgements

Research reported herein was supported in part by funds the National Institute of Allergy and Infectious Diseases of the National Institutes of Health under Award Number R01AI063517 to RAB and SH. The content is solely the responsibility of the authors and does not necessarily represent the official views of the Department of Veterans Affairs or the National Institutes of Health.

Additional information

Funding

National Institute of Allergy and Infectious Diseases (R01AI063517)