key: cord-0904653-nzm9058p
authors: Hall-Swan, Sarah; Devaurs, Didier; Rigo, Mauricio M.; Antunes, Dinler A.; Kavraki, Lydia E.; Zanatta, Geancarlo
title: DINC-COVID: A webserver for ensemble docking with flexible SARS-CoV-2 proteins
date: 2021-10-15
journal: Comput Biol Med
DOI: 10.1016/j.compbiomed.2021.104943
sha: 42e914ae87f9d787fce7c9ee56834b2e67f64c26
doc_id: 904653
cord_uid: nzm9058p

An unprecedented research effort has been undertaken in response to the ongoing COVID-19 pandemic. This has included the determination of hundreds of crystallographic structures of SARS-CoV-2 proteins, and numerous virtual screening projects searching large compound libraries for potential drug inhibitors. Unfortunately, these initiatives have had very limited success in producing effective inhibitors against SARS-CoV-2 proteins. A reason might be an often overlooked factor in these computational efforts: receptor flexibility. To address this issue we have implemented a computational tool for ensemble docking with SARS-CoV-2 proteins. We have extracted representative ensembles of protein conformations from the Protein Data Bank and from in silico molecular dynamics simulations. Twelve pre-computed ensembles of SARS-CoV-2 protein conformations have now been made available for ensemble docking via a user-friendly webserver called DINC-COVID (dinc-covid.kavrakilab.org). We have validated DINC-COVID using data on tested inhibitors of two SARS-CoV-2 proteins, obtaining good correlations between docking-derived binding energies and experimentally-determined binding affinities. Some of the best results have been obtained on a dataset of large ligands resolved via room temperature crystallography, and therefore capturing alternative receptor conformations. In addition, we have shown that the ensembles available in DINC-COVID capture different ranges of receptor flexibility, and that this diversity is useful in finding alternative binding modes of ligands. Overall, our work highlights the importance of accounting for receptor flexibility in docking studies, and provides a platform for the identification of new inhibitors against SARS-CoV-2 proteins.

[1]organization=Department of Computer Science, Rice University, city=Houston, postcode=77005, state=Texas, coun-try=United States [2] organization=MRC Institute of Genetics and Cancer, University of Edinburgh, city=Edinburgh, postcode=EH4 2XU, country=United Kingdom [3] organization=Department of Biology and Biochemistry, University of Houston, city=Houston, postcode=77005, state=Texas, country=United States [4] organization=Department of Physics, Federal University of Ceará, city=Fortaleza-CE, country=Brazil

The respiratory disease COVID-19, caused by the novel coronavirus SARS-CoV-2, went from an outbreak to a pandemic in just a few months [1] . In response, there have been unprecedented global efforts to develop vaccines and effective treatments. Several effective vaccines have been validated and used in massive vaccination campaigns [2] , particularly in developed countries. However, the need for pharmacological treatments for infected patients persists due to unequal vaccination coverage across the globe [3, 4] and to the rise of more virulent variants that can cause symptoms even in fully-vaccinated individuals [5, 6] .

Among potential pharmacological targets, proteins involved in the viral replication have been used in several computational studies focused on drug design, drug repurposing and virtual screening [7, 8] . An impressive number of SARS-CoV-2 protein structures have been solved in a short period of time, but computational studies targeting these proteins have been mostly limited to the use of a single experimental structure [9] [10] [11] [12] [13] [14] [15] [16] . Although this approach is expected to correctly reproduce similar experimental structural data, it tends to fail in exploring the binding of ligands with diverse chemical features (e.g., large ligands or ligands with alterna- 1 These authors contributed equally to this work tive binding modes), as it does not account for the inherent flexibility of proteins in solution [17] [18] [19] .

One example of the importance of protein malleability on SARS-CoV-2 drug screening involves the Main protease (Mpro). Although Mpro has been, so far, one of the most explored SARS-CoV-2 targets in computational studies, there remain numerous open questions for the design of effective inhibitors. Indeed, the malleability of the Mpro active-site cavity remains the greatest challenge in the development of effective inhibitors. By collecting X-ray data from ligands bound to Mpro at room temperature, Kneller et al. have demonstrated experimentally that Mpro has the ability to substantially distort its shape in response to binding [18] . Such malleability allows Mpro to accommodate a larger diversity of physico-chemical features (e.g., chemical groups, size, charge distribution) in ligands.

Accounting for protein flexibility adds complexity and increases the computational cost of molecular docking studies. Therefore, explicitly sampling full receptor flexibility while also sampling alternative ligand conformations is generally unfeasible. Various strategies have been proposed to make this problem computationally tractable [20] . A popular solution involves explicitly sampling only selected parts of the receptor, as in selective docking with flexible bindingsite residues. Another one, known as ensemble docking, only implicitly accounts for full receptor flexibility [21] : a separate sampling method is used to explore protein flexibility, and a set of representative conformations is extracted to create an ensemble for docking [22] ; ligand sampling is then conducted as in a regular molecular docking job. By docking the ligand against all conformations in the ensemble, rather than a single conformation, ensemble docking implicitly accounts for receptor flexibility [23, 24] .

Note that each of these strategies suffers from its own limitations and that successful predictions usually require significant knowledge about the system of interest and the required methodologies. On one hand, selective docking requires knowledge of which residues should be prioritized for explicit sampling during docking. In addition, execution time grows exponentially with the number of flexible bonds, which might prevent the exploration of complex binding sites or large ligands. On the other hand, ensemble docking requires expertise on tools that can sample protein conformations (e.g., molecular dynamics, Langevin dynamics, Monte Carlo simulations, or coarse-grained simulations [25, 26] ). These methods can produce a huge number of conformations, in turn requiring user expertise on methods that can select representative conformations (e.g., dimensionality reduction, clustering, or free energy estimation) [22, 27] . All these tasks are very time-consuming, taking days to weeks to implement and execute, depending on available computational resources [22, 23] .

