Molecular docking is a powerful method for high-throughput virtual screening of small-molecule inhibitors and predicting the ligand-protein interactions in computer-assisted drug design. It consists of three main types: rigid docking, semi-flexible docking, and flexible docking, which are categorized in ascending order of protein flexibility. Semi-flexible docking was chosen for its balance between computational efficiency and accuracy, allowing modifications in the conformations of ligands and residues within a defined area.

The crystal structure of the CA II (PDB: 12CA) was sourced from the RCSB library and does not include post-translational modifications (Berman et al., 2000). Both PFOS-A and PFOS-B were fully optimized utilizing the Gaussian09 software (Frisch MJT et al., 2009) at the density functional theory level, employing the B3LYP function with the London-dispersion correction by Grimme and the def2-SVP basis set (Schäfer et al., 1994; Grimme et al., 2011).

Autodock Vina was utilized for the docking of PFOS with CA II (Trott and Olson, 2010). A 10 × 10 × 10 nm box, set up using default parameters, encompassed almost the entire active surface of CA II, centered on the typical binding pocket identified in previous research. This pocket, delineated by hydrophobic residues (Ile91, Val121, Phe131, Leu141, Val143, and Leu198) and hydrophilic residues (Asn62, His64, Asn67, and Gln92), is crucial for rapid catalytic conversion of CO₂ to bicarbonate. The energy range was set to 5.0 kcal/mol, and the exhaustiveness parameter was configured as 22 (Agarwal and Smith, 2023). The optimal pose with the best scores within multiple iterations of the docking process was selected as the initial structure for the subsequent MD simulation, according to the affinity assessment through scoring function and binding mode analysis (Tanchuk et al., 2016).

The frontier molecular orbital analysis was conducted to assess the docking reactivity of these molecules due to the significance of HOMO and LUMO (E_HOMO and E_LUMO) energy gaps in chemical reactivity. Various physical and chemical parameters, including chemical hardness (η), softness (s), and electrophilicity index (w), were calculated and analyzed.

Based on the initial structures obtained from the molecular docking results, MD simulations were conducted using the GROMACS 2016 software package, (Abraham et al., 2015), and the proteins and small molecules were parameterized by the all-atom CHARMM36 and CHARMM general force field (CGenFF), respectively.

The complex was put in a 10 × 10 × 10 nm periodic cubic box and then filled with the TIP3P water molecules. (Mark and Nilsson, 2001). Then, the chloride ions were added to neutralize the system charge. During the energy minimization process, the unreasonable interactions and atomic collisions were removed using the steepest decent and conjugated gradient methods. Next, the temperature of each system was incrementally increased from 0 K to 300 K during the 1 ns simulation using an NVT ensemble. The temperature of the protein-ligand and solvent-ion groups was controlled separately using a Berendsen thermostat to obtain more reliable results. After that, a Parrinello–Rahman barostat was used to maintain the pressure of the systems at 0.1 MPa in NPT equilibration for 1 ns. Finally, MD simulations were conducted for 100 ns, and the last 10 ns trajectory was taken every 10 ps for subsequent analysis.

During the simulation, the length of the hydrogen bonds was constrained by LINear Constraint Solver algorithm. The cut-off for short-range nonbonded interactions was set at 1.0 nm, while for long-range electrostatic interactions, the Particle Mesh Ewald algorithm was employed without an explicit cut-off. The determination of cut-offs was guided by empirical practices widely adopted in previous studies (Essmann et al., 1995; Kabedev et al., 2023; Sun et al., 2023). Newton’s classical equations of motion were integrated using the Verlet leapfrog algorithm. (Hockney et al., 1974). The stability of the systems was checked by the root mean square deviation (RMSD) of the protein backbone.

To further assess the interaction between PFOS and CAII, the binding energy was calculated and partitioned into different parts according to the type of interaction by the MM-PBSA method implemented in the gmx_MMPBSA software package (Miller et al., 2012; Suri and Naik, 2015). It should be noted that the binding free energies calculated by the MM-PBSA method is system-dependent and sensitive to the chosen parameters, which is somewhat less accurate compared to the free energy perturbation or thermodynamic integration (Breznik et al., 2023). However, the computationally efficient MM-PBSA method has proven to be reliable in the field of computational drug design and evaluation (Genheden and Ryde, 2015; Wen et al., 2019).

