Pharmacoinformatics-Based Approach for Uncovering the Quorum-Quenching Activity of Phytocompounds against the Oral Pathogen, Streptococcus mutans

Streptococcus mutans, a gram-positive oral pathogen, is the primary causative agent of dental caries. Biofilm formation, a critical characteristic of S. mutans, is regulated by quorum sensing (QS). This study aimed to utilize pharmacoinformatics techniques to screen and identify effective phytochemicals that can target specific proteins involved in the quorum sensing pathway of S. mutans. A computational approach involving homology modeling, model validation, molecular docking, and molecular dynamics (MD) simulation was employed. The 3D structures of the quorum sensing target proteins, namely SecA, SMU1784c, OppC, YidC2, CiaR, SpaR, and LepC, were modeled using SWISS-MODEL and validated using a Ramachandran plot. Metabolites from Azadirachta indica (Neem), Morinda citrifolia (Noni), and Salvadora persica (Miswak) were docked against these proteins using AutoDockTools. MD simulations were conducted to assess stable interactions between the highest-scoring ligands and the target proteins. Additionally, the ADMET properties of the ligands were evaluated using SwissADME and pkCSM tools. The results demonstrated that campesterol, meliantrol, stigmasterol, isofucosterol, and ursolic acid exhibited the strongest binding affinity for CiaR, LepC, OppC, SpaR, and Yidc2, respectively. Furthermore, citrostadienol showed the highest binding affinity for both SMU1784c and SecA. Notably, specific amino acid residues, including ASP86, ARG182, ILE179, GLU143, ASP237, PRO101, and VAL84 from CiaR, LepC, OppC, SecA, SMU1784c, SpaR, and YidC2, respectively, exhibited significant interactions with their respective ligands. While the docking study indicated favorable binding energies, the MD simulations and ADMET studies underscored the substantial binding affinity and stability of the ligands with the target proteins. However, further in vitro studies are necessary to validate the efficacy of these top hits against S. mutans.


Introduction
Streptococcus mutans is a gram-positive bacterium commonly found in the human oral cavity. While it is considered part of the normal microbial flora in the oral cavity, it is also the primary causative agent of dental caries [1,2]. Despite numerous studies that

Homology Modeling and Model Validation of the Target Proteins
SWISS-MODEL performs BLAST for each target protein against PDB and provides a list of templates. In the present study, for all the seven target proteins, model structures were built using SWISS-MODEL since none of them had a crystal structure available. For each of the selected target proteins, five templates were selected based on the identity and sequence coverage, and the 3D models were built.
The modeled structures that were built were assessed and validated using a Ramachandran plot using the SWISS-MODEL structure assessment tool, and one out of five models was selected for further studies (Supplementary Figure S1). The predicted structures of all the target proteins are depicted in Supplementary Figure S2.

Molecular Docking
Among the 110 ligands tested, 17 compounds were found to bind effectively with all the target proteins (Table 1). Among the 17 compounds that exhibited higher binding towards the selected targets of S. mutans, 15 compounds were from A. indica, and 1 each from M. citrifolia and S. persica. Among the 17 lead compounds, the top hits that exhibited lower binding energy and, hence, high binding affinity against each target were then selected for further simulation analysis ( Table 2). The 3D and 2D interactions between amino acid residues of the target proteins and their respective top hit compounds were also plotted using Discovery Studio. The list of amino acids present in each of the targets that interact with the ligands is presented in Table 3. The 3D and 2D plots depicting the interactions between target proteins and their ligands are presented in Figures 1 and 2, respectively.