Despite these challenges, recent studies have highlighted the importance of considering multiple receptor conformations in molecular docking studies [22, 24, 28] , and even presented preliminary results on the use of ensemble docking to screen for drug inhibitors against SARS-CoV-2 proteins [23, 29] . However, the resources and expertise required for conducting an ensemble docking study still prevents its use by most researchers interested in testing new drugs or natural products against SARS-CoV-2 proteins.

To bridge this gap, we present DINC-COVID, a userfriendly webserver for ensemble docking of small molecules and peptides to SARS-CoV-2 proteins. The sampling of binding modes is performed with DINC, a parallelized metadocking approach that has been shown to outperform conventional docking tools on several challenging datasets, especially for large or flexible ligands (e.g., peptides or peptidomimetics) [30] . The proteins currently available for docking through DINC-COVID are the Main protease (Mpro) [9] [10] [11] [12] [13] [14] [15] [16] , the Papain-like protease (PLpro) [31] , and the RNAdependent RNA polymerase (RdRp) [32] [33] [34] . For Mpro, in addition to the catalytic site, we allow for predictions targeting an allosteric site [35] . Note that the ensembles we have pre-computed cover different ranges of protein flexibility, for each targeted binding site. This initial selection of SARS-CoV-2 proteins and binding sites was performed based on the availability of crystallographic structures in the Protein Data Bank (PDB) and the potential importance for drug discovery. Nonetheless, our team is continuously creating additional ensembles, covering other targets of interest, which will be later added to the DINC-COVID webserver.

We have carried out a proof-of-concept validation of DINC-COVID using datasets of experimentally-determined inhibitors for SARS-CoV-2 proteins. As such data are still scarce, we have focused our validation on a small dataset of PLpro inhibitors [31] and three small datasets of Mpro inhibitors [18, 36, 37] . For all datasets, we have obtained good correlations between predicted binding energies and experimentally-determined 50 values. We have also compared DINC-COVID against two other online servers targeting SARS-CoV-2 proteins, namely the COVID-19 Docking Server [38] , and the DockThor-VS webserver [39] . Note that DINC-COVID is the only webserver that performs the simultaneous docking of a ligand against multiple receptor conformations and automatically aggregates the results (as suggested by the term ensemble docking), while the two other webservers provide more docking targets. We provide a detailed discussion on how our results are affected by choices related to method settings, dataset composition and statistical analysis. Finally, we provide further insight on the range of receptor conformational flexibility captured by DINC-COVID predictions, depending on the data/methodology used to generate the input ensemble and on the scoring function used to rank the output binding modes.

To enable fast predictions of protein-ligand binding modes while accounting for receptor flexibility, we have decoupled the steps required for ensemble generation from those directly related to the ensemble docking procedure. For each binding target of interest, three ensembles have been precomputed and stored on our webserver (see Figure 1 ). Therefore, all ensembles mentioned in this manuscript are available for docking. After accessing the DINC-COVID webserver, users can choose a binding target, select one of the available ensembles, and upload a ligand of interest (e.g., a drug-like ligand or a peptide). These files are then used as input for ensemble docking, which involves the parallelized meta-docking approach DINC [30] . All generated binding modes are rescored and ranked using three scoring functions: AutoDock Vina [40] (a.k.a. Vina), AutoDock 4 [41] (a.k.a., AD4) and Vinardo [42] . Finally, the top-scoring binding modes (for each scoring function) are returned to the user, reflecting the flexibility of both the ligand and the receptor.

Four target binding sites of SARS-CoV-2 proteins are currently available through the DINC-COVID webserver: catalytic site and allosteric site of Main protease (Mpro), catalytic site of Papain-like protease (PLpro), and catalytic site of RNA-dependent RNA polymerase (RdRp) (see Figure 2 ). For each target binding site, we first compiled three ensembles of conformations (see Appendix A for more information): one ensemble containing crystal structures, and two ensembles containing conformations extracted from in silico molecular dynamics (MD) simulations.

We performed the MD simulations with the GROMACS 2019 package [43] , using two different force fields: CHARMM36 [44] and GROMOS53a6 [45] . With each force field we ran five independent MD simulations of 200 ns, for a total of 1 s. All simulations were performed with full-atom representation and explicit water molecules. Note that CHARMM36 simulations used TIP3 water models, while GROMOS53a6 used SPC water models. In addition, ions Na+ and Cl-were added at concentration of 0.15 M to neutralize the net charge of the system. Each system was minimized through the steepest descent algorithm, followed by two equilibration steps of 1 ns each, using NVT and NPT ensembles, respectively. Positions of protein atoms were restrained during equilibration. During the simulations, long-range electrostatics were modelled with the Particle Mesh Ewald (PME) method [46] , temperature coupling was set at 310.15 K using the V-rescale thermostat [47] and the Parrinello-Rahman barostat [48] with a reference pressure of 1 bar; a compressibility of 4.5 × 10 −5 bar −1 was applied for pressure control. Covalent bonds were constrained to their equilibrium length by the LINCS algorithm [49] . The integration steps of all simulations were set to 2 fs. The stability of each structure during the simulation time was assessed through root-mean-square deviation (RMSD) calculations using the "gmx rms" algorithm. Eventually, we extracted a set of roughly 100,000 conformations from each MD simulation with MDtraj [50] . These come in addition to the 179 crystallographic structures (156 for Mpro, 12 for PLpro, and 11 for RdRp) we extracted from the Protein Data Bank (PDB).

To create the final ensembles we implemented a datadriven protocol that could extract representative conformations using algorithms from the scikit-learn package [51] . Briefly, minimal distances between all pairs of amino acids in the targeted binding site (i.e., distances between the closest pair of heavy atoms belonging to two residues) were used as features for a Principal Component Analysis (PCA). The first principal components accounting for around 80% of variance in the data were used to project each set of conformations into a lower-dimensional space. In this space, the Elbow method of scikit-learn was used to determine the ideal number of clusters (K). Representative members of each cluster were identified with K-means, and extracted to build the final ensembles. Although other strategies could have been adopted, this is a reasonable protocol to build a representative ensemble of conformations, which has produced good results in our validation experiments (see Results section). Finally, all selected structures were protonated at pH 7.0 using the PROPKA algorithm [52] at the PDB2PQR server [53] . In summary, for each target binding site, we produced 