This calculation was conducted using a set of snapshots extracted from the final 10 ns trajectory of the MD simulations, sampled at 10 ps intervals. The binding free energy calculations in the MM-PBSA method can be concisely represented by the following equations: Δ G = G c o m p l e x − G r e c e p t o r + G l i g a n d = Δ H − T Δ S (1) Δ G = Δ E g a s + Δ G s o l − T Δ S (2) where, Δ E g a s = Δ E b i n d e d + Δ E n o n b i n d e d (3) Δ E b i n d e d = Δ E b i n d + Δ E a n g l e + Δ E d i h e d r a l (4) Δ E n o n b i n d e d = Δ E e l e + Δ E v d w (5) and, Δ G s o l = Δ G p o l a r + Δ G n o n p o l a r = Δ G P B + Δ G S A (6) where, Δ G S A = N P T E N S I O N × Δ S A S A + N P O F F S E T = γ ⋅ Δ S A S A (7)

Within these equations, G_complex, G_receptor, and G_ligand are the free energy of the complex system, CA II and small molecule inhibitors, respectively, collectively contributing to the binding free energy (∆G). This free energy can be further partitioned into the sum of the binding enthalpy (∆H) and the conformational entropy change upon inhibitor binding (-T∆S), as depicted in Eq. 1.

The enthalpic term, denoted as ∆H, can be formulated as the sum of ∆E_gas and ∆G_sol, as elaborated in Eq. 2. ∆E_gas accounts for the energy changes computed through molecular mechanics in the gas phase, while ∆G_sol pertains to solvent-related energy contributions.

In Eqs 3–5, the term ∆E_gas is divided into two separate contributions: ∆E_binded and ∆E_nonbinded. ∆E_binded encompasses energy contributions arising from binding, angle, and dihedral alterations, while ∆E_nonbinded can be further partitioned into the ∆E_ele and the ∆E_vdw.

The computation of the solvation-free energy involves two components: ∆G_polar, assessed via the Poisson–Boltzmann approach, and the ∆G_nonpolar, estimated through a SASA-only model (Eq. 6), assuming proportionality to the solvent accessible surface (∆SASA) outlined in Eq. 7.

Notably, the entropy change (∆S) effect is ignored due to computational constraints and the substantial uncertainty it introduces into the results, stemming from considerable errors during the calculation process. Several studies have shown that the MM-PBSA method can give high calculation accuracy without accounting for ∆S (Qin et al., 2010; Hou et al., 2011; Guo et al., 2021).

The trajectory obtained from the MD simulation was analyzed using auxiliary tools provided within the GROMACS 2016 package. The interaction between small molecules and CA II receptors is analyzed by Protein-Ligand Interaction Profiler (PLIP). Furthermore, the occupancies of the intermolecular hydrogen bonds within the PFOS systems were calculated using the Visual Molecular Dynamics (VMD) 1.9.3 software. With the focus on direct interaction between PFOS and CAII, the trajectory used in hydrogen binding analysis have been pre-processed to remove the waters. Hydrogen bonds were identified based on established criteria: an acceptor-hydrogen-donor angle exceeding 135° and an acceptor-hydrogen atom distance less than 3.5 Å.

In addition, a principal component analysis was conducted to delineate the dominant motion of the CA II and the correlated motions among the atoms extracted from the MD trajectories within the PFOSs system (Ma et al., 2017; John et al., 2020). Notably, an intermolecular motion was excluded by eliminating protein translation and rotation via the auxiliary tools available in the GROMACS suite. The eigenvectors, which signify the directions of the atomic motions, along with their corresponding eigenvalues that indicated their magnitudes, were calculated by evaluating the covariance matrix of the α-carbon atoms of the CA II protein. This process aimed to elucidate the impact of ligands on the conformational distribution of the protein (Mu et al., 2005; Papaleo et al., 2009). Subsequently, the first two principal components (PC1 and PC2) were utilized as reaction coordinates to construct the free energy landscape, and the free energy values were calculated with Eq. 8. G = − k B T × ⁡ ln ⁡ P (8)

Where k_B represents Boltzmann constant, T stands for the simulation temperature, and P is the probability of the conformational distribution.

