Challenging Dissimilatory Adenosine 5’-Phosphosulfate (APS) Reductase with Assimilatory APS Reductase Inhibitors by Docking and Molecular Dynamics Simulations: A Computational Study

Abstract

Dissimilatory adenosine 5’-phosphosulfate reductase (APSrAB) is a metalloenzyme of the metabolic pathway of sulfate-reducing microorganisms, which generates H₂S in oil production wells, causing losses in the oil industry. A set of 32 compounds derivatives of adenosine 5’-phosphosulfate (APS), previously evaluated as inhibitors of the assimilatory adenosine 5’-phosphosulfate reductase (APSr), were challenged in this current work by in silico molecular docking and molecular dynamics (MD) techniques as potential inhibitors of the APSrAB (APS is substrate of both enzymes). From this set, 20 compounds showed the highest affinity by the APSrAB (tendency to remain in the binding site pointed out by the docking experiments). The residues of the active site (71-398 region) had hydrogen bonds with a lifetime of more than 10.00% (mainly Arg265). They were responsible for the binding of these ligands over time, while the binding energy (ΔG_binding) values showed the energetic contribution of these residues to the stabilization of the APSrAB-ligand complex (mainly Arg265). Thus, compounds 13d, 14a, 14d, 16c, and 16d showed ΔG_binding values more favorable than -30.00 kcal mol^-1 and had more affinity for the enzyme than the APS substrate, especially 16d, which can be pointed as an APSrAB potential inhibitor.

Keywords:
sulfate-reducing microorganism; H₂S generating; dissimilatory adenosine 5’-phosphosulfate reductase inhibitors; molecular docking; molecular dynamics simulations

Introduction

In the sulfur biogeochemical cycle, the sulfate-reduction pathway can be assimilatory (ASR) and dissimilatory (DSR), constituting one of the oldest and most studied metabolic pathways on Earth.¹ The ASR pathway is mainly related to generating sulfur organic compounds (e.g., cysteine and methionine amino acids). In contrast, in the DSR pathway, sulfate is the final electron acceptor, generating H₂S as a final product.²

Sulfate-reducing microorganisms (SRM), like bacteria and archaea, are found in natural (e.g., groundwater) and industrial (e.g., exploration oil wells) environments and metabolize sulfate in the DSR pathway.³ In industry, H₂S generation can cause corrosion in oil extraction structures, reducing oil quality.⁴

In both reduction pathways, ASR and DSR, the catalysis begins with the reaction between SO₄²⁻ and adenosine triphosphate (ATP) catalyzed by ATP-sulfurylase (EC 2.7.7.4), generating adenosine 5’-phosphosulfate (APS).⁵ APS is the substrate of two enzymes, adenosine 5’-phosphosulfate reductase, assimilatory (APSr, EC: 1.8.4.10) and dissimilatory (APSrAB, EC 1.8.99.2), being converted in SO₃²⁻ and adenosine monophosphate (AMP) by both.⁶ In the final, SO₃²⁻ is reduced to S²⁻ by dissimilatory sulfite reductase enzyme (DSrAB, EC 1.8.99.3) or organic sulfur.⁷

Although APSr and APSrAB metabolize the same substrate (APS) and contain Fe-S clusters (i.e., [Fe₄S₄]²⁺) in their structures, they have different molecular architectures and cofactors (Figure 1).⁸ Pseudomonas aeruginosa APSr (thioredoxin dependent) (PDB ID: 2GOY), for example, is a homotetramer with one Fe-S cluster in the active site of each monomer, which acts as a cofactor.⁶ On the other hand, Archaeoglobus fulgidus APSrAB (PDB ID: 2FJA)⁹ is a heterodimer (AB). It contains two [Fe₄S₄]²⁺ groups in chain-B, far from the active site, but responsible for directing the electron flow to the FAD cofactor (flavinadenine dinucleotide) that reacts with the APS substrate.³

Figure 1
Ribbon diagram of two adenosine 5’-phosphosulfate reductases. (a) APSrAB, dissimilatory, from Archaeoglobus fulgidus (PDB ID: 2FJA); the active site contains substrate (APS) and cofactor (FAD), while the [Fe₄S₄]²⁺ cluster is far away from the site.⁹ (b) APSr, assimilatory, from Pseudomonas aeruginosa (PDB ID: 2GOY); the active site contains the same substrate (APS) and the cofactor is [Fe₄S₄]²⁺.⁶ The APS, FAD, and [Fe₄S₄]²⁺ molecules are in stick model and colored by elements.

The reduction of two electrons from APS to AMP is a fundamental step in both pathways for both eukaryotic and prokaryotic cells.¹ In the ASR pathway, thioredoxin (Trx) is the electron donor, while in the DSR pathway, this role is played by the QmoABC protein complex (quinone-interacting membrane-bound oxidoreductase).¹⁰,¹¹ Thus, both enzymes are pivotal objects of study.

APSr is an essential enzyme in the production of organic sulfur compounds, and the inhibition of the APSr of Mycobacterium tuberculosis in the sulfate assimilation pathway in mycobacteria is a research strategy involving the fight against tuberculosis.¹² The similar binding modes of APS and AMP to the enzyme encourage the search for new inhibitors based on the structures of these compounds.¹³,¹⁴

