AI-guided pipeline for protein-protein interaction drug discovery identifies a SARS-CoV-2 inhibitor

Abstract

Protein-protein interactions (PPIs) offer great opportunities to expand the druggable proteome and therapeutically tackle various diseases, but remain challenging targets for drug discovery. Here, we provide a comprehensive pipeline that combines experimental and computational tools to identify and validate PPI targets and perform early-stage drug discovery. We have developed a machine learning approach that prioritizes interactions by analyzing quantitative data from binary PPI assays or AlphaFold-Multimer predictions. Using the quantitative assay LuTHy together with our machine learning algorithm, we identified high-confidence interactions among SARS-CoV-2 proteins for which we predicted three-dimensional structures using AlphaFold-Multimer. We employed VirtualFlow to target the contact interface of the NSP10-NSP16 SARS-CoV-2 methyltransferase complex by ultra-large virtual drug screening. Thereby, we identified a compound that binds to NSP10 and inhibits its interaction with NSP16, while also disrupting the methyltransferase activity of the complex, and SARS-CoV-2 replication. Overall, this pipeline will help to prioritize PPI targets to accelerate the discovery of early-stage drug candidates targeting protein complexes and pathways.

Keywords: AlphaFold; Machine Learning; Protein–Protein Interactions; SARS-CoV-2; VirtualFlow.

Figure 1. Developing a maSVM algorithm to classify protein pairs from hsPRS-v2 and hsRRS-v2 using the LuTHy assay.

(A) Schematic overview of the maSVM learning algorithm. Step 1: assembly of reference set; Step 2: feature selection and data normalization for training and test set; Step 3: assembly of ‘e’ training sets (ensembles) by unweighted sampling ‘j’ protein pairs from the reference set to train ‘e’ maSVM models, where the training classifier labels are reclassified in ‘i’ iterations; Step 4: prediction of test set protein pairs excluding training set pairs using the respective maSVM model. Scatter plot showing (B) log-transformed and normalized in-cell mCitrine expression (x axis) against normalized cBRET ratios (y axis limited to ‘>-10’) or (C) the number of proteins pairs (x axis) against log-transformed and normalized cLuC ratios (y axis) for all hsPRS-v2 (blue) and hsRRS-v2 (magenta) protein pairs from all eight tagging configurations. Average classifier probabilities from the 50 maSVM models are displayed as the size of the data points and as a colored grid in the background. Scatter plot showing (D) normalized cBRET ratios (x axis) or (E) normalized cLuC ratios (x axis) against classifier probability (y axis) for all hsPRS-v2 (blue) and hsRRS-v2 (magenta) protein pairs from all eight tagging configurations. Receiver characteristic analysis comparing sensitivity and specificity between (F) cBRET ratios or (G) cLuC ratios and classifier probabilities. The calculated areas under the curve are displayed. Bar plots showing the fraction of hsPRS-v2 and hsRRS-v2 protein pairs that scored above classifier probabilities of 50%, 75%, or 95% with (H) LuTHy-BRET or (I) LuTHy-LuC. Only the highest classifier probability per tested tagging configuration is considered. (J) Heatmaps showing the highest classifier probabilities for the hsPRS-v2 (top) and hsRRS-v2 (bottom) protein pairs per tested tagging configuration. Due to different reference set interactions, heatmaps for SIMPL data from Yao et al (Yao et al, 2020) are shown in Appendix Fig. S3 and Appendix Fig. S4. (K) Bar plots showing the fraction of hsPRS-v2 and hsRRS-v2 protein pairs that scored above classifier probabilities of 50%, 75%, or 95% or above-fixed cutoffs at maximum specificity, mean or median plus one standard deviation for seven binary PPI assays. Only the highest classifier probability per tested tagging configuration is considered. All LuTHy experiments from this study were repeated twice with n = 2, biological replicates, each containing n = 3 technical replicates; all other from Choi et al (Choi et al, 2019) and Yao et al (Yao et al, 2020). Bars and error bars in this figure represent mean values and standard error of the proportion, respectively. Source data are available online for this figure.

