CACHE Challenge #2: Targeting the RNA Site of the SARS-CoV-2 Helicase Nsp13

Computational workflows were diverse

The computational workflows represented diverse design strategies, techniques and tools (Figure 2). Out of 23 teams, ten used neural networks to generate or evaluate compounds, eight used crystallized fragments in the PDB to guide their design, seven used molecular dynamics simulations to account for protein flexibility, four used free energy calculation and two quantum mechanics to refine their prediction.

For example, the Poda-Hoffer team (workflow 1448 - WF1448), adopted a conservative, purely physics-based but well-established screening pipeline where Glide (Schrodinger, New York, Inc.) was used to screen a large and diverse library, with pharmacophoric constraints, against a conformational ensemble extracted from fragment-bound Nsp13 crystal structures in the PDB, along with a few conformationally refined snapshots from quick molecular dynamics simulations. The output was refined with another scoring function (HYDE, BioSolveIT) after considering crystallographic water molecules from the system. Both computational and medicinal chemists visually inspected the top-scoring molecules to finalize the selection.

The Moretti-Meiler team (WF1414) implemented the challenge on Drug-it within the Fold-it platform¹², where citizen scientists use an online gaming interface to grow fragments bound to Nsp13 available in the PDB. After multiple rounds of chemical modification, the closest commercial analogs were re-docked with RosettaLigand¹³ and ranked based on neural network-predicted absolute binding free energies¹⁴. Interestingly, these widely divergent workflows ended-up producing the two best Nsp13 binders.

Compounds were drug-like and chemically diverse