Molecular Dynamics (MD) Simulation
MD simulation assessed the inter-molecular protein-ligand interactions under an artificial environment with specified thermodynamical conditions such as temperature, volume, density, and pressure for 100 ns time duration. Further, a final production step aided in the exploration of the structural modification of the complex. The analysis of unique parameters such as root mean square deviation (RMSD) (Figure 3), simulation interactions diagram (Figure 4), protein-ligand contact ( Figure 5), and timeline representation of interaction and contacts ( Figure 6) between the target protein and their respective ligands, aided in the analysis of structural changes at each level.    The hydrogen bonding interaction between amino acid residues of the target protein and their respective ligand analyzed in MD simulation was found to correlate with that of AutoDock results (Table 3). The simulation study provided insight into the stability of the interaction between the top hits and their respective target proteins.
In order to perform the post-MM-GBSA (molecular mechanics with generalized Born and surface area solvation) analysis of the free binding energy calculation, 0-2002 frames with a sampling size of roughly 10 steps were produced. During the MM-GBSA calculation of the 100 ns MD data of the target proteins with their respective ligands, a total of 201 frames were processed and analyzed. All the complexes showed good binding affinity, thereby validating the docking and MD results (Table 4).

ADMET Analysis
Analysis of ADMET properties of the top hits was performed using the SwissADME online tool and the pkCSM tool, and the results are presented in Tables 5 and 6, respectively. All six compounds had molecular weights greater than 400 g/mol. Except for meliantrol, which was predicted to be moderately soluble in water, all other selected bioactive compounds were predicted to be poorly soluble in water. Solubility of ligands aid in the proper solvation and absorption into the host body. This also aids in the formulation of a solvent or carrier for delivery. Due to their poor solubility, all the hits violated one out of five of Lipinski's rules. The RADAR plot ( Figure 7) shows the overall drug-likeness of the top hit ligands. Meliantrol passes all the criteria for drug-likeness; however, other ligands violate the drug-likeness due to their poor solubility in water. The BOILED-Egg plot ( Figure 8) shows all the compounds, except citrostadienol, in the area of intestinal absorption. However, only meliantrol was predicted to exhibit good GI absorption. None of the hits cross the blood-brain barrier. The results were correlated with that of pkCSM.