Figure 2. Benchmarking AFM using well-established positive and random reference sets.

(A) Schematic overview of AlphaFold-multimer (AFM) benchmarking. First, the hsPRS-v2 and hsRRS-v2 were filtered for protein pairs with less than 1400 amino acids combined, resulting in 51 positive (hsPRS-AF) and 67 random reference set pairs (hsRRS-AF). For these 118 protein pairs, five structural models were predicted each using ColabFold through the AFM algorithm (590 total structures). Following, PAE and pLDDT values were extracted from the AFM-predicted structures, and inter-subunit amino acids were filtered for pLDDT >50. If >10 inter-subunit amino acids remained, PAE values were k-means clustered. If clustering failed, the mean PAE of the unclustered amino acids was calculated, else the average PAE for each of the eight clusters were calculated and the cluster with the lowest average PAE was selected as the amino acid region with the minimal distance between the two proteins. In addition, PDBePISA was used to determine the solvation-free energy (ΔG) and the area (iA) of the interface region for 521 of the 590 structures. For the remaining 69 structures PDBePISA could not identify an interface. Finally, a multi-adaptive maSVM learning algorithm was trained on the PAE and iA features of the hsPRS-AF and hsRRS-AF as outlined in Fig. 1A. (B) Heatmap of the PAEs, ΔGs, iAs, and predicted probabilities for protein pairs of the hsPRS-AF and hsRRS-AF. Shown are the minimum PAE values after k-means clustering. If <10 amino acids had a pLDDT >50, the PAE values were not used and are shown in black. Protein pairs where no interaction interface was detected by PDBePISA are shown in gray. (C) Representative example for the k-means clustering strategy of AFM reported PAE values. Heatmap shows the PAEs for the protein pair BAD + BCL2L1 (hsPRS-AF) rank 1 model. The intra-molecular PAEs are shown with 50% opacity. The predicted local distance difference test (pLDDT) for all five predicted models (rank 1–5) are shown as line graphs on top and on the right of the heatmap. Inter-molecular PAE regions with pLDDT >50 that were used for k-means clustering are highlighted with arrows. (D) Clustering results of regions highlighted in (C). Cluster numbers are indicated. (E) Average PAE values for the eight clusters from (C,D). The arrow indicates the cluster with the lowest average PAE value. (F) Scatter plot showing inter-PAE (x axis) against interface area (y axis) for all models of the hsPRS-AF (blue) and hsRRS-AF (magenta) protein pairs. Average classifier probability from the 100 maSVM models are displayed as the size of the data points and as a colored grid in the background. (G) Scatter plots showing PAE (x axis, left panel) or interface area (x axis, right panel) against classifier probability (y axis) for all hsPRS-AF (blue) and hsRRS-AF (magenta) protein pairs. (H) Bar plots showing the fraction of hsPRS-AF and hsRRS-AF protein pairs that scored above classifier probabilities of 50%, 75% and 95%. Bars and error bars represent mean values and standard error of the proportion, respectively, with n = 5 structural models predicted. Note that this analysis includes interactions in the hsPRS-AF that have experimentally solved structures. A comparison between interactions with and without structural information can be found in Fig. EV2B. Source data are available online for this figure.

Figure 3. Validating the maSVM algorithm by mapping interactions within multiprotein complexes using the LuTHy and mN2H assays.

