Identification of New Potential Inhibitors of Quorum Sensing through a Specialized Multi-Level Computational Approach

Biofilms are aggregates of microorganisms anchored to a surface and embedded in a self-produced matrix of extracellular polymeric substances and have been associated with 80% of all bacterial infections in humans. Because bacteria in biofilms are less amenable to antibiotic treatment, biofilms have been associated with developing antibiotic resistance, a problem that urges developing new therapeutic options and approaches. Interfering with quorum-sensing (QS), an important process of cell-to-cell communication by bacteria in biofilms is a promising strategy to inhibit biofilm formation and development. Here we describe and apply an in silico computational protocol for identifying novel potential inhibitors of quorum-sensing, using CviR—the quorum-sensing receptor from Chromobacterium violaceum—as a model target. This in silico approach combines protein-ligand docking (with 7 different docking programs/scoring functions), receptor-based virtual screening, molecular dynamic simulations, and free energy calculations. Particular emphasis was dedicated to optimizing the discrimination ability between active/inactive molecules in virtual screening tests using a target-specific training set. Overall, the optimized protocol was used to evaluate 66,461 molecules, including those on the ZINC/FDA-Approved database and to the Mu.Ta.Lig Virtual Chemotheca. Multiple promising compounds were identified, yielding good prospects for future experimental validation and for drug repurposing towards QS inhibition.


Introduction
Bacterial biofilms are aggregates of microorganisms anchored to a surface and embedded in a self-produced matrix of extracellular polymeric substances [1]. Due to the properties of the matrix and to the intercellular interactions between the bacteria within, the biofilm develops, thus becoming increasingly sophisticated and offering the constituting bacteria several advantages over their planktonic counterparts [2,3]. These include an increased antibiotic resistance, elevated levels of lateral gene transfer, higher stress resistance and subversion of host defense mechanisms [4][5][6]. Quorum sensing (QS), a process of cell-to-cell communication, is important for the control of several virulence factors and plays a key role during biofilm formation. In this process, cells communicate using auto-inducer signals. QS allows bacteria to synchronously adjust gene expression to alter their behavior in response to changes in population density and the surrounding bacterial community [7].
Several clinically important pathogenic bacteria, such as cystic fibrosis-associated Pseudomonas aeruginosa, surgical sites associated with Staphylococcus aureus and many others, cause infection through biofilm formation [8]. This can have devastating consequences. Microbes that reside in biofilms may not be eliminated by traditional antibiotics because of metabolic dormancy or molecular resistance mechanisms [9]. The US National Institutes of Health estimate that 80% of all bacterial infections occurring in the human body are

PDB Structures Identification
The structures were identified using the Biofilms Structural Database [38] and then downloaded from the Protein Data Bank [39]. These X-ray structures and the corresponding resolution value, ligand and bacterial strain are summarized in Table 1. The CviR structures belong to two different strains. Additionally, different ligands with different activities are complexed with each structure. 3QP1 is complexed with its native ligand C6-HSL, which is a full agonist in strain 31532. 3QP2 and 3QP4 are bonded to ligands with longer acyl chains, which fail to fully activate CviR. C8-HSL leads to 40% of the original activity, and C10-HSL elicits only 6%. These ligands also work as a partial antagonist in the presence of C6-HSL. 3QP5 is bonded to chlorolacetone, an even stronger antagonist than C10-HSL. 3QP6 is bonded to C6-HSL, which functions as an antagonist on strain 12472 [20]. Finally, 3QP8 is bonded to C10-HSL. This is an agonist, which is located closer to this strain's native ligand (3-hydroxy-C10-HSL). The proteins from each strain are 87% identical in the amino acid sequence, with its difference in autoinducer being partially due to a naturally occurring serine instead of a methionine at position 89. This is a key residue that occupies the opening of the ligand-binding pocket. The importance of Met 89 is reinforced by the fact that, when using antagonists, the side-chain swings away. The extension of this movement increases with the strength of the antagonist. It was observed that Met89 s side-chain swings away from the ligand-binding pocket along with an increase in length of the ligand's acyl chain. The extent of this variation relates to the antagonist capabilities of the ligand. 3QP1 is bonded to the native ligand, C6-HSL, and Met89 is in its original pose. In 3QP2, there are two conformations: one similar to 3QP1 and the other in which there is a small, intermediate variation on the side-chain orientation. In both 3QP4 and 3QP5, the side-chain has fully changed its position, swinging far away from the center of the binding pocket [20]. The receptors were then isolated from the ligands using PyMOL [40]. In the case of 3QP5 and 3QP8, it was also necessary to isolate one of the four chains that are in the PDB file. In both cases, chain A was the one selected.

