Evaluation of propolis activity as sucrose-dependent and sucrose-independent Streptococcus mutans inhibitors to treat dental caries using an in silico approach

Propolis screening

A propolis sample was provided by PT. Nano Herbaltama Internasional, Indonesia. Methanol 9 ml and propolis 1 ml were added to a 10 ml volumetric. The Minisart SRP15 Syringe Filter PTFE 0.2 μm (Sartorius, Germany) was used to filter the solution. LC-MS/MS analyses were done by adding a 2 ml vial of the filtered solution in the LC-MS/MS Xevo G2-XS Quadrupole Time-of-Flight Mass Spectrometer (Waters, New Zealand). The separation process was figured out using the ACQUITY UPLC HSS-T3 column (100 × 2.1 mm, 1.7 μm). The temperature was set at 40 and 25?C to maintain the column and sample, respectively. Formic acid of 0.1% in water (A) and formic acid of 0.1% in acetonitrile (B) were used as the mobile phase at a flow rate of 0.6 ml/minute. The preparation of the elution gradient was done using the following protocol: starts at 1%, at 0.5 minutes 1% A, at 16 minutes 35% A, at 18 minutes 100%, and 20 minutes 1% A. Flavonoids screening was done by injecting 2 μl of the test solution, and the chromatograph was recorded for 20 minutes. Electrospray Ionization positive ion mode was set with the range scan of mass at 50–1,200 Da. The source and desolvation temperatures were 120^?C and 550^?C, respectively. The low collision energy was 6 eV, and the high collision energy was 15–40 eV (ramp). The mass correction was done using a 1 ng/ml solution of leucine enkephalin as a lock spray during acquisition. Finally, a flow rate of 10 μl/minute was used to inject the lock spray, and the scan time was 0.1 seconds with an interval of 30 seconds.

Receptor and ligand preparation