DINC-COVID uses parallelization to speed up both the sampling and the scoring of protein-ligand binding modes (see Figure 3 ). Once a ligand structure and an ensemble of receptor conformations have been selected as input, the algorithm triggers multiple parallel docking jobs with the fast docking method Vina. Each independent job starts with a distinct, randomized conformation of the ligand and one of the receptor conformations from the selected ensemble. This first batch of docking jobs produces a diverse set of binding modes for the ligand against a single receptor conformation. This process is then repeated for each of the receptor conformations available in the ensemble. Once all docking jobs are completed, the rescoring phase starts. In this phase, all binding modes predicted for each receptor conformation are rescored using three popular scoring functions: Vina, Vinardo and AD4. At the end of the whole procedure, a userspecified number of top scoring conformations is provided as output, for each scoring function. These results include both alternative conformations of the ligand and alternative conformations of the receptor.

The DINC-COVID ensemble docking procedure is made available through an intuitive webserver interface (see Figure 4) . It allows users to upload a ligand of interest and select the targeted receptor ensemble. Optional ligand preparation (e.g., adding hydrogens and charges) is also available. Note that we plan to allow for batch submission in future work. After submission and execution, users can visualize the selected binding modes through the viewer embedded in the "Results" page. If they want to perform offline analyses, users can download all the results, which include top-scoring binding modes in pdb format, as well as a text file listing the selected conformations and their binding energies with respect to all scoring functions. More details are available at dinc-covid.kavrakilab.org/method/.

The DINC-COVID webserver is implemented using Docker [54] , allowing it to run with its own dependencies on any machine with Docker installed. The main backend is implemented with Django [55] . Submitted jobs are managed by Celery [56], a distributed task queue. The webserver is currently hosted on a 16-cores virtual machine in the Owl Research Infrastructure Open Nebula (ORION) VM Pool on Rice University Campus. DINC-COVID currently uses 16 parallel threads, and our virtual infrastructure allows for future expansion.

As in any virtual screening study, identifying better inhibitors for SARS-CoV-2 proteins requires ranking sets of ligands with respect to their estimated binding affinity with a given receptor. To validate the ability of our method to accurately rank ligands, we computed correlations between docking-derived binding energies and the natural logarithm of experimentally-determined half-maximal inhibitory concentration, ( 50 ); these correlations were evaluated using Pearson's and Spearman's . As experimental data, we used ( 50 ) instead of 50 because it is expected to better correlate with predicted binding energies [57] [58] [59] [60] . As DINC-COVID-predicted data, for each ligand and each scoring function, we used the binding energy predicted for the top-scoring binding mode of this ligand according to this scoring function. Note that our goal was to assess the overall ranking accuracy of our ensemble docking results, as opposed to assessing the accuracy of individual binding energy predictions.

To compare DINC-COVID with other methods, we performed the same experiments with the COVID-19 Docking Server (hereafter referred to as COVID-19-DS) [38] and the DockThor-VS webserver [39, 61] . For COVID-19-DS we selected Mpro or PLpro as the "nCoV Protein Target", with the recommended exhaustiveness of 8. Although the top 10 models are returned with a score value (in kcal/mol), we used only the score value of the top 1 model to perform our correlation analysis.

For DockThor-VS, we selected either the Nsp5-Main protease (wild type) or the Nsp3-PLpro (wild type) as the target protein. Since two crystal structures of the Mpro receptor (with PDB codes 6LU7 and 6W63) are made available, we used both of them as input receptors for independent docking runs. Similarly, for PLpro, the structures with PDB codes 6W9C and 6WX4 were used as input receptors for independent docking runs. As docking parameters, we manually defined the catalytic binding site of both proteins and used the standard algorithm precision (i.e., 1,000,000 evaluations, population size of 750, and 24 runs). DockThor-VS provides a single output binding mode with a corresponding binding energy, which we used to compute correlations with experimental data.

We used four datasets of ligands with experimentallydetermined binding affinities in our validation experiments: one for the PLpro receptor and three for the Mpro receptor. More specifically, to create the first dataset (Table 1) , we selected seven PLpro inhibitors with known 50 values published by Osipiuk et al. [31] : The first Mpro dataset (Table 2 ) comprises 14 compounds identified by Li et al. through free energy perturbation-based virtual screening and validated experimentally through an enzymatic assay [36] . In their study, 50 values for all compounds were measured through a fluorescence assay using an inhibitory curve with 500 nM of enzyme, 20 M of substrate and six different concentrations of each compound. The second Mpro dataset (Table 3 ) was initially gathered by Ngo et al. to validate a free energy perturbation protocol to assess Mpro inhibitors during virtual screening [37] . It comprises 11 Mpro inhibitors whose experimental 50 values were retrieved from the literature and therefore obtained in different experiments. The third Mpro dataset (Table 4) was obtained from a study demonstrating the capacity of leupeptin and three hepatitis C clinical protease inhibitors to bind and inhibit SARS-CoV-2 Mpro [18] . Kneller et al. characterized these four ligands by X-ray crystallography at near-physiological room temperature, thus capturing Mpro motions that would not be observable at lower temperature. 

The accuracy of molecular docking tools is usually assessed through self-docking experiments, where protein-ligand complexes previously determined by experimental methods are used as references to be reproduced. In a self-docking experiment, the ligand conformation is randomized before docking, and the quality of predicted binding modes is quantified by computing their deviation from the reference experimental structure. These experiments are useful to validate the geometries of predicted binding modes, which is usually summarized in terms of RMSD values. However, it is important to highlight that reproducing the bound geometries of available crystal structures is not the intended application of our webserver. DINC-COVID's goal is to efficiently account for receptor flexibility during the docking of a ligand. This allows for the sampling of alternative low-energy binding modes that are not captured by rigid crystal structures or by conventional molecular docking techniques. Therefore, by design, the geometries of output complex conformations can significantly differ from those observed in available crystal structures, especially when using an ensemble of receptor conformations derived from molecular dynamics. Our underlying hypothesis is that these alternative lowerenergy binding modes can better approximate the population of receptor-bound ligand conformations that exist in solution (in vitro/in vivo).