Paritala et al.¹¹ synthesized a series of 32 adenosine and AMP analogs (Figures 2 and 3), also assayed as potential inhibitors of the APSr enzyme from M. tuberculosis. The set of 32 compounds is composed by eight adenosine derivatives (Figure 2), with R² = H and R¹ = a (4a), b (4b), c (4c), d (4d), e (4e), f (4f), g (4g) and h (4h); eight AMP derivatives (Figure 2), with R² = PO₃²⁻ and substituted at the terminal amine with R¹ = a (8a), b (8b), c (8c), d (8d), e (8e), f (8f), g (8g) and h (8h); and 16 AMP derivatives (Figure 3), with R² = H and substituted at the PO₃(CH₂)_nCH₂R₃ group, with R³ = HS⁻ (compounds 13a, 13b, 13c and 13d), 1,2,3-triazol (compounds 14a, 14b, 14c and 14d), COOH (compounds 15a, 15b, 15c and 15d) and CONHOH (compounds 16a, 16b, 16c and 16d), with variation in the n value (1, 2, 3 and 4). The experimental enzyme inhibition data indicates that the substituted AMP derivatives at the terminal amine of the adenine group showed better equilibrium binding constant (K_d) values than the adenosine derivatives. In the best docking pose obtained using the AutoDock software, the most potent inhibitor (compound 8f, K_d = 4 μM; R = f, Figure 2) occupied the substrate-binding pocket, doing favorable interactions with key residues in the active site, in addition to a favorable interaction between S-atom of the ethyl-thiol of 8f with the S-atom of the Fe–S cluster of the APSr of Pseudomonas aeruginosa (PDB ID: 2GOY).¹¹

Figure 2
The general structure of the adenosine analogs (compounds 4a-4h, R¹ = a-h, R² = H) and substituted AMP analogs at the terminal amine (compounds 8a-8h, R¹ = a-h, R² = PO₃H₂).¹¹

Figure 3
The general structure of the substituted AMP analogs at the phosphate group.¹¹

Recently, our group¹⁵ developed parameters compatible with the CHARMM36 force field for Fe–S systems, which was used in a molecular dynamics (MD) simulation study of a series of chelating agents as potential APSrAB inhibitors using the GROMACS 2019 package, evaluating the interaction of these compounds with both, the APS substrate binding pocket (i.e., active site) and the binding pocket containing the [Fe₄S₄]²⁺ clusters.¹⁶ This forcefield was also used in other research groups,¹⁷ such as in the modeling of base excision repair (BER) enzymes.

Since the compounds synthesized by Paritala et al.¹¹ as APSr inhibitors are structurally similar to the APS substrate of the sulfate-reduction enzymes of both pathways (ASR and DSR), these compounds may interact with the catalytic site of the APSrAB, which can be investigated by theoretical methods, such as molecular docking and MD simulations,¹⁸,¹⁹ as done in our previous work involving APSrAB inhibitors.¹⁶ Then, in the present work, the compounds of Paritala et al.¹¹ will be evaluated as potential inhibitors of the APSrB enzyme using in silico docking and MD simulation techniques to obtain putative binding modes of these compounds with the APSrAB active site.

Methodology

Ligands and protein preparation

The structures of all ligands were previously submitted to energy minimization by the density functional theory (DFT) using the B3LYP exchange-correlation hybrid functional²⁰,²¹ and the 6-311G(d,p) basis set;²²,²³ the absence of imaginary frequencies indicated all structures corresponded to energy minima. The polarizable continuum model using the integral equation formalism variant (IEFPCM)²⁴ was included in the calculations to model the solvent effect (water). The ionization state of all ligands was adjusted to a pH = 7.4 (physiological pH) using Avogadro software.²⁵

The crystal structure of the Archaeoglobus fulgidus APSrAB protein constituted by two α:β heterodimers (chains A:B and C:D) and containing as ligands the FAD cofactor and APS substrate, resolved by X-ray diffraction, was obtained in the Protein Data Bank (PDB) server with the 2FJA code (resolution: 2.00 Å).⁹ The protein was prepared according to our previous work:¹⁶ in short, the structures of the C:D heterodimer (including Fe–S clusters on chain D) and FAD (on chain C) of the 2FJA.pdb file was maintained while the structures of the A:B heterodimer, APS (on chain C), and all water molecules were removed. The C:D heterodimer was chosen because it contained the APS substrate coordinates.

All the pdbqt files were generating using the AutoDockTools (ADT) graphical user interface.²⁶

Molecular docking protocol

The docking protocol of our previous work¹⁶ was used to evaluate the compounds of Paritala et al.¹¹ (Figures 2 and 3). The box with dimensions 40 × 40 × 40 ų was centered at the x, y, z coordinates of the APS substrate located in the active site on chain-C of the protein (x = 49.819, y = −12.188, z = 86.640), using a grid step size equal to 0.375 Å. The protein was maintained rigid, and the ligands and FAD were flexible. APS was removed during the experiments. The docking runs were executed with AutoDock Vina (gap energy = 4 kcal mol⁻¹; exhaustiveness = 8).²⁷

Our previous work¹⁶ selected the best APS pose based on its low energy (a maximum difference of 1.0 kcal mol⁻¹ from the lowest energy pose) and the lowest root mean square deviation (RMSD) to the crystal. Thus, for the other ligands, the pose with the shortest FAD-ligand distance was selected as the best pose, respecting the energy range of 1.0 kcal mol⁻¹. The FAD-ligand distance was considered as the distance between the N8 atom of the FAD (N atom that does a nucleophilic attack on the S-atom of the APS) and the X-atom close to N8, or N8(FAD)-X(ligand).

Molecular dynamics simulations

The MD simulations procedure followed the methodology adopted in our previous work:¹⁶ the best docking pose of each ligand in the protein structure (i.e., C:D heterodimer, two [Fe₄S₄]²⁺ clusters, and FAD cofactor) was considered as the initial protein-ligand complex; all energy minimization steps and MD simulations were executed with GROMACS 2019 package,²⁸ using the CHARMM36 force field²⁹ containing parameters for [Fe₄S₄(Cys)₄]²⁻.¹⁵ The topologies and parameters of the ligands were obtained in the CGenff server.³⁰