A three-dimensional crystal protein structure of glucosyltransferase [Protein data bank (PDB) ID: 3AIC] and Ag I/II (PDB ID: 3IPK) from S. mutans were retrieved from Protein Data Bank PDB (www.rcsb.org). These proteins were analyzed for their similarity to the human protein, also Lactobacillus reuteri and Bifidobacterium lactis as oral probiotics using the Protein Basic Local Alignment Search Tool (BLAST-P) on the NCBI website (https://www.ncbi.nlm.nih.gov/). After conducting BLAST-P, all proteins were prepared by the AutoDockTools version 4.2 program (Scripps Research Institute, La Jolla, CA) to remove the cocrystallized ligand and water while missing polar hydrogen was added and merged with nonpolar hydrogen.

The structures of up to a total of 93 ligands from propolis screening compounds were retrieved from PubChem (https://pubchem.ncbi.nlm.nih.gov/) in the 3D structure data format form. Some ligands that did not have a 3D structure from PubChem were generated using Discovery Studio Visualizer 2021. All ligands were converted into the PDB, partial charge (Q), and atom type (T) form using PyRx 0.8 (Dallakyan and Olson, 2015). Thus, these were prepared to add the Gasteiger charges and merged with the nonpolar hydrogen molecules using AutoDockTools. Rotatable torsions of the ligands were activated. The grid box coordinate was determined using the center grid box option of the receptor-ligand complex.

Molecular docking validation and visualization

AutoDockTools was used to validate the molecular docking method. Grid coordination of the native ligand and protein was saved as a grid parameter file to create a grid log file (GLG) file using AutoGrid. The GLG file was used to build docking parameter file and docking log file (DLG) files while running AutoDock. The maximum number of evals while running was at a medium level using a Lamarckian genetic algorithm as the output function, and the rest of the parameters were set to default. The root mean square deviation (RMSD) value was analyzed from the DLG file and visualized using Discovery Studio Visualizer 2021.

Molecular docking and visualization

Molecular docking was carried out using the PyRx 0.9 version, which used AutoDock 4.2 as the docking software (Dallakyan and Olson, 2015) underperformed on Windows Pro 10, Intel Core^® 7, and RAM 16 GB. All 93 ligands were loaded into PyRx and the grid box coordination that was obtained from AutoDockTools during preparation was set (193.74; 49.075; 195.987 and 1743.498; 44.166; 198.644 of x, y, z coordination 3AIC and 3IPK, respectively, with spacing in angstrom of 0.375). The docking process was done in PyRx using the Lamarckian genetic algorithm. All the other parameters were kept as default. Discovery Studio Visualizer 2021 was used to visualize the docking result.

ADMET analysis

Absorption, Distribution, Metabolism, and Excretion (ADME) analyses were done using the SwissADME web server (https://www.swissadme.ch/) (Daina et al., 2017). The top 10 ligands after conducting molecular docking were moved forward for ADMET analysis. The ligands of this study were conducted for dental oral local treatment; therefore, some parameters such as the maximum two violation of Lipinski’s rule, water-soluble drug, and both negatives of the blood-brain barrier and P-gp substrate were considered (Al-Khafaji et al., 2021; Sharif et al., 2013). Moreover, for toxicity for lethal dose 50 (LD₅₀ mg/kg) prediction, the ProTox-II program was used (https://tox-new.charite.de) (Banerjee et al., 2018). A ligand that has been classified as having above class 3 toxicity was selected.

Molecular protein interaction

The best two ligands after molecular docking and ADMET screening were selected for molecular protein interaction analysis. Discovery Studio Visualizer 2021 was used to analyze and visualize some parameters such as pose conformation, amino acid residue, and chemical interaction between ligand and protein.

Conserved amino acid analysis

The ConSurf web server was used (https://consurf.tau.ac.il/) to analyze the conserved amino acid sequence. This server used the Bayesian method to identify the residues in the conserved region (Ashkenazy et al., 2016). PDB IDs of 3AIC and 3IPK were uploaded to the web server, and the identification region was limited to chain A only. Multiple sequence alignment (MSA) was automatically made from ConSurf, and the rest of the parameters were kept at default. Scaling conserved amino acid result was displayed in the gradation color starting from 1 as the variable until 9 as the conserved region.

Molecular dynamics

The molecular dynamics simulation studies were executed using GROMACS 2021.5 (Berendsen et al., 1995). A unit cell with a rhombic dodecahedral shape was utilized for the protein-ligand complex and positioned at its center at least 10 Å from the edge of the box. The system was neutralized by replacing some water molecules with Na+ and Cl−. The energy minimization was performed with the steepest descent algorithm until the maximum force reached less than 1,000 kJ mol⁻¹ nm⁻¹. Equilibration was conducted for 100 ps at 300 K and separated into two stages, the isothermal-isochoric ensemble and the isothermal-isobaric ensemble. The particle-mesh Ewald was used for electrostatic interaction calculation, and the short-range cut-off of Van der Walls was 12 Å. The 50 ns molecular dynamics production was carried out with a modified Berendsen thermostat as the temperature coupling and a Parrinello-Rahman barostat as the pressure coupling. RMSD and root mean square fluctuation (RMSF) were calculated from the molecular dynamics trajectories.

Receptor and ligand preparation

To address the potential comprehensive inhibition of S. mutans as an important cariogenic bacterium, we used glucosyltransferase protein (Gtf) (PDB ID: 3AIC) as a sucrose-dependent agent and Ag I/II as a sucrose-independent agent (PDB ID: 3IPK). From the BLAST-P (NCBI) result, the 3IPK and 3AIC proteins were less identical to L. reuteri protein sequences. Moreover, 3AIC was also less identical to human protein sequences (Supplementary Figs. S1–S3).

As many as 93 compounds were identified in propolis screening by the LC/MS-MS method. Most of these compounds were determined as being in the flavonoids group. Flavonoids are built by a 15-carbon skeleton consisting of two rings of benzene that are linked by a heterocyclic pyrene ring (Khameneh et al., 2019; Wagh, 2013). Among 93 flavonoid ligands, the most prominent flavonoid subgroups were flavones, isoflavones, and flavonols. Most flavone groups have a hydroxyl group in the primary ring of position 5 while flavonols have the hydroxyl group in the 3 positions of ring C. Furthermore, the isoflavones subgroup has a substituent benzene ring at 3 positions of ring C (Mutha et al., 2021; Panche et al., 2016). Flavonoids are usually linked with health-promoting effects, and therefore they are common as components in nutraceutical, pharmaceutical, medicinal, and cosmetic products (Wagh, 2013). Their broad activities are structure dependent. For instance, the hydroxyl functional group in flavonoids induces an antioxidant effect by scavenging free radicals or chelating metal ions (Ko?evar et al., 2007). In addition, the antimicrobial effect of flavonoids was also recorded, which is dependent on the hydrophobic substituents’ structure such as the prenyl groups (Xie et al., 2014). There are several mechanisms in the flavonoid group to combat bacteria, such as inhibition of nucleic acid synthesis, inhibition of cytoplasmic membrane function by disrupting biofilm formation, and some crucial enzymes (Barbieri et al., 2017; Górniak et al., 2019; Khameneh et al., 2019). Furthermore, Yuan et al. (2021) revealed that flavonoids mainly disrupt the cell membrane of Gram-positive bacteria such as S. mutans as the major target. To examine the flavonoid content as a potential anticaries agent, an in silico method was conducted. All ligands’ details are shown in Supplementary Table S3.

Before conducting the molecular docking process between the ligand and protein, docking validation was carried out. Acarbose as a native ligand from the 3IAC protein and phenylmethanesulfonic acid (PMS) as a native ligand from the 3IPK protein were used for this validation. From validation of the molecular docking method, the lowest binding affinity of acarbose was −7.16 kcal/mol and −4.75 kcal/mol of PMS with the RMSD values of 1 and 1.24, respectively (Supplementary Tables S1 and S2). Supplementary Figure S4 showed the result from the docking validation and the binding site of the 3AIC and 3IPK proteins.

Molecular docking and ADMET analysis

The top 10 docking results of the 3AIC and 3IPK proteins were summarized in Tables 1 and 2, with the range of binding affinity from all ligands of 3AIC being −7.53 to −8.44 kcal/mol and the range of 3IPK from −6.5 to −7.32 kcal/mol. On the other side, acarbose and PMS as native ligands were −5.78 kcal/mol and −4.64 kcal/mol, respectively.

During the in silico process, molecular docking is important to conduct drug ability screening. The molecular docking process shows the list of interactions between two molecules to build a stable complex (Trezza et al., 2020). The lower binding affinity indicates the more stable and closer native state of protein-ligand efficiency. It exhibits a good potential drug (Fu et al., 2018). According to the docking results (Tables 1 and 2), it seemed that lower binding affinity results of some ligands after running molecular docking from the 3AIC and 3IPK proteins were affected by the more complex structure. This can be seen in the best ligand of both proteins, especially diosmetin, cosmosiin, genistin, and 3′-methoxy-5′-hydroxyisoflavone-7-O-β-D glucoside. The interactions among these proteins as the best ligand of the 3AIC and 3IPK proteins were visualized (Supplementary Fig. S5). Hydrogen, hydrophobic, Van der Waals, and electrostatic bonds are the main interactions between ligand-protein complexes, and these are very important to creating stable bonding with the receptor (Freitas and Schapira, 2017; Shawon et al., 2021). In contrast, acarbose and PMS as native ligands generally only had hydrogen and Van der Waals bonds without other interactions, such as a hydrophobic interaction, except for an electrostatic bond in the PMS-3IPK complex (Supplementary Fig. S10). It may result in the lower binding affinity of these native ligands than diosmetin, cosmosiin, genistin, and 3′-methoxy-5′-hydroxyisoflavone-7-O-β-D glucoside.

Table 1. Top 10 binding affinity of 3AIC protein. [Click here to view]

Table 2. Top 10 binding affinity of 3IPK protein. [Click here to view]

The top 10 best ligands from molecular docking were moved forward for ADMET analysis. Diosmetin and cosmosiin were the best ligand interactions with the 3AIC receptor that met all drug-like requirements besides their binding affinity score (Table 3). Furthermore, in the 3IPK receptor, genistin and 3′-methoxy-5′-hydroxyisoflavone-7-O-β-D glucoside were the best ADME (Table 4) and docking results. To assess toxicology, the ProTox-II web server showed the classification of toxic doses. Diosmetin and cosmosiin were classified as having class 5 toxicity, which is known as a safe compound (Supplementary Fig. S5). Moreover, genistin and 3′-methoxy-5′-hydroxyisoflavone-7-O-β-D glucoside were classified as having classes 5 and 6 toxicity (Supplementary Fig. S6), respectively.

Molecular protein-ligand interaction

According to the binding affinity and ADMET result, the two best ligands of each protein receptor were chosen to be visualized (Fig. 1). Diosmetin and cosmosiin were selected from the 3AIC protein (Fig. 1A–D) while genistin and 3′-methoxy-5′-hydroxyisoflavone-7-O-β-D glucoside were selected from the 3IPK protein (Fig. 1E–H). Diosmetin formed three classical hydrogen bonds with ASP 477, TRP 517, and GLU 515. Besides classical hydrogen, it also formed a weak hydrogen bond as a carbon-hydrogen bond with ASP 909 and ASN 481. There were four amino acids from diosmetin that built a hydrophobic interaction (TYR 916 as pi-pi interaction, HIS 587, ALA 478, and LEU 433 as pi-alkyl interaction). Another noncovalent bonding of a small ligand was the Van der Waals interaction. It was built by diosmetin to interact with the 3AIC protein with the residues ASN 914, ARG 475, ASN 862, and ASP 480. Furthermore, an additional interaction involved pi-anion bonding from ASP 588 also existing between diosmetin and the 3AIC protein. On the other side, cosmosiin did not build any electrostatic interactions like diosmetin. Cosmosiin built classical hydrogen bonds with ASP 477, ASP 480, ASP 588, GLU 515, and GLN 960. Weak hydrogen bonds only interacted with ASN 481 and ASP 909. Furthermore, hydrophobic interactions existed among TRP 517 as pi-pi interaction, LEU 382, LEU 433 as pi-sigma interaction, and ALA 478 as pi-alkyl interaction. The last interaction in cosmosiin with the 3AIC protein was the Van der Waals interaction that involved TYR 430, LEU 434, ARG 475, HIS 587, GLN 592, PHE 907, and TYR 916.

In the 3IPK protein, genistin built a hydrogen bonding interaction with ASP 512, ARG 824, ASP 760, TRP 816, and THR 586 as conventional hydrogen and SER 697 as a weak hydrogen bond. A hydrophobic interaction also existed in genistin interaction LYS 822 as pi-alkyl and TRP 816 as pi-pi stacked. The Van der Waals interaction involves LEU 653, ASN 820, SER 762, and VAL 587. Furthermore, 3′-methoxy-5′-hydroxyisoflavone-7-O-β-D glucoside-built hydrogen bonds with THR 652, ASN 820, ARG 824, TRP 816, and THR 586 as conventional hydrogen bonds and ASP 909 and ASN 481 as weak hydrogen bonds. The hydrophobic interaction involved TRP 816 as pi-pi stacked, LEU 653 as alkyl, and LYS 822 as pi-alkyl. The Van der Waals interaction involved ASP 512, SER 762, and VAL 587. Detailed information on molecular interaction was displayed in Tables S4 and S5.

Table 3. Ligand ADME analyses result in the 10 best potential ligands for 3AIC protein. [Click here to view]

Table 4. Ligand ADME analyses result in the 10 best potential ligands for 3IPK protein. [Click here to view]

Freitas and Schapira (2017) showed that the hydrophobic interaction was the most prominent among the high-efficiency ligands. It is different from protein-fragment complexes that are more stable using hydrogen bond interaction. In the hydrophobic interaction, the aliphatic carbon of the receptor and aromatic carbon of the ligand are most prevalent. Regarding the aromatic carbon, a small ligand as an inhibitor generally contains aromatic carbon such as a benzene ring to build the interaction with its receptor, and it is owned by the flavonoids group, as seen in Supplementary Table S3. Moreover, leucine and valine are frequently involved in the hydrophobic interaction (Freitas and Schapira, 2017), and they also exist in the ligand-protein complex of both the 3AIC and 3IPK proteins, which can be seen in Supplementary Tables S4 and S5.

Analyzing the data from acarbose, diosmetin, and cosmosiin, amino acids Glu 515, Leu 433, ASN 481, and ASP 477 were the residues that were always involved in the interaction between the ligand and 3AIC receptor shown by their existence in the interaction (Supplementary Table S4). Moreover, ARG 824 and TRP 816 were amino acid residues that existed in PMS as the native ligand of 3IPK also in both genistin and 3′-methoxy-5′-hydroxyisoflavone-7-O-β-D glucoside as the best ligand compound (Supplementary Table S5). It can be known that these amino acids acted as key residues in the 3AIC and 3IPK interactions, respectively. The key residue that acts as a functional amino acid is only distributed in the active binding site of the protein (Chen et al., 2006). Amitai et al. (2004) and Ma et al. (2003) stated that the key residue is generally well conserved in the protein structure, which has an important role in protein interaction and drug development. It can be considered because of its highest binding energy related to the others. As displayed in Supplementary Figures S7 and S8, the key residues in the 3AIC and 3IPK proteins were well conserved compared to the other organisms that have high similarity to these proteins. Furthermore, a mutation in the key residue region will most affect the protein stability and substrate specificity (Yang et al., 2021).

Molecular dynamics

After conducting the molecular dynamics (MD) process within 50 ns, the RMSD and RMSF values from each protein-ligand complex were obtained (Fig. 2). From Figure 2B, it can be observed that the 3AIC protein with diosmetin reached stability after 10 ns with the structure overall remaining at 3.5 Å while cosmosiin kept stable at 4 Å (Fig. 2C) after 10 ns. Furthermore, the RMSF value from both ligands almost had a similar pattern (Fig. 2D).

Figure 1. Protein-ligand interaction results from potential ligands of the 3AIC (A–D) and 3IPK (E–H) proteins. Molecular interaction (A) and hydrophobicity area (B) of diosmetin. Molecular interaction (C) and hydrophobicity area (D) of cosmosiin. Molecular interaction(E) and hydrophobicity area (F) of genistin. Molecular interaction (G) and hydrophobicity area (H) of 3′-methoxy-5′-hydroxyisoflavone-7-O-β-D glucoside in the 3IPK protein. [Click here to view]

Figure 2. Molecular dynamics results from potential ligands of the 3AIC (A–F) and 3IPK (G-L) proteins. RMSD (A–C) and RMSF (D–F) of diosmetin and cosmosiin, RMSD (G–I) and RMSF (J–L) of genistin and 3′-methoxy-5′-hydroxyisoflavone-7-O-β-D glucoside. The residue number displayed within the RMSF result showed the binding site residue from both ligands (7E, 7F, 7K, and 7L). [Click here to view]

Similar to the ligand of the 3AIC protein, genistin and 3′-methoxy-5′-hydroxyisoflavone-7-O-β-D glucoside as the best candidate ligands of the 3IPK protein also were examined using MD. The RMSD value showed that genistin reached stabilization at around 3 Å (Fig. 2H) while 3′-methoxy-5′-hydroxyisoflavone-7-O-β-D glucoside was stable at around 4 Å (Fig. 2I) with more fluctuations than the 3AIC complexes. On the other side, the RMSF value showed that both ligands had a similar pattern (Fig. 2J).

In contrast to molecular docking, molecular dynamics can show the flexibility of the ligand and receptor within a certain solvent, therefore allowing the receptor to change its conformation (Prasasty et al., 2020). The two most common methods to calculate the conformation change are RMSD and RMSF. RMSD shows the average atomic displacement from the entire protein against its reference structure in a time-dependent manner (Martínez, 2015; Ramírez and Caballero, 2018). Overall, both diosmetin and cosmosiin almost had a similar pattern of RMSD with a value maximum at around 5 Å (Fig. 2A). It showed that both ligands might have the same influence on the 3AIC protein. Conversely, the RMSD values from genistin and 3′-methoxy-5′-hydroxyisoflavone-7-O-β-D glucoside were slightly different from each other (Fig. 2G). Ligand 3′-methoxy-5′-hydroxyisoflavone-7-O-β-D glucoside showed more fluctuation than genistin where the maximum fluctuation of both ligands was at around 12 Å. Some fluctuations in both complexes indicated that the protein searches for other conformations to stabilize (Prasasty et al., 2020).

Alongside the RMSD value, we also examined the RMSF value of the protein-ligand complex. RMSF is used to examine the fluctuation of each residue in the entire protein. Overall, RMSF from the 3AIC and 3IPK complexes almost showed a similar pattern (Fig. 2D and J). Yang et al. (2017) and Baweja et al. (2017) stated that the loop and surface region display higher RMSF values. Based on the RMSF value, the residues of 173 until 182 and 276 until 283 from diosmetin and cosmosiin (Fig. 2D) that exhibited high fluctuation were located in the loop and surface region of the 3AIC protein (Supplementary Fig. S9a). Similar to the 3AIC-ligand complexes, residues from 1 until 23 and 482 until 489 were highly fluctuating in both ligands of 3IPK (Fig. 2J) settled in the loop and surface region of the 3IPK protein (Supplementary Fig. S9b).

The binding site residues in both diosmetin and cosmosiin were LEU 433, ARG 475, ASP 477, ALA 478, ASN 481, TRP 517, GLU 515, ASP 588, and HIS 587 (Supplementary Table S4) while in the 3IPK-ligand complex with genistin and 3′-methoxy-5′-hydroxyisoflavone-7-O-β-D glucoside these were ARG 824, TRP 816, SER 762, ASP 512, THR 586, LYS 822, ASN 820, and VAL 587 (Table S5). Regarding the RMSF result (Fig. 2E, F, K, and _L), all binding sites almost did not show any high fluctuations. Some research reported that binding residues within the structural protein have a lower RMSF value (Freitas and Schapira, 2017; Shawon et al., 2021). Moreover, De Vita et al. (2021) stated that a lower RMSF value in the binding site indicates strong interaction of the ligand to the protein. Therefore, these ligands from both proteins were able to stabilize strongly to the protein regarding the low fluctuation in the binding site residues.