Discussion
Streptococcus mutans is one of the oral commensals that opportunistically cause dental caries. One of the pathogenic mechanisms adopted by S. mutans to cause dental caries is the formation of biofilm on the tooth surface which results in the erosion of tooth enamel through acid production. Biofilm formation requires quorum sensing, a densitydependent communication mechanism. The development of antimicrobial resistance in these pathogens can be avoided by using newer approaches rather than conventional antibiotics. These approaches should control the virulence of the pathogen instead of killing them. Hence, the objective of this study was to develop a compound that can control the quorum sensing in S. mutans rather than killing it. Pharmacoinformatics-based approaches such as drug target identification, molecular docking, and molecular dynamics simulations have accelerated the process of drug discovery. In the present study, we used bioinformatics tools and software to screen phytochemicals from A. indica, M. citrifolia, and S. persica for their potential as anti-quorum sensing agents against the selected target proteins of S. mutans.
The target proteins involved in quorum sensing were selected from our previous study [21]. These seven targets included CiaR (putative response regulator CiaR), LepC (signal peptidase I), OppC (putative transmembrane protein, permease OppC), SecA (protein translocase subunit SecA), SMU1784c (putative Eep protein-like protein), SMU_659 (putative response regulator SpaR), and YidC2 (membrane protein insertase YidC2). The target protein, CiaR, is a response regulator in a two-component signal transport system that controls several virulent characteristics of S. mutans. These virulent characteristics include mutacin I activity, oxidative stress tolerance, acid tolerance, and biofilm formation [22,23]. Studies on S. sanguinis showed the development of a fragile biofilm as a result of the CiaR mutation, which resulted in decreased polysaccharide synthesis [24]. The product of the lepC target gene is a signal peptidase that aids in the export of several virulent proteins. It has also been used as a housekeeping gene for several studies [10,25]. The target protein, OppC, is an oligopeptide permease of the ABC transporter family. It helps the bacteria to regulate XIP (sig X inducing peptide) production and in competence development [9,26]. SecA is a membrane protein translocase that helps the bacteria export proteins, leading to an increase in the virulence of the organism [10,27,28]. Another target, SMU1784c, plays an important role in the management of oxidative and acid stress, EPS production, and biofilm formation [11]. The SpaR protein is a response regulator of the spa (surface protein antigen) family, which is one of the virulence factors of S. mutans [28,29]. The membrane protein insertase of S. mutans, Yidc2, helps in EPS production and biofilm formation [10,30,31].
The resolved structures of the target proteins were not available in PDB or any other structural databases. Hence, the SWISS-MODEL online tool was used to predict the 3D structure of target molecules. Once the protein sequence is submitted as a query, SWISS-MODEL performs BLAST against PDB and gives a list of templates. Based on high identity and sequence coverage, five templates were selected for each protein, and the 3D models were built. A similar approach has been followed by researchers who have used SWISS-MODEL to model NOX 2 of S. mutans by using the crystal structure of NADH oxidase from S. pyogenes as a template [32,33]. Various proteins of S. mutans that were modeled using SWISS-MODEL include domain V of glucosyltransferase (GTF-SI) [34]; SMU.63, an amyloid-like secretory protein of S. mutans [35]; the Spase I protein and β-sheet-rich Nterminal collagen-binding domain (CBD) of Cnm, a collagen-and laminin-binding surface adhesin protein of S. mutans [36,37]. Since the models are predicted in silico, it requires validation before further processing. Hence, the Ramachandran plot from the structure assessment tool of SWISS-MODEL was used for validation of the predicted structures. All five models had more than 90% residues in the allowed region. The models from other templates with less than 90% residues in the allowed region were rejected. The remaining selected models were taken for further docking studies. The modeled structures of fibronectin/fibrinogen binding protein (FBP) from S. mutans, phospholipase D (F13) protein of monkeypox virus, 3-chymotrypsin and papain-like proteases of SARS-CoV2, and U box domain-containing protein gene (GsPUB8) from Glycine soja were all validated using the Ramachandran plot [38][39][40][41].
Computational drug discovery studies evaluate the binding of ligands with the target protein, but the agonist or antagonistic effect of ligands on the target protein is only validated through in vitro and in vivo experiments [42,43]. In the current study, based on an initial screening of 110 compounds using AutoDock, 17 high binders were selected that were found to bind with all the target proteins efficiently. Among the 17 compounds, 15 were from A. indica, and 1 each from M. citrifolia and S. persica. Campesterol, meliantrol, stigmasterol, isofucosterol, and ursolic acid were the top binders specific for the target proteins, CiaR, LepC, OppC, SpaR, and Yidc2, with a binding energy of −8.76, −10.16, −6.75, −9.1, and −9.73 kcal/mol, respectively. Citrostadienol was the high binder against two of the selected targets, SMU1784c and SecA, with a binding score of −8.45 and −9.88 kcal/mol, respectively. The source of ursolic acid is M. citrifolia, whereas all other leads are constituents of A. indica. Molecular docking is a versatile tool that is very useful in screening hundreds of compounds before testing the effective ones using in vitro studies. Molecular docking using AutoDock tools has been used previously to study both agonist and antagonist activity of various natural and synthetic ligands. In silico and in vitro agonistic activity of ligands have been studied for the treatment of diabetes [44,45], Parkinson's disease [46], and cardiac diseases [47]. The antagonistic activity of ligands against S. mutans [48,49], Leishmania donovani [50], Helicobacter pylori [51], and SARS-CoV-2 [52] has also been studied. All these research works corroborate the necessity of in vitro experiments in the validation of computational analysis. However, they also demonstrate that phytocompounds and synthetic compounds can cause competitive or non-competitive inhibition of target proteins involved in microbial diseases.
The RMSD of the protein-ligand complex is plotted to evaluate the stability of the interaction between the protein and the ligand. The RMSD plot of target protein CiaR in complex with its top binding ligand, campesterol, displayed a fluctuation in RMSD up to 7 Å to 13 Å (Figure 5a). Ligand RMSD was stable, and fluctuations were between 9 and 10.5 Å. The complex of LepC and its top hit ligand meliantrol shows RMSD fluctuation up to 8 and 12 Å (Figure 5b). Ligand RMSD was stable, and fluctuations were between 4 and 10.6 Å. The simulation of target protein OppC in complex with stigmasterol shows the stability of protein at 9 to 16 Å (Figure 5c). Ligand RMSD was stable, and fluctuations were between 9 and 1.6 Å. In the SecA-citrostadienol complex, the protein remains stable at 3 to 4 Å (Figure 5d), whereas the ligand fluctuation is at 14-16 Å and, hence, shows stable interaction between the protein and the ligand. The SMU1784c-citrostadienol complex shows fluctuation up to 2 to 4 Å (Figure 5e), whereas the ligand fluctuation is at 3-6 Å and, hence, shows stable interaction between the protein and the ligand. The Yidc2 and ursolic acid complex shows good RMSD results. The protein fluctuation is up to 6 Å and the ligand up to 9 Å. The RMSD plot converges from 20 ns till the end of the simulation at 100 ns, thus showing a stable interaction between the protein and the ligand. The hydrogen bond and other interactions plotted in the simulation interactions diagram correlate with the docking results. The hydrogen bond interaction between residues ASP86, LEU155, and PRO101 in CiaR, OppC, and SpaR with their respective ligands was observed both in AutoDock and MD simulations. Similarly, in the LepC protein, amino acid residues ILE69 and GLU178 were observed. In SecA, SMU1784c, and Yidc2 though there were no common residues binding through a hydrogen bond, the same residues bind with ligands through other types of bonds. The highly interacting amino acids were ASP86 of CiaR, ARG182 of LepC, ILE179 of OppC, GLU143 of SecA, ASP237 of SMU_1784c, PRO101 of SpaR, and VAL84 of YidC2. Similarly, molecular dynamic simulations and an energy calculation method have been used by researchers to study the LPXTG sequence in the C-terminus of surface proteins, the substrate of the cysteine transpeptidase sortase A (SrtA) enzyme, to better understand how leucine residue affects the dynamics of the enzyme-substrate complex structure. According to the findings, the substrate's 'Leu' residue appears to be essential for anchoring and guiding the conformational shift of the enzyme SrtA [53]. Molecular docking and dynamics simulation studies have been exploited to study the inhibition of glucan sucrase-mediated biofilm formation of S. mutans by thiosemicarbazide derivatives [54]. Similar techniques have also been used to investigate the stability of phosphodiesterase type 5 (PDE5) in complexes with bioactive compounds from Mimosa pudica to understand their aphrodisiac performance [55]. Similarly, a pharmacoinformatics-based molecular docking and dynamics simulation analysis of bioactive components from Indian cuisine, rasam, was conducted against MAPK6 (mitogen-activated protein kinase 6), a family of serine/threonine protein kinases that is crucial in regulating extracellular signaling into a variety of cellular functions, including ROS production [56].
MM-GBSA also validates the molecular docking and MD simulation studies as it shows binding energy ranging from −49 to −79 kcal/mol in all the protein-ligand complexes studied. MM-GBSA has previously been used to study and validate in silico interaction of ligands with SARS-CoV-2 protease [57]. A similar MD simulation approach has also been utilized for screening substrate analog inhibitors of L-Ornithine-N5-monooxygenase (PvdA) to control Pseudomonas aeruginosa infections [58].
Analysis of ligands for ADMET properties using SwissADME and pkCSM shows that the molecular weight of all six hits was larger than 400g/mol. Only meliantrol was expected to have moderate water solubility. All other hits had low water solubility. The right solvation and absorption into the host body are made possible by the ligands' solubility. Additionally, it helps with the formulation of the solvent or delivery vehicle. Except for meliantrol, none of the selected compounds exhibited good GI absorption. All the hits broke one out of five of Lipinski's rules because of their poor solubility. The Bioavailability RADAR plot depicts the overall drug-likeness of the top binding ligands. Meliantrol fulfills all criteria for drug-likeness, while other ligands fail due to their low solubility in water. The BOILED-Egg plot displays all of the hits in the area of intestinal absorption except citrostadienol. None of the ligands were predicted to penetrate the blood-brain barrier. Analysis of ADMET properties of an array of ligands using Swiss ADME and pkCSM tools has been previously followed by many researchers to control S. mutans biofilm [59,60].