Protein-Ligand Docking Protocol Validation
For the validation of the molecular docking protocol, redocking and cross-docking studies are performed. The goal of the redocking studies is to see if the molecular docking protocol can accurately reproduce the experimental complexed pose of the known ligands. Several parameters, such as the size and coordinates of the search area, are evaluated, and the results were compared. The RMSD between the re-docked and the original poses was calculated to better evaluate the results. In cross-docking studies, each ligand from each structure is docked into all the other structures. Considering that all structures are of the same protein, this is an assessment of how accurately each structure can accommodate in its pocket the ligands from the other structures, i.e., its general usefulness [41]. As the performance of different docking programs/scoring functions is known to depend among other features on the specific type of target under consideration [42], in this study, four popular protein-ligand docking programs [43,44] and seven different scoring functions were taken into consideration: AutoDock 4 [45]; AutoDock Vina (Vina) [46]; LeDock [47]; and GOLD (with the scoring functions CHEMPLP, ASP, CHEMSCORE and GOLDSCORE) [48]. The RMSD calculations were done using DockRMSD [49]. For AutoDock and Vina, the search area was a box designed using the AutoDock/Vina plugin for PyMol by Daniel Seeliger [50]. This search area was based on the residues known to be involved in the interaction between the ligand and CviR. These are Tyr80, Trp84, Asp97 and Ser155. For LeDock, the same box dimensions were considered. With the different GOLD scoring functions, the search area considered involved spheres with the same center coordinates and volume as in AutoDock, Vina and LeDock. With each docking program/scoring function, the protocol was optimized for each protein target to minimize the RMSD in the redocking predictions of the reference ligand by comparison with the crystallographic structure of the corresponding complex. The optimized parameters for each program/scoring function were: AutoDock-docking box position, docking box dimensions; Vina-docking box position, docking box dimension, exhaustiveness; LeDock-docking box position, docking box dimension; GOLD (CHEMPLP, ASP, CHEMSCORE, and GOLDSCORE)-binding pocket center, docking region radius, search efficiency, number of runs.
It should be noted that the different docking program results have different values and scales. The score for all the GOLD scoring functions is dimensionless, and the higher the score, the better the binding affinity. AutoDock, Vina and LeDock scoring functions, on the other hand, use a metric that resembles binding free energy (in kcal/mol), so a more negative value means a better-predicted affinity.
The best combinations of protein-ligand/scoring functions were used in the next stage.