(A) Structures of the protein complexes analyzed in this study: LAMTOR (PDB: 6EHR), BRISC (PDB: 6H3C) and MIS12 (PDB: 5LSK). (B) Binary interaction approach to systematically map PPIs within distinct complexes. Every protein subunit from each complex was screened against every other one (all-by-all, 16 × 16 matrix). (C–E) Scatter plot showing (C) log-transformed and normalized in-cell mCitrine expression (x axis) against normalized cBRET ratios (y axis), (D) number of protein pairs (x axis) against log-transformed and normalized cLuC ratios (y axis) or (E) the number of protein pairs (x axis) against the log-transformed and normalized mN2H ratios (y axis) for all protein pairs of the LAMTOR (yellow), BRISC (blue) and MIS12 (green) complexes and inter-complex (magenta) protein pairs from all eight tagging configurations. Average classifier probabilities from the 50 maSVM models for LuTHy-BRET (C) and LuTHy-LuC (D) or 100 maSVM models for mN2H (E) are displayed as the size of the data points and as a colored grid in the background. (F–H) Scatter plot showing on the x axis the normalized (F) cBRET ratios, (G) cLuC ratios, or (H) mN2H ratios against classifier probabilities (y axis) for protein pairs of the BRISC (blue), LAMTOR (yellow) and MIS12 (green) complexes and inter-complex (magenta) protein pairs from all eight tagging configurations. (I–K) Bar plots showing the fraction of protein pairs of the LAMTOR (yellow), BRISC (blue) and MIS12 (green) complexes and inter-complex protein pairs that scored above the classifier probabilities of 50%, 75% or 95% by (I) LuTHy-BRET, (J) LuTHy-LuC, and (K) mN2H. Only the highest classifier probability per tested tagging configuration is considered. (L) Tile plots showing the classifier probabilities for the Donor/F1 protein pairs (x axis) against the Acceptor/F2 protein pairs (y axis) for LuTHy-BRET (orange, left), LuTHy-LuC (purple, middle) and mN2H (green, right) for protein pairs above 75% or 95%. Only the highest classifier probability per tested tagging configuration is shown. LuTHy experiments were performed in HEK293 cells two times with n = 2 biological replicates, each containing n = 3 technical replicates. mN2H experiments were performed in HEK293T cells four times with n = 4 biological replicates and n = 1 technical replicate. Tiles of not expressed constructs are filled black and respective protein names are colored in red. Bars and error bars in this figure represent mean values and standard error of the proportion, respectively. Source data are available online for this figure.

Figure 4. Mapping binary interactions between SARS-CoV-2 proteins.

(A) Search space between SARS-CoV-2 proteins tested by LuTHy. (B,E) Scatter plots showing (B) log-transformed and normalized in-cell mCitrine expression (x axis) against normalized cBRET ratios (y axis) or (E) number of protein pairs (x axis) against log-transformed and normalized cLuC ratios (y axis) for SARS-CoV-2 (orange) protein pairs from all eight tagging configurations. Average classifier probability from the 50 maSVM models is displayed as the size of the data points and as a colored grid in the background. (C, F) Bar plots showing the fraction of hsPRS-v2 and hsRRS-v2 protein pairs that scored above classifier probabilities of 50%, 75%, or 95% by (C) LuTHy-BRET or (F) LuTHy-LuC. Only the highest classifier probability per tested tagging configuration is considered. (D, G) Tile plots showing SARS-CoV-2 protein pairs with >95% and >75% classifier probability detected with (D) LuTHy-BRET and (G) LuTHy-LuC. Only the highest classifier probability per tested tagging configuration is shown. All LuTHy experiments were performed in HEK293 cells two times, with n = 2 biological replicates, each containing n = 3 technical replicates. (H) Scatter plot showing PAE (x axis) against interface area (y axis) for 23 SARS-CoV-2 protein pair structures predicted with AlphaFold-Multimer, with n = 5 models each. Average classifier probability from the 100 maSVM models trained on the hsPRS-AF and hsRRS-AF (Fig. 2F) is displayed as the size of the data points and as a colored grid in the background. (I) Bar plots showing the number of AlphaFold-Multimer predicted SARS-CoV-2 protein pair structures that scored above classifier probabilities of 50%, 75%, and 95%. (J) Heatmaps showing the classifier probabilities for the AFM-predicted SARS-CoV-2 protein pair structures. Source data are available online for this figure.

Figure 5. Predicting the NSP10-NSP16 PPI complex with AFM to target the interaction interface by ultra-large virtual drug screening.