Target Proteins in S. mutans
The target proteins were selected based on our earlier study, which demonstrated an in silico subtractive proteomics approach for screening drug targets in S. mutans UA159. Among the 13 novel drug targets that were identified, 7 proteins were found to be involved in quorum sensing [21]. Hence, these target proteins have been subsequently used in the current study. The target proteins include two putative response regulators, a signal peptidase, a putative EEP protein-like protein, a putative transmembrane protein-permease, a protein translocase subunit, and a membrane protein insertase (Table 7).

Selection of Ligands
The bioactive compounds from three Indian medicinal plants, Azadirachta indica, Morinda citrifolia, and Salvodora persica, were selected as ligands from the literature, as well as various databases such as Dr. Duke's Phytochemical and Ethnobotanical Databases (https: //phytochem.nal.usda.gov/; accessed on 11 April 2022), IMPPAT (https://cb.imsc.res. in/imppat/; accessed on 13 April 2022), PubChem (https://pubchem.ncbi.nlm.nih.gov/; accessed on 15 April 2022), and Drug Bank (https://go.drugbank.com/; accessed on 18 April 2022). However, among the hundreds of compounds, the ones with a background of antimicrobial activity in the literature were selected. The screening process involved manually searching each compound in the PubMed database for its antimicrobial activity against S. mutans or other oral pathogens. Finally, 110 compounds were shortlisted for analysis. The 3D structures of these ligands were either collected from PubChem or drawn using ChemSketch Version 12.00 (Table 8).

