Insight into the PmrB structures of colistin-resistant Gram-negative bacteria through the multi-template ligand-guided homology modeling and in silico mutagenesis

Retrieval of amino acid sequence and template identification

The amino acid sequence of the PmrB histidine kinase originating from the selected Gram-negative pathogens was extracted from the GenBank database (Sayers et al., 2020). A comprehensive search was conducted using specified keywords for bacterial strain in combination with “histidine kinase” and “PmrB”. The most relevant sequences (as shown in Table S1) were then selected for downstream analysis.

In order to identify appropriate templates for homology modeling, a Position-Specific Iterated BLAST (PSI-BLAST) search was conducted on the Protein Data Bank (PDB) (Berman et al., 2000). The investigation was performed iteratively twice to minimize incorrect alignments that occurred from randomization. To be considered as templates for homology modeling, each structure must meet specific criteria, including containing the ATPase domain in holo-form with a relevant ligand in the binding pocket, which is essential for maintaining a proper binding site architecture of the protein-ligand complex during model construction. Moreover, we also extended the template searching to the AlphaFold database (Varadi et al., 2022) for use as complementary templates to enhance the predictive performance of homology modeling. The templates that satisfied these requirements are summarized in Table 1.

Multi-template homology modeling with ligand-specific restraints

To address the challenges posed by a low sequence identity of less than 30 percent on average and a limited amino acid sequence coverage of currently available templates, a multi-template homology modeling strategy was employed using MODELLER (version 10.4) (Webb & Sali, 2016). The approach involved the utilization of multiple templates (Table 1) incorporated with ligand information to construct more precise and dependable models. Prior to performing homology modeling, the target-template sequence alignment was carried out using ClustalOmega (Sievers et al., 2011) and prepared as inputs for homology modeling. Consequently, we generated an ensemble of 1,000 homology models with customized parameters by maximizing the structure refinement and optimization to ensure the accuracy of protein models. FoldX (version 5.0) (Delgado et al., 2019) was employed to provide additional optimization by fixing torsional clashes of amino acid side chains. Furthermore, the amino acid sequences were submitted to the SWISS-MODEL (Waterhouse et al., 2018), I-TASSER (Yang & Zhang, 2015), and AlphaFold (Varadi et al., 2022) to generate homology models for comparison.

Molecular docking and model selection

To select the appropriate protein conformation that matches ATP structure, a natural ligand of histidine kinase, molecular docking was conducted using AutodockFR (Ravindranath et al., 2015). This process was assisted by in-house developed Python scripts and the ChemAxon Marvin Suite for drawing and characterizing chemical structures (ChemAxon, 2023). The genetic algorithm parameters were customized as follows:

• Number of Genetic Algorithm evolutions (nb_runs): 20

• Maximum evaluations of the scoring function per Genetic Algorithm run (max_evals): 4,000,000

• Termination of Genetic Algorithm evaluations after consecutive generations with no improvement (no_improve_stop): 5

• Maximum number of generations (max_gens): 10,000,000

• Seed for the random number generator to ensure reproducibility (seed_value): 1.

