The following four isoforms of BabA complexed with Le^b were used in the current study: 1) Iso1 (PDB ID 5F7M), 2) Iso2 (PDB ID 5F9A), 3) Iso3 (PDB ID 5F8R), and 4) Iso4 (PDB ID 5F93) (Moonens et al., 2016). We also simulated apo-BabA (without Le^b) and Le^b in an aqueous solution. The AMBER ff14SB force field (Maier et al., 2015) was used for the protein modeling, and glycam06j-1 (Kirschner et al., 2008) was employed to describe Le^b. The protonation state of each residue was estimated using the Propka 3.1 webserver (Olsson et al., 2011), corresponding to pH seven for all four different Isoforms of BabA. Each system was solvated in a truncated octahedral box using the TIP3P (Price and Brooks, 2004) water molecules. The distance between the wall and solute was kept as 10 Å. Each solvated system was neutralized by adding an appropriate number of counterions. All bonds, including hydrogen atoms, were restrained using the SHAKE algorithm (Kräutler et al., 2001). The long-range electrostatic interactions were treated using the particle mesh-Ewald (PME) method (Darden et al., 1993). A cut-off distance of 10 Å was set for non-bonded interactions. The temperature was maintained using the Langevin thermostat with a collision frequency of 2 ps⁻¹ and the pressure using the Berendsen barostat (Berendsen et al., 1984). Two stages of minimization were performed using the steepest descent followed by the conjugate gradient algorithm for 500 steps with a harmonic restraint of 2 kcal mol⁻¹ Å⁻². It was then subjected to gradual heating from 0 to 300 K with a restraint force of 2 kcal mol⁻¹ Å⁻². Before the production run, we conducted a 1 ns simulation at the NPT ensemble without restraint on the solute for equilibration. Finally, the production simulations for complexes and apo-BabA were performed in the NPT ensemble with a timestep of 2 fs for 1 µs using the pmemd.cuda module of AMBER18(Case et al., 2018). We recorded the coordinates every 10 ps, resulting in 100,000 configurations. To ensure better sampling, we conducted two more replica simulations for apo and complex systems. We also conducted 2 × 3 µs long simulations for the free glycan in the aqueous solution using a similar protocol as discussed in our previous work (Roy et al., 2022).

In addition, to investigate the role of an acidic environment on the complexation of BabA/Le^b, we conducted constant pH molecular dynamics (CpHMD) simulation at pH 4.5 for 2 × 250 ns for Iso1 by titrating acidic amino acids in the crown region (Donnini et al., 2011). In the selected region, six ionizable residues were found (E80, E158, D161, D172, D199, and D213), which we attempted to titrate during the simulation. After every 180 fs of explicit MD simulation, the protonation change event was attempted for all mentioned residues in the GB implicit model (Mongan et al., 2004).

All the trajectories were analyzed using the cpptraj module of AmberTools19 (Roe and Cheatham, 2013). The first 200 ns data from each run were discarded for the trajectory analysis to minimize initial noises. After that, the combined trajectories for three replicas were prepared and subjected to various calculations related to the dynamics of the systems. The initial stability and flexibility were estimated by calculating the root mean squared deviations (RMSD) and the root mean squared fluctuations (RMSF) concerning the well-equilibrated conformations. Conformational space of the sugar molecules in the free and bound state was also estimated with the help of glycosidic dihedral angle, Cremer–Pople (CP) (Cremer and Pople, 1975) parameters of ring puckering. The glycosidic dihedral angles for the α/β-n linkages are defined as φ = O 5 − C 1 − O − C n and ψ = C 1 − O − C n − C ( n − 1 ) . The puckering state of any sugar as per the Cremer–Pople convention was estimated using the spherical coordinate set, Q, θ, and φ. The spherical space was further subdivided into 38 IUPAC recognized regions, where each space is dedicated to a specific ring conformation. To estimate the global and local conformations of several loops related to binding, the Cartesian principal component analysis (PCA) and dihedral principal component analysis (dPCA) were performed, which were discussed in detail in our previous study (Roy et al., 2021). The coordinate of C α atom of each amino acid and dihedral angle of peptide bonds of the selected loop region were used for PCA and dPCA, respectively. Furthermore, the free energy surface (FES) was constructed using the Boltzmann equation: ΔG = −k _B T ln(ρ) - k, where k _B is the Boltzmann constant, T is the temperature, and ρ is the probability density of the geometric coordinate x. The parameter (k) was chosen such that the global minimum was located at 0 kcal/mol. In addition, a distance cut-off of ≤3.5 Å and an angle cut-off of ≥120° were used to evaluate the hydrogen bonding (Jeffrey, 1997).