(A) Heatmap showing the predicted alignment error (PAE) of the AlphaFold-Multimer predicted NSP10-NSP16 complex for the rank 1 model. The intra-molecular PAEs are shown with 50% opacity. The predicted local distance difference test (pLDDT) for all five predicted models (rank 1–5) are shown as line graphs on top and on the right of the heatmap. (B) The five models of the AlphaFold-Multimer predicted NSP10-NSP16 complex and the published crystal structure (PDB: 6W4H) are shown. Structures were overlaid using the “matchmaker” tool of ChimeraX. (C) Scatter plot showing for each amino acid (x axis) the solvation-free energy (ΔG, y axis, fill color) upon formation of the interface, in kcal/mol, as determined by PDBePISA. Dots represent the mean ΔG for the five predicted models and error bars correspond to the standard deviation from n = 5 AFM-predicted structural models. The x axis indicates the amino acid positions of the whole complex structure starting from NSP10’s N-terminus and ending with NSP16’s C-terminus. Lysine 93 (Lys93) of NSP10 and aspartate 106 (Asp106) of NSP16, which showed the strongest solvation-free energy gain upon complex formation, are indicated, respectively. (D) Zoom-in into the NSP10-NSP16 complex showing the contacts of NSP10’s Lys93 and NSP16’s Asp106 as determined using ChimeraX, using the Contacts tool with the parameters, “VDW overlap ≥ -0.40 Å”, “Limited by selection: with at least one selected” of NSP10 Lys93 and NSP16 Asp106; “Include intramodel”; “Display as pseudobonds”. (E) LuTHy-BRET donor saturation assay, where constant amounts of NSP10-NL WT or Lys93Glu are co-expressed with increasing amounts of mCitrine-NSP16 WT or mCitrine-NSP16 Asp106Lys. Nonlinear regression was fitted through the data using the “One-Site – Total” equation of GraphPad Prism. Data points represent mean values from two n = 2 (NSP10 + NSP16, NSP10 Lys93Glu + NSP16) or n = 4 (NSP10 + NSP16 Asp106Lys) biological replicates each containing n = 2 technical duplicates. (F) Docking box on the NSP10 structure (PDB: 6W4H) used for the ultra-large virtual screen. (G) Schematic overview of the workflow of the virtual docking screen using VirtualFlow. (H) Docking scores of the top 100 molecules identified by virtual screening. The horizontal line indicates mean docking score and error bars the standard deviation, with the virtual screen performed once (n = 1). Source data are available online for this figure.

Figure 6. Compound 459 inhibits the NSP10-NSP16 interaction and reduces SARS-CoV-2 replication.

(A) Schematic overview of the NSP10-NSP16 methyltransferase (MTase) assay. (B) Heatmap showing the result of the MTase activity of the NSP10-NSP16 complex in the absence or presence of 100 µM of the top 15 compounds. Statistical significance was calculated with a kruskal-wallis test (P value = 9.7e-5, chi-squared = 47.656, df = 17, the experiment was performed once with n = 3, technical replicates), followed by a post hoc Dunn test and adjusted p-values are shown. (C) Compound 459 docked onto the NSP10 structure (PDB: 6W4H). (D) Chemical structure of compound 459. (E) Assay principle of the microscale thermophoresis (MST) assay. The fluorescence intensity change of the labeled molecule (purple) after temperature change induced by an infrared laser (red) is measured. The binding of a non-fluorescent molecule (blue) can influence the movement of the labeled molecule. (F) Representative MST traces of labeled NSP10 and different concentrations of unlabeled compound 459. The bound fraction is calculated from the ratio between the fluorescence after heating (F₁) and before heating (F₀). (G) Scatter plot showing the 459 concentration (x axis) against the fraction of 459 bound to NSP10 (y axis). Nonlinear regression was fitted through the data using the “One-Site –  Total” equation of GraphPad Prism (the experiment was repeated three times with n = 3, biological replicates). (H) Scatter plot showing the 459 concentration (x axis) against the normalized BRET ratio (nBRET ratio) for the interaction between NSP10-NL and mCit-NSP16 measured in HEK293 cells. Nonlinear regression was fitted through the data using the “log(inhibitor) vs. response–Variable slope (four parameters)” equation of GraphPad Prism (the experiment was repeated four times with n = 4, biological replicates, each containing n = 3 technical replicates). (I) Scatter plot showing the 459 concentration (x axis) against the relative luminescence measured from icSARS-CoV-2-nanoluciferase in HEK293-ACE2 cells. Nonlinear regression was fitted through the data using the “log(inhibitor) vs. normalized response” equation of GraphPad Prism (the experiment was repeated three times, with n = 3 for 0.1, 100, 200 µM; n = 6 for 0.2, 0.4 µM; all other n = 9; all biological replicates; error bars represent the standard deviation). (J) Barplot showing the relative luminescence measured from icSARS-CoV-2-nanoluciferase in HEK293-ACE2 cells upon incubation with 0, 25, 50, or 100 µM of compound 459 together with 2.5 µM AZ1 or without AZ1 (0.0 µM). Statistical significance was calculated in GraphPad Prism by a “two-way ANOVA”, where each cell mean was compared to the other cell mean in that row using “Bonferroni’s multiple comparisons test” (the experiment was repeated three times, with n = 3, biological replicates; error pars represent the standard deviation; source of variation: 57.91% 459 concentration, P < 0.0001; 28.33% AZ1 concentration, P < 0.0001; 11.40% 459/AZ1 interaction, P < 0.0001).