To evaluate DINC-COVID predictions, we used a recently published dataset of seven PLpro inhibitors (see Table 1 ) [31] . This dataset is mostly composed of small ligands (i.e., up to seven flexible bonds) that are all strong binders (i.e., 50 values ranging from 2.3 M to 43.2 M). It is therefore quite challenging with respect to the task of reproducing a ranking of binding affinities.

In spite of that, we achieved high correlations between ( 50 ) values and binding energies predicted by DINC-COVID using the Crystal ensemble (i.e., Pearson's = 0.9 and Spearman's = 0.89), far above correlations achieved by binding energies predicted by other ensembles and by other servers (see Table 1 ). In addition, when comparing predicted binding modes with crystal structures available for three of these compounds, it appears that DINC-COVID produced binding modes that are very similar to the corresponding crystal structures (see Figure 5 ).

The best results were obtained with the Vina scoring function and the Crystal ensemble, which might be due to the nature of this dataset (i.e., small ligands and limited receptor flexibility). However, good correlations were also observed in other settings, for example using Vinardo or the Charmm ensemble (see Table 1 ). Note that the highest correlation with ( 50 ) values is not a spurious one (see Appendix B, Fig. 11 ). Although the small size of this dataset is likely to have contributed to this high correlation, the structural agreement between predicted binding modes and ex-perimental structures corroborates the validity of this correlation. In addition, the structural examination of binding modes provides a rationale for the lower correlations observed for other servers in this experiment.

To further validate DINC-COVID predictions, we used three datasets of ligands that were tested in inhibitory assays with Mpro, and for which 50 values are available.

The first dataset contains 14 drug-like ligands that are experimentally characterized SARS-CoV-2 Mpro inhibitors (see Table 2 ) [36] . Performing a docking experiment with this dataset produced rather low correlations with experimental ( 50 ) values for all tested settings. According to Pearson's , the best correlation is observed for DockThor-VS using the 6LU7 receptor ( = 0.36), but according to Spearman's , the best correlation is observed for DINC-COVID using Vinardo and the Charmm ensemble ( = 0.55). This experiment highlights the challenges of dealing with outliers (see Appendix B, Fig. 12 ). Indeed, two ligands have a particularly strong impact on the obtained correlations (see Table 2 ). Removing Indinavir from the correlation analysis produces better correlations for all settings, the best results being achieved by DockThor-VS with the 6LU7 receptor (Pearson's = 0.65) and DINC-COVID with Vinardo and the Charmm ensemble (Spearman's = 0.75). Although removing Dipyridamole alone from the correlation analysis does not produce substantial changes, removing both Dipyridamole and Indinavir substantially improves Pearson's correlations, with the best results being achieved by DockThor-VS with the 6LU7 receptor ( = 0.74) and DINC-COVID with Vinardo and the Charmm ensemble ( = 0.65). These results show that, when considering only the 12 remaining ligands, good agreement is reached between predicted binding energies and experimental binding affinities.