The initial binding structures of CA II with each PFOS compound were obtained through molecular docking. Firstly, to validate the reliability of the docking method, redocking experiments were performed on the complex of CA II with N-[(4-methylphenyl)methyl]-thiophene-2,5-disulfonamid (PDB: 1BNA, inhibitor: AL5) (Boriack-Sjodin et al., 1998). The results demonstrate that the conformation of the complex obtained by docking method matches with the native crystal structure located in the active pocket of the CA II, as shown in Supplementary Figure S1, with an RMSD of 0.78 Å. Furthermore, the interactions of the docked pose, assessed using the PLIP program, (Adasme et al., 2021), are consistent with the research conducted by Sjodin, (Boriack-Sjodin et al., 1998), further validating the accuracy of the docking method. The correlative analysis is shown in Supplementary Figure S2. Having validated the docking method, the two PFOS molecules were then docked to the CA II using the same docking method (Figure 2). The conformations with the best score for both structures are in the pocket adjacent to the zinc ions.

The docking results reveal that both PFOS structures effectively bind to the active site on CA II as classical inhibitor AL5 used for docking verification. The chemical softness of the small molecules underscores their binding reactivity, with values of 7.58 for PFOS-A and 7.96 for PFOS-B, since a higher value of chemical softness typically correlates with increased molecular reactivity. The additional details on other chemical reactivity parameters are provided in Supplementary Table S1.

Analysis based on PLIP reveals that PFOS-A predominantly binds to residues ASN-62, HIS-64, ASN-67, ASP-72, PHE-131, GLY-132, THR-200, and PRO-201, while PFOS-B exhibits binding interactions with TRP-5, ASN-62, HIS-64, ASN-67, GLN-92, HIS-94 and THR-200. These residues are primarily located on the hydrophilic surface of the protein, attributed to the presence of multiple hydroxyl groups in these glycosyl molecules. However, the neglect of the structural flexibility inherent in the majority of protein structures in the docking method is a formidable obstacle to accurate binding analysis. Nevertheless, the docking conformations can still serve as valuable initial input structures for subsequent MD simulations and as a reference for the binding analysis. The off-target effects of PFOS have also been preliminarily assessed by the molecular docking method (Supplementary Figure S3). It should be noted that the activity of the non-target protein CA I may be slightly affected.

With the initial structures obtained by molecular docking, 100 ns MD simulations of the systems were performed to get a more detailed and accurate understanding of their behavior. The RMSD was calculated for the resides within 10 Å of the ligand to assess the structural stability throughout the 100 ns simulations. As illustrated in Figure 3A, both PFOS complexes reached equilibrium after 50 ns, exhibiting the RMSD values of less than 2.8 Å for the majority of the equilibration stage, which are more stable than the CA II without inhibitor (averaging 3.0 Å in the last 50 ns). This observation suggests that the presence of the PFOS small molecules stabilized the core backbone of the CA II. Notably, during the equilibration period, the PFOS-B system demonstrated slightly higher stability with lower RMSD values than the PFOS-A system, while the opposite trend was observed during the early stage. This discrepancy may be attributed to the disruption or instability of the hydrogen bonds present in the initial docking structure of the PFOS-A system. Nevertheless, since all three systems reached equilibrium, further analysis can be conducted based on the obtained MD simulation trajectory.

The Root mean square fluctuation (RMSF) was further calculated to measure the per-residue flexibility of the CA II, PFOS-A, and PFOS-B complex systems. As illustrated in Figure 3B, the closely overlapping profiles of PFOS complexes and CA II suggest that the inhibitor has little effect on the regional flexibility of the protein. The stable regions mainly occurred at the β-sheets with the RMSF less than 0.8 Å, while the helix domain showed higher RMSF values up to 2.5 Å.

The secondary structure of CA II is shown in Figure 3C, which consists of 10 stranded twisted β-sheet and 7 helices and the zinc ions presented in the center bottom of the protein chain surrounded by the β-sheet regions. According to the secondary structure, all helix regions with high RMSF values were far away from the binding pocket and had little effect on the binding process. On the contrary, the β-sheet, especially the C-sheet, D-sheet, and G-sheet segments, showed relatively lower RMSF values in complex with PFOS as inhibitors compared to the CA II system. These segments are identified as the key sites for the combination of PFOS. In addition, systems with PFOS-B as the inhibitor are found to have higher residue stability than PFOS-A, which is in accordance with the previous RMSD analysis.

Further RMSF analysis focused on residues within the active site of CA II (Figure 3D). Upon binding with PFOS, the flexibility of these residues generally decreased, particularly for residues 131PHE and 135VAL. These two residues have also been identified as key contributors to the protein-ligand interaction in subsequent energy analysis to accommodate the hydroxyl group of PFOS.