A triclinic simulation box with periodic boundary conditions (spc216.gro) with the TIP3P water model³¹ was used for the solvation of the protein-ligand complex. Before each simulation, each system was energy-minimized sequentially using the steepest-descent and conjugate gradient algorithms, and then equilibrated sequentially with the NVT and NPT ensembles (N = number of particles, V = volume, T = temperatures, P = pressure), considering a temperature of 300 K, pressure of 1 bar, and position restraint of the whole system, except for ions and water molecules. The total simulation time (production step) was 200 ns, with an integration time of 2 fs and a cutoff radius of 10 Å for long-range interactions.

All systems were evaluated based on RMSD (gmx rms module), FAD-ligand distance (gmx distance), and hydrogen bonds (gmx hbond, cut-off radius = 3.5 Å; cut-off angle = 30.00°). The frequency of hydrogen bond interactions was determined using HbMap2Grace³² from the 41 frames, with the most representative frame selected based on the clustering method.

Protein-ligand binding energy estimation

The trajectories of the systems that showed the best protein-ligand interaction (indicated by RMSD, FAD-ligand distance, and hydrogen bond analyses) during the MD simulations in GROMACS were submitted to the gmx_mmpbsa module³³ (igb = 5, saltcon = 0.150, idecomp = 2, dec_verbose = 3, print_res = “within 4”) to perform end-state free energy calculations: ΔG_binding = ΔG_complex − ΔG_receptor − ΔG_ligand. This module uses the Molecular Mechanics/Generalized Born Surface Area (MMGBSA) method, which is one of the most popular methods to estimate binding free energies because it balances accuracy and computational speed³⁴ and has been used in several recent MD studies, such as studies about anticancer agents,³⁵,³⁶ and neurodegenerative diseases,³⁷,³⁸,³⁹ for example. Additionally, the GBSA method was chosen because it is faster than the PBSA (Poisson-Boltzmann surface area) method, considering the same parameters.³³

ΔG_binding is the summation of the energies associated with the interactions in the vacuum (ΔG_gas) and solvent (ΔG_solv): ΔG_binding = E_bonded + E_el + E_vdW + ΔG_gb + ΔG_surf = ΔG_Gas + ΔG_Solv, where ΔG_gas is associated with the bonded (E_bonded), van der Waals (E_vdw) and electrostatic (E_el) interactions, while ΔG_solv represents the polar (ΔG_gb) and non-polar (ΔG_surf) interactions with the solvent.

The analysis was carried out using representative frames during the stability period of each complex (the central pose, obtained by the clustering method gmx_cluster, was the most representative).⁴⁰ The contribution of each residue to total energy was also evaluated.¹⁶

Results and Discussion

Molecular docking

The adjustment of the protonation state of the ligands presented the phosphate group completely ionized in all compounds. The carboxylic acid groups in compounds 4a-4b, 8a-8b, and 15a-15d were also ionized. The other groups remained not ionized.

The choice of the best pose in a redocking protocol should consider criteria such as low energy and RMSD values, indicating high affinity for the binding site and low deviation from the crystallographic structure, respectively.¹⁸,¹⁹ Thus, the docking protocol used in our previous work¹⁶ to reproduce the experimental binding mode of APS in the active site of the C-chain of the heterodimer (PDB ID: 2FJA) indicated pose No. 4 as having more favorable energy (binding affinity = -8.8 kcal mol⁻¹) and low RMSD about the crystallographic pose (RMSD = 1.8 Å).

Therefore, the poses of each compound with the shortest N8(FAD)-X(ligand) distance (Å) were selected (Table S1, Supplementary Information (SI) section). Overall, the affinity order with the active site was 8a-8h (average = −8.9 kcal mol⁻¹) > 13a-16d (mean = −8.6 kcal mol⁻¹) > 4a-4h (mean = −7.9 kcal mol⁻¹). The calculated N8(FAD)-S(APS) distance (3.6 Å) corresponds to the experimental value,⁹ while the N8(FAD)-X(ligand) distances were calculated in the range 2.1-4.7 Å (Table S1), indicating that all ligands presented a distance to the FAD close to the experimental N8(FAD)-S(APS) distance with a maximum deviation of around +2.1 Å.

The best poses of each compound obtained by molecular docking were used as starting structures for the MD simulations, a strategy employed in many recent works involving biological targets.⁴¹,⁴²,⁴³,⁴⁴,⁴⁵,⁴⁶,⁴⁷

Molecular dynamics simulations

The stability of the studied APSrAB system (PDB code 2FJA) was evaluated over 200 ns (Figure S1, SI section). The RMSD values for the protein, FAD cofactor, and [Fe₄S₄]²⁺ groups showed that the forcefield used adequately reproduced the evaluated system.

The RMSD values indicated the affinity of APS and the other ligands by docking site during the simulation (Table S2, SI section). APS showed affinity by the docking region (Figure S2, SI section), with an RMSD value of 2.93 ± 0.37 Å, in agreement with our previous works.¹⁶,⁴⁸ For the adenosine derivatives, 4b and 4c showed RMSD values compatible with their affinity in the docking region throughout the entire simulation (200 ns). In contrast, the other derivatives tended to move away from the docking site (Figure 4). For AMP analogs, only compounds 8a, 8c, and 8e showed affinity in the docking region. Although compounds 8f and 8h have presented stable RMSD values, these ligands showed affinities in areas far away from the active site (Figure 4).

Figure 4
RMSD values of the substituted adenosine and AMP analogs compared to APS (black) during 200 ns: compounds 4a-4d (A), 4e-4h (B), 8a-8d (C), and 8e-8h (D). Indexes: a = red, b = green, c = blue, d = yellow, e = orange, f = cyan, g = magenta, h = violet.