The second one is a mixed dataset composed of 11 ligands of various types that are experimentally characterized SARS-CoV-2 Mpro inhibitors (see Table 3 ) [37] . This dataset is quite challenging for the task of reproducing a ranking of binding affinities because all ligands are very strong binders, with 50 values ranging from 0.04 M to 21.39 M. After performing a docking experiment with this dataset, reasonably good correlations between predicted binding energies and ( 50 ) values were obtained across most settings (see Table 3 ). The best correlations were achieved by COVID-19-DS (with Pearson's = 0.74 and Spearman's = 0.73). The second best correlations were achieved by DINC-COVID using the Crystal ensemble, with Vinardo according to Pearson's (= 0.7), but with Vina according to Spearman's (= 0.68). Results obtained on this dataset are also strongly influenced by two outliers: 13a and Shikonin (see Appendix B, Fig. 13 ). Removing these outliers from the correlation analysis leads to improved correlations for all methods. The best correlations are then achieved by DINC-COVID using Vina and the Gromos ensemble (Pearson's = 0.88 and Spear- Table 1 Correlations between docking-derived binding energies and experimentally-determined binding affinities of PLpro drug-like inhibitors [31] . Predicted binding energies are reported for DINC-COVID (using three ensembles and three scoring functions), COVID-19-DS and DockThor-VS (using two receptor conformations). Correlation coefficients between binding energies and ( 50 ) values are reported below; in each row, the two highest coefficients are highlighted. man's = 0.95). This represents a very good agreement between predicted binding energies and binding affinities. The third dataset contains four peptidomimetics that are known antiviral compounds: Leupeptin, Telaprevir, Narlaprevir and Boceprevir (see Table 4 ) [18] . At the time the Mpro Crystal ensemble for DINC-COVID was built, no experimental structure of Mpro bound to these compounds was available. Indeed, it was only recently that the crystallographic structures of these complexes were solved at room temperature (and deposited under PDB codes 6XCH, 6XQS, 6XQT, and 6XQU). Working at room temperature allowed for the induced-fit process to occur without the constraints of crystallographic packing. Even without access to these additional structures, DINC-COVID was able to reproduce the experimentally-determined binding affinities with good correlation (i.e., Pearson's = 0.86), using Vina and the Charmm ensemble. Unfortunately, due to the very small number of ligands, this correlation is not of high quality (see Appendix B, Fig. 14) . However, other correlations (of higher quality) achieved by DINC-COVID are superior to those achieved by DockThor-VS and COVID-19-DS (see Table 4 ). Interestingly, results produced by the Vinardo scoring present the worst correlations for all three DINC-COVID ensembles in this experiment. This is in contrast with good results produced by Vinardo in previous experiments, and might indicate a limitation of this function when ranking peptidomimetic ligands.

These three datasets contain Mpro inhibitors with experimentally determined 50 values. However, these values were produced by different groups using different protocols, and might not be fully comparable. That is why we initially decided to use these datasets separately. On the other hand, working with small datasets introduces other issues regarding the reliability of the obtained ligand rankings. Therefore, we performed an additional correlation analysis in which the three Mpro datasets were combined into a larger one. Although the observed correlations are not very high, they are positive for all methods (see Table 5 ). The best correlations were achieved by DINC-COVID using Vinardo and the Crystal ensemble (with Pearson's = 0.41 and Spearman's = 0.42). Again, Indinavir was clearly an outlier (see Appendix B, Fig. 15 ), and removing it produced better correlations. Without it, the best results are achieved by DINC-COVID using Vinardo and the Crystal ensemble according to Pearson's (= 0.49), or by DockThor-VS using the 6LU7 receptor according to Spearman's (= 0.52). Finally, using a larger dataset allowed us to compute another statistic for all methods: the area under the receiver operating characteristic (ROC) curve (see Appendix B, Fig. 16 ). This analysis confirmed the good performance of DINC-COVID when using Vinardo and the Crystal ensemble, but also indicated a competitive performance of DockThor-VS when using the 6LU7 receptor.

Taken together, these results constitute a proof-of-concept validation of the predictive capabilities of DINC-COVID, as it consistently produced predictions in good agreement with experimental data across all studied datasets. The results also provide insights on the differences between docking outputs produced by different ensembles, as well as a comparison with webservers that do not perform ensemble docking. For instance, when making predictions for ligands that do not require significant adjustments in the receptor conformation, good results can be obtained with an ensemble of crystal structures (i.e., with only limited receptor flexibility) and even without an ensemble docking protocol, as illustrated by results obtained with the other webservers (see Tables 2  and 3 ). On the other hand, MD-derived conformations have greater chances of successfully accommodating ligands that require greater receptor flexibility for binding (see Table 4 ).

Overall, our results highlight the benefits of using ensemble docking as an efficient way to account for receptor flexibility in molecular docking targeting SARS-CoV-2 proteins. In DINC-COVID, receptor flexibility is not adjusted on-the-fly during docking, as it is accounted for implicitly through the use of pre-computed ensembles [20] . However, we tried to compensate for this limitation by generating different ensembles for each target binding site. Having several pre-computed ensembles at their disposal allows users to explore different ranges of receptor flexibility, according to their needs.

The differences in range between the provided ensembles are real and can be quantified (see Figure 6 ). For instance, the Crystal ensemble created for Mpro captures mostly sidechain rearrangements in the receptor, and only subtle variations of the backbone: the all-heavy-atom RMSD (for binding site residues) between members of the Crystal ensemble is lower than 0.5 Å. On the other hand, the Gromos ensemble created for Mpro captures much greater conformational changes, with up to 5 Å all-heavy-atom RMSD in relation to the structures in the Crystal ensemble. Interestingly, the Charmm ensemble created for Mpro includes conformations in an intermediate range of flexibility.

Most importantly, these differences between ensembles are useful, since the best results are not always provided by the same ensemble. As an extension, the diversity within ensembles is also useful, since the best ligand binding modes are not always obtained with the same receptor conformation from a given ensemble (see Figure 7 ). In summary, our analysis shows that both levels of diversity -i.e., within an ensemble, and between ensembles -contribute to the good quality of DINC-COVID predictions.

A third level of diversity in DINC-COVID output is introduced by the scoring function used to rank all sampled Table 2 Correlations between docking-derived binding energies and experimentally-determined binding affinities of Mpro drug-like inhibitors [36] . Predicted binding energies are reported for DINC-COVID (using three ensembles and three scoring functions), COVID-19-DS and DockThor-VS (using two receptor conformations). Correlation coefficients between binding energies and ( 50 ) values are reported below; in each row, the two highest coefficients are highlighted. Results for correlation analyses without outliers (i.e., Dipyridamole and Indinavir) are also reported. Figure 7) . Some of the starkest differences are observed when using the Charmm ensemble: the greatest diversity in receptor conformations included in the top conformations is observed with the AD4 scoring function, while Vina and Vinardo produce lower levels of diversity. It is worth emphasizing that in its current implementation DINC-COVID is not automatically selecting one of the available scoring functions, and is not providing a consensus ranking across all three functions. Instead, each output produced by the server includes the three alternative rankings of all sampled conformations, according to the three available scoring functions. This approach provides users with the flexibility to either follow one of the rankings or to analyze the results further, considering expert knowledge about the scoring functions or the ligand/receptor of interest.

DINC-COVID offers a ready-to-use solution for researchers to account for protein flexibility while testing compounds against SARS-CoV-2 proteins. It allows users to run ensemble docking experiments without the additional burden of time-consuming simulations required for ensemble generation and docking preparation. Since sampling conformational flexibility for large proteins and multimeric complexes is a challenging and computationally intensive task, we decoupled it from the docking procedure itself. By making the pre-computed ensembles readily available in our webserver, we aim to facilitate and speed up the use of ensemble docking by a broader audience of users searching for new SARS-CoV-2 inhibitors.

Users of our webserver can choose between four target binding sites: the catalytic site of Main protease (Mpro), Papain-like protease (PLpro) and RNA-dependent RNA polymerase (RdRp), as well as the allosteric binding site of Mpro. For each target binding site, three distinct ensembles are provided, including experimental data (i.e., structures from X- Table 3 Correlations between binding energies and experimental binding affinities of Mpro inhibitors [37] . Binding energies are reported for DINC-COVID (using three ensembles and three scoring functions), COVID-19-DS and DockThor-VS (using two receptor conformations). Correlation coefficients between binding energies and ( 50 ) values are reported below; in each row, the two highest coefficients are highlighted. Results for correlation analyses without outliers (i.e., 13a and Shikonin) are also reported. ray crystallography) or simulated data (i.e., structures from molecular dynamics with different force fields). The choice of the force field has significant impacts on the results of MD simulations, since it encompasses a series of empiricallydetermined parameters [62, 63] . By relying on two force fields that use distinct representations for hydrogen atoms and different water models, we aim to limit the impact of force field bias on conformational sampling. It is important to recognize that the conformational ensembles available on the DINC-COVID webserver reflect different scales of protein flexibility, from subtle side-chain rearrangements (e.g., in the Crystal ensemble) to larger backbone motions (e.g., in the Gromos ensemble). Note that we have carefully prepared the structures provided in the ensembles, which includes verifying the protonation state of His, Glu and Asp residues. We have also provided users with the possibility of automatically preprocessing the ligand before docking. Altogether, these features make our server of potential interest even for advanced users, adding robustness to docking analyses with SARS-CoV-2 proteins.