Figure EV1. (related to Fig. 1). Effect of different scoring approaches on recovery rates.

(A) Schematic overview of the LuTHy-BRET and LuTHy-LuC assays. X: Protein X, Y: Protein Y, D: NanoLuc donor, A: mCitrine acceptor, AB: antibody. (B) With the LuTHy assay, each protein pair X–Y can be tested in eight possible configurations (N- vs. C-terminal fusion for each protein), and proteins can be swapped from one tag to the other resulting in 16 quantitative scores for each protein pair, i.e., eight for LuTHy-BRET and eight for LuTHy-LuC. (C) Line plots showing the fraction of protein pairs that scored positive (y axis) dependent on the quantitative interaction scores (x axis) for 10 binary PPI assay versions. For each tested protein pair, the tagging configuration with the highest interaction score is used. For LuTHy all eight tagging configurations were tested, whereas for MN2H, VN2H, YN2H, GPCA, NanoBi four and for KISS, MAPPIT and SIMPL two tagging configurations were tested. Recovery rates at maximum specificity, i.e., where none of the protein pairs in the RRS scored positive (0%), are indicated. Note that in Choi et al (Choi et al, 2019) recovery rates at maximum specificity were calculated by using distinct cutoffs for each tagging configuration. (D) Line plots showing the fraction of protein pairs that scored positive (y axis) dependent on the distribution of interaction scores, i.e., the mean of all interaction scores + n*(sd) (x axis) for 10 binary PPI assay versions. Recovery rates at mean + 1 standard deviation are indicated. (E) Line plots showing the fraction of protein pairs that scored positive (y axis) dependent on the distribution of interaction scores, i.e., the median of all interaction scores + n*(sd) (x axis) for 10 binary PPI assays. Recovery rates at median + 1 standard deviation are indicated. LuTHy experiments from this study were repeated twice with n = 2, biological replicates, each containing n = 3 technical replicates; SIMPL from Yao et al (Yao et al, 2020); all other from Choi et al (Choi et al, 2019). Note that the SIMPL assay was benchmarked by Yao et al (Yao et al, 2020) against 88 positive proteins pairs derived from the hsPRS-v1 (Venkatesan et al, 2009) and as a random reference set against “88 protein pairs with baits and preys selected from the PRS but used in combinations determined computationally to have low probability of interaction” (Yao et al, 2020).

Figure EV2. (related to Fig. 2). Training a maSVM algorithm to classify AFM-predicted structures.