Among the phosphate analogs, only compounds 13a, 13d, 14a, 14c, 14d, 15b, 15c, 16b, 16c, and 16d tended to remain in the docking site as can be concluded by analysis of the RMSD values. Although compounds 14b and 15d presented stability during the simulation, their affinity was to a region far away from the enzyme’s active site (Figure 5, Table S2). These results pointed to compounds 4b, 4d, 8a, 8c, 8e, 13a, 13d, 14a, 14b, 14c, 14d, 15b, 15c, 16b, 16c, and 16d showed the highest affinity by the APSrAB active site.

Figure 5
RMSD values of the substituted phosphate analogs compared to APS (black) during 200 ns: compounds 13a-13d (A), 14a-14b (B), 15a-15d (C), and 16a-16b (D). Indexes: a = red, b = green, c = blue, d = yellow.

The N8(FAD)-X(ligand) distance was evaluated over time (Table S2). This distance was taken as the distance between N8(FAD) and atom X of the adenosine group (involving ribofuranose and the nitrogenous base). For the APS, the initial distance (t = 0) of 7.400 Å increased to an average distance of 10.500 ± 1.241 Å, which is a variation of about 3.00 Å. This value better correlates with the variation shown in Figure S3 (SI section).

It was observed that average N8(FAD)-X(ligand) distances (200 ns) in the following ranges 3.824-40.213 Å (4a-4h), 10.379-21.247 Å (8a-8h), and 5.076-38.563 Å (13a-16d) (Table S2). In general, the results agreed with the tendency observed in RMSD values: the ligands 4b, 4d, 8a, 8c, 8e, 13a, 13d, 14a, 14c, 14d, 15b, 15c, 16b, 16c, and 16d remained in the active site of the enzyme. An exception was observed for compound 14c, with variation exceeding 10.00 Å, although it can be related to the greater flexibility of the ligand. Anyway, 14c and the other ligands that tended to remain in the docking region showed great flexibility from the docking region, as observed in Figures 6 and S4 (SI section).

Figure 6
Frames plotting of t = 0 ns (carbons colored in gray) and structure generated by the gmx cluster module (carbons colored in magenta) for ligands.

The replicates of the systems APS (reference), compound 4a (tendency to leave the active site), and compound 16d (tendency to remain in the active site) were evaluated to determine if the observed tendencies were maintained. There was no significant difference between the original values (APS, 2.93 ± 0.37 Å; 4a, 11.39 ± 8.62 Å; and 16d, 3.74 ± 0.82 Å) and the respective average value of the triplicates (APS, 2.75 ± 0.57 Å; 4a, 7.78 ± 4.17 Å; and 16d, 4.66 ± 1.12 Å), reinforcing the observed tendencies.

The hydrogen bond interactions of the ligands 4b, 4d, 8a, 8c, 8e, 13a, 13d, 14a, 14c, 14d, 15b, 15c, 16b, 16c, and 16d were analyzed (considering only the interactions with lifetime ≥ 10.00%). The results showed that the residues in the 71-398 region, such as Asn74, Arg317, Asn318, His398, and, especially, Arg265, were mainly responsible for the docking of the ligands in the binding pocket over the simulation (Table S3, SI section, which contains the figures of the interactions). For APS, the hydrogen bonds involve mainly the residues Arg265, Val273 (92.68%), and Gly274 with the OPO₂OSO₃²⁻ terminal group, which participates directly in the catalysis.⁹

As for the inhibitor candidates, compound 4b presented interactions between the D-ribofuranose group of the nitrogenous base (Asn74) and, mainly, the terminal carboxylate in R¹ (Arg265 and Ser71, with a lifetime of about 100.00%); while compound 4d, showed interactions between Arg265 and the D-ribofuranose group, and between Asn74, Gl145, Trp144, and His446 with the adenosine (nitrogenous base) and R¹. The absence of the OPO₃³⁻ group in these adenosine derivatives favored the interactions involving the R¹ group, mainly the oxygen atoms.

For the AMP derivatives, compound 8a presented interactions Gln145---ribofuranose and Arg317---R¹; compound 8c, interactions between the PO₃²⁻ terminal and Arg265 (until 100.00%) and Ar317, beyond interactions between R¹ with Ala284 and Tyr292; for compound 8e, Arg265---PO₃²⁻ (about 95.00%). The presence of the OPO₃³⁻ group in these three compounds favored the interactions involving the oxygen atoms of this group, competing with the interactions involving R² group, mainly in 8c.

For the AMP derivatives substituted in the phosphate group, compound 13a had interactions of PO₃²⁻ with Arg265 and adenosine; for compound 13d, interactions of Arg265 with adenosine, ribofuranose and PO₃²⁻, between adenosine and Trp234 and His398, beyond PO₃²⁻---Val273 and PO₃²⁻---Gly274 interactions; for compound 14a, interactions of the adenosine group with Asn74, Trp234, Arg265, and His398 (100.00%), beyond the Arg265---PO₃²⁻, Glu141-sugar, and Asn74---ribofuranose interactions; compound 14c showed interactions of R³ with Phe264 and Arg265, beyond the Thr314---ribofuranose and PO₃²⁻---Arg317 interactions; for 14d, interactions of adenosine with FAD, Asn74, and Gln145, interactions between ribofuranose and Asn74 and Tyr95, beyond interactions of R³ with Arg265, Phe264 and Asn318. In general, the substituent R³ group in OPO₃³⁻ favored the interactions involving the oxygen atoms of this group and, mainly, the nitrogen atoms of the adenosine group. There was a significant contribution of the R³ group in compounds 14a-14c (nitrogen atoms of the triazole ring).