The molecular mechanics/generalized Born surface area (MM/GBSA) scheme (Onufriev et al., 2004) was employed to estimate binding free energy. For protein-glycan systems, it has been reported that MM/GBSA outperforms the computationally expensive MM/PBSA scheme (Mishra and Koča, 2018; Sood et al., 2018; Roy et al., 2020a). It offers a good trade-off between speed and accuracy compared to other computationally expensive methods, such as thermodynamics integration (TI) and free energy perturbation (FEP). The methodological details of the MM/GBSA protocol were discussed in our previous works (Kar et al., 2013; Roy et al., 2020b; Jonniya et al., 2021). The binding free energy ( Δ G b i n d i n g ) was determined from the free energies of the receptor/BabA (protein), the ligand (carbohydrate), and the complex (complex) according to the following equation (Kollman et al., 2000): Δ G b i n d i n g =   G c o m p l e x − ( G c a r b o h y d r a t e +   G p r o t e i n ) (1)

The free energy of each species (protein, carbohydrate, and complex) was estimated using the following formula: G = E v d W + E e l e + G G B + G S A (2) where E v d W and E e l e represent the van der Waals and electrostatic interaction energy, respectively. The polar solvation free energy ( G G B ) was estimated using the generalized Born equation, while the non-polar solvation free energy ( G S A )   was determined using the following equation: G S A =   γ A + b (3) with γ = 0.00542   kcal . mol − 1 . Å − 2 and b = 0.92  kcal . mol − 1 . The symbol A denotes the solvent-accessible surface area (SASA). The GB model (igb = 2) developed by Onufriev et al. was used in this study (Onufriev et al., 2004). The configurational entropy was not considered in the current study due to the high computational cost. We considered 25000 frames from the last 500 ns of each production simulation for the free energy calculation. The overall binding free energy was estimated by averaging the binding free energies obtained from the three independent runs. In addition, we performed a per-residue decomposition of the binding free energy by the MM/GBSA scheme. On the other hand, 7500 frames from the last 150 ns of the trajectory were employed for the binding free energy calculation at pH 4.5. The final values were determined by averaging both replica runs.

First, we investigated the structural stability of each system by analyzing the root mean squared deviations (RMSD) with respect to the corresponding crystal structure. The time evolution of RMSD for multiple runs is presented in Supplementary Figures S2-S3. It demonstrates that simulations are stable for each BabA isoform in its apo (see Supplementary Figure S2) and glycan-bound forms (Supplementary Figure S3). The handle region of BabA was more flexible as it was not fully crystallized, which led to an increase in the overall RMSD. Next, we determined the RMSD of each BabA–Le^b complex isoform by excluding the handle region and computed the corresponding probability distribution by combining all three replica runs (see Figure 2). In addition, each complex was compared with the respective free form. The average RMSD values for multiple runs were computed and listed in Table 1. It is evident from Table 1 that the average value varies between 2.96 Å and 3.60 Å for all apo isoforms, while the RMSD value decreases upon complexation, with the values ranging from 2.27 to 3.16 Å for all four complexes. It is evident from Figure 2A-D and Table 1 that compared to apo, BabA–Le^b complexes got relatively stabilized. However, for Iso4, a slightly higher RMSD was observed. These results indicate that the interaction between Le^b and BabA leads to overall structural stability. All complexes except Iso3 showed similar distribution for apo and bound states. In Iso3, the apo form showed a stable conformation across all three independent runs, which got relatively flexible after the binding of Le^b. In the apo case, the distribution was narrow with a sharp peak, while a broader distribution was observed when a complex was formed. However, the peak shifted to a lower RMSD value compared to apo Iso3.