(A) Receiver characteristic analysis comparing sensitivity and specificity between the five AFM-predicted structural models for PAE, ΔG and iA of the hsPRS-AF and hsRRS-AF. (B) Bar plots showing the fraction of hsPRS-AF and hsRRS-AF interactions with structures deposited in PDB that scored above classifier probabilities of 50%, 75% and 95% by AlphaFold-Multimer (i) by LuTHy (ii) or the mean recovery of N2H (MN2H, VN2H, YN2H), GPCA, KISS, MAPPIT and NanoBiT (iii). Data for the SIMPL assay was excluded for this analysis due to the different composition of the reference sets. LuTHy experiments from this study were repeated two times with n = 2, biological replicates, each containing n = 3 technical replicates; AFM was used to predict n = 5 structural models; all other from Choi et al (Choi et al, 2019). Bars and error bars in this figure represent mean values and standard error of the proportion, respectively.

Figure EV3. (related to Figs. 4 and 5). Validating SARS-CoV-2 protein interactions using the mN2H assay and predicting SARS-CoV-2 protein complexes structures using AlphaFold-Multimer.

(A) Venn diagrams showing the overlap between interactions recovered by LuTHy at >50%, >75% and >95% probabilities and interactions deposited in the IntAct database (Orchard et al, 2014). (B) Scatter plot showing normalized mN2H ratios (y axis) of each of the eight SARS-CoV-2 interactions newly identified with LuTHy (x axis). Average classifier probabilities obtained from the hsPRS-v2/hsRRS-v2 mN2H models are displayed as the size of the data points and as a colored grid in the background. (C) Scatter plot showing normalized mN2H ratios (x axis) against classifier probabilities (y axis) for the newly identified SARS-CoV-2 interactions selected for validation. (D) Bar plots showing the fraction (left y axis) and number (right y axis) of newly identified SARS-CoV-2 interactions selected for validation that scored above classifier probabilities of 50%, 75% or 95% with mN2H. Bars and error bars represent mean values and standard error of the proportion, respectively, with n = 3 biological replicates. (E) Heatmaps showing the mN2H classifier probabilities for the newly identified SARS-CoV-2 interactions selected for validation. (F) Boxplots showing predicted alignment error (PAE), solvation-free energy (ΔG) and interface area (iA) from AlphaFold-Multimer (AFM) predicted SARS-CoV-2-AF structures. Boxplots display the median, lower and upper hinges of the 25th and 75th percentiles and lower and upper whiskers extending from the hinges with 1.5× the interquartile range. Each dot represents one predicted structural model. (G) Scatter plot showing PAE (x axis) against interface area (y axis) for all SARS-CoV-2-AF (orange) protein pairs. Average classifier probability predicted by the 100 maSVM models trained by the hsPRS-AF and hsRRS-AF set (see Fig. 2F), is displayed as the size of the data points. Each point in the colored grid in the background displays the average classifier probabilities from the 100 maSVM models. (H) Scatter plot showing the ΔG (x axis) for all five AFM-predicted structural SARS-CoV-2-AF models against the LuTHy-BRET determined binding strengths (BRET₅₀, see Appendix Fig. S7A,B). The respective log-transformed interface areas are indicated by the fill color of the data points. A linear regression fit through the data is shown and the Spearman correlation coefficient (R) and P value are indicated. (I) Barplot showing the fraction of AFM-predicted structures with 0–75%, 75–95% and >95% classification probability that have an experimentally reported structure deposited to the PDB (Berman et al, 2000) database. (J,K) Luminescence (J) and fluorescence (K) values from LuTHy-BRET donor saturation experiments, where constant amounts of NSP10-NL WT or K93E (Lys93Glu) are co-expressed with increasing amounts of mCitrine-NSP16 WT or D106K (Asp106Lys). Experiments with NSP10-NL WT and K93E were repeated two times, with n = 2, biological replicates, and each with n = 2 technical replicates; experiments with mCit-NSP16 WT and D106K were repeated four times, with n = 4, biological replicates, and each with two technical replicates, n = 2. Bars and error bars represent the mean and standard deviation, respectively.