The main contribution of our webserver is to account for receptor flexibility while sampling ligand binding modes. This is made possible through a parallelized ensemble dock-ing protocol using DINC. DINC has been extensively validated in previous work, through both self-docking and crossdocking experiments [30, 64, 65] . On the other hand, the DINC-COVID webserver cannot be properly assessed through conventional self-docking experiments, since the goal of ensemble docking is precisely to find alternative low-energy binding modes. For instance, when using an MD-derived ensemble, the best ranked receptor conformations will usually not be the ones with the lowest RMSD to a reference crystal structure. Our approach relies on the fact that a set of diverse receptor conformations can accommodate a much broader range of ligand conformations and sizes than a single receptor conformation [22] . To sum up, the goal of ensemble docking is not to reproduce binding geometries observed in crystallographic structures, but to produce conformations that better reflect the possible binding modes of a proteinligand complex in solution.

In this context, we searched for experimental binding affinities of inhibitors for SARS-CoV-2 proteins, to evaluate whether DINC-COVID could properly rank these ligands. Good correlations between predicted binding energies and experimental binding affinities were obtained with both PLpro inhibithors (e.g., Pearson's = 0.9, see Table 1 ) and Table 4 Correlations between docking-derived binding energies and experimentally-determined binding affinities of Mpro peptidomimetic inhibitors [18] . Predicted binding energies are reported for DINC-COVID (using three ensembles and three scoring functions), COVID-19-DS and DockThor-VS (using two receptor conformations). Pearson's correlation coefficient between binding energies and ( 50 ) Mpro inhibithors. Interestingly, we obtained good correlations between predicted binding energies and experimental binding affinities with a recently published set of antiviral peptidomimetic inhibitors for Mpro (e.g., Pearson's = 0.86, see Table 4 ). As noted, the binding of these peptidomimetics requires malleability of Mpro catalytic site [18] , and as such could not be predicted by docking methods using a single receptor conformation (i.e., COVID-19-DS and DockThor-VS). These results highlight the potential of DINC-COVID to identify completely novel binding modes. Such binding modes could involve conformational changes of the receptor's binding site, and could potentially be predicted regardless of the availability of experimental data on such receptor's conformation. As a note of caution, we should reiterate that the analyzed datasets were small, which can introduce biases regarding their composition (e.g., impact of outliers) and the choice of statistical analysis. We tried to address these limitations by discussing the impact of different scenarios in the analysis of these datasets. In the case of the Mpro datasets, we also performed an alternative analysis in which we aggregated all results into a single dataset. Despite all the observed variability across different scenarios, the ensemble docking strategy of DINC-COVID produced results that were consistently in agreement with experimental data, outperforming other webservers in many of these experiments.

It is important to acknowledge that accurate scoring remains an open challenge in the field of molecular docking [66] , especially when considering large and flexible ligands [30] . In the context of ensemble docking, accurately scoring and ranking sampled binding modes becomes even more challenging. To try and compensate for the limitations of relying on a single scoring function, our meta-docking approach provides scoring of the results with three alternative functions. Users can then decide if they want to rely on a specific scoring function, or build a consensus across the three functions. In fact, our results show that the scoring function can affect the diversity of binding modes produced by DINC-COVID, which is a factor that users might want to explore in different ways depending on their specific application. Most importantly, despite the challenges and known limitations of available scoring functions, we were able to obtain good correlation with experimental data in our validation experiments.

There are countless scoring functions proposed in previous studies and implemented in other tools. DINC-COVID predictions could certainly be improved by using scoring functions that might be tailored to a particular type of ligand. In order to explore this possibility, users can request to download a larger number of output conformations, and rescore them locally with the most appropriate scoring function. Future work on DINC-COVID will include the implementation of additional scoring functions during the sampling phase, and consensus ranking during the rescoring phase.

The quality of results produced by ensemble docking is directly influenced by the quality of the conformations included in the ensemble [22] . For instance, it was shown that including too many conformations in an ensemble can deteriorate docking results [22, 67] , since it exacerbates existing limitations of both sampling and scoring. As briefly discussed in the Methods section, we used a reasonable and efficient strategy to extract conformations from MD simulations. By combining standard procedures for dimensionality reduction and clustering, we tried to maximize coverage of the conformational landscape, while minimizing the number of conformations included in each ensemble. In fact, we used the same strategy to generate the Crystal ensem-bles, as opposed to simply using all available crystal structures. Again, using a different strategy for selecting conformations could potentially improve DINC-COVID results [29, 68] . However, demonstrating the optimality of our selection strategy goes beyond the scope of this study. It is worth noting that our proof-of-concept validation indicates the success of decisions made for conformation generation and conformation selection, providing us with a useful protocol that can be further improved upon.

Finally, the choice of making DINC-COVID available through a webserver also comes with trade-offs. On one hand, it takes away most of the burden associated with software installation and the preparation of input files. We hope that this format will enable users who are not familiar with molecular docking to go ahead and run their own tests, potentially evaluating new compounds that are not yet available in public databases. On the other hand, it constrains the use of our ensemble docking strategy to the specific algorithms and parameters we considered. Fortunately, our implementation offers the possibility to make this pipeline available for local customization and execution. In particular, we are already working on leveraging the use of Docker containerization and widely-used python packages, to develop a future version of DINC-COVID that could be deployed on local computational resources. This will enable a full customization of DINC-COVID, and foster its use for large-scale virtual screening of viral variants. 

The file ligands-mol2.zip containing the ligands that have been evaluated in this study is available at http://dinc-covid.kavrakilab.org/static/pdb/datasets/ligands-mol2.zip. All ligands have been prepared for docking; they are organized based on the datasets/tables presented here.

This work was funded in part by the National Science Allosteric binding site: The selected Mpro crystal structures for the allosteric site are shown in Table 7 and Figure 8 .B. The number of structures in the Crystal, Charmm, and Gromos ensembles are 24, 10 and 9, respectively. The scoring box center was set to -39.23, 7.85 and 57.56 (X, Y and Z, respectively) with dimensions 82×78×82 Å for the Crystal ensemble; -36.23, 7.09 and 85.56 (X, Y and Z, respectively) with dimensions 82×76×82 Å for the Charmm ensemble; and -37.23, 6.09 and 58.56 (X, Y and Z, respectively) with dimensions 84×76×82 Å for the Gromos ensemble. Again, the box used for sampling was tuned to the binding site of each receptor conformation.

Three ensembles were built for PLpro, based on conformational changes in the catalytic binding site. To build the Crystal ensemble, monomeric forms of PLpro were extracted from 12 crystals representing wild type structures (see Table 8 and Figure 8 .C). Crystals with double occupancy were split into two files to ensure each file would represent one distinct orientation. In addition, two new ensembles were produced by sampling conformations from MD simulations of PLpro in the apo state. For these simula- As for Mpro, we have performed an RMSD analysis (for binding site residues only) to quantify differences between the three receptor ensembles of PLpro (see Figure 9 ). Again, 

Three ensembles were built for RdRp, based on conformational changes in the catalytic binding site. The position of remdesivir in the crystallographic structure with PDB code 7BV2 was taken as a reference for the targeted site. In total, 7 crystallographic structures, out of 11 analyzed, were included into the RdRp Crystal ensemble (see Table 9 and Figure 8 .D). In addition, two new ensembles were produced by sampling conformations from MD simulations of RdRp in the apo state. In total 100,000 snapshots were extracted for each force field. The number of structures in the Crystal, Charmm, and Gromos ensembles are 7, 14 and 13, respectively. The scoring box center was set to 90.02, 90. 17 We have also performed an RMSD analysis (for binding site residues only) to quantify differences between the three receptor ensembles of RdRp (see Figure 10 ). The allheavy-atom RMSD between members of the Crystal ensemble is lower than 0.89 Å. The all-heavy-atom RMSD between members of the Crystal ensemble and of the Charmm ensemble ranges from 1.36 to 2.43 Å. The all-heavy-atom RMSD between members of the Crystal ensemble and of the Gromos ensemble ranges from 1.89 to 3.58 Å.

In addition to reporting correlation coefficients in the main text, in this appendix we present scatter plots illustrating these correlations. This allowed us to determine which correlations were of high quality and to identify outliers in the datasets. This appendix includes the correlation analysis of the PLpro dataset (see Figure 11 ), of the three Mpro datasets (see Figures 12, 13 and 14) , and of the large Mpro dataset obtained by combining these three small datasets together (see Figure 15 ).

In addition, we present a receiver operating characteristic (ROC) analysis performed on the large, combined Mpro dataset (see Figure 16 ). This dataset is the only one that is large enough to allow for this analysis, although a typical ROC analysis would include many more ligands. ROC analyses are traditionally performed to evaluate the ability of a docking method at distinguishing binding from non-binding ligands. However, in our case, all ligands are Mpro binders. Therefore, we performed this analysis by arbitrarily dividing the dataset between "strong binders" and "weak binders" using an IC 50 threshold of 10 M. Figure 11 : Correlation analysis for seven PLpro drug-like inhibitors [31] . Experimentally-determined binding affinities are plotted against docking-derived binding energies obtained with DINC-COVID (using three ensembles and three scoring functions), COVID-19-DS and DockThor-VS (using two receptor conformations). The best linear fit and Pearson's are also reported on each plot. [36] . Experimentally-determined binding affinities are plotted against docking-derived binding energies obtained with DINC-COVID (using three ensembles and three scoring functions), COVID-19-DS and DockThor-VS (using two receptor conformations). The best linear fit and Pearson's are also reported on each plot. Two outliers are highlighted: Dipyridamole (in blue) and Indinavir (in red). [37] . Experimentally-determined binding affinities are plotted against docking-derived binding energies obtained with DINC-COVID (using three ensembles and three scoring functions), COVID-19-DS and DockThor-VS (using two receptor conformations). The best linear fit and Pearson's are also reported on each plot. Two outliers are highlighted: 13a (in blue) and Shikonin (in red). [18] . Experimentally-determined binding affinities are plotted against docking-derived binding energies obtained with DINC-COVID (using three ensembles and three scoring functions), COVID-19-DS and DockThor-VS (using two receptor conformations). The best linear fit and Pearson's are also reported on each plot. 

SARS-CoV-2/COVID-19 and advances in developing potential therapeutics and vaccines to counter this emerging pandemic

Landscape and progress of global COVID-19 vaccine development

World Bank grants for global vaccination -why so slow?

The unequal scramble for coronavirus vaccines -by the numbers

New infections by SARS-CoV-2 variants of concern after natural infections and post-vaccination in Rio de Janeiro, Brazil

Transmission of SARS-CoV-2 variant B.1.1.7 among vaccinated health care workers

A critical overview of computational approaches employed for COVID-19 drug discovery

Computational drug discovery and repurposing for the treatment of COVID-19: A systematic review

Putative inhibitors of SARS-CoV-2 main protease from a library of marine natural products: a virtual screening and molecular modeling study

Drug repurposing studies targeting SARS-CoV-2: an ensemble docking approach on drug target 3C-like protease (3CLpro)

Potential Inhibitors for Novel Coronavirus Protease Identified by Virtual Screening of 606 Million Compounds

Rapid Identification of Potential Inhibitors of SARS-CoV-2 Main Protease by Deep Docking of 1.3 Billion Compounds

Characterization and Computational Study on Potential Inhibitory Action of Novel Azo Imidazole Derivatives against COVID-19 Main Protease (Mpro: 6LU7)

Statins and the COVID-19 main protease: in silico evidence on direct interaction

Interactions Between Remdesivir, Ribavirin, Favipiravir, Galidesivir, Hydroxychloroquine and Chloroquine with Fragment Molecular of the COVID-19 Main Protease with Inhibitor N3 Complex (PDB ID:6LU7) Using Molecular Docking

Fragment Molecular Orbital Based Interaction Analyses on COVID-19 Main Protease -Inhibitor N3 Complex (PDB ID: 6LU7)

Structural and Evolutionary Analysis Indicate That the SARS-CoV-2 Mpro Is a Challenging Target for Small-Molecule Inhibitor Design

Malleability of the SARS-CoV-2 3CL Mpro Active-Site Cavity Facilitates Binding of Clinical Antivirals

Flexibility and mobility of SARS-CoV-2-related protein structures

Understanding the challenges of protein flexibility in drug design

Ensemble Docking in Drug Discovery

Ensemble Docking in Drug Discovery: How Many Protein Configurations from Molecular Dynamics Simulations are Needed To Reproduce Known Ligand Binding?

Supercomputer-Based Ensemble Docking Drug Discovery Pipeline with Application to Covid-19

Drug design and repurposing with DockThor-VS web server focusing on SARS-CoV-2 therapeutic targets and their non-synonym variants

Molecular dynamics, monte carlo simulations, and langevin dynamics: a computational review

Coarse-grained conformational sampling of protein structure improves the fit to experimental hydrogenexchange data

Extracting representative structures from protein conformational ensembles

Multi-target, ensemble-based virtual screening yields novel allosteric KRAS inhibitors at high success rate

Targeting multiple conformations of SARS-CoV2 Papain-Like Protease for drug repositioning: An in-silico study

Using parallelized incremental meta-docking can solve the conformational sampling issue when docking large ligands to proteins

Structure of papain-like protease from SARS-CoV-2 and its complexes with non-covalent inhibitors

RNA-Dependent RNA Polymerase as a target for COVID-19 drug discovery

Analysis of SARS-CoV-2 RNA-dependent RNA polymerase as a potential therapeutic drug target using a computational approach

Druggable targets of SARS-CoV-2 and treatment opportunities for COVID-19

Molecular characterization of ebselen binding activity to SARS-CoV-2 main protease

Identify potent SARS-CoV-2 main protease inhibitors via accelerated free energy perturbation-based virtual screening of existing drugs

Assessing potential inhibitors of SARS-CoV-2 main protease from available drugs using free energy perturbation simulations

COVID-19 Docking Server: a meta server for docking small molecules, peptides and antibodies against potential targets of COVID-19

Drug design and repurposing with DockThor-VS web server focusing on SARS-CoV-2 therapeutic targets and their non-synonym variants

AutoDock Vina: improving the speed and accuracy of docking with a new scoring function, efficient optimization, and multithreading

A semiempirical free energy force field with charge-based desolvation

Vinardo: A Scoring Function Based on Autodock Vina Improves Scoring, Docking, and Virtual Screening

GROMACS: fast, flexible, and free

Optimization of the additive CHARMM all-atom protein force field targeting improved sampling of the backbone , and side-chain (1) and (2) dihedral angles

A biomolecular force field based on the free enthalpy of hydration and solvation: the GROMOS force-field parameter sets 53A5 and 53A6

A smooth particle mesh Ewald method

Canonical sampling through velocity rescaling

Polymorphic transitions in single crystals: A new molecular dynamics method

Lincs: A linear constraint solver for molecular simulations

MDTraj: A Modern Open Library for the Analysis of Molecular Dynamics Trajectories

Duchesnay, Scikit-learn: Machine learning in Python

PROPKA3: Consistent Treatment of Internal and Surface Residues in Empirical pKa Predictions

PDB2PQR: an automated pipeline for the setup of Poisson-Boltzmann electrostatics calculations

Empowering App Development for Developers | Docker

The Web framework for perfectionists with deadlines | Django, ???? [56] Celery -Distributed Task Queue âĂŤ Celery 5.0.5 documentation

Crystal structure of SARS-CoV-2 main protease provides a basis for design of improved α-ketoamide inhibitors

Structure of Mpro from SARS-CoV-2 and discovery of its inhibitors

Structural basis for the inhibition of SARS-CoV-2 main protease by antineoplastic drug carmofur

Structure-based design of antiviral drug candidates targeting the SARS-CoV-2 main protease

New machine learning and physics-based scoring functions for drug discovery

Current status of protein force fields for molecular dynamics simulations

All-Hydrocarbon Staples and Their Effect over Peptide Conformation under Different Force Fields

General Prediction of Peptide-MHC Binding Modes Using Incremental Docking: A Proof of Concept

HLA-Arena: A Customizable Environment for the Structural Modeling and Analysis of Peptide-HLA Complexes for Cancer Immunotherapy

True accuracy of fast scoring functions to predict high-throughput screening data from docking poses: The simpler the better

Recipes for the selection of experimental protein conformations for virtual screening

A selective method for optimizing ensemble docking-based experiments on an InhA Fully-Flexible receptor model

We thank the Center for Research Computing (CRC) at Rice University for supporting our use of ORION VM Pool. Use of CRC resources is supported by the Data Analysis and Visualization Cyberinfrastructure funded by NSF (OCI-0959097) and by Rice University. We also thank the Centro Nacional de SupercomputaÃğÃčo (CESUP/UFRGS, Brazil), whose resources were used to perform our MD simulations.

The authors declare no conflict of interest. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript, or in the decision to publish the results. Neither the manuscript nor any parts of its content are currently under consideration or published in another journal.

To build the Crystal ensembles, the dimeric state of 156 crystallographic structures of Mpro was considered, based on conformational changes in the catalytic or allosteric binding pocket. For the MD simulations, the dimeric form of Mpro was built from the structure with PDB code 6LU7.Catalytic binding site: The selected Mpro crystal structures for the catalytic site are shown in Table 6