Next, we calculated the RMSD of the ring atoms of Le^b in the free and protein-bound states concerning the initial coordinates. The corresponding RMSD distribution is displayed in Figure 2E. A broad distribution was observed for the free glycan with a mean value of ∼0.9 Å. In the case of the Iso3 complex, a very sharp unimodal distribution with a peak at ∼0.5 Å was observed, indicating a firm binding. A similar distribution was observed for Iso2/Le^b. On the other hand, moderate flexibility was observed for the glycan when bound to Iso1 and Iso4. In the case of Iso1, the most probable peak is located at 1.0 Å, with a secondary peak at ∼2.6 Å. In the case of Iso4, the RMSD distribution is characterized by two peaks, situated at ∼0.5 Å and ∼1.1 Å. The degree of flexibility of the glycan in the protein-bound form may reflect its binding affinity, as discussed in the later section.

Next, the average values of the radius of gyration (RoG) were computed to check the compactness of the protein and are listed in Table 1. The time evolution of RoG of each isoform in free and ligand-bound states is shown in Supplementary Figures S4, S5. It is evident from Supplementary Figure S5 that in the case of Iso4, run2 showed distinct time series compared to other runs. It is worth mentioning that the deviation in RoG was mostly because of the flexible handle region, which is located far from the binding region. Overall, the calculated average RoG values indicated no significant changes for each isoform. Next, we computed the SASA and found that the SASA value of apo variants was more significant than their respective complex form, as seen in Table 1 and Supplementary Figures S6, S7. The distribution of RoG and SASA for apo and complexes are shown in Supplementary Figure S8. These observations agree with the RMSD profile, suggesting that non-polar amino acid residues form interactions during complexation with the glycan. In addition, it suggests that the overall stability of BabA is enhanced upon complexation with the glycan, indicating that the hydrophobic core of BabA is more protected from the external environment upon binding with Le^b (see Supplementary Figure S8).

Finally, to provide further insights into the flexibility of different regions of BabA, the root means square fluctuations (RMSF) of C_α atoms were computed. The difference RMSF (∆RMSF) plot of the complex from the respective apo conformation was computed and shown in Figure 3. The N-terminal domain of BabA was relatively more flexible, and it corresponds to the handle region. However, this region is located far from the binding site and can be ignored for our current study. Furthermore, it is evident from Figure 3 that among different significant loop regions of BabA, DL2 of apo Iso3 exhibited large fluctuations compared to the corresponding complex. This indicates that the DL2 loop underwent structural changes upon interaction with Le^b and could play a crucial role in the molecular specificity of Le^b across isoforms of the H. pylori adhesin BabA.

Overall, Iso3 showed differences in dynamics and flexibility upon complexation, especially near the binding site, compared to the other three isoforms. This structural change may influence the Le^b recognition process.

The flexible glycosidic linkage of Le^b yields several conformations, which may play a critical role in protein-glycan recognition. The flexibility was investigated by calculating the free energy surface (FES) using φ/ψ angles of each glycosidic linkage and shown in Figure 4. The flexibility of each linkage stays similar for free and bound states. Although the conformational sampling across the dihedral space remained almost unchanged after Le^b binding, a minor variation for specific linkages, such as [Gal (5)β1-4Glc (6)] and [GlcNAc(3)β1-3Gal (5)], was observed. Both linkages displayed higher flexibility compared to others. Furthermore, linkages such as Fuc (4)α1-4GlcNAc(3), Fuc (1)α1-2Gal (2), and Gal (2)β1-3GlcNAc(3) stayed at a single minimum along with a very narrow conformational basin irrespective of free or bound states. The GlcNAc(3)β1-3Gal (5) linkage showed two equiprobable minima around (−80, −150) and (−80, −80), which got diminished after ligand binding. In the bound state, the flexibility of Le^b remained almost similar across all isoforms, indicating a stable conformation around the binding pocket. However, in the case of Iso4, the [Gal (5)β1-4Glc (6)] linkage showed two closely separated minima regions, indicating a slightly different conformation.

Another important parameter for characterizing the carbohydrate conformation is the ring geometry, for which we evaluated the puckering profile for each monosaccharide. We used the Cremer–Pople parameter (θ) to construct the conformational free energy to show the transition from ⁴C₁ to ¹C₄. In our study, no monosaccharide was able to complete the chair (⁴C₁)-to-chair (¹C₄) transition, except for glucose (6) and galactose (5) in the free system (see Figure 5). However, other galactose, such as Gal (2), stayed around the ⁴C₁ conformer, indicating a similar rigidity as other studies (Alibay et al., 2018). In the case of β-glucose, the complete transition was observed in the free simulation where the ¹C₄ state is energetically much higher than the other chair form ⁴C₁ (∼2 kcal/mol). On the contrary, in the case of Iso4, both the pucker states were sampled in an equiprobable way. The flexible nature of this ring conformation was also reflected in the dihedral space, as shown in the FES of the [Gal (5)β1-4Glc (6)] linkage (Figure 4). An earlier study also reported a change in the puckering state for glucose in a microsecond timescale (Alibay and Bryce, 2019). Fucose residues stayed only in the ¹C₄ state, as observed in another protein–carbohydrate recognition study (Roy et al., 2020a). Overall, Le^b displayed similar flexibility across all four isoforms.