The ligand 15b showed interactions R³---Arg317 (100.00%), PO₃²⁻---Gln145, and interactions of adenosine with Asn74, Ser71, and FAD; 15c, interactions of R³ with Arg265 (92.68%) and Arg317 (100.00%), NH₂(terminal)---FAD, and interactions of ribofuranose with Gln145, and Arg265; 16b, interactions of adenosine with Ser71 and Gln145, PO₃²⁻---Arg265 and PO₃²⁻---Arg317 interactions; for 16c, interactions PO₃²⁻---Arg265 and between adenosine and Ser71 and Asn74; finally, for 16d, interactions of adenosine with Trp234 (100.00%) and Arg265, of ribofuranose with FAD, His398 (100.00%) and Asn74, beyond of R³ with Asn318, Tyr292, and Arg317. In general, for 15b-15c, the main interactions involved the R³ group (oxygen atoms of terminal carboxylate), while for 16b-16d, the oxygen and nitrogen atoms of phosphate and adenosine groups.

Taken together, analyses of hydrogen bonds interactions (Table S3) revealed that residues in the 71-398 region interact primarily with the adenosine group (nitrogenous base and R¹, present in three classes evaluated) and also with the PO₃²⁻ group (present in classes 8a-8h and 13a-16c), with Arg265 showing the most significant participation in the APSrAB-ligand interaction, with interactions of at least 80.00% lifetime for 4b, 8c, 8e, 13a, 13d, 14a, 14d, 15b, 15c, 16b, 16c, and 16d. The comparison of the results obtained for the ligands and the APS indicated that the studied ligands could make a more significant number of longer-lasting hydrogen bonds with the active site, competing with APS for binding in this region.

Ligands containing a terminal carboxylate or hydroxamate group, or a PO₃²⁻ group, showed a more effective interaction with the enzyme. All ligands evaluated at this step (4b, 4d, 8a, 8c, 8e, 13a, 13d, 14a, 14c, 14d, 15b, 15c, 16b, 16c, and 16d) showed a similar interaction with the enzyme (interaction involving the nitrogenous base and the PO₃²⁻ and the terminal carboxylate/hydroxamate groups if they were present in the molecule). Even those that showed few hydrogen bond interactions with a lifetime greater than 10.00% (8e, 13a, and 14c) cannot be discarded, as the great flexibility of the ligands may be responsible for the low number of interactions, as observed in the RMSD analysis. Thus, all these 14 ligands and APS were evaluated in terms of binding energy (ΔG_binding).

Protein-ligand binding energy estimation

The APS substrate presented ΔG_binding of −3.74 ± 4.85 kcal mol⁻¹ (Table 1), and this value was used as a reference to evaluate the affinity of other ligands with the binding site.

Table 1
ΔG_binding values obtained with the MMPBSA method

No.	ΔGbinding / (kcal mol−1)
4b	−24.27 ± 4.02
4d	−11.48 ± 3.12
8a	4.63 ± 7.15
8c	−26.44 ± 9.15
8e	−0.60 ± 5.33
13a	−15.55 ± 2.82
13d	−30.05 ± 4.00
14a	−37.76 ± 6.51
14c	−23.23 ± 4.33
14d	−32.41 ± 4.52
15b	−13.41 ± 7.21
15c	−16.64 ± 3.32
16b	−26.98 ± 4.56
16c	−30.56 ± 2.48
16d	−49.45 ± 3.05

• The individual contribution of each energy term for the ΔG_binding is in Table S4 (SI section). For the APS (substrate), a ΔG_binding = −3.74 ± 4.85 was calculated. The energy components for the ΔG_binding are listed in the Table S5 (SI section).

The ligands 4b, 4d, 8a, 8c, 8e, 13a, 13d, 14a, 14c, 14d, 15b, 15c, 16b, 16c, and 16d were selected for this step due to their tendency to remain in the binding site (Table 1). The ΔG_binding values revealed a binding affinity to protein in the following order: 16d > 14a > 14d > 16c > 13d > 16b > 8c > 4b > 15c > 13a > 15b > 4d > 8e > 8a. The enzyme had the most significant contribution to the stability of the APSrAB-ligand complex (receptor, Table S4, SI section). Except for the ligands 4b, 4d, 13a, 13d, 16b, 16c, and 16d, all the other ligands contributed to stabilizing the protein-ligand complex (Table S4). The AMP derivatives (8a and 8e) had the most favorable contribution (−226.94 ± 11.30 and −122.85 ± 9.95 kcal mol⁻¹, respectively).

According to the ΔG_binding values, all ligands, except 8a and 8e, bind more favorably to APSrAB than APS (from the thermodynamic point of view). The ligands 13d, 14a, 14d, 16c, and 16d showed ΔG_binding values more favorable than −30.00 kcal mol⁻¹, and 16d was the most favorable one (−49.45 ± 3.05 kcal mol⁻¹). In general, there was highest contribution for the ligand-APSrAB binding from the term related to solvation (ΔG_Solv) than the term of the gas phase (ΔG_Gas), mainly from the polar interactions term (ΔG_GB) (Table S5, SI section).

Taking into account the ΔG_binding and H-bond values, it was observed that the analogs substituted in the OPO₃³⁻ group (compounds 13-16) showed the highest affinity, followed by the adenosine derivatives (compounds 4), suggesting that the ligand-APSrAB interaction increases when the OPO₃³⁻ is absent or substituted. For these compounds, the FAD--- OPO₃³⁻ approximation, which is essential for catalysis involving APS, is compromised, especially in compounds 13-16, where the R³ groups cause an FAD---adenosine approximation, taking the OPO₃³⁻ far from the catalysis region. Although the absolute number of H-bonds did not bring a clear correlation with ΔG_binding values, the H-bonds involving the group OPO₃-R³ caused the FAD---adenosine interactions instead of the expected FAD---OPO₃³⁻ interactions. This reinforces the idea that the insertion of a R group of high size in OPO₃³⁻ (mainly a group containing nitrogen and oxygen atoms) increases the affinity of the ligand by the APSrAB.

