High-throughput matrix screening identifies synergistic and antagonistic antimalarial drug combinations

Drug resistance in Plasmodium parasites is a constant threat. Novel therapeutics, especially new drug combinations, must be identified at a faster rate. In response to the urgent need for new antimalarial drug combinations we screened a large collection of approved and investigational drugs, tested 13,910 drug pairs, and identified many promising antimalarial drug combinations. The activity of known antimalarial drug regimens was confirmed and a myriad of new classes of positively interacting drug pairings were discovered. Network and clustering analyses reinforced established mechanistic relationships for known drug combinations and identified several novel mechanistic hypotheses. From eleven screens comprising >4,600 combinations per parasite strain (including duplicates) we further investigated interactions between approved antimalarials, calcium homeostasis modulators, and inhibitors of phosphatidylinositide 3-kinases (PI3K) and the mammalian target of rapamycin (mTOR). These studies highlight important targets and pathways and provide promising leads for clinically actionable antimalarial therapy.

. Single agent data for compounds screened against 3d7, Hb3 and Dd2 P. falciparum isolates. Table S2. Summary of matrix combination screening data against P. falciparum Table S3. Complete analysis of each combination screened against P. falciparum Table S6. Combinations for minimum spanning tree analysis Table S14. P. falciparum combination high-throughput screen data quality control filtered Table S15. P. falciparum self-cross interactions Table S16. P. falciparum 3D7 cHTS data quality control filtered Table S17. P. falciparum Dd2 cHTS data quality control filtered Table S18. P. falciparum HB3 cHTS data quality control filtered Table S19. P. falciparum 3D7 interactions with a DBSum negative value less than 3 Table S20. P. falciparum Dd2 interactions with a DBSum negative value less than 3 Table S21. P. falciparum HB3 interactions with a DBSum negative value less than 3 Table S22. P. falciparum 3D7 interactions with a DBSum positive value greater than 3 Table S23. P. falciparum Dd2 interactions with a DBSum positive value greater than 3 Table S24. P. falciparum HB3 interactions with a DBSum positive value greater than 3 Table S25. Aggregated P. falciparum interactions with artemether with a DBSum negative value less than 3 Table S26. Aggregated P. falciparum interactions with artesunate with a DBSum negative value less than 3 Table S27. Aggregated P. falciparum interactions with dihydroartemisinin with a DBSum negative value less than 3 Table S28. Aggregated P. falciparum interactions with lumefantrine with a DBSum negative value less than 3 Table S29. Aggregated P. falciparum interactions with halofantrine with a DBSum negative value less than 3 Table S30. Aggregated P. falciparum interactions with mefloquine with a DBSum negative value less than 3 Table S31. Aggregated P. falciparum interactions with chloroquine with a DBSum negative value less than 3 Table S32. Aggregated P. falciparum interactions with atovaquone with a DBSum negative value less than 3 Table S33. Aggregated P. falciparum interactions with pyronaridine with a DBSum negative value less than 3 Table S34. Aggregated P. falciparum interactions with piperaquine with a DBSum negative value less than 3 Table S35. Aggregated P. falciparum interactions with tafenoquine with a DBSum negative value less than 3 Table S36. Aggregated P. falciparum interactions with amodiaquine with a DBSum negative value less than 3 Table S37. Aggregated P. falciparum interactions with sulfadoxine with a DBSum negative value less than 3 Table S38. Aggregated P. falciparum interactions with NITD609 with a DBSum negative value less than 3 Table S39. Aggregated P. falciparum interactions with NITD609 with a DBSum negative value less than 3 QC criteria The quality control score is a numerical characterization of the quality of a combination that attempts to take into account single agent performance and the presence of noise in the dose combination region of the matrix. It is composed of a number of heuristics, developed by examination of a series of matrix screening runs. The calculation is performed on the bounded (0-100) response matrix and considers the following conditions: 1. DMSO response lies between 80 and 100.
2. Both single agents should display a valid dose response and provide a valid IC 50 value. 3. If single agents have a valid curve fit, the responses (across the whole concentration range as well when the top concentration is excluded) should have a relative standard deviation greater than 20 (35) 4. The relative standard deviation of the dose combination sub-matrix should be greater than 25 5. The responses in the dose combination sub-matrix should exhibit spatial autocorrelation (i.e., combination responses should not be randomly distributed within the dose matrix). This is measured by Morans I (39). We consider a combination to exhibit spatial autocorrelation if the p-value is less than 0.05 Each combination is tested against the five conditions and the result is expressed as a vector, defined as = { }, = 1 … 5 where = 1 if the 'th condition is true, 0 otherwise.
The final score is obtained by computing the dot product of the binary vector and a weight vector, of the same length as . The latter allows us to assign an importance to each of the five conditions noted above. The final score is given by Currently, the weight vector is defined as = {2,5,3,5,3}, allowing for a minimum score of 0 (the combination passes all QC criterion) and maximum score of 18 (combination fails all QC criterion).
A combination with a score of 0 is deemed as good quality and suitable for further consideration. Increasing values of the QCScore indicate poorer quality combination responses ( Fig S22).
It should be noted that the score can identify false positivecombinations that are scored well, but on visual inspection turn out to be poor quality. The score is meant to provide an initial ranking and should not preclude manual inspection. Each matrix block was ultimately manually classified as synergistic, additive, antagonistic, or inconclusive.
Additional matrix metrics In addition to the previously described metrics (15) we expanded our analysis to include two new quantifiers of synergy -DBSumNeg and DBSumPos. Considering all dose combinations tested, these are defined as the sum of positive deviations from the Bliss model and the sum of the negative deviations from the Bliss model, respectively. In contrast to simply summing all deviations from the Bliss model, these two variables characterize the extent of synergy and antagonism, respectively, within a set of dose combinations.
Expanded isobologram analyses Many studies evaluating drug synergies rely upon static values such as isobolographic analyses or combination indices (CI) derived from the multiple drug effect equation developed by . The CI method is based upon discreet concentration responses and reports generally rely upon using data thought to reflect the overall combination response of the two agents. We chose not to list CI values within our web-based ranking options as many combinations revealed overlapping combination effects (synergy and additivityorsynergy and antagonism). Although isobolograms act an additional means to visualize synergy between two agents we were not able to automate the drawing of isobolograms for each data set. Further, our dosing was reliant on standard dilutions from a starting value and do not necessarily provide large numbers of fixed ratio dosing. Isobolograms for selected combinations were generated using Compusyn (Combosyn, Inc, version 1.0) and are presented in Fig 2A (panel 4). The 10×10 data sets used to derive them are shown in Fig S9.  D-F (for ATM + reserpine). From these, we derived Fa-CI plots that demonstrate these combinations are synergistic at relevant doses (Fa>0.5 or Fa>0.75) at all ratios tested.
Clustering analyses Based on the QCScore values we selected a subset of 1760 combinations that were deemed to be of good to medium quality (QCScore < 5). We then characterized each combination in terms of three parameters We then converted these labels to a binary fingerprint using the following scheme. The potency class for each single agent is coded as a 2-bit vector where the first bit is set to 1 if the single agent is labeled as MP and the second bit is 1 if it is labeled as HP. Similarly, the synergy class is encoded as a 4-bit vector, where a single bit is set to 1 depending on the synergy class label assigned to the combination. Finally, the MOA class for each single agent is represented as a 6-bit vector.
Thus for a given combination, we end up with five bit vectors and the final fingerprint for each combination is simply the concatenation of the individual bit vectors. Thus the fingerprint is a 20-bit vector, whose bit positions are labeled as shown in Fig. S23. We compute the fingerprint for all 1760 combinations and then evaluated a pairwise similarity matrix using the Tanimoto metric. We investigated the use other similarity metrics (cosine and Tversky) and observed no differences in the resultant analyses.
The similarity matrix was converted to a dissimilarity matrix (by subtracting each element from 1.0), which was then used to compute a hierarchical clustering. While there is much literature on identifying clusters, there is no consensus on how one might identify the optimal number of clusters.
In this study we were interested in identifying clusters that were statistically enriched (see below for details) in one or more MOA class pairs (i.e., considering the MOA of both single agents, table S40). As a result we identified the largest number of clusters such that every cluster was enriched in at least one MOA class pair. This procedure identified 6 clusters. The complete dendrogram, color coded according to arbitrary cluster number is shown in Fig. 1D.
It is useful to summarize the combinations belonging to a given cluster. We considered two approaches. First, we examined the occurrence of specific fingerprint bits within clusters. Second, we computed enrichments of MOA combinations within each clusters.