It has been shown that the PFOS effectively stabilize the CA II protein via binding with some key residues situated within the β-sheet region. However, the lack of a zinc-binding group distinguishes them from most existing CA II-targeting inhibitors. It is difficult for PFOS structures to bind directly to the catalytic zinc ion, which plays a critical role in the catalytic function of the enzyme. Hence, it is crucial to validate the exact binding affinities of these two small molecules with CA II, wherein the MM-PBSA method was employed. This calculation relied on the snapshots extracted from the final 10 ns of the equilibrium stage in the MD simulation of each respective system. The same methodology was also applied to the experimentally tested inhibitor AL5, employed in the docking validation, serving as the benchmark (Table 1).

The two PFOS structures show comparable binding affinity to the CA II protein (Table 1). The PFOS-B (−8.09 kcal/mol) shows a slightly lower ΔG compared to PFOS-A (−7.27 kcal/mol), a trend consistent with the RMSD and RMSF (Figure 3) results. From a thermodynamic perspective, these negative values corroborate the spontaneity of the binding process and the stability of the combination in the solvent. This level of binding stability aligns with findings in similar research, where a small glycosyl molecule interacts with the human sweet taste receptor, exhibiting a ΔG of −8.53 ± 1.78 kcal/mol, as validated through MD simulations and experimental data. (Mahalapbutr et al., 2020). However, the biological activity of PFOS-A and PFOS-B is relatively lower when compared to typical CA II inhibitors exemplified by the benchmark compound AL5. This difference in activity can be attributed to the absence of the zinc-binding group in the PFOS structure, which can form much more favorable electronic interactions to maintain the combination of the complex.

To further investigate the potential inhibition mechanism of PFOS on CA II, common features of the two structures are first identified by energy decomposition analysis. The total binding energy can be decomposed into different components including electrostatic energy, van der Waals energy, polar solvation-free energy, and non-polar solvation energy. Thus, relative importance of each interaction within the binding process can be assessed.

The dominant propellant for binding for both the PFOS-A and PFOS-B systems is the total nonpolar binding energy, which can be further separated into van der Waals energy and non-polar solvation-free energy (Table 1). Both these components contribute positively to the binding process. In addition, negative electrostatic energy was observed in all systems, with significantly lower values in both PFOS systems compared to AL5. This discrepancy can be attributed to the abundant hydroxyl groups in PFOS-A and PFOS-B (Figure 1), which exhibit substantial potential for establishing hydrogen bonding interactions. However, a significant unfavorable influence was detected concerning the total polar free energy, primarily attributed to the desolvation process, chiefly driven by the considerable polarity exhibited by the PFOS molecules. The formation of intramolecular hydrogen bonds by PFOS molecules may also impose limitations on their interactions with CA II molecules. The slightly weaker binding free energy of PFOS-A observed during dynamic simulations will be explained by following energy decomposition.

In order to elucidate the detailed interaction between PFOS molecules and the CA II protein, the identification of key residues critical for the binding process was pursued through energy decomposition analysis. The analysis encompassed residues within a 10 Å radius of the PFOS molecules. The ten residues exhibiting the most substantial contributions to the binding energy were ranked and presented in Figure 4. Key residues identified by contributions exceeding the threshold of −0.5 kcal/mol were summarized in Table 2 and Table 3. The detailed energy decomposition result is shown in Supplementary Tables S2, S3.

PFOS-A structure interacted with 5 key residues: TRP-5, PRO-202, PHE-20, PHE-131, and VAL-135 (Figure 4). Meanwhile, the PFOS-B structure interacted with 9 pivotal residues: PHE-131, VAL-135, LEU-198, PRO-202, LEU-204, TYR-7, THR-200, LEU-141, and TRP-5, which is nearly twice the key residues that participated in the combination of PFOS-A system. Based on the analysis of the docking conformation of the two systems, the reason for the fewer key residues in PFOS-A may be attributed to its stick line backbone structure in the combination process, which makes it less possible to access the residues located at the entrance of the docking pocket to form van der Waals interaction with them. The fusiform structure of PFOS-B makes it more suitable to block the docking pocket with more accessible residues. In addition, the van der Waals energy contribution of the key residues is significantly larger compared to the electrostatic energy, which is consistent with the previous analysis in Table 1. The TRP-5 in the PFOS-A system and the PHE-131 in the PFOS-B system exhibited the most favorable ΔG_total values. For both PFOS structures, residues PHE-131, VAL-13 were found to play a significant role in the van der Waals interactions at the hydrophobic side of the binding pocket entrance of the CA II protein. These findings are consistent with the prior studies suggesting the importance of these residues in the binding of the hydrophobic tail in Zinc binding group inhibitors. (Alterio et al., 2012).