A thorough structural investigation was conducted to gain insights into the dynamicity of the crown region of BabA, which mainly interacts with the glycan molecule. This region is comprised of several key loops, such as CL1-3 and DL1-2, which are the key functional regions for the glycan binding. The principal component analysis (PCA) of BabA was conducted, and the free energy landscape was generated using the first two principal components (see Figure 6). The apo structures of all four BabA isoforms showed several co-existing minima connected via shallow energy barriers, suggesting a very dynamical system in solution. However, upon Le^b binding, these regions shrunk to a certain extent, indicating the stability of the protein–carbohydrate complexes, especially for Iso1 and Iso2. In addition, for Iso4, both in apo and bound states, several minima regions were separated by a high energy barrier (more than 4 kcal/mol). However, the Cartesian coordinate-based PCA fails to provide a correct separation of internal and overall motions. The projection on the low-dimensional space offers only limited information about the multidimensional energy landscape (Krivov and Karplus, 2004).

We also generated the porcupine plot corresponding to PC1 displaying the directions of movements (see Figure 7). In Figure 7, the arrow corresponds to the direction of the movement, while the length of the arrow represents the movement’s strength. Except for the terminal region, most regions showed stability or minor movements. However, the crown region or the carbohydrate-binding region displayed the least movement for apo and complexes except for Iso1. In the case of apo Iso1, DL2 showed a movement toward the binding region. However, in the case of the bound state, the direction of the movement was completely reversed. A similar trend was observed for apo Iso3. Apart from the DL2 loop, some parts of the CL3 loop showed upward movements or toward the crown region. However, this cysteine–cysteine loop region is far from the binding region. So, the flexible movement of this region may not affect the glycan binding.

Further, the structural stability of the crown region was investigated by estimating the number of intra-residue hydrogen bonds (H-bonds). We evaluated the occupancy of H-bonds of the crown region using the Bridge2 software (Siemers and Bondar, 2021). The H-bonds heatmap corresponding to the crown region is shown in Figure 8, where only H-bonds with more than 70% occupancy are listed. For all cases, two to three stable H-bonds were observed, which may determine the structural integrity of the region. The key residues obtained from this analysis are primarily located in the DL2 loop and its surrounding region.

Furthermore, we also estimated the joint occupancy of H-bonds, which reports true when all the estimated H-bonds are present and false otherwise (Siemers and Bondar, 2021). The joint occupancy for Iso2 was 83.9%, the highest among all isoforms (Iso1 ∼61.5%, Iso3 ∼56.3%, and Iso4 ∼57.3%). Despite the high occupancy of different hydrogen bond pairs, the DL2 region was more flexible across several isoforms of BabA.

Furthermore, to provide a deeper understanding of the dynamicity of significant loop regions, an inspection of loop regions was conducted by estimating the RMSD. The time evolution of all three cysteine–cysteine loops (CL1, CL2, and CL3) RMSD was estimated over all three MD runs and shown in Supplementary Figures S9–S11. In all three replicas, no major conformational changes were observed. This behavior was also true for the diversity loop 1 (DL1), as shown in Supplementary Figure S12. We performed the backbone dihedral angles principal component analysis (dPCA) to extract more information about the conformational changes. CL1 and CL3 loops, away from the binding region, showed a moderate conformational basin across all four isoforms (see Supplementary Figures S13 and 14). For apo structures, the principal minima of CL1 were separated from the rest of the conformational spaces by a high energy barrier. However, after the complex formation, these regions were connected by the shallow free energy barrier. On the other hand, for CL3, the conformation variability remains unchanged before and after the complex formation. As these regions are far apart from the carbohydrate-binding site, the dynamics of these regions may not interfere with Le^b binding.