In Round 1, each team was asked to select up to 100 in-stock or make-on-demand compounds from the Enamine catalog, leading to a collection of 1957 compounds quite evenly distributed between participants (61 to 97 compounds each, Figure 3a). Participants were also encouraged to use badapple (https://datascience.unm.edu/badapple/) to filter out promiscuous compounds¹⁵, though doing so was not mandatory. Overall, compounds displayed drug-like properties, as reflected by the distribution of their Lipinski descriptors¹⁶ (Figure 3b). While three of the four fragments crystallized in the RNA site of Nsp13 included a carboxylic acid attached to a ring, compounds were diverse, as illustrated by a pairwise distance matrix of Tanimoto distances based on ECFP4 Morgan fingerprints calculated with RDKit (Figure 3c). Chemical diversity was also observed within selections from each team, with rare exceptions, outlined by darker squares along the diagonal of the distance matrix.

Only 20 compound pairs selected from different participants had a Tanimoto distance of 0.3 or lower, based on ECFP4 fingerprints. Not surprisingly, all closest analogs selected by different participants (Figure 3d) were also close analogs of the crystallized fragments found in the PDB (Figure 1), however none of these were ultimately confirmed experimentally. Indeed, in the previously reported fragment screen by crystallography, Nsp13 crystals were soaked in 50 mM fragments solutions⁵, which can lead to the capture of fragments that are too weak to be detected by SPR (maximum concentration of 200 μM). Yet, crystallographically captured fragments were successfully grown into 20–40 μM hits, as detailed below.

Experimental testing of Round 1 compounds

Helicases are complex and structurally dynamic enzymes that couple ATP (or other nucleotides) hydrolysis at one site with RNA or DNA duplex unwinding at another. Given that the fragments in the targeted Nsp13 structure (PDB codes 5RLH, 5RLZ, 5RML, and 5RMM) bound to full-length Nsp13 in the absence of ATP or RNA⁵, a similar form of the protein was used in a surface plasmon resonance (SPR) assay to measure the direct binding of the 1957 Round 1 compounds to the full-length protein (Table S2). Nsp13 is a core component of the replication-transcription complex that also includes the viral RNA-dependent RNA polymerase (RdRp)¹⁰, but the isolated protein was used in the assay for two reasons: first, fragments in the PDB were bound to the isolated monomer, and second, binding to RdRp would have obscured the results. All compounds were also tested in an ATPase assay (Table S2), but we saw no correlation between SPR and ATPase assays and decided to rely on direct binding (SPR) to advance compounds to Round 2. Indeed, false positives in the ATPase assay that may bind to other assay-specific molecular components should be true negative in SPR, while true positives binding the RNA site in the SPR assay may not inhibit the ATPase activity. We also cannot discount the possibility that SPR hits may bind at unexpected and functionally neutral sites.

All compounds were tested at 50 μM in both assays. 300 compounds had acceptable SPR sensorgram profiles with a binding signal above 50% of the expected signal (based on the amount of protein captured on the SPR streptavidin chip), and were advanced to dose-response by SPR. Another 54 compounds that inhibited the ATPase activity by 40% or more at 50 μM were selected for SPR dose-response. Dose response measurements were conducted on the resulting 354 compounds by SPR, as well as on 96 compounds in the ATPase assay. Binding to WDR5, an unrelated protein, was also measured by SPR for selected hits to flag non-specific binders. Compounds of interest with signs of poor solubility or aggregation (<80% detected laser power at 100 μM) as measured by dynamic light scattering (DLS) were also flagged but were not dismissed to avoid false negatives (Table S2). Indeed, unlike a typical drug discovery program, no active compound should be left behind in CACHE, as this would defeat the purpose of evaluating the efficiency of computational predictions. In the end, 46 compounds selected by 18 teams had a K_D < 150 μM, a binding signal between 30% and 150% of the expected signal, and were advanced to Round 2 (Figure 4, Table S2). While most hit rates were between 0 and 3%, workflows WF1454, WF1418 and WF1456 had significantly higher hit rates (9%, 8% and 7% respectively). The overall Round 1 hit rate was 2.3%.

Selection and Experimental Testing of Round 2 Compounds

The goal of the second round was to build confidence in advanced hits by experimentally verifying that their chemical analogs were also binding to the target. Compounds associated with experimental orange flags, such as signs of aggregation or poor solubility, were advanced to Round 2 to avoid false negatives and unfairly discounting computational methods. Seventeen teams selected up to 50 analogs of their Round 1 compounds of interest (compounds showing a binding signal by SPR), leading to 618 Round 2 molecules that were screened at 50 μM in an SPR binding assay, followed by dose-response and measurements of aggregation and solubility, as in Round 1 (Tables S3, S4). Compounds were also tested in an ATPase assay (Tables S5 and S6 respectively), and no correlation was observed with SPR data, as in Round 1. ¹⁹F-NMR was used as an orthogonal binding assay for fluorinated molecules.

Multiple chemical series emerged from this exercise (Table 1, Figure 5, and Supplementary material).

Table 1:

13 Compounds with a score greater than 10

Round 1 hit	Score	Most potent analog	Structure	KD(μM)
CACHE_1413_19	12.5	CACHE2-HO_1413_3		20
CACHE_1414_40	20.2	CACHE2-H0_1414_20		27
CACHE_1419_42	17.3	CACHE2-HO_1419_28		17
CACHE_1421_21	13.3	CACHE2-H0_1421_16		30
CACHE_1421_62	12.7	CACHE2-HO_1421_29		18
CACHE_1422_15	10.6	CACHE2-HO_1422_23		19
CACHE_1430_25	14.2	CACHE2-HO_1430_33		36
CACHE_1438_39	15.6	CACHE2-HO_1438_19		80
CACHE_1448_70	18.5	CACHE2-HO_1448_8		26
CACHE_1454_91	12.7	CACHE2-HO_1454_36		88
CACHE_1454_98	16.3	CACHE2-HO_1454_45		50
CACHE_1456_73	13.3	CACHE2-HO_1456_42		29
CACHE_1456_42	14.4	CACHE2-HO_1456_26		86

Open in a new tab

Thirteen of the high ranked compounds as well as compound derivatives that scored lower but exhibited high binding affinity by SPR were further assayed for inhibition of double-stranded RNA unwinding activity in a FRET-based assay. Two compounds were found to potently inhibit Nsp13 helicase activity (CACHE2-HO_1431_6: K_D 770 nM ± 180, unwinding IC₅₀ 8.6 μM ± 1.7; CACHE2-HO_1454_15: K_D 31 μM ± 0.7, unwinding IC₅₀ 57 μM ± 2). Inhibition of dsRNA unwinding by CACHE2-HO_1431_6 was also confirmed in a gel-based unwinding assay with CACHE2-HO_1454_15 partially inhibitory, consistent with CACHE2-HO_1431_6 having a more potent unwinding activity and stronger binding affinity and lower IC₅₀ value. (Figure 6). Note that many compounds in this series have an ester group linker that is likely to be hydrolyzed in cells and represents a serious medicinal chemistry liability, which penalized the final score of this chemical series. However, modifying the ester linker to a more stable group could easily address this liability while conserving potency. But this medicinal chemistry work is beyond the scope of the CACHE study.

Evaluation of experimental data and computational workflows

The biophysical data and structure-activity relationship (SAR) of Round 1 hits and their Round 2 analogs were evaluated by an independent Hit Evaluation Committee composed of industry experts in biophysics, medicinal chemistry and computational chemistry (Table S1), leading to a final score assigned to each Round 1 hit (Table S8). Overall, 13 compounds had a score greater than 10 (Table 1), reflecting robust experimental confirmation, which corresponds to a hit rate of 0.7%.

The computational workflows of CACHE #2 participants were then evaluated based on the aggregated score of Round 1 compounds, and based on the best scoring Round 1 molecule (Figure 7a,b, Table S9). In a separate evaluation scheme, all participants were asked to predict Nsp13 ligands from the merged collection of 1957 Round 1 compounds before the experimental data were generated. The aggregated score of predicted hits, normalized based on the number of hits predicted, was used to rate the computational workflows (Figure 7c). This scheme is complementary as here, all teams predicted hits from the same library, while in Rounds 1 and 2, participants screened compound collections from the Enamine catalog that may vary widely in size to best align with their computational methods and resources.

While these combined metrics provide a complete evaluation of computational workflows used in CACHE #2, a list of six well-performing workflows was compiled for further analysis, including WF1454 and WF1456, which had the best two aggregated scores, WF1414, WF1448, and WF1419 which predicted the three best scoring chemical series, and WF1438 that did best in predicting hits out of the 1957 Round 1 compounds (Figure 8). Importantly, absence from this selection focused on top-performing computational pipelines does not imply that a workflow failed.

Trends and strategies from the best performing computational workflows

Most of the best scoring compounds were docked to the RNA binding groove, at the site occupied by fragments found in PDB structures 5RMM, 5RLZ and 5RLH, the target site defined for this CACHE challenge (Figure 9). An exception is CACHE _1454_98, which is predicted to occupy an unrelated binding pocket. In this workflow (WF1454), compounds were docked onto a receptor grid spanning most of the target protein. Six of the eight Round 1 compounds from WF1454 that advanced to Round 2 occupied the RNA-binding groove, one the ATP site and one (CACHE _1454_98) an unrelated site. While top hits from other workflows occupy the expected site, they do not share pharmacophoric features or conserved interactions. CACHE_1414_40 was obtained from growing the crystallized fragment found in the PDB structure 5RMM (Figure 1) and is predicted to loosely overlap with the bound fragment.

The six best-performing workflows (Figure 8a) can be divided into three groups (Figure 10). WF1414 and WF1438 both adopted strategies where fragments from the PDB were gradually grown and commercial analogs identified along multiple iterative cycles but their implementations were drastically divergent: WF1414 relied on citizen scientists and the gaming interface provided by Foldit to grow fragments, followed by RosettaLigand¹³, a physics-based docking tool, and BCL-AffinityNet, a feed-forward deep neural network, for final scoring¹⁴; WF1438 used FEgrow¹⁷ to enumerate fragments in the binding pocket based on a hybrid machine learning (ML) / molecular mechanics energy function leveraging the ANI neural network potential¹⁸ for ligand energetics, and final evaluation with the convolutional neural network scoring function GNINA¹⁹.

Another selection strategy adopted in workflows WF1456 and WF1454 was to dock a small and diverse library with GNINA¹⁹ or Vina^20, respectively, to initiate iterative active learning cycles where a ML model is trained on a small set of docking scores to predict ML-scores for billions of commercial compounds, and where ML-scores are used to select the next small subset for docking and refinement of the ML model. In WF1454, the selection was further refined with a round of consensus scoring.

Finally, WF1448 and WF1419 implemented a more direct approach where a large and diverse library was docked with pharmacophore constraints followed by orthogonal re-scoring. WF1448 used purely physics-based approaches for docking (Glide) and Scoring (HYDE), followed by visual inspection and selection of top compounds by both computational and medicinal chemists. WF1419 used the popular open-source software Vina²⁰ for docking combined with ML/deep learning scoring functions RF-Score-VS²¹ and SCORCHs²².

Overall, five of the top six performing workflows combined physics-based and ML techniques. All five workflows used ML to score docked poses, and two (WF1454 and WF1456) used ML to accelerate screening within active learning cycles. A more conventional, purely physics-based approach (WF1448) also proved successful, demonstrating that well-established physics-based virtual screening techniques remain competitive when deployed by experienced computational chemists. While only 22% of the workflows (five out of 23) used in Round 1 explicitly accounted for protein flexibility using conformational ensembles (WF1419, WF1422, WF1447, WF1448) or flexible docking (WF1414) (Figure 2), they represented 50% of the most successful workflows (three out of six) (Figure 10). Considering the well-known conformational dynamics of helicases²³, including Nsp13⁵, accounting for receptor flexibility may indeed have increased chances of success.

Computational workflows

Computational methods are available from https://cache-challenge.org/results-cache-challenge-2

Surface Plasmon Resonance (SPR)

The assay was conducted using a Biacore^™ 8K (Cytiva) at 20 °C. Biotinylated Nsp13_SARS2, with approximately 4900–5100 response units (RU), was immobilized onto the flow cell two of a streptavidin-conjugated streptavidin chip following the manufacturer’s protocol. The flow cell one served as a reference for subtraction for each channel. Compounds were initially dissolved in 100% DMSO to create 10 mM stock solutions, which were subsequently serially diluted (factor: 0.5) to obtain six concentration points in 100% DMSO. For the SPR run, these serially titrated compound stocks were diluted at the ratio 1:50 in HBS buffer, containing Mg2+ (10 mM HEPES pH 7.4, 150 mM NaCl, 5 mM MgCl2, 0.03% (v/v) Tween 20) to achieve a final DMSO concentration 2%. Binding experiments used multi-cycle kinetics with a contact time of 60 seconds and a dissociation time of 180 seconds at a flow rate of 40 μL/min at 20 °C. The dissociation constant (KD) values were determined using steady-state affinity 1:1 binding with the Biacore^™ Insight Evaluation software (Cytiva).

Dynamic Light Scattering

The solubility of compounds was estimated by DLS that directly measures compound aggregates and laser power in solution. Compounds were serially diluted directly from DMSO stocks, then diluted 50x into filtered 10 mM Hepes pH7.4, 150 mM NaCl, 5 mM MgCl2, 0.03% Tween20 (2% DMSO final). The resulting samples were then distributed into 384-well plates (black with a clear bottom, Corning 3540), with 20 μL in each well. The sample plate was centrifuged at 3500 rpm for 5 min before loading into DynaPro DLS Plate Reader III (Wyatt Technology).

ATPase activity

The level of ATP consumed by Nsp13 was quantified by measuring the amount of remaining ATP using a luciferase-based assay as previously described³⁶. The inhibitory effects of compounds were assessed in 384-well format (14 μL final volume) using reactions composed of 50 mM HEPES, pH 7.5, 5% Glycerol, 5 mM magnesium acetate, 5 mM DTT, 0.01% Triton X-100, 0.01% BSA, 0.1 nM Nsp13, 3.5 nM 30b PolyT ssDNA, 2.5 μM ATP, and 2% DMSO. Samples containing DMSO only (no compounds) were used as a control. Reactions were started by the addition of substrate and incubated for 60 min at room temperature. Then, 10 μL of the reactions were transferred into 384-well white plates containing 10 μL luciferase reagent (Cat# V6712; Promega, Madison, WI, USA) and incubated for another 20 min at room temperature. Compounds that were followed up for dose-response experiments were tested using the same luciferase reagent, and the data were analyzed using GraphPad Prism 9.

¹⁹F-NMR Spectroscopy

The binding of fluorinated compounds was assayed by assessing the broadening and/or perturbation of ¹⁹F resonances upon addition of Nsp13 (at protein to compound ratios of 2:1 to 3:1) in PBS buffer (pH 7.4, 137 mM NaCl, 2.7 mM KCl, 10 mM Na2HPO4, 1.8 mM KH2PO4, and with 5% D2O). 1D-¹⁹F spectra were collected at 298 K on a Bruker AvanceIII spectrometer, operating at 600 MHz, and equipped with a QCI probe. Two to four thousand transients were collected with an acquisition period of 0.2 s, over a sweep width of 150 ppm, a relaxation delay of 1.5 s, and using 90° pulses centered at −120 ppm. The concentration of the compounds in both reference and protein-compound mixtures was10 μM. TFA (20 μM) was added as an internal standard for referencing. Prior to Fourier transformation, an exponential window function was applied (lb = 1 to 3) to the FID. All processing was performed at the workstation using the software Topspin 3.5.

Unwinding assays

The FRET-based dsRNA unwinding assay using recombinant nsp13 purified from E. coli and gel-based dsRNA unwinding assay using mammalian cell expressed proteins were conducted as previously described using custom dsRNA templates obtained from International DNA Technologies³⁷.