The models were ranked based on their binding energy, and the top-ranked models were chosen to compare with the model generated from other approaches. The Protein-ligand Interaction Profiler (PLIP) (Adasme et al., 2021) was also used to analyze the interactions between ATP and PmrB obtained from the best-performing model. The results were visualized in PyMOL software (Schrödinger, 2024) for a comprehensive structural analysis. The transmembrane regions were identified using the Phobius web server (Käll, Krogh & Sonnhammer, 2007). The structural validation was then performed to assess the model quality with multiple validation matrices, including Discrete Optimized Protein Energy (DOPE) score (Shen & Sali, 2006), MolProbity (Williams et al., 2018), QMEANDisCo (Studer et al., 2020), and Ramachandran plots using the SWISS Structure Assessment web server (Waterhouse et al., 2024). The detailed Python scripts implementing this multi-template ligand-guided homology modeling workflow are publicly available in our GitHub repository (https://github.com/nuttapat/Modeller-AutodockFR).

Molecular dynamics simulation

Molecular dynamics (MD) simulations were performed using GROMACS version 2024.3 (Abraham et al., 2024). The topology and parameters for the ATP ligand were obtained from the CGenFF webserver (Vanommeslaeghe, Raman & MacKerell Jr, 2012), ensuring compatibility with the CHARMM36 force field (Vanommeslaeghe et al., 2010), employed for all simulations. Each PmrB-ATP complex was carefully prepared for simulation by solvating it in a cubic box with extended simple point charge (SPC/E) water molecules, ensuring a minimum distance of 10 Å between the protein and the box edges. The protonation states of the N- and C-termini of the protein were fixed to match the requirements of the CHARMM36 force field, and counterions were added to maintain the overall charge neutrality of the system.

At the pre-MD simulations step, the system underwent a two-stage equilibration process, consisting of a 100 ps canonical ensemble (number of particles, volume, and temperature, NVT) equilibration at 300 K followed by a 100 ps isothermal-isobaric ensemble (number of particles, pressure, and temperature, NPT) equilibration at 300 K and 1 bar pressure. Subsequently, two independent MD simulations were then conducted for 200 ns at 300 K and one bar pressure with random seed number. The trajectories of all atoms in the system were recorded and analyzed using Python scripts with the Matplotlib library to calculate various parameters, including the root mean square deviation (RMSD) of the protein backbone and the ATP ligand, the root mean square fluctuation (RMSF) of individual residues, the radius of gyration (Rg) of the protein, the solvent accessible surface area (SASA), and the number of hydrogen bonds formed between PmrB and ATP. Furthermore, the average atomic distance between key interacting residues and functional moieties of ATP during 200 ns simulations was calculated to emphasize the dynamic binding behavior at the binding sites. These investigations provided valuable insights into the stability, flexibility, and interactions of the PmrB-ATP complex.

Binding free energy calculation

Binding free energies between PmrB and ATP were estimated using the Molecular Mechanics Generalized Born Surface Area (MM/GBSA) method implemented in the gmx_MMPBSA tool (Valdés-Tresanco et al., 2021). The generalized Born (GB) model was employed to estimate the polar solvation energy efficiently. Per-residue energy decomposition analysis was performed using the default parameters in gmx_MMPBSA to identify key residues contributing to the PmrB-ATP interaction. The information obtained from this analysis is essential for further understanding the main driving force of the complex formation.

Structural analysis of PmrB mutation

In order to investigate the effects of PmrB mutation on ligand binding, in silico mutagenesis was conducted on the PmrB model of K. pneumoniae due to the availability of a large amount of genetic characterization data that have been previously identified in genomic studies of colistin-resistant K. pneumoniae clinical isolates. This was based on prior research indicating amino acid substitutions that lead to an increase in vitro colistin resistance (Elias, Duarte & Perdigão, 2021). The BuidModel function of FoldX was utilized to assign amino acid substitutions with a particular focus on the histidine kinase domain, including T246A, R256G, K281L, L339C, H340I, N341T, R342D, Q343S, L344P, and P346Q, to the PmrB model and evaluated for internal energy contribution. The molecular docking was performed to calculate the binding energy and predict the binding mode of ATP to the mutants. Finally, the structural comparison between wild-type and mutants was carried out in PyMOL.

Ligand-guided homology modeling

Multi-template homology modeling was performed to create PmrB histidine kinase homology models for four bacterial species: A. baumannii, E. coli, K. pneumoniae, and P. aeruginosa. Suitable templates were identified using PSI-BLAST searches and supplemented with structures from AlphaFold, as shown in Table 1. The search for templates revealed that three histidine kinase crystal structures in complex with their ligand (4BIW, 5C93, and 4KP4) were closely related to our query amino acid sequence, based on the highest alignment score (MaxScore). Although there was limited sequence coverage for PmrB in the selected bacterial species, with a percentage of sequence coverage of 65%, 60%, 44%, and 59%, respectively, the amino acid sequence of all selected templates covered the histidine-kinase ATPase domain. This domain is crucial for studying protein-ligand interaction.

Despite the limitations of low sequence identity and coverage, we utilized multi-template homology modeling in MODELLER to improve the accuracy of protein models, especially at the ATP-binding site, which is covered by a highly flexible lid loop. This involved generating 1,000 protein models for each deployed template (both singles and multiple templates) and refining them using FoldX for side-chain optimization. The models were ranked by calculating the DOPE score, a statistical potential that estimates the energy of the protein structure and distinguishes between proper and incorrect folds. Lower DOPE scores indicate more favorable energy states and may represent more accurate models (Khare et al., 2019). It can be noted that the PmrB models generated using the AlphaFold template exhibit lower DOPE scores compared to the PDB templates identified from the PSI-BLAST, as shown in Table 2 and Table S3. Furthermore, the multi-template models also provide a superior DOPE score than those constructed from a single PDB template. This may result from the inaccuracy of modeling missing regions in sequence alignment. Therefore, utilizing multi-template and ensemble protein modeling may improve the structural prediction of proteins with low sequence identity and coverage to the template. These findings demonstrate the potential of using the AlphaFold as a complementary template, which can be beneficial in the modeling of PmrB structures.

Although the DOPE score can be used as a default scoring for structural assessment, it did not provide information on the arrangement of protein structural motifs and binding pockets. To improve the model selection with precise binding pocket architecture, molecular docking was additionally performed to assess the binding of ATP, the natural ligand of histidine kinases, to all generated PmrB models. This served to select the representative protein conformations based on the binding fitness to the endogenous ligand, thus ensuring the correction of the ligand binding site (Schaller et al., 2019). The best docking score of all PmrB models generated by ligand-guided homology modeling with different types of templates for all bacterial species exhibited superior binding affinity as indicated by docking score in the range of −10.00 kcal/mol to −12.80 kcal/mol, while the average docking scores ranged from −7.11 kcal/mol to −9.49 kcal/mol suggested lower binding affinity to the targets, as shown in Table 2. Interestingly, the docking score of ATP against all generated models from the ligand-guided homology modeling also significantly outperformed the models generated by the AlphaFold, I-TASSER, and SWISS-MODEL. These results highlight the synergistic advantage of combining ensemble model generation with molecular docking for model selection, enabling the identification of PmrB conformations that accurately represent the stable, ligand-bound state. This integrated approach demonstrates its superiority over conventional homology modeling methods by consistently yielding models with improved ATP-binding site architecture, as evidenced by the superior docking scores.

Analysis of the top 10 models, ranked by binding energy from the ligand-guided homology modeling, consistently revealed similar ATP binding sites. In contrast, the top 10 models generated from conventional homology models, which were ranked by DOPE score, exhibited substantial deviations in the ATP binding pocket. This discrepancy arises from variations in ATPase domain orientation and, most notably, the conformation of the flexible lid loop region within the binding pocket (Fig. S1). These findings underscore the importance of utilizing a ligand-guided approach for model selection, prioritizing binding energy as a key metric to ensure an accurate representation of the ligand-bound state. The final refined PmrB models and their binding modality with ATP for A. baumannii, E. coli, K. pneumoniae, and P. aeruginosa are illustrated in Figs. 2A–2D.

Structural assessment of PmrB models

The structural integrity and reliability of the final PmrB homology models were rigorously evaluated using multiple assessment metrics, as shown in Fig. 3. MolProbity scores for all models ranged from 2.58 to 2.82, indicating good stereochemical quality. Ramachandran plot analysis further confirmed this, with over 95% of residues in favored and allowed regions with minimal outliers. The overall model quality, as determined by QMEANDisCo global scores, demonstrated consistently high values, ranging from 0.55 to 0.57. This data further validates the structural reliability of the models. Importantly, all models exhibited normalized QMEAN scores (Z-scores) less than 1, indicating quality comparable to experimentally determined structures in the Protein Data Bank.

Additionally, the DOPE scores of our ligand-guided homology models consistently outperformed those generated using SWISS-MODEL and I-TASSER alone, as summarized in Table 2. This underscores the effectiveness of our modeling strategy, especially for low-sequence identity targets like PmrB. This comprehensive assessment provides strong evidence for the accuracy and suitability of our models for further investigations into PmrB structure, function, and interactions with ATP, paving the way for potential therapeutic interventions against colistin resistance.

Analysis of MD trajectory

To assess the stability and dynamic behavior of the constructed PmrB homology models, we conducted two independent all-atom MD simulations for 200 ns. Details of the dynamic behavior of all MD trajectories were summarized and illustrated in Figs. 4A–4F, Figs. S2, and S3. The RMSD of the protein backbone was analyzed to evaluate the structural deviation over time. As shown in Fig. 4A, all four PmrB models reached equilibrium after approximately 125 ns, indicating that the structures stabilized during the simulation. The average RMSD values for A. baumannii, E. coli, K. pneumoniae, and P. aeruginosa were 1.73 ± 0.26 nm, 0.89 ± 0.14 nm, 1.16 ± 0.15 nm, and 1.89 ± 0.37 nm, respectively. Minor fluctuations were observed in the A. baumannii and P. aeruginosa models, suggesting some degree of flexibility in these proteins. The RMSD of the ATP ligand within the ATP-binding site (Fig. 4B) also showed stabilization with some minor fluctuations later, with average RMSD values of 0.23 ± 0.02 nm, 0.24 ± 0.03 nm, 0.21 ± 0.02 nm, and 0.17 ± 0.03 nm for A. baumannii, E. coli, K. pneumoniae, and P. aeruginosa, respectively. This indicates a lower degree of structural deviation of the ligand inside the binding pocket of PmrB, conferring stability of binding from all constructed models.

The flexibility of different regions of the PmrB protein was assessed by analyzing the RMSF of individual residues, as shown in Fig. 4C. It can be noted that loop regions and terminal ends exhibited higher RMSF values, indicating greater flexibility during simulation. The average RMSF values for A. baumannii, E. coli, K. pneumoniae, and P. aeruginosa were 0.58 ± 0.17 nm, 0.49 ± 0.20 nm, 0.54 ± 0.36 nm, and 0.91 ± 0.27 nm, respectively. In contrast, the core regions of the protein showed lower RMSF values, suggesting that these regions maintain a relatively rigid and stabilized structure. Interestingly, similar fluctuation patterns were observed amongst the catalytic domain of all four bacterial species, suggesting the conserved structural dynamics behavior of PmrB.

The overall compactness of the PmrB models was evaluated using the radius of gyration (Rg) as shown in Fig. 4D. The average Rg values for A. baumannii, E. coli, K. pneumoniae, and P. aeruginosa were 4.21 ± 0.18 nm, 3.77 ± 0.09 nm, 3.62 ± 0.18 nm, and 4.54 ± 0.30 nm, respectively. The Rg values for all four species remained relatively stable after an initial simulation period, further supporting the structural stability of the models.

The solvent-accessible surface area (SASA) was analyzed to assess the extent of the protein surface exposed to the solvent. The average SASA values for A. baumannii, E. coli, K. pneumoniae, and P. aeruginosa were 279.99 ± 3.59 nm², 230.23 ± 2.90 nm², 236.76 ± 9.93 nm², and 292.02 ± 10.00 nm², respectively. While some initial fluctuations were observed in Fig. 4E, particularly for P.aeruginosa, the SASA values generally stabilized over time, suggesting that the overall solvent exposure remained relatively consistent.

We also analyzed the number of hydrogen bonds (Hbond) formed between PmrB and the ATP during the simulation, as demonstrated in Fig. 4F. The average numbers of hydrogen bonds for A. baumannii, E. coli, K. pneumoniae, and P. aeruginosa were 5.35 ± 1.23, 4.58 ± 1.30, 8.09 ± 1.16, and 5.82 ± 1.16, respectively. The plots revealed dynamic fluctuations in the number of hydrogen bonds, which is expected due to the dynamic nature of these interactions. It is notable that the Hbond of K. pneumoniae was significantly greater than that of the others, indicating higher stabilization of the complex during MD simulation. The results regarding the average Hbond offered insights into the overall stability of the PmrB-ATP complex.

MM/GBSA and per-residue decomposition analysis

To further investigate the PmrB-ATP interaction, we performed MM/GBSA analysis to calculate the energy change upon binding (ΔG) for each PmrB-ATP complex. The ΔG values, along with their standard deviations (SD) and standard errors of the mean (SEM), are summarized in Table 3, while the details of energy contributions and per-residue energy decomposition from each MD trajectory are presented in Figs. S4 and S5. The results showed that all four PmrB-ATP complexes exhibited negative ΔG values, indicating that the formation of the complex is energetically favorable from an enthalpic perspective. The individual energy components contributing to ΔG, namely van der Waals (VDWAALS), electrostatic (EEL), polar solvation (EGB), and nonpolar solvation (ESURF) energies, are shown in Figs. 5A–5D. It can be noted that VDWAALS and EEL interactions were the major driving forces for binding in all complexes, contributing significantly to the negative ΔG values. In contrast, the EGB and ESURF terms represented unfavorable contributions due to the desolvation penalty upon complex formation.

To identify the key residues involved in PmrB-ATP binding, we conducted per-residue energy decomposition analysis. The results, presented in Figs. 5E–5H, and Table S3, revealed that several amino acid residues made substantial contributions to the ΔG obtained from the PmrB-ATP complex formation. In A. baumannii, GLY399, ASN341, and GLY401 exhibited the most substantial contributions, with average interaction energies contribution of −2.63, and −2.57 kcal/mol, respectively. Similarly, in E. coli, ASN264, TYR268, and LEU322 were the major contributors, with interaction energies ranging from −2.68 to −2.05 kcal/mol. In K. pneumoniae, GLY320, TYR269, GLY322, and ASN265 played a crucial role, with interaction energies contribution ranging from −4.63 to −2.42 kcal/mol. Finally, in P. aeruginosa, ASN363, GLY423, LEU424, and GLY421 were the major contributors, with interaction energies ranging from −6.10 to −2.63 kcal/mol. These residues, along with other key contributors, form a network of interactions that anchor ATP within the binding site. The details of energy contributions and per-residue energy decomposition from each MD trajectory are presented in Figs. S4 and S5.

Interestingly, several residues exhibiting strong favorable interactions with ATP, such as ASN265 (−2.42 kcal/mol) in K. pneumoniae, are not directly involved in binding to the ATP but are located adjacent to the residues in the binding pocket. This suggests that these residues may play a crucial role in supporting the alignment of residues within the binding site, facilitating an induced fit mechanism that enhances the interaction with ATP. This observation underscores the importance of considering not only the direct interactions between the protein and ligand but also the contributions of neighboring residues in mediating the binding process. These findings are consistent with our observation from the in silico mutagenesis study, where mutations in residues adjacent to the binding pocket, such as T246A and L344P in K. pneumoniae, led to significant changes in ATP binding affinity. This further highlights the importance of these neighboring residues in modulating the PmrB-ATP interaction. The structural superimposition between the initial structures and the equilibrated models possessing the lowest ΔG were depicted in Fig. S6.

Structural analysis of PmrB models

The ligand-guided homology models of PmrB across the four Gram-negative pathogens reveal the conserved 3D architecture of this histidine kinase while also highlighting some species-specific variations, as presented in Fig. 2. Each PmrB model consists of a signal peptide, two transmembrane helices (TM1 and TM2), an extracellular sensor domain, a HAMP domain, a DHp domain, and a catalytic domain that includes the ATP-binding site (Wang, 2012). The ATPase lid loop, a flexible region critical for ATP binding, is also located within the catalytic domain. The summary of the amino acid position of each substructural domain from all four bacterial strains was summarized in Table S2.

Across all species, the transmembrane helices (TM1 and TM2) are predicted to be conserved in terms of length and sequence, consistent with their role in anchoring the protein to the cell membrane and mediating signal transduction. In contrast, the extracellular sensor domain, responsible for detecting environmental stimuli, exhibits notable differences in length between the Enterobacterales family (E. coli and K. pneumoniae) and non-fermentative bacteria (P. aeruginosa and A. baumannii). In Enterobacterales, this domain consists of approximately 30 amino acids, whereas in non-fermentative bacteria, it is considerably larger, with 116 and 127 amino acids for A. baumannii and P. aeruginosa, respectively. This variability likely reflects adaptations to the distinct environmental niches and signaling mechanisms employed by these two groups of bacteria.

The HAMP and DHp domains, crucial for signal transduction and autophosphorylation, respectively, are relatively conserved in length and sequence across the four species (Table S2). This conservation underscores their essential role in PmrB function and colistin resistance. The catalytic domain, which contains the ATP-binding site, is the largest domain in all models and shows moderate variations in length. Notably, the ATPase lid loop, a flexible region that regulates access to the ATP-binding site, is highly conserved in length and sequence across all four species, with the exception of P. aeruginosa, which has a slightly longer loop (18 amino acids compared to 16 in the other species). This high degree of conservation suggests a critical role for this loop in ATP binding and PmrB function (Marina et al., 2001; Marina, Waldburger & Hendrickson, 2005; Albanesi et al., 2009; Trajtenberg et al., 2010).

Analysis of protein-ligand interaction profiles of all generated PmrB models in complex with ATP was performed to investigate the key interacting residues and their moieties of preference among PmrB from all four bacterial species. The results reveal a conserved pattern of ATP binding across the four bacterial species, with some notable species-specific variations, as shown in Table S2. All models exhibit hydrogen bonding between the adenine moiety of ATP and residues in the catalytic domain, highlighting the importance of these interactions for ATP binding. Notably, direct π-π stacking interactions were observed between adenine and Y268 in E. coli and Y367 in P. aeruginosa with average atomic distances of 0.25 ± 0.06 and 0.28 ± 0.05 nm, respectively (Fig. 6). In K. pneumoniae and A. baumannii, the corresponding tyrosine residues (Y269 and Y345, respectively) are located close to the adenine moiety, suggesting a strong potential for π − π stacking interactions. Although the distance exceeds the strict criteria for classification by PLIP, the average atomic distances calculated from two independent MD simulations are 0.33 ± 0.12 and 0.30 ± 0.05 nm, respectively, indicating a significant likelihood of binding throughout the simulation. Additionally, all models display hydrogen bonds between the ribose moiety of ATP and residues in the lid loop or catalytic domain, further stabilizing the ligand-protein complex.

While the overall pattern of ATP binding is conserved, species-specific interactions are observed, particularly in the residues involved in hydrogen bonding and salt bridge formation with the phosphate groups of the ATP (Table S4). For example, A. baumannii and E. coli exhibit salt bridges with residues K344 (0.35 ± 0.10 nm) and R310 (0.44 ± 0.18 nm), respectively, while the other species primarily interact with the phosphate groups through hydrogen bonds. Additionally, K. pneumoniae and P. aeruginosa might exhibit hydrogen bonds between the phosphate groups and residues in the DHp domain (T189 and Q193 in K. pneumoniae, Q294 in P. aeruginosa), suggesting a potential role for this domain in stabilizing ATP binding in these species.

The catalytic domain and the lid loop emerge as the primary regions involved in ATP binding across all four species, consistent with their known function in ATPase activity. However, the observed interactions in the DHp domain of K. pneumoniae and P. aeruginosa suggest a more complex interplay between domains in these species. These findings underscore the importance of considering both conserved and species-specific interactions when investigating PmrB function and designing potential inhibitors.

These structural insights, derived from the ligand-guided homology models, provide a foundation for understanding the molecular basis of PmrB function and colistin resistance. The conserved architecture across species suggests a common mechanism of action, while the species-specific variations highlight potential targets for the development of tailored therapeutic interventions. Further investigation of these structural differences may reveal insights into the varying colistin susceptibility observed among different Gram-negative pathogens.

Impact of PmrB mutation upon ATP-binding activity

Despite the major functionality of PmrB being regulated by ATPase activity, there are many reports about mutations located outside the ATPase binding pocket that can increase the level of colistin resistance. In silico mutagenesis of the K. pneumoniae PmrB model was conducted to gain further understanding of the effects of amino acid substitution upon ATP binding. The initial analysis of PmrB models revealed the presence of T246A/L344P double mutations in the original amino acid sequence obtained from the Genbank. Therefore, the wild-type PmrB model was subsequently generated by reverting these mutations to their original residues (A246T, P344L) and will be used as a baseline model for mutational analysis. All mutations that deployed in this study included L17Q, L82R, S85R, T140P, D150H, T157P, S205P, S208N, T246A, R256G, K281L, L339C, H340I, N341T, R342D, Q343S, and P346Q, all of which have been previously associated with colistin resistance. Herein, a schematic summary of all amino acid substitutions and their locations relative to each functional domain on the PmrB structure of K. pneumoniae was depicted in Fig. 7.

Interestingly, the results revealed a diverse range of effects on ATP binding affinity and protein stability, as summarized in Table 4 and Table S5. It can be noted that K281L and L339C mutations within the catalytic domain resulted in a decrease in binding affinity of 0.6 and 0.8 kcal/mol, respectively, compared to the wild type. These residues may likely play crucial roles in stabilizing the ATP-bound conformation, and their substitution may disrupt key interactions within the binding pocket. Conversely, the T246A and Q343S mutation in the catalytic domain led to an increase in ATP binding affinity of −1.0 and −0.8 kcal/mol, respectively. The major notable increase in the docking score was observed in T246A and L344P double mutations presented in our constructed PmrB models, which increase the binding affinity of −1.7 kcal/mol, suggesting the synergistic effects to enhance the ATP-binding.

In terms of protein stability, mutations T140P in the HAMP domain, R342D in the catalytic domain, S205P in the DHp domain, as well as L82R in the transmembrane region exhibited the most significant destabilizing effects over other mutations in their corresponding domain, with increases in total energy of 4.177, 2.989, 1.329, and 0.509 kcal/mol, respectively, compared to the wild-type baseline. The notable increase of contributed energy in T140P mutation was found in Van der Waals clashes, Backbone HBond, and Sidechain HBond with values of 4.133, 1.873, and 0.589 kcal/mol, respectively, as shown in Table S5. A similar pattern was found in the S205P mutation with the value of 2.560, 1.354, and 0.632 kcal/mol, respectively. The results suggested that these mutations were likely to increase inter-atomic forces and disrupt essential intramolecular interactions, such as hydrogen bonds or salt bridges, leading to a less stable protein than the wild type. In the meanwhile, the total energy change of R342D mutation was driven by Backbone HBond, Sidechain Hbond, Van der Waals, Electrostatics, Solvation polar, and Solvation hydrophobic with the value of 0.451, 0.643, 1.062, 1.6, 0.153, and 0.957 kcal/mol, respectively. Given that the values for Van der Waals, Electrostatics, and Solvation Hydrophobic are relatively high, these interactions are likely the main contributors to the destabilization of the protein upon mutation. A different energy pattern was observed in L82R, which functions to maintain the structural integrity inside the cell membrane. The total energy of this mutation was contributed by Backbone HBond, Electrostatics, Solvation polar, Solvation hydrophobic, Van der Waals clashes, Entropy sidechain, Entropy mainchain, and Torsional clash with the values of 0.02, 0.05, 0.079, 0.121, 0.009, 0.096, 0.122, and 0.064 kcal/mol, respectively. However, these values are all relatively small, so the degree of destabilization in transmembrane mutation may not be very significant as compared to the aforementioned mutations.

The mutations in L344P and T246A in the catalytic domain, T157P in the DHp domain, and S85R in the transmembrane region conferred increased stability. This was shown by a reduction in total energy of −1.381, −0.934, −1.113, and −0.738 kcal/mol, respectively, compared to the wild type. The L344P mutation displays a notable increase in stability, primarily attributed to a substantial favorable change in mainchain entropy, reduced backbone, and torsional clash with the values of −0.955, −0.571, and −0.222 kcal/mol. In contrast, the results of T246A mutation exhibited different energetic changing patterns in which Van der Waals clashes, solvation polar, and entropy sidechain gain more flavor with the values of −1.473, −1.232, and −0.551 kcal/mol, respectively. The DHp domain mutant of T157P achieved more stability, predominantly due to favorable changes in the entropy of the main chain, reduction of backbone clash, and more flavor solvation hydrophobic with the value of −1.007, −0.957, and −0.735 kcal/mol, respectively. Another mutation in the transmembrane region was observed in S85R, which increases stability slightly. The favorable increase of total energy change of this mutant was contributed by a reduction in solvation hydrophobic, Van der Waals, and backbone HBond with the values of −0.826, −0.503, and −0.407 kcal/mol, respectively.

Interestingly, the T246A/L344P double mutant, representing our PmrB models from the obtained K. pneumoniae sequence, presents an intriguing case due to a greater decrease in total energy of −4.3 kcal/mol, which indicates a potential stabilizing effect over other mutants. The favorable energy change of this combination was primarily driven by a significant reduction in torsional clash (−3.717 kcal/mol) and solvation polar (−1.897 kcal/mol), suggesting the positive effects of mutation towards the torsion angles in the protein backbone and interaction between polar groups in the protein against polar solvent.

The results of in silico mutagenesis revealed a diverse range of effects on protein stability and their energy contribution. However, this experiment did not provide a direct explanation of their effects on the ATP binding affinity. Thus, we conducted a molecular docking study between all PmrB mutants with ATP to investigate the effects of this substitution upon binding. The docking scores of all mutants and wild-type PmrB models were summarized in Table 4. It can be noted that the best docking score was observed in T246A/L344P double mutation with the value of −10.8 kcal/mol followed by −10.1, −9.9, −9.9, and −9.8 kcal/mol of T246A, Q343S, S205P, and L344P, respectively. While the wild-type model exhibited −9.1 kcal/mol as a reference point, all aforementioned mutants achieved a significant increase in ATP binding affinity with over 0.5 kcal/mol from the wild type. These findings suggested that the mutation in these positions may confer the intricate relationship between PmrB structure and their ATPase function, which could potentially lead to PmrB signaling and the response to colistin resistance.