Parameter distributions
To examine how individual bits (say, the bit corresponding to synergy class of synergistic) are distributed within a cluster we use the concept of a bit spectrum (44). Given combinations each represented using an -bit fingerprint, we evaluate the fraction of combinations for which the 'th bit is set to one. This is repeated for all bit positions, resulting in a floating-point vector that is termed the bit spectrum for that set of combinations. Formally, the 'th value of the bit spectrum is defined as where , represents the 'th bit in the fingerprint for the 'th combination and 1 < < .
Since the values range from 0 to 1, this allows us to compare distributions of bit positions between different sets of combinations. We evaluated bit spectra for the 6 clusters obtained above and these are summarized in Fig. S24.
For example, the 886 combinations comprising cluster 1 (red) have single agents that are primarily anti-malarial and ion channel inhibitors. In contrast the 14 combinations comprising cluster 6 have single agents that are nearly all signaling/transport modulators. Single cell photometry Live parasites within iRBC were imaged under constant perfusion using a custom singlecell photometry apparatus described previously (48). Coverslips were prepared by adding 500 µL of 0.1% (w/v) poly-L-lysine and incubating for 10 min to coat the coverslip. They were then dried and stored at 4º C until use. Trophozoite-infected RBCs were resuspended at 0.5% hematocrit in incomplete media with 25 mM HEPES, pH 7.4, and 200 µL was placed on a poly-L-lysine coated coverslip. The cells were incubated for 3 min under standard cell culture atmosphere. Non-adherent cells were washed off and the coverslip was mounted on a custom-designed perfusion chamber for microscopy work. For Ca 2+ imaging, cells were incubated with Fura-2 AM at a final concentration of 5 µM. 0.1% v/v Pluronic F-127 was added to improve dye loading into iRBCs. For perfusion experiments, the parasites were maintained under constant perfusion at 1 mL/min of physiologic buffer (HBSS balanced with 5% O 2 / 5%CO 2 / 90% N 2 at 37 ºC). Changes in cytosolic Ca 2+ upon drug addition were monitored by switching to identical perfusate containing drug at the indicated concentrations. Fura-2 ratiometric values vs time were collected for individual cells, data were converted to [Ca 2+ ] as described below, and Ca 2+ -transients for ≥ 15 parasites were then averaged at each drug concentration.
The single cell photometry (SCP) apparatus includes a custom-built Nikon Diaphot epifluorescence microscope equipped with a 100X oil immersion objective capable of UV transmission (Fluoro, N.A. 1.25, 160/017), and a 16-bit Sensys CCD camera (Tucson, AR) attached to the side port of the microscope. Excitation light was provided by a computer controlled xenon arc lamp (LAMBDA LS, Novato, CA). 2 bandpass filters at 340 nm and 380 nm filtered UV light for ratiometric illumination of Fura-2 (Asahi Spectra Co., Ltd., San Jose, CA), and were housed in a Lambda-10 filter wheel controlled via acquisition software (Imaging Workbench). The excitation light was transported by a liquid light guide (Novato, CA), collimated, and passed through the microscope optics. A filter cube housing a 400 nm dichroic mirror combined with a 410 nm long-pass filter separated excitation from emission. Light power before the objective was monitored with a near UV power meter (Metrologic Model No. 45-545). Exposure time was typically set at 500 ms followed by a 10 sec recovery in total darkness. For ratiometric measurements, data at alternate excitation wavelengths (340 and 380 nm) were collected from a cytosolic region of interest (ROI) within the parasite, while a second region of interest of the same dimensions was drawn outside the iRBC and imaged simultaneously background subtraction.
In situ calcium calibration of Fura-2 For calibration, parasites were loaded with Ca 2+ dye as above, and perfusates were prepared from 2 stocks: 1) HBSS + 25 mM HEPES + 10 mM EGTA (pH 7.4); 2) HBSS + 10 mM EGTA +25 mM HEPES + 10 mM CaCl 2 (pH 7.4). Solutions with specific concentrations of free Ca 2+ were prepared by mixing the stocks at different ratios. The free Ca 2+ concentration was calculated using the "MAXCHELATOR" program (provided by Professor Chris Patton from Stanford University http://maxchelator.stanford.edu/index.html); parameters were set at pH 7.4, 37 ºC, EGTA, ionic strength 150 mM. Before calibration experiments, 10 µM ionomycin was added in order to clamp the internal free calcium concentration to the concentration in the perfusate. Intact trophozoite-stage parasites were first perfused with Ca 2+ -free buffer at a rate of 1 mL/min, and data collected for ≥ 15 parasites. This process was repeated for each perfusate containing ionomycin and specific [Ca 2+ ] (54, 127, 217, 488, 1,000 and 23,700 nM); ratiometric data were collected for ≥ 15 cells at each concentration and averaged. Data were processed and analyzed with Microsoft Excel and SigmaPlot 11.0. Calibration curves (e.g. Fig S26) were fitted to sigmoidal curves and Fura-2 ratiometric data obtained under perfusion with physiological perfusate +/-drug were extrapolated to determine Ca 2+ concentrations (shown on the Y axis for plots of Ca 2+ vs time, see example Fig S27).
Reactive oxidative stress assays Measurement of reactive oxidative species (ROS) production was assessed in synchronized trophozoite parasites as previously described, with minor modifications (50). Briefly, parasites were incubated at room temperature with 10 µM H 2 DCFDA (2',7'dichlorodihydrofluorescein diacetate; Invitrogen, Carlsbad, CA) in the dark. The sample was then centrifuged at 713 × g, washed once with complete media and resuspended at 4% hematocrit. 100 µl was pipetted to each well of a 96-well plate that contained the preserially diluted compound of interest and wells without drug as negative controls assays (49). Plates were then incubated in the dark in 5% O 2 /5% CO 2 /90% N 2 gassed chambers for 3 hrs at 37 °C. After incubation, wells were resuspended by pipetting, and 30 µl was transferred to a replicate 96-well plate and 100 µl of 0.9% NaCl/0.2% Dextrose (Baxter Healthcare, Deerfield, IL) with 500 nM Syto61 (Invitrogen) was added. The sample was incubated again in the gassed chamber at 37 °C for an additional 30 min. Subsequently, the plate was centrifuged at 500 × g for 1 min. The supernatant was removed and the cells were washed once with 0.9% NaCl/0.2% Dextrose solution. The plate was centrifuged again, supernatant removed, and the pellet was resuspended in 100 µl of 0.9% NaCl/0.2% Dextrose solution. The plate was then analyzed by flow cytometry on an Accuri C6 flow cytometer. With an initial gate for events containing DNA (Syto61) and then for modulation of H 2 DCFDA fluorescence by drug induced ROS, normalized to the no drug control (media only) wells, as previously described (50).
Homology model The homology model of PfVps34 was constructed by using X-ray crystal structure of Drosophila melanogaster Vps34 (PDB ID: 2X6H)* as the template. PfVps34 shares over 36% identical residues with its template in the kinase domain. The homology modeling was carried out using the molecular modeling software MOE (v. 2012.10)(Chemical Computing Group, Montreal, Canada). The force field of Amber12:EHT was employed in homology modeling and energy minimization.
Manual drug combination assays using the Chou-Talalay method The effect of combining two drugs together was assayed through use of the Chou-Talalay method of fixed-ratio analysis (40, 53). Briefly, the two compounds were initially screened for their individual cytostatic and cytocidal activities, as described above. Combination analyses were initially run at fixed ratios of these determined values (1:1 unless noted). This results in a set of concentrations all at a multiple or fraction of the compounds' individual IC 50 s or LD 50 s. A combination stock, set at four times the IC50 or LD50, is serially diluted seven times to yield a range of concentrations that can generate a combination growth or survival curve, see Table S41 for a diagram of the plate layout.
These samples are processed as described above, using conditions of either the cytostatic or cytocidal assay. These data are then plotted versus each drug's individual concentrations. For each drug in each assay, there is a 50% point, termed the "pseudo"-IC 50 or "pseudo"-LD 50 , that represents the activity of one drug in the combination. This value is then used to generate a fractional inhibitory concentration (FIC) or fractional lethal dose (FLD) using the following equations for compounds A and B (40, 53): (1) The FIC index of a combination is then used to assign synergy, additivity, or antagonism. The cut-off values for these designations vary (54), but those selected for this work were (15): FICindex, FLDindex and FAIVCindex values were averaged from at least 2 independent trials, each trial done in triplicate and are shown +/-standard error of the mean.
In vivo manual drug combination assays using the Chou-Talalay method In vivo ratios were estimated based on previously determined in vivo levels for ATM (55), LUM (56) and NVP-BGT226 (57).
NITD-609 Synthesis NITD-609 was synthesized in house according to literature procedure 17  *Note: representative heatmaps shown below will be given the notation of assay ID and block ID in the format of "(AID/BID)" which corresponds to the representative line in Table S3.       Compound effect on P. falciparum Dd2 mitochondrial transmembrane potential as assessed by JC1 staining. Evaluation of each compounds ability to perturb the mitochondrial membrane potential in a dose dependent manner. Graphs indicate the dose dependent ΔΨ dissipation expressed as FL2/FL1 (Red/green fluorescence ratio). Each value has been normalized versus the FL2/FL1 ratio of the vehicle treated control to which an arbitrary value of 1 has been assigned. Each point is the mean of three independent experiments performed in duplicate, error bars represent SEM.  Compound Interaction with Artemether on P. falciparum mitochondrial transmembrane potential as assessed by JC1 staining. Evaluation of each compounds interaction with Artemether in modulating the mitochondrial membrane potential in a fixed ratio assay. Graphs indicate the fixed ratio dependent ΔΨ dissipation expressed as FL2/FL1 (Red/green fluorescence ratio). Each value has been normalized versus the FL2/FL1 ratio of the vehicle treated control to which an arbitrary value of 1 has been assigned. Each point is the mean of three independent experiments performed, error bars represent SEM.   Synchronized trophozoite P. falciparum 3D7 parasites labeled with H 2 DCFDA were incubated for 3 h in the presence of serially diluted compound. Cells were then washed and the mean DCF signal of the parasite infected erythrocyte population was analyzed by flow cytometry. DCF signal was normalized to untreated controls; error bars represent SEM of duplicate measurements from at least three independent assays. Synchronized trophozoite P. falciparum Dd2 parasites labeled with H 2 DCFDA were incubated for 3 h in the presence of serially diluted compound. Cells were then washed and the mean DCF signal of the parasite infected erythrocyte population was analyzed by flow cytometry. DCF signal was normalized to untreated controls; error bars represent SEM of duplicate measurements from at least three independent assays.

Fig. S16
PfATG8 positive puncta, quantified as described previously (25) were measured upon treatment with LD 50doses of ATM, LF, ATM+LF (CoArtem) as well as either ATM or LF in the presence of two PI3K inhibitors (GSK2126458 or NVP226). The first two bars (left) and the last two bars (right) in each plot are not statistically different (P value > 0.05), but either of the first two bars are statistically different from either of the second two bars in both plots (P value < 0.05).

Fig. S19
A modified Peter's Suppressive Test was used to initially assess drug-drug interactions in the mouse malaria model. Female CD1, used in panels A and B, or BALB/c mice were infected i.p. with 10 6 P. berghei N parasites. Between two and three hours post infection drug treatment, as indicated above, was initiated by oral gavage of compound resuspended in standard suspension vehicle (0.5% hydroxyethyl cellulose and 0.1% Tween-80; Sigma, St. Louis, MO). Compound administration was continued for an additional two days (3 consecutive days of compound administration in total). Compound dosage was adjusted by weight. Parasitemias were determined by microscopic examination of Giemsa-stained blood films taken daily, beginning on day two. ATM, artemether. Note: The propafenone and GSK-2126458 experiments were performed at the same time, and used the same vehicle and ATM single agent control mice, the graphs have been separated for ease of presentation.

Fig. S20
A modified Thompson's 30-day Test was used to assess the curative activity of NVP-BGT226 (BGT) in combination with artemether (ATM) or lumefantrine (LUM). A) Single agent rate of action, compared to vehicle alone, was assessed by comparison of the natural log transformed number of infected erythrocytes (RBCs) per µl. B) Alteration in mouse mass associated with NVP-BGT226 treatment in addition to P. berghei infection.
In addition to the transient drug-induced loss of weight mice also experienced lethargy, hunched posture and decreased grooming. All drug-induced effects were temporary and decreased after cessation of drug treatment.

Fig. S21
The minimum spanning tree derived from the full network of 2,134 combinations, with edges colored by manually assigned combination class. This network serves to highlight that manual inspection must be performed due to the lack of robustness in the DBSumNeg metric, which may flag a combination as synergistic whereas visual inspection of the response matrix would suggest otherwise.

Fig. S23
The layout of the fingerprint representation used to characterize the potency, synergy and MOA of each combination.

Fig. S24
A bit spectrum for each of the six clusters, highlighting distribution of parameter values within each cluster.

Zaldaride
Compounds listed passed our quality control analysis (see Methods) with at least one assay producing a DBSumNeg value ≤ -3 with the compound listed in the heading. Compounds producing a DBSumNeg ≤ -3 with more than one of the indicated compounds are shown in bold, those present in more than 3 are also highlighted in yellow.