Among all the cysteine-cysteine and diversity loops, the major loops are CL2, DL1, and DL2, which play a crucial role in ligand binding (Moonens et al., 2016). For CL2, the free energy surface was constructed using the first two dPCs, as shown in Figure 9. For all the cases, the free energy surface was characterized by two equiprobable minima, separated by a high energy barrier (Figure 9). However, the minima regions shrunk in the bound state compared to the apo form. The key change was observed in Iso3, where the apo BabA was sampled over four different areas separated by a high energy barrier. However, this intrinsic flexibility diminished after binding with Le^b, indicating a better ligand binding phenomenon. Furthermore, to investigate the conformational dynamics of DL1 and DL2 in detail, dPCA was conducted, and the resultant FES is shown in Supplementary Figures S15 and S26. The DL1 loop showed a very stable profile throughout the replica runs, as is evident from the time series of RMSD (Supplementary Figure S12). However, flexibility was observed for all apo and complexes. Among all isoforms, Iso3 showed a more stable conformational space, and the corresponding free energy surface was characterized by two principal minima in the apo and complex forms. It should be noted here that the location of the minima region shifted after the complexation.

In DL2, the flexibility was observed only in cases of apo forms of Iso1 and Iso3, which decreased after the binding of Le^b, as evident from Supplementary Figure S16A and C (for apo cases) and Supplementary Figure S16E and G (for complexes). Similar behavior was observed from the RMSD plot of the DL2 loop in apo and complex forms (see Figure 10A). In apo cases, all three replica runs showed flexibility over the entire 1.0 µs for Iso1 and Iso3. The RMSF analysis (Figure 3) also supports the same observation.

To investigate the reason behind the shift in RMSD, we conducted the DSSP analysis using the cpptraj module of AmberTools19 (Roe and Cheatham, 2013), where specifically a change in β-sheet and turn was observed. The time evolution of the secondary structure in Iso1 for three replica runs is shown in Figure 10B. The flexible behavior of the DL2 loop in apo cases resonates with the percentage change in the turn of this segment. However, DL2 did not form any β-sheet in the apo simulation of BabA isoforms. In contrast, after the complexation with Le^b, a change in coordination of sheet and turn was observed, agreeing with the RMSD profile. Specifically, when the RMSD was stable, a steady percentage of turn was observed, which diminished with the fluctuation. Contrarily, we observed the appearance of β-sheet with the increase in RMSD. Therefore the occurrence of β-sheet was inversely proportional to the decrease in the percentage of turn in the loop structure. The location of the β-sheet could be seen in the representative structure from the respective segment of RMSD (Figure 10C, Iso1). Similarly, the 3D conformation was also shown for Iso3 to indicate the change in the structure of the DL2 loop region (Figure 10D).

Overall, our study revealed that the intrinsic flexible nature of Iso1 persists after the glycan binding. In the case of apo Iso3, the flexibility of the crown region was relatively higher than the other three isoforms, which got drastically reduced after the complex formation. It indicates a stable Le^b binding for Iso3 compared to the other isoforms.

To elucidate the binding mechanism of Le^b against variants of BabA, we calculated the binding free energy for each complex and its different components (see Table 2). We calculated the average binding free energy over the three replica runs and found −37.70, −34.96, −40.26, and −38.28 kcal/mol for Iso1, Iso2, Iso3, and Iso4, respectively. Overall, it suggests that Le^b has a stronger affinity toward Iso3 of BabA than the other three isoforms. In contrast, Iso2 showed the lowest affinity for Le^b (see Figure 11A). It is to be noted that the electrostatic contribution (ΔE_elec) was more favorable than the van der Waals (ΔE_vdW) interactions for all complexes. ΔE_elec varied between −46.72 kcal/mol and −76.25 kcal/mol, while ΔE_vdW varied between −41.77 kcal/mol and −44.71 kcal/mol. The total favorable interaction energy (ΔE_MM), which comprises ΔE_vdW and ΔE_elec, was found to be more favorable for Iso3 (ΔE_MM = −120.96 kcal/mol) than other complexes. The lowest value was obtained for Iso2 (ΔE_MM = -89.65 kcal/mol). In the case of Iso3, both ΔE_vdW (−44.71 kcal/mol) and ΔE_elec (−76.25 kcal/mol) were more favorable than other complexes. Although ΔE_elec was more favorable to the complexation compared to ΔE_vdW, the net electrostatic component, i.e., the sum of intermolecular electrostatic and polar solvation free energy (ΔE_elec+ ΔG_pol), was unfavorable to Le^b binding for all complexes. It suggests that the intermolecular van der Waals interactions mainly governed the binding of Le^b. Overall, for Iso3, both components (ΔE_elec and ΔE_vdW) favor the BabA–Le^b complexation the most compared to other isoforms, resulting in a better affinity for Le^b.