Regarding the residues present in the active site, Arg317 plays a fundamental role in maintaining the APS conformation, while Asn74, Arg265, and His398 are important during catalysis.³ In our previous work,¹⁶ we observed the contribution of these residues to the APSrAB-APS binding, in addition to the contribution of the FAD cofactor, indicating that these species are essential for anchoring the substrate in the active site. The contribution of APS to binding destabilization is in agreement with the proposed mechanism, where the substrate is converted into the products AMP and SO₃²⁻, which are released from the active site.

The analysis of the individual energy contribution (kcal mol⁻¹) of the protein residues, FAD, and APS to the formation of the APSrAB-APS complex revealed that the most favorable contributions are from Arg265 (−12.318 ± 9.522), Pro272, Val273, Gly274, Phe277, Leu278, Arg317, and His398 (until about −5.00), while APS had an unfavorable contribution (26.278 ± 12.804). These results align with available experimental data, indicating the involvement of residues in the 71-398 region of the α-subunit, particularly Arg265, Arg317, and His398.³,⁹

Since 13d, 14a, 14d, 16c, and 16d showed more significant affinity to APSrAB, the individual contribution of protein residues, FAD, and these ligands to the energy of each APSrAB-ligand complex was evaluated (Table S6, SI section), revealing, as for APS, the most favorable contribution of residues in region 71-398 of the C-chain (α-subunit) in the binding of ligands, such as Asn74, Gln145, Trp234, Arg265, Phe277, Thr314, Arg317, Ans318, and His398, residues present in the active site of the enzyme.⁹ Residue Arg265 made the most significant contribution to binding, with values ranging from −4.817 and −12.897 kcal mol⁻¹.

Figure 7 showed that, for ligand 13d, it was the region containing residues 234-398, with the main contribution coming from Arg265 (−8.833 ± 8.608) followed by Pro272, Val273 and the 13d ligand itself (until about −5.00). For 14a, there was a majority contribution from Arg265 (−12.897 ± 10.809), with significant contributions from FAD and the ligand (until −6.00); for 14d, residues along region 74-318, mainly Arg317 (−7.675 ± 8.218), followed by Arg265 (−4,817 ± 7,719), Asn318, and Gln145; for 16c, the most significant contributions came from Arg265 (−10.154 ± 10.405), Arg317 (−2.639 ± 10.552) and the ligand itself (−6.031 ± 11.541). At the same time, for 16d, there were high contributions from Arg265 (−11.267 ± 10.555), FAD (−4.568 ± 20.286), and the ligand itself (−14.919 ± 11.332), in addition to residues Tyr292 and Asn318.

Figure 7
Individual energy contribution (kcal mol⁻¹) of each protein residue, FAD cofactor, and APS substrate for ΔG_binding of 13d, 14a, 14d, 16c, and 16d.

The 13d, 14a, 14d, 16c, and 16d ligands have, in common, the presence of a substituent of high size in OPO₃³⁻, which makes the FAD---PO₃³⁻ approximation difficult, suggesting that the steric obstruction of OPO₃³⁻ should play an essential role in their affinity to the enzyme. In general, it was observed that residues in regions 74-398 contributed in some way to the stabilization of the APSrAB-ligand complex, where the most significant contribution came from residue Arg265.

The results of the APSrAB-ligand binding energy for the ligands set 4b, 4d, 8c, 13a, 13d, 14a, 14c, 14d, 15b, 16b, 16c, and 16d are consistent with the observed results for the APS molecule in our previous work,¹⁶ indicating that all these ligands demonstrate affinity to the enzyme active site, especially 13d, 14a, 14d, 16c, and 16d due to the ΔG_binding values found to these six ligands.

The best result was observed for ligand 16d, which should have the most remarkable ability to compete with the APS substrate for the binding site on APSrAB.

Conclusions

The results obtained in this study identified the ligands 4b, 4d, 8c, 13a, 13d, 14a, 14c, 14d, 15b, 15c, 16b, 16c, and 16d as potential inhibitors of APSrAB. All presented values of variation in the FAD-ligand distance of up to 5.00 Å, indicating a tendency for these compounds to remain in the binding site predicted by molecular docking.

The analysis of hydrogen bond interactions showed that residues in the active site of the enzyme, all located in the 71-398 region, had lifetimes of more than 10.00%, especially Arg265, with values of up to 100.00% for some ligands. These residues in the 71-398 region made the main energetic contributions to stabilize the APSrAB-ligand complex, mainly Arg265, which contributed values for ΔG_binding more favorable than −10.00 kcal mol⁻¹ for 14a, 16c, and 16d. Ligands 13d, 14a, 14d, 16c, and 16d showed ΔG_binding values more favorable than −30.00 kcal mol⁻¹, indicating a high affinity for the enzyme, possibly, due to the steric obstruction caused by the R³ substituent in the POO₃³⁻ group. The best result was obtained for 16d, which presented the most favorable binding energy with the enzyme, indicating that this ligand is the most promising potential inhibitor of APSrAB in SRM.

Data Availability Statement

The authors confirm that all data generated and analyzed during this study are available in the article and the supplementary information.

Supplementary Information

Supplementary information (Figures S1-S4 and Tables S1-S6) is available free of charge at http://jbcs.sbq.org.br as PDFfile.

Acknowledgments