Computational Analysis
Computational analysis was performed using Debian Linux Operating System running on a 3.10 GHz Intel ® Core™ i5-4440 CPU (Acer, Bengaluru, India) with 8 GB RAM. For the analysis of docking results, the Accelrys ® Discovery Studio Visualizer (Accelrys Software Inc., San Diego, CA, USA) that runs on the Windows operating system was used.

Homology Modeling and Model Validation
The 3D structures of target proteins were modeled using the SWISS-MODEL online tool (https://swissmodel.expasy.org/ accessed on 20 April 2022) [61]. The modeled structures were assessed and validated using the Ramachandran plot analysis.

Molecular Docking
The target proteins were prepared by removing the water molecules using Discovery Studio. The parameters for docking were set up by using the graphical user interface (GUI) of AutoDock Tools [62]. The preparation of the target protein structure involves the addition of hydrogen atoms and formal charges. The number of rotatable bonds and determination of root for the ligand were set to default. Later, molecular docking was performed using AutoDock 4.2. The AutoGrid4 program of AutoDock allowed the generation of grid maps for target proteins embedded in a three-dimensional grid of manually set parameters (Table 9). The binding energies between the target protein and ligand were calculated by running the AutoDock4 program using pre-set grid maps. The result analysis was performed using the GUI of AutoDock Tools and Discovery Studio. The binding energies for various conformations of the ligand with the target proteins were determined, and the best conformation was chosen based on the binding energy and the number of hydrogen bonds that they formed with the protein.

Molecular Dynamics
The protein-ligand complexes (PLCs) were prepared using PyMOL [63], and these PLCs were taken as input for Molecular Dynamics (MD) simulations (Table 10). The MD simulation of PLCs aids in the visualization of target binding sites and also provides information concerning the binding stability of the PLCs under physiological conditions. The Desmond module of Schrödinger was used to perform MD simulations (academic license, Version 2020-1) [64]. Initially, an explicit water model was prepared as an orthorhombic simulation box with Simple Point-Charge (SPC). This was designed with the builder panel of the system in such a way that a minimum distance of 10 Å is maintained between the protein surface and the solvent surface. Then, the PLCs with receptors were solvated with the orthorhombic TIP3P water model [65]. The addition of counter ions and limitation of salt concentration of the physiological system to 0.15M was performed for the neutralization of the solvated system (Table 11).  Then, the PLC system was designated with the OPLS AA force field [66]. For the RESPA (reversible reference system propagator algorithms) integrator [67], Nose-Hoover chain thermostat [68], and Martyna-Tobias-Klein barostat, a relaxation time of two seconds was used. The equilibrated system was used to perform the final production process of MD simulations. In this study, the default parameters set for relaxation before MD simulations include a time duration of 100 ns, temperature of 310 K, a pressure of 1.0 bar with NPT (Isothermal-Isobaric ensemble, constant temperature, constant pressure, constant number of particles) ensemble [69] (Table 12). The MD simulation tool was used to run the simulation, and the output file was used to retrieve trajectories and create the movie. The output file in .cms format was imported into the software, and the created simulation movie was exported at a higher resolution of 1280 × 1024 px for better quality. The trajectory of the overall MD simulation process was written in 1000 frames. The frames of the protein backbone were aligned to the backbone of the first frame. This provides a better understanding of the stability of the PLCs during MD simulations. Finally, the Root Mean Square Deviation (RMSD) and Root Mean Square Fluctuation (RMSF) of PLCs and the simulation interaction diagrams were analyzed [57,70]. For each protein-ligand complex, three simulations studies were run, and the best run was taken for analysis.
The free binding energies of the protein and ligand complexes were also examined using Molecular Mechanics, the Generalized Born model, and Solvent Accessibility (MM-GBSA). The "thermal_MMGBSA.py" script and "MM-GBSA ∆G Bind:" parameter from the Prime/Desmond module of the Schrodinger suite were employed to carry out the post-MM-GBSA analysis [71].

ADMET Analysis
The ADMET properties of the high binding ligands were analyzed using the SwissADME online tool (http://www.swissadme.ch; accessed on 10 May 2022) [72] and pkCSM online tools (http://biosig.unimelb.edu.au/pkcsm/prediction; accessed on 10 May 2022) [73]. The input was given as SMILES notation of the ligands, and the output was retrieved as a ".csv" file. In SwissADME, the BOILED-Egg plot and radar plot for all the leads were also retrieved as a single image [74]. The ADMET analysis helps to predict the drug-likeness of the top hits.

Conclusions
S. mutans is a major contributor to dental caries. In this study, an in silico approach was employed to identify a plant metabolite with potential efficacy against specific target proteins of S. mutans involved in quorum sensing. The selected target proteins included: Membrane protein insertase YidC 2, Permease OppC, Putative Eep protein-like protein, Putative transmembrane protein, Putative response regulator CiaR, Putative response regulator SpaR, and Signal peptidase I LepC. The three-dimensional models of these target proteins were generated using SWISS-MODEL and validated using the Ramachandran plot. Plant metabolites derived from A. indica (Neem), M. citrifolia (Noni), and S. persica (Miswak) were evaluated for their potential binding affinity to the selected target proteins. Molecular docking studies were performed using AutoDock Tools. From a total of 110 ligands, 6 top hits were identified for each target protein. These ligands, namely campesterol, meliantrol, citrostadienol, stigmasterol, isofucosterol, and ursolic acid, were subjected to molecular simulation analysis to assess their stability and interaction patterns. The highly interacted amino acid residues identified were ASP86, ARG182, ILE179, GLU143, ASP237, PRO101, and VAL84, corresponding to the proteins CiaR, LepC, OppC, SecA, SMU1784c, SpaR, YidC2, respectively. Furthermore, the ADME characteristics of the identified ligands were evaluated using the SwissADME program. Collectively, the results suggest that these phytosterols have the potential to serve as quorum-quenching agents. However, further in vitro analysis is required to confirm their efficacy against S. mutans and to enable their application in the treatment of dental caries.