It is worth noting here that the binding free energy across the triplicate runs varied to a certain extent for all except Iso2. To investigate the origin of variation, we calculated the RMSF of the bound Le^b (Supplementary Figure S17). It is evident from Supplementary Figure S17 that in the case of Iso2, Le^b showed the least flexibility compared to the other three isoforms. Consequently, Iso2 showed almost similar binding free energy in all three replicas. However, Iso2 showed the lowest binding free energy among all four cases. It may indicate the fluctuation of a few monosaccharides can influence the protein-glycan binding free energy. However, the key residue, such as fucose, displayed the least flexibility across all four isoforms that govern the protein-glycan recognition.

The percentage contribution of each monosaccharide to the total binding was plotted and shown in Figure 11B. Fuc (1) and GlcNAc(3) contributed more favorably to binding than other monosaccharides for all complexes. Furthermore, our study revealed that Fuc (1) contributed ∼40–53% to the total binding free energy among monosaccharides in each isoform, and GlcNAc(3) contributed ∼15–27%. In Iso2, Fuc (1) contributed the highest among other isoforms, accounting for 53%, while the contribution of GlcNAc(3) was the lowest (∼15%) compared to other isoforms. It agrees with a previous study where it was reported that Fuc and GlcNAc contributed more toward the total binding free energy that favors the cholera toxin (CT)/LewisY complexation (Roy et al., 2020a).

Next, we computed the contribution of each amino acid of BabA to the total binding free energy to find out the hotspot residues. All amino acids contributing more than −1.5 kcal/mol of energy are listed in Supplementary Table S1 and graphically shown in Figure 12. We noted that the highest contributing residues toward Le^b binding (>−3.5 kcal/mol) come from Iso3 and Iso4 in agreement with the binding energetics given previously, especially the residue S440 from Iso3 and S247 from Iso4. Similarly, S244 of Iso1 contributed high in the binding with a -3.10 kcal/mol value. In contrast, in the case of Iso2, the highest contributing residue, A240, showed a contribution of -2.52 kcal/mol, which also agrees with the lesser binding affinity than others toward Le^b. The other contributing residues involved in the binding included T246, G191, V243, N194 (Iso1), A240, Q235, A241, T243, G191 (Iso2), V243, G191, T246, N194 (Iso3), and G192, T249, Q246, N195 (Iso4). Overall, these results suggest that primary interacting residues for recognizing Le^b come from CL2 (G191 and N194) and DL2 (V243, Q235, and S244) loops, which agrees with the above mentioned RMSF and free energy surface analysis. Moreover, these results are consistent with the experimental findings in explaining the significant role of CL2 and DL2 loops in recognizing glycans for the attachment of BabA of H. pylori (Moonens et al., 2016).

In addition, we calculated the contributions of each monosaccharide to ΔG_bind for all complexes and listed them in Supplementary Table S2. Fuc (1) contributed more favorably to the protein-glycan binding than other monosaccharides for all complexes. The contribution of Fuc (1) ranged from -7.65 to -8.33 kcal/mol and was more favorable for Iso3 and Iso4 than Iso1 and Iso2. In addition, GlcNAc also contributed more favorably to the binding. The contribution of GlcNAc was found to be the highest for Iso3 (-5.33 kcal/mol), followed by Iso4 (-5.11 kcal/mol), Iso1 (−4.77 kcal/mol), and Iso2 (−2.25 kcal/mol). However, in Iso2, Gal (1) contributed more significantly (−3.19 kcal/mol) than other monosaccharides. Therefore, it suggests that GlcNAc is the next significant monosaccharide favoring the association between Le^b and BabA. Furthermore, it is revealed from Supplementary Table S2 that the electrostatic contribution (T _ele) of Fuc (1) to the total binding was more favorable than T _vdW for each complex ranging from -12.04 kcal/mol to -19.83 kcal/mol. However, T _ele was overcompensated by the disfavourable polar solvation of energy (T _GB) in each complex. Therefore, the key contributing component of Fuc (1) was vdW (T _vdW) which favored the binding.