The authors thank the following Brazilian governmental agencies: CAPES (Coordenação de Aperfeiçoamento de Pessoal de Nível Superior), FAPERJ (Fundação de Amparo à Pesquisa do Estado do Rio de Janeiro), and CNPq (Conselho Nacional de Desenvolvimento Científico e Tecnológico). This research was developed with the help of CENAPAD-SP (Centro Nacional de Processamento de Alto Desempenho - São Paulo), UNICAMP/FINEP-MCT (Universidade Estadual de Campinas / Financiadora de Estudos e Projetos - Ministério da Ciência, Tecnologia e Inovação) Project (proj 835). This study was financed by CAPES, FAPERJ (scholarship E-26/010.210.513/2019 for Sergio P. Machado; SEI-260003/007043/2022 and SEI-260003/003788/2022 for Camilo H. S. Lima), and CNPq (grant 304402/2017 and grant 304105/2021-0 for Sergio P. Machado; SI scholarship 147933/2024-2 for Henrique N. Lança).

References

• ¹ Kushkevych, I.; Abdulina, D.; Kováč, J.; Dordević, D.; Vítězová, M.; Iutynska, G.; Simon, K.-M. R.; Rittmann, S. K. M. R.; Biomolecules 2020, 10, 921. [Crossref]
  » Crossref
• ² Wu, B.; Liu, F.; Fang, W.; Yang, T.; Chen, G. H.; He, Z.; Wang, S.; Sci. Total Environ. 2021, 778, 146085. [Crossref]
  » Crossref
• ³ Wójcik-Augustyn, A.; Johansson, A. J.; Borowski, T.; Biochim. Biophys. Acta, Bioenerg. 2021, 1862, 148333. [Crossref]
  » Crossref
• ⁴ de Jesus, E. B.; de Andrade Lima, L. R. P.; Bernardez, L. A.; Almeida, P. F.; Braz. J. Pet. Gas 2015, 9, 95. [Crossref]
  » Crossref
• ⁵ Grein, F.; Ramos, A. R.; Venceslau, S. S.; Pereira, I. A. C.; Biochim. Biophys. Acta, Bioenerg. 2013, 1827, 145. [Crossref]
  » Crossref
• ⁶ Chartron, J.; Carroll, K. S.; Shiau, C.; Gao, H.; Leary, J. A.; Bertozzi, C. R.; Stout, C. D.; J. Mol. Biol. 2006, 364, 152. [Crossref]
  » Crossref
• ⁷ Oliveira, T. F.; Vonrhein, C.; Matias, P. M.; Venceslau, S. S.; Pereira, I. A. C.; Archer, M.; J. Biol. Chem. 2008, 283, 34141. [Crossref]
  » Crossref
• ⁸ Palde, P. B.; Bhaskar, A.; Pedró Rosa, L. E.; Madoux, F.; Chase, P.; Gupta, V.; Spicer, T.; Scampavia, L.; Singh, A.; Carroll, K. S.; ACS Chem. Biol. 2016, 11, 172. [Crossref]
  » Crossref
• ⁹ Schiffer, A.; Fritz, G.; Kroneck, P. M. H.; Ermler, U.; Biochemistry 2006, 45, 2960. [Crossref]
  » Crossref
• ¹⁰ Duarte, A. G.; Santos, A. A.; Pereira, I. A. C.; Biochim. Biophys. Acta, Bioenerg. 2016, 1857, 380. [Crossref]
  » Crossref
• ¹¹ Paritala, H.; Suzuki, Y.; Carroll, K. S.; Nucleosides, Nucleotides Nucleic Acids 2015, 34, 199. [Crossref]
  » Crossref
• ¹² Bhave, D. P.; Muse III, W. B.; Carroll, K. S.; Infect. Disord.: Drug Targets 2007, 7, 140. [Crossref]
  » Crossref
• ¹³ Feliciano, P. R.; Carroll, K. S.; Drennan, C. L.; ACS Omega 2021, 6, 13756. [Crossref]
  » Crossref
• ¹⁴ Eric, S.; Cvijetic, I.; Zloh, M.; J. Serbian Chem. Soc. 2021, 86, 561. [Crossref]
  » Crossref
• ¹⁵ da Silva, T. U.; Pougy, K. C.; Albuquerque, M. G.; Lima, C. H. S.; Machado, S. P.; J. Biomol. Struct. Dyn. 2022, 40, 3481. [Crossref]
  » Crossref
• ¹⁶ da Silva, T. U.; Pougy, K. C.; Albuquerque, M. G.; Lima, C. H. S.; Machado, S. P.; J. Biomol. Struct. Dyn. 2023, 41, 2466. [Crossref]
  » Crossref
• ¹⁷ Olivieri, A.; Nardi, A. N.; D’Abramo, M.; Biochemistry 2024, 63, 1991. [Crossref]
  » Crossref
• ¹⁸ dos Santos, E. S.; de Souza, L. C. V.; de Assis, P. N.; de Almeida, P. F.; Ramos-De-Souza, E.; J. Biomol. Struct. Dyn. 2015, 33, 1176. [Crossref]
  » Crossref
• ¹⁹ dos Santos, E. S.; de Souza, L. C. V.; de Assis, P. N.; Almeida, P. F.; Ramos-De-Souza, E.; J. Biomol. Struct. Dyn. 2014, 32, 1780. [Crossref]
  » Crossref
• ²⁰ Becke, A. D.; J. Chem. Phys. 1993, 98, 5648. [Crossref]
  » Crossref
• ²¹ Lee, C.; Yang, W.; Parr, R. G.; Phys. Rev. B 1988, 37, 785. [Crossref]
  » Crossref
• ²² Krishnan, R.; Binkley, J. S.; Seeger, R.; Pople, J. A.; J. Chem. Phys. 1980, 72, 650. [Crossref]
  » Crossref
• ²³ McLean, A. D.; Chandler, G. S.; J. Chem. Phys. 1980, 72, 5639. [Crossref]
  » Crossref
• ²⁴ Cancès, E.; Mennucci, B.; Tomasi, J.; J. Chem. Phys. 1997, 107, 3032. [Crossref]
  » Crossref