Virtual Screening Protocol Optimization
The optimization of the virtual screening protocol was performed via a virtual screening benchmark test on the CviR system, considering a library of experimentally confirmed active molecules for CviR and randomly generated decoys [42], i.e., molecules with similar 1-D physicochemical properties to known actives (e.g., molecular weight, calculated LogP), but dissimilar 2-D topology to be likely non-binders.
A total of 46 actives ( Figure S1 in the supplementary materials) were used, 23 obtained from ChEMBL [51] and 23 found in the literature [28,[52][53][54][55][56][57][58][59][60][61][62]. Structures were downloaded from PubChem [63] or modeled using the Avogadro molecular editing program [64]. These are listed in the supplementary materials. To obtain the decoys, the 46 actives were submitted to the Database of Useful Decoys-Enhanced (DUD-E) server [65]. DUD-E generated 50 decoys per each active. This led to a final library of 2346 molecules, constituting a CviR-specifically designed and challenging test set to evaluate and optimize the performance of different virtual screening protocols in their ability to discriminate between CviR binders and non-binders.
Different metrics are routinely used to evaluate the performance of VS protocols, including receiver operating characteristics (ROC) curves and the respective area under the curve (AUC), enrichment factor (EF), robust initial enhancement (RIE), Boltzmann enhanced discrimination of ROC (BEDROC), predictive curves (PC) and total gain (TG) [42,66]. The calculation of the metrics to evaluate the performance of VS protocols was performed using two different applications. The ROC and enrichment factors were calculated on Microsoft Excel. The web-based application Screening Explorer [66] was used for calculating Total Gain and BEDROC.
In this step, Vina, all scoring functions of GOLD and LeDock were used, with Vina, GOLD with the CHEMPLP scoring function (GOLD/CHEMPLP) and LeDock being the best performing VSs. The optimized parameters for Vina were an exhaustiveness parameter of 8 and a search box with the following dimensions: 15.0 × 15.0 × 14.0 Å. For GOLD, the binding pocket radius was set to 9.1 Å, and the number of independent genetic algorithms (GA) runs set to 10. For each GA run, the population size was 100 with a selection pressure of 1.1. The number of operations was 100,000, with the crossover, mutation and migration frequency being 95, 95 and 10, respectively. For LeDock, a clustering by RMSD was performed, with a value of 1 Å set as the threshold, and the number of selected binding poses was 20.

Virtual Screening for the Identification of New Potential CviR Binders
Three optimized virtual screening protocols selected in the previous stage and based on Vina, GOLD/CHEMPLP and LeDock were used for identifying new potential CviR binders. Two virtual databases were considered: the US Food and Drug Administration (USFDA) approved database available on ZINC [67], and the Mu.Ta.Lig. Virtual Chemotheca [68]. ZINC is a public access database, which is used to obtain compounds for several uses. These include virtual screening, ligand discovery and force field development. It was created in the Department of Pharmaceutical Chemistry at the University of California (UCSF) and included 220 million molecules [67]. In the present work, 1657 FDA-approved molecules were available in this database and were assessed in this work to evaluate the possibility of drug repurposing. The Mu.Ta.Lig. (Multi-Target Ligand) Chemotheca database was created from COST Action CA15135-Multi-target paradigm for innovative ligand identification in the drug discovery process and was developed with the goal of identifying multi-target agents and repurposing known compounds. A large number of molecules with promising pharmaceutical relevance are developed every year and are forgotten when they fail to have the desired effect. However, these molecules can have a positive effect on other targets. At the time of this work, Mu.Ta.Lig Virtual Chemotheca features 64,804 compounds [68].
From the virtual screening results, a selection of 10 molecules from the FDA-approved database and 10 molecules from the Mu.Ta.Lig Virtual Chemotheca was done among the top results in all VSs, taking into consideration the chemical diversity and predicted physical-chemical properties of the molecules.

Molecular Dynamics (MD) Simulations
Molecular dynamics (MD) simulations were carried out using the Amber18 software package [69]. Molecular mechanics parameters for the selected molecules were assigned using ANTECHAMBER and considering the general amber force field [70], with RESP charges calculated using HF/6-31G(d) with Gaussian16 [71]. The protein was described through the amber14sb force field [72]. The protein-ligand complexes were placed into a box of TIP3P water molecules, whose edges are placed at least 12 Å away from each atom of the complex. All the MD were performed using periodic boundary conditions. Long-range electrostatic interactions were calculated using the particle mesh Ewald summation method. The cutoff value for the short-range electrostatic and Lennard-Jones interactions was set at 10.0 Å. All bonds involving hydrogen atoms were constrained using the SHAKE algorithm, enabling applying a 2 fs time step.
All ligand-CviR solvated complexes went through 4 minimization steps, with a maximum of 2500 cycles each. After 1250 cycles, the minimization method was switched from steepest descent to conjugate gradient. In the first minimization, all molecules except water molecules were restrained. In the second minimization, only hydrogen atoms were not restrained. During the third minimization, only the protein backbone was restrained. Finally, for the last minimization step, there were no restraints. Following the minimization phase, all solvated complexes were heated from 0 to 310.15 K over 50 ps. They were further equilibrated at 310.15 K during 50 ps to stabilize the density. Finally, for each ligand-CviR complex, the production phase was run for a total of 100 ns in an NPT ensemble at a pressure of 1 bar and a temperature of 310.15 K.
The analysis of the trajectories was carried out using the cpptraj tool [73] and VMD [74].

MM/PB(GB)SA Calculations and Replica MDs
The MM/PB(GB)SA calculations were performed using the MM/PBSA.py script available in AMBER [75]. The calculations considered the last 40 ns of the MD simulation of every complex, using an interval of 100 ps. This means that the program used every 10th frame of the simulation. In the MM/PBSA calculations, the following constants were used: ionic strength of 0.100 mol dm −3 , the external dielectric constant of 80.0 and the internal dielectric constant of 4. In the MM/GBSA calculations, a salt concentration of 0.100 mol dm −3 was used.
The Free energy decomposition option was used to obtain information about the local interactions of the complex. Using per-residue decomposition, the contribution of each residue to the total free energy was estimated.
For the 4 ligands with the strongest binding free energy identified, 3 additional 100 ns replica molecular dynamics simulations were performed, using the protocol and approximation described in 2.5. Table 2 compares the best redocking performance obtained with Autodock4, Vina, CHEMPLP, CHEMSCORE, GOLDSCORE, ASP and LeDock in terms of RMSD between the predicted pose and the X-ray crystallographic pose. One important detail is that 3QP2 was separated into two different files, one for each conformation of the important residue Met89. 3QP2a has the intermediate conformation, while 3QP2b has a similar conformation to 3QP1. The RMSD values show that the best docking program for reproducing experimental poses of the ligand is LeDock, with an average RMSD of 0.74 Å. All scoring functions from GOLD also have a good performance. On the other hand, the programs, which performed worse were Vina and AutoDock 4, with an average RMSD of 1.44 Å and 1.13 Å, respectively. Vina has the worse overall RMSD value for the redocking of 3QP6. In this case, the program placed the molecule the opposite way when compared to the crystallographic structure. However, this could be due to the fact that a smaller molecule is being docked in a binding pocket previously accommodating a larger ligand. In reality, when C6-HSL binds to this structure (in this strain, C6-HSL functions as a partial antagonist), there is a shift in the DNA-binding domain, which will help fill the binding pocket [20]. Nevertheless, most RMSD values are below 2 Å, which indicates a good overall performance from all the different molecular docking software/scoring functions. This good performance is mainly seen with the lactone head group, with most programs struggling with the acyl chain. The scores from all redocking procedures are available in Table S1 in the supplementary materials.

Results and Discussion
The cross-docking studies (details on Tables S2 and S3 in the supplementary materials) show that, on average, the best results are seen on structures 3QP6 and 3QP8, both from strain 12472. From strain 31532, the structure, which generates the best scores, is 3QP4. In contrast, the worst results in all ligands and software are seen on structure 3QP5. This may be due to the worse resolution of this structure. Similar to the redocking studies, C6-HSL shows less affinity when compared to all other ligands, while C10-HSL and CL show the best results. This behavior is shown in most structures and software. One major difference from the redocking and cross-docking studies is the behavior of 3QP6. As expected, the bad performance during the redocking studies was due to the ligand, C6-HSL. When docking other ligands, 3QP6 is frequently the best structure that leads to the best results.
Globally, the redocking and cross-docking studies suggest that most combinations of docking program/scoring function, when optimized, perform similarly. For these structures, the poorest results were obtained with AutoDock 4 that is also the most timeconsuming program.
Regarding the structures, the main conclusions emerging from these results are that both structures from strain 12472 (3QP6 and 3QP8) displayed the highest scores with all scoring functions, from which 3QP6 is the structure that generated the best results. From strain 31532, 3QP4 displayed the higher and most consistent scores for all scoring functions, while 3QP1, 3QP2a and 3QP2b presented more variable scores. Finally, 3QP5 was consistently the worst performing of all available structures of CviR.
After evaluating and optimizing the docking of known CviR ligands with wellcharacterized and experimentally defined binding poses, the next step was to evaluate the ability of this protocol in distinguishing between known CviR binders and non-binders.
Tables 3 and 4 compare the ability of Vina, CHEMPLP, CHEMSCORE, GOLDSCORE, ASP, LeDock, as well as all the different PDB structures, in distinguishing experimentally confirmed CviR binders (actives) and non-binders (decoys) through the usage of several different evaluation metrics. The best-performing programs were consistent across most of the metrics, with Vina and GOLD with the CHEMPLP scoring function (GOLD/CHEMPLP) being more capable of distinguishing the experimentally confirmed binders from the decoys, generally followed by LeDock. With this in mind, the programs used in the virtual screening stage of this work were Vina, GOLD/CHEMPLP and LeDock (see Figure 1).
Finally, it is important to select the receptor structures. Across nearly all metrics, the best discrimination ability between binders and non-binders was obtained when using 3QP6 and 3QP8. However, since they are both from strain 12472, the decision was made to use the best performing structure from each strain. The structure from strain 31532 that generally generated the best results was 3QP4. Therefore, the two CviR structures used in the virtual screening step of this work were 3QP4 and 3QP6.
The optimized protein-ligand docking and virtual screening protocols were then applied to two libraries of compounds: the US Food and Drug Administration (USFDA) approved database available on ZINC (ZINC/FDA) to consider possible drug repurposing alternatives, and the more general Mu.Ta.Lig. Virtual Chemotheca. The results of all virtual screenings are in Figure 2. In this study, different molecular docking programs and different PDB structures were used in a consensus approach. This approach can explore many potential ligands with greater accuracy, contributing to a more rigorous selection of the molecules for testing. From the top solutions in the ranked lists obtained from the VS, 10 molecules from each library were selected for the MD simulations. For each library, the selection was based on the position of the molecules on the different VS protocols. The choice took into consideration the following criteria: (i) molecules among the top 25 solutions for the two structures considered in the VS; (ii) molecules among the top 25 solutions with more than one scoring function; (iii) molecular diversity taking into consideration the physical-chemical properties ( Figure 3). The molecules from the ZINC/FDA database were atovaquone, famotidine, iloprost, mebendazole, mirabegron, montelukast, paliperidone, pimozide, glycerol phenylbutyrate, and sulfasalazine ( Figure 4). The molecules from Chemotheca were CMLDID2574, CML-DID5450, CMLDID17434, CMLDID18049, CMLDID23812, CMLDID35542, CMLDID38590, CMLDID40723, CMLDID50121 and CMLDID60399 ( Figure 4). As a display of the variability of the selected molecules, multiple important properties (cLogP, molecular mass, total surface area, hydrogen donors and acceptors) of the 20 selected molecules are displayed in Figure 3. These properties and graphs were calculated using the data visualization and analysis program, DataWarrior [76]. Comparisons between the score of the 20 selected molecules and the scores obtained by the known active molecules is available in Figures S2-S7 in the supplementary materials.
Based on the EF 1% for the actives/decoys training set, we estimated the minimum number of actives expected to be found on the 25 first positions of each VS ( Table 5).
The next step of this work involved 100 ns of molecular dynamics simulation for each of the above-mentioned molecules in complex with 3QP6. This was done to validate the protein-ligand docking results, to evaluate the structural stability of the protein-ligand complex and to carry out the MM/PB(GB)SA calculation. For reference, MD simulations of 3-hydroxy-C10-HSL and C10-HSL in complex with the protein were also performed. While the native ligand of the strain 12472 of C. violaceum is 3-hydroxy-C10-HSL, C10-HSL is an agonist for this protein [20,77]. Additionally, this ligand was the agonist that generated the highest scores on the cross-docking studies with 3QP6. Consequently, C10-HSL, together with 3-hydroxy-C10-HSL, was selected as the reference ligands. Since these MD simulations were used for the refinement of the virtual screening results, only the ligand-binding domain was considered. In the future, we intend to perform further studies for a better understanding of how these ligands could affect the DNA-binding domain.
To assess the structural stability of the complexes, RMSD calculations were performed for the Cα atoms of each complex and for the ligands.  Considering the selected ligands from the ZINC/FDA database and from the Chemotheca database, all complexes exhibit low RMSD values through the simulation. Most ligands also display low RMSD values than the docking prediction (Table 6). However, multiple molecules display higher values, suggesting an induced-fit adjustment to the CviR-binding pocket during the MD simulation. The low standard deviation confirms that after the initial 40 ns of simulation, the ligand positions are well stabilized. This is further confirmed by the solvent-accessible surface area analysis. An increase in SASA from the initial pose predicted from docking would imply that the ligand was exiting the binding pocket. Fortunately, all ligands display a stable SASA value along the simulation. This, together with the RMSD results, confirms that all the selected ligands form stable complexes with CviR. Further information on the MD simulations is available in the supplementary materials ( Figures S8-S11 and Table S4). The final step of this study was to perform a MM/PB(GB)SA analysis for each ligandreceptor complex to determine the binding free energy of each ligand. For these calculations, only the last 40 ns of MD simulation were considered. Results are presented in Table 7 and Figure 5.
The results demonstrate that most ligands present a lower affinity towards CviR than the reference ligands (MM/GBSA predicted values between −22.7 ± 0.1 and −67.2 ± 0.2, against −49.0 ± 0.2 and −51.0 ± 0.2 kcal/mol and MM/PBSA predicted values between −13.0 ± 0.1 and −89.6 ± 0.2, against −51.7 ± 0.2 and 59.6 ± 0.2 kcal/mol) The only ligand from the ZINC/FDA database, which displays higher affinity towards CviR is pimozide (−67.2 ± 0.2 and −89.7 ± 0.2 kcal/mol, respectively). Pimozide is used as an antipsychotic agent and for the suppression of vocal and motor tics in patients with Tourette syndrome. Although its exact mechanism of action is unknown, it is thought that it inhibits the dopamine D2 receptor [78].   To provide additional confirmation of these results, and to better assess if the best performing ligands are indeed true binders, replicate MD runs were performed. For each of the four best-performing ligands, three additional MD simulations were performed, and the average RMSD, RMSF, SASA and number of H bonds were calculated. These results are shown in Table 8.
Across all evaluated metrics, the replicas display similar results to the original runs. Furthermore, the percentage of SASA buried on all replicas gives additional strength to the prediction that these ligands are true binders and remain tightly bound.
To further analyze the affinity between the ligands and the receptor, the overall Gibbs energy of association was decomposed into the contribution associated with each residue (Figures 6 and 7). On 3-hydroxy-C10-HSL, the residues with the higher contribution to the Gibbs energy of association are Tyr80, Trp84, Tyr88, Asp97 and Ser155. As for C10-HSL, using both methods, the amino acids that have a bigger impact are Tyr80, Trp84, Tyr88 and Ser155. When comparing with the already known interactions map obtained with LigPlot+ [79] (Figure 8), these results show great similarity. The only difference is the contribution of Asp97 in the interaction with C10-HSL. While this residue is considered to have an important role in the binding of ligand to CviR, the calculations show a less favorable result for Asp97.
Overall, the residues with the greatest contribution for the affinity of the ligands towards CviR are Met72, Tyr80, Trp84, Leu85, Tyr88 and Ser155. For pimozide, the bigger contribution is obtained from Tyr80 (MM/PBSA), Tyr88 (MM/GBSA), and Asp97. Whereas on most ligands, the calculations show a less favorable contribution from this residue, with pimozide, it is the major contribution for the affinity between the ligand and the protein.
For CMLDID17434, the residues with the highest contribution to the overall result using both methods are Leu85, Tyr88 and Ser155. In the case of CMLDID50121, Tyr80, Leu85 and Tyr88 are the residues with the most impact in the predicted affinity towards CviR. Lastly, for CMLDID60399, the amino acids with the greater impact are Met72, Tyr80, Tyr88.  Figure 6. Per-residue decomposition of the free energy calculations using MM/GBSA for the best performing ligands. Table 3. hydroxy-C10-HSL; the green dotted line represents C10-HSL.  The molecular dynamics simulations and the MM/PB(GB)SA resulted in four molecules with higher or comparable binding affinities to the reference ligands, 3-hydroxy-C10-HSL and C10-HSL. These molecules were pimozide (Figure 9) from the ZINC/FDA database and CMLDID17434 (Figure 10), CMLDID50121 ( Figure 11) and CMLDID60399 ( Figure 12) from Chemotheca.

Conclusions
Because microorganisms embedded in biofilms have several advantages, infections associated with biofilms have been accepted as a significant danger to our society. The recalcitrance of these structures towards existing antimicrobial approaches made necessary the discovery of novel methods to inhibit their mechanisms of formation. Inhibiting the formation of biofilms by disrupting quorum-sensing is the most promising strategy. With that in mind, this work was to use computer-assisted drug design to find promising molecules to block quorum-sensing and, therefore, prevent biofilm formation. This was achieved by optimizing molecular docking, and virtual screening protocols focused on the quorum-sensing receptor from Chromobacterium violaceum, CviR.
The optimized molecular docking and virtual screening protocols were applied to two databases, which generated multiple promising molecules that were then analyzed via molecular dynamics simulations and MM/PB(GB)SA calculations, finally resulting in four compounds with higher potential for blocking quorum-sensing and 16 other promising candidates. More experimental studies are required to validate these promising molecules as CviR binders and to confirm them as antagonists.
Supplementary Materials: The following are available online, Figure S1. Known actives for CviR, Figure S2. Comparison of the Vina scores of the 10 chosen molecules from the ZINC/FDA library and the known actives, Figure S3. Comparison of the GOLD/CHEMPLP scores of the 10 chosen molecules from the ZINC/FDA library and the known actives, Figure S4. Comparison of the LeDock scores of the 10 chosen molecules from the ZINC/FDA library and the known actives, Figure S5. Comparison of the Vina scores of the 10 chosen molecules from the Chemotheca and the known actives, Figure S6. Comparison of the GOLD/CHEMPLP scores of the 10 chosen molecules from the Chemotheca and the known actives, Figure S7. Comparing the LeDock scores of the 10 chosen molecules from the Chemotheca and the known actives, Figure S8. Root-mean-square deviation plots of the Cα atoms for the CviR-ligand complexes for the selected molecules from the ZINC/FDA database, Figure S9. Root mean square deviation plots of the ligands on the CviR-ligand complexes for the selected molecules from the ZINC/FDA database, Figure S10. Root mean square deviation plots of the Cα atoms for the CviR-ligand complexes for the selected molecules from the Chemotheca database, Figure S11. Root-mean-square deviation plots of the ligands on the CviR-ligand complexes for the selected molecules from the Chemotheca database, Table S1. Redocking scores for all available structures of CviR and all molecular docking programs used, Table S2. Cross-docking Results for every structure with all programs used, Table S3. RMSD for the cross-docking procedure, Table S4. Average RMSD values (Å) of the last 40ns of the simulation for the CviR-ligand complexes for the selected molecules. Funding: This work was supported by the Applied Molecular Biosciences Unit-UCIBIO, which is financed by national funds from FCT (UIDP/04378/2020 and UIDB/04378/2020). This work was also supported by national funds from Fundação para a Ciência e a Tecnologia (grant numbers: IF/00052/2014, UID/QUI/50006/2019, and 2020.01423.CEECIND). Some of the calculations were produced with the support of INCD funded by FCT and FEDER under project 01/SAICT/2016 number 022153 and projects CPCA/A00/7140/2020 and CPCA/A00/7145/2020.

Conflicts of Interest:
The authors declare no conflict of interest.