Subsequently, we computed the H-bond occupancy formed between Le^b and BabA for all isoforms, as listed in Supplementary Table S3. Supplementary Figure S18 shows that isoforms 1, 3, and 4 showed high H-bond frequency compared to Iso2. The critical BabA–Le^b H-bond interactions are noted in Supplementary Table S3. It suggests that Fuc (1) formed multiple H-bonds with BabA residues (G191, T246, N194, and C189), with the occupancy varying between ∼67 and ∼90% for Iso1, while it ranged between ∼62 and ∼88% for Iso3. Similarly, T249, C190, G192, and N195 from Iso4 interacted with Fuc (1) via H-bond, with the occupancy varying between ∼56 and ∼86%. In the case of Iso2, Fuc (1) formed H-bond with C156, N161, G158, and T210, with an occupancy ranging from ∼40 to ∼87%. Compared to other isoforms, the number of residues forming H-bonds with Le^b was less in Iso2. It is evident from Supplementary Figure S18A that the average number of H-bonds in isoforms 1, 2, and 3 was estimated as 8, while on average, 6 H-bonds were observed for Iso2.

Further, we computed the hydrophobic contacts between BabA–Le^b as shown in Supplementary Figure S18B. The average value of the hydrophobic contact was the highest for Iso2 (∼12), while for other isoforms, on average, five to six contacts were formed. Altogether, it suggests that the electrostatic interactions favored the association between Le^b and BabA for isoforms 1, 3, and 4, while the hydrophobic interaction favors the complexation in Iso2.

The abovementioned analysis was complemented by plotting the interaction diagram of BabA/Le^b for all four complexes (see Supplementary Figure S19). As seen from Supplementary Figure S19C, in Iso3, Fuc (1) formed H-bond interactions with G191, C189, T246, and N194. Similarly, GluNAc(3) formed H-bonds with D233, Gal (5) with S244, and Fuc (4) with N208. The spike representation was used to depict the hydrophobic interactions. Overall, Supplementary Figure S19A–D indicate fewer interactions in Iso2 than other isoforms, which agrees with the aforementioned binding affinity. Overall, the interaction profile suggests a strong hydrogen bonding and hydrophobic interactions between BabA and Le^b for Iso3 compared to other isoforms of the BabA family, which could unravel the molecular specificity of hexasaccharide Le^b against isoforms of H. pylori adhesin BabA.

Finally, we estimated the binding affinity of Iso1 (17875 strain) at acidic pH (4.5) to elucidate the effect of pH on glycan binding. The RMSD of BabA for both runs was estimated and shown in Supplementary Figure S20. Both runs attained equilibrium after 100 ns, and the last 150 ns were used to calculate the binding free energy using the MM/GBSA protocol. The binding affinity of Iso1 against Le^b in the acidic condition was estimated as −34.20 kcal/mol, which is ∼4 kcal/mol lower than what was evaluated in neutral solution (−37.70 kcal/mol). This finding was also supported by the recent experimental results where BabA senses the change in pH in the environment (Bugaytsova et al., 2017). In acidic conditions, the binding affinity of BabA against Le^b was dramatically reduced for the strain 17875, agreeing with our MM/GBSA calculation.

Details of the binding free energy and its components are listed in Supplementary Table S4. An increase in the polar solvation energy and decrease in the van der Waals component mainly lowered the binding free energy. However, Iso1 lacks the key salt bridge interaction between D198 and R207, which plays a significant role in acid sensitivity; still, it senses the change in environment as expected by the experimental findings (Hatakeyama, 2017). The short α helix, or the key-coil (residue 199-202) region located in the DL1 loop, is present in Iso1. However, the central residue in the key-coil region for acid sensitivity and binding affinity, i.e., P199, was not present in the 17875 strain, which may reduce the sensitivity against the pH change (Bugaytsova et al., 2017). The critical residues involved in the binding are shown in Supplementary Figure S21, indicating a similar pattern relative to the neutral condition. Irrespective of pH, the key residues are the same, with a minor variation determining the strength of the binding. In acidic conditions, direct contributions from the amino acids of the key-coil region were not observed.