• ²⁵ Hanwell, M. D.; Curtis, D. E.; Lonie, D. C.; Vandermeersch, T.; Zurek, E.; Hutchison, G. R.; J. Cheminform. 2012, 4, 17. [Crossref]
  » Crossref
• ²⁶ Morris, G. M.; Huey, R.; Lindstrom, W.; Sanner, M. F.; Belew, R. K.; Goodsell, D. S.; Olson, A. J.; J. Comput. Chem. 2009, 30, 2785. [Crossref]
  » Crossref
• ²⁷ Trott, O.; Olson, A. J.; J. Comput. Chem. 2009, 31, 455. [Crossref]
  » Crossref
• ²⁸ Abraham, M. J.; Murtola, T.; Schulz, R.; Páll, S.; Smith, J. C.; Hess, B.; Lindahl, E.; SoftwareX 2015, 1-2, 19. [Crossref]
  » Crossref
• ²⁹ Huang, J.; MacKerell, A. D.; J. Comput. Chem. 2013, 34, 2135. [Crossref]
  » Crossref
• ³⁰ Vanommeslaeghe, K.; Hatcher, E.; Acharya, C.; Kundu, S.; Zhong, S.; Shim, J.; Darian, E.; Guvench, O.; Lopes, P.; Vorobyov, I.; Mackerell, Jr. A. D.; J. Comput. Chem. 2010, 31, 671. [Crossref]
  » Crossref
• ³¹ Mark, P.; Nilsson, L.; J. Phys. Chem. A 2001, 105, 9954. [Crossref]
  » Crossref
• ³² Gomes, D. E. B.; Silva, A. W.; Lins, R. D.; Pascutti, P. G.; Soares, T. A.; HbMap2Grace, v7.f90, Universidade Federal do Rio de Janeiro, RJ, Brazil, 2009. [Crossref]
  » Crossref
• ³³ Valdés-Tresanco, M. S.; Valdés-Tresanco, M. E.; Valiente, P. A.; Moreno, E.; J. Chem. Theory Comput. 2021, 17, 6281. [Crossref]
  » Crossref
• ³⁴ Demir, S.; Alparslan, G. T.; J. Mol. Graph. Model. 2025, 137, 108994. [Crossref]
  » Crossref
• ³⁵ Albassam, H.; Almutairi, O.; Alnasser, M.; Altowairqi, F.; Almutairi, F.; Alobid, S.; J. Enzyme Inhib. Med. Chem. 2025, 40, 2468852. [Crossref]
  » Crossref
• ³⁶ Azam, F.; Khan, A.; Ahsan, M. J.; J. Mol. Struct. 2025, 1331, 141574. [Crossref]
  » Crossref
• ³⁷ Bhagat, K.; Yadav, A. J.; Padhi, A. K.; J. Phys. Chem. B 2025, 129, 3366. [Crossref]
  » Crossref
• ³⁸ Liu, B.; Zhao, L.; Tan, Y.; Yao, X.; Liu, H.; Zhang, Q.; ACS Chem. Neurosci. 2025, 16, 1617. [Crossref]
  » Crossref
• ³⁹ An, H.; Shao, C.; He, Y.; Zhou, H.; Wang, T.; Xu, G.; Yang, J.; Wan, H.; ACS Chem. Neurosci. 2025, 16, 1550. [Crossref]
  » Crossref
• ⁴⁰ Camargo, P. G.; Suzukawa, H. T.; Pereira, P. M. L.; Silva, M. L.; Macedo, Jr. F.; Albuquerque, M. G.; Rodrigues, C. R.; Yamada-Ogatta, S. F.; Lima, C. H. S.; Bispo, M. L. F.; Sci. Rep. 2025, 15, 465. [Crossref]
  » Crossref
• ⁴¹ Camargo, P. G.; dos Santos, C. R.; Albuquerque, M. G.; Rodrigues, C. R.; Lima, C. H. S.; Sci. Rep. 2024, 14, 11575. [Crossref]
  » Crossref
• ⁴² Chen, R.; Liu, H.; Meng, W.; Sun, J.; Sci. Rep. 2024, 14, 21043. [Crossref]
  » Crossref
• ⁴³ Shamsi, A.; Khan, M. S.; Yadav, D. K.; Shahwan, M.; Furkan, M.; Khan, R. H.; Sci. Rep. 2024, 14, 19439. [Crossref]
  » Crossref
• ⁴⁴ Shamsi, A.; Khan, M. S.; Yadav, D. K.; Shahwan, M.; Furkan, M.; Khan, R. H.; Sci. Rep. 2024, 14, 21778. [Crossref]
  » Crossref
• ⁴⁵ Das, R.; Bhattarai, A.; Karn, R.; Tamang, B.; Sci. Rep. 2024, 14, 19585. [Crossref]
  » Crossref
• ⁴⁶ Moussaoui, M.; Baammi, S.; Soufi, H.; Baassi, M.; Allali, A. E.; Belghiti, M. E.; Daoud, R.; Belaaouad, S.; Sci. Rep. 2024, 14, 16418. [Crossref]
  » Crossref
• ⁴⁷ Li, C.; Zhuo, C.; Ma, X.; Li, R.; Chen, X.; Li, Y.; Zhang, Q.; Yang, L.; Wang, L.; Schizophrenia 2024, 10, 80. [Crossref]
  » Crossref
• ⁴⁸ da Silva, T. U.; da Silva, E. T.; Pougy, K. C.; Lima, C. H. S.; Machado, S. P.; Inorg. Chem. Commun. 2022, 135, 109120. [Crossref]
  » Crossref

Edited by

• Editor handled this article:
  Giovanni Wilson Amarante (Associate)

Publication Dates

• Publication in this collection
  15 Sept 2025
• Date of issue
  2025

History

• Received
  19 June 2025
• Accepted
  11 Aug 2025