The electrostatic interaction also plays an important role in inducing the complexation of PFOS molecules and CA II, in addition to van der Waals and the non-polar solvation effect. Further characterization of this kind of interaction anchoring the PFOS small molecules into the binding pocket was conducted. The occupancy of the direct intermolecular hydrogen bonds was statistically analyzed over the last ten MD simulations using the PLIP and Visual Molecular Dynamics (VMD) program. (Humphrey et al., 1996). The percentage occupancy of hydrogen bonds in the PFOS-A and PFOS-B systems are presented in Supplementary Tables S4, S5. The stable hydrogen bonding has been summarized according to the standard for stable hydrogen bonding, that the hydrogen bonding donor-acceptor interaction occupancy is greater than 10% (Table 4; Table 5). As expected, both PFOS systems showed the formation of hydrogen bonds with polar residues within the C-sheet, D-sheet and G-sheet of the CA II protein. However, it can be seen that the number and stability of these hydrogen bond formations are different for the two systems.

In PFOS-A system, three main hydrogen bonds formed, including H-N2 of ASN-67 with O4 and H-N2 of TRP-5 with adjacent O72 and O70, respectively (Figure 5). This observation is corroborated by the aforementioned energy decomposition analysis, which revealed the negative electrostatic energy contributions from these residues, particularly notable for 5TRP with a value of −1.94 kcal/mol. Conversely, in the PFOS-B and CA II complex, the formation of hydrogen bonds was chiefly facilitated by PRO-201, GLN-92, and GLN-136. Notably, more hydrogen bonds are present in system PFOS-B compared to PFOS-A, suggesting that PFOS-B has a stronger electrostatic interaction, in agreement with free energy calculations. However, the most stable H-bond was detected in the PFOS-A system, which persisted for 39.6% of the simulation time. While, the majority of the other hydrogen bonds endured for less than 20% of the simulation duration. In addition, the overall analysis of the H-bonds occupation data in the table indicates general consistency with the results of molecular docking, which also shows the H-bonds formation of HIS64, ASN67, THR200, and PRO201 in PFOS-A system and ASN62, ASN67, GLN92 and GLN 136 in PFOS-B system.

Considering the presence of multiple hydroxyl groups in PFOS, the formation of hydrogen bond largely depends on the structure and conformation of PFOS, which may be dynamic. Whereas, the ASN67, GLN92, GLN135 and PRO201 showed as the most potential residue to form hydrogen bond with PFOS.

Based on the above analysis, it can be concluded that the van der Waals interaction domain governs the binding of PFOS molecules, while the electrostatic energy contribution resulting from the formation of hydrogen bonds also has a crucial influence on the binding of PFOS molecules to the hydrophilic residue.

By integrating the comprehensive analyses of both van der Waals and electrostatic interactions, the elucidation of the binding site and interaction mechanism between PFOS small molecules and the CA II protein has been elucidated in Figure 6A. Both PFOS structures effectively bind to the entrance of the pocket that accommodates the Zn ion within CA II. While, the double-branched fusiform PFOS-B structure is better at blocking the pocket. The small molecules bind to the hydrophobic surface of the protein through van der Waals interactions and form hydrogen bonds with the hydrophilic surface through electrostatic interactions. This dual interaction strategy effectively contributes to the binding of PFOS, which blocks the carbon dioxide conversion pathway and then results in the inhibition of CA II activity.

Based on the characterization of the binding pocket, the free energy landscapes (FEL) of the three systems, CA II, PFOS-A, and PFOS-B, were analyzed to further investigate the impact of the small molecules on the overall conformational dynamics of the CA II protein and also furnish a qualitative description of their inhibitory potential. The FEL was constructed by projecting the 1000 snapshots from the MD simulation onto the first two principal components, which were employed as the coordinates, as illustrated in Figure 6B.

In the FEL of the PFOS-B system, a distinct blue region indicates the high probability of the existence of a stable conformation. In contrast, the two prominent conformational regions were identified for CA II in the absence of an inhibitor molecule, indicating the stabilizing effect of PFOS-B binding. In the case of the PFOS-A system, its conformational space is characterized by a large region containing two distinct energy traps. This means that the PFOS-A system has a comparatively lower stability. In addition, mirroring the results of the RMSF analysis, the FEL of the PFOS-A system shows numerous smaller basins that are absent from the FEL graph of the PFOS-B system, implying a higher flexibility of the complex binding with PFOS-A.