Functionally defined therapeutic targets in diffuse intrinsic pontine glioma.

Diffuse intrinsic pontine glioma (DIPG) is a fatal childhood cancer. We performed a chemical screen in patient-derived DIPG cultures along with RNA-seq analyses and integrated computational modeling to identify potentially effective therapeutic strategies. The multi-histone deacetylase inhibitor panobinostat demonstrated therapeutic efficacy both in vitro and in DIPG orthotopic xenograft models. Combination testing of panobinostat and the histone demethylase inhibitor GSK-J4 revealed that the two had synergistic effects. Together, these data suggest a promising therapeutic strategy for DIPG.

Introduction

Diffuse Intrinsic Pontine Glioma (DIPG) is the leading cause of pediatric brain tumor death. Median survival is only 9 months, and 5-year survival is less than 1%[1]. Historically, DIPG research has been limited by a dearth of tumor tissue available for study and lack of experimental model systems. Recently, patient-derived DIPG cell cultures and orthotopic xenograft models have been established[2,3]. Genomic studies identified recurrent K27M mutations in genes encoding histone H3[4-6], resulting in inhibition of the Polycomb Repressive Complex-2 (PCR2) methyltransferase and global hypomethylation of H3K27[7,8]. Armed with these important insights and new tools to design clinical trials based on preclinical data, an international collaboration has formed to pool DIPG tissue resources and evaluate functional targets for DIPG therapy.

Results

We assembled a panel of 14 patient-derived DIPG cell cultures, representing neurosphere and adherent models, obtained from both biopsy or autopsy samples (methods and [9-11]; Fig. 1a), representing the breadth of DIPG cell cultures available worldwide at the initiation of the study. Two additional cell cultures (methods and [3]) were subsequently added (Fig. 1). Clinical, histopathological and molecular characteristics of the cultures are summarized in Supplementary Tables 1–2. We performed whole exome sequencing (WES) of culture models used for drug screening when primary tumor DNA was available (WES for 7/9 of tumors not previously reported[10] are described here); additional DIPG tumors not associated with cultures (n = 15) were also sequenced to add to publically available DIPG genomic data (n = 22 DIPG tumors with WES data; Supplementary Fig. 1 and Supplementary Table 3).

Figure 1

Functionally-defined targets in DIPG therapy

(a) Chemical screen of 14 patient-derived DIPG cell cultures. Heatmap demonstrating DIPG cell line sensitivity to each of the 83 agents tested. The values shown are the absolute IC50 divided by the maximum dose. All drugs had a maximum dose of 10 µM except for vismodegib, sodium butyrate, pazopanib, alisertib, and vemurafenib, which had a maximum dose of 100 µM. Values are shown as gradations of red to white, with red representing submicromolar IC50 values, white indicating IC50 greater than the maximum dose for that drug (i.e. 10 µM or 100 µM), and pink showing the range in between. Grey boxes indicate those drugs not included in the screen for that cell line. Numbers corresponding to the drug names in the key below the heatmap are listed on the horizontal axis, culture IDs listed on the vertical axis. The histone status of each culture used in the screen is indicated by green (wild type, WT), yellow (H3.3K27M, H3F3A-K27M) or blue (H3.1K27M, HIST1H3B-K27M); see also Supplementary Table 1. Recurrent “hits” are visualized as a column of red or pink. (b) Dose-response curves: Patient-derived DIPG lines (SU-DIPG IV, SU-DIPG-VI, SU-DIPG-XIII, JHH-DIPG1, SF77613) were treated with the indicated drugs at 0.001/0.01/0.1/1/10 µM or 0.1% DMSO control in at least triplicate (n = 3 wells) and cell viabilities were assessed at 72 hr. Data are expressed as relative to the 0.1% DMSO control values. A pediatric cortical GBM line (SU-pcGBM2; histone WT; orange curves) was treated in parallel for comparison in a subset. Data are shown as mean ± SD. (c) Panobinostat time course: DIPG cells were treated with panobinostat in quadruplicate (n = 4 wells) at indicated concentrations (25 nM – 500 nM) or 0.1% DMSO vehicle control. Cell viabilities were assessed at 0, 24, 48 and 72 hrs of treatment. Data are shown as mean ± SD. **P < 0.01, ***P < 0.001 (two-tailed t test results shown for the lowest concentration to reveal a significant difference at 48 or 72 hours).

Using this DIPG culture panel, we screened 83 drugs selected by pediatric neurooncologists as promising small molecule compounds or traditional chemotherapeutic agents. Readout was cell viability at 72 hours. We found sensitivity to a limited number of drugs (14/83 compounds demonstrated activity against 3 or more DIPG cultures; Fig. 1a, Supplementary Fig. 2). The chemical screen and subsequent validation experiments revealed DIPG cell sensitivity to histone deacetylase (HDAC) inhibitors, as well as a limited number of other agents (Fig. 1a–c, Supplementary Fig. 2). Notable “misses” that highlight the resistance of this tumor to traditional chemotherapies include temozolomide, carboplatin and vincristine (Fig. 1a).

To examine the molecular context of functional drug screening results, RNA deep sequencing (RNAseq) was performed to analyze expression in primary DIPG tumors versus normal, paired cerebral cortex tissue (Supplementary Fig. 3). Screening and RNAseq data were then used in an integrative predictive modeling approach to tumor drug sensitivity using Target Inhibitor Maps (TIMs)[12,13]; the DIPG TIM predictions and in vitro validations are presented in Supplementary Figures 3–4, Supplementary Table 4.

The multi-HDAC inhibitor panobinostat[14], among the most effective agents (Fig. 1, Supplementary Fig. 2), was selected for further exploration. In all, 12/16 DIPG cultures demonstrated sensitivity (Fig. 1). Compared to vorinostat, panobinostat illustrated substantially greater potency (Fig. 1b, Supplementary Fig. 5). Testing cell viability at 24, 48 and 72 hours at various dose levels of panobinostat demonstrated decreased viability by 48 hours (Fig. 1c). Expression of proliferation-associated genes MKI67 and CCND1 decreased with time after panobinostat exposure (Supplementary Fig. 6a). FACS analysis of EdU incorporation and Annexin V staining illustrated dose-dependent decreased cell proliferation and increased cell death (Fig. 2a, Supplementary Fig. 6b). shRNA-mediated knock-down of HDAC1 or HDAC2 decreased cell viability (Fig. 2b, Supplementary Fig. 6c – e). Western blot analyses demonstrated a dose-dependent increase in H3 acetylation and H3K27-trimethylation following panobinostat (Fig. 2c), suggesting partial rescue of the H3K27M-induced global hypotrimethylation phenotype. Increased K27 trimethylation is an unexpected effect of the drug, but is consistent with recent findings that acetylated H3K27 can “detoxify” K27M-induced inhibition of PCR2[15]. RNAseq performed on panobinostat- or vehicle-treated DIPG cells revealed sweeping changes in gene expression, including normalization of the K27M gene expression signature[16] and decrease in the oncogenic MYC target gene signature (Supplementary Figs. 7–8).

Figure 2

Panobinostat is a promising therapy for DIPG

(a) FACS analysis of DIPG tumor cell proliferation and cell death: Top row: Overlapping histogram plots of EdU FACS analyses are shown on the left; quantifications of EdU+ cell population levels from each condition are shown in bar plots on the right for DIPG cell cultures SU-DIPG-VI and SU-DIPG-XIII (both H3.3K27M mutant cell lines). Bottom row: Left, Overlapping plots of Annexin V, DAPI FACS analyses; Right, bar plots show early apoptotic (AV+DAPI+) or late apoptotic (AV+DAPI+) cell population levels from each condition for each cell line as above. (b) HDAC1 and HDAC2 knock-down in DIPG cells using shRNA verify panobinostat mechanism of action in four DIPG cell lines. Cell viability assays at each time point for each cell line were performed in triplicate (n = 3 wells); data are expressed relative to Day 0 and are shown as mean ± SD. Note the varying growth rates of cell cultures result in varying y-axes depicting relative change in cell viability. **P < 0.01; ***P<0.001 (Two-way ANOVA). (c) Panobinostat increases histone-3 acetylation and restores H3K27 trimethylation. Western blot analyses of histone-3 acetylation and H3K27 trimethylation (H3K27me3) in H3K27M mutant DIPG cell lines SU-DIPG-VI and SU-DIPG-XIII (left blots) and in 293T cells expressing a mutant H3.3K27M-HA tagged construct (293-H3.3-K27M-FH[8]; right blots). Controls included total protein levels of H3, HDAC1, HDAC2 and EZH2. Expression of the HA tag in the 293T cells confirms expression of the H3.3K27M-FH construct. (d) Schematic illustrating convection enhanced delivery strategy to infuse drug into brainstem. Blue illustrates approximate distribution of the infused solution. (e) Distribution of infusate illustrated by delivering blue dye to the brainstem by CED. Ventral side of a mouse brain is shown immediately following CED delivery of Coomassie Blue dye. Scale bar = 3 mm (f) in vivo bioluminescent imaging of DIPG xenografts 7 days following CED delivery of panobinostat (T = treated with panobinostat) or vehicle control (C = control). The heat map superimposed over the mouse head represents the degree of photon emission by DIPG cells expressing firefly luciferase. Scale bar = 3.5 cm. (g) in vivo DIPG xenograft tumor growth as measured by change in bioluminescent photon emission over the seven days following (g) CED delivery of panobinostat. panobinostat = red squares (n = 5 mice) and vehicle control = blue circles (n = 4 mice). Data points represent the change in maximum photon flux (percent of baseline) between Day 0 and Day 7 for each mouse. (h) As in (g), with systemic administration of panobibostat. Three systemic dose levels were used, 1 mg/kg (n = 6 control, 8 treated mice) or 10 mg/kg (n = 5 mice per group) delivered IP on M,W,F or 20 mg/kg (n = 7 mice per group) delivered once per week. Error bars, s.e.m. *P < 0. 0.5; **P < 0.01; N.S. indicates P > 0.05 (two-tailed t test). (i) Systemic delivery of panobinostat prolongs survival in a histone H3 wild type DIPG orthotopic xenograft model IBs-W0128DIPG. Panobinostat 10 mg/kg I.P. doses given as indicated by arrows. n = 10 per group; P = 0.0179 (log-rank analysis).

We next used convection-enhanced delivery (CED)[17] to directly infuse panobinostat into the pons of a DIPG orthotopic xenograft model engineered to express firefly luciferase (SU-DIPG-VI-luc; H3.3K27M subtype). In vivo bioluminescent imaging was performed immediately prior to and 7 days after panobinostat (or vehicle control) to assess tumor growth (Fig. 2d – g). We confirmed that panobinostat does not down-regulate the luciferase transgene (Supplementary Fig. 9a). A marked effect of panobinostat on the rate of DIPG xenograft growth was observed: vehicle-treated control mice exhibited xenograft growth ∼6.5 fold greater than mice treated with a single dose of panobinostat by CED (n = 4 vehicle controls, 5 treated mice; P < 0.05 by two-tailed t test; Fig. 2g and Supplementary Fig. 9b). We next investigated pontine penetration of panobinostat following systemic administration using liquid chromatography-tandem mass spectrometry (LC-MS/MS). Following a single 20 mg/kg IP dose, pontine panobinostat levels were 0.068 ng/mg, or ∼200 nM (Supplementary Table 5). As the DIPG IC50 for panobinostat was found to be ∼100 nM (Fig. 1c), we reasoned that systemic delivery may prove effective. Mice with SU-DIPG-VI-luc orthotopic xenografts were treated with 1, 10 or 20 mg/kg IP panobinostat. At one week, a significant reduction in tumor growth in panobinostat-treated animals in 10 mg/kg and 20 mg/kg treatment groups occurred compared to vehicle-treated controls (n = 5 per group and 7 per group; P < 0.05 and 0.01 by two-tailed t test, respectively; Fig. 2h). At 4 weeks, tumor growth was reduced only in the 20 mg/kg group (Supplementary Fig. 9c). We tested a second patient-derived DIPG brainstem orthotopic xenograft model (LI-F; H3WT subtype) using an alternative panobinostat dose schedule (10 mg/kg IP, 5 days on, 5 days off). Mice treated with panobinostat exhibited significantly prolonged survival compared to vehicle-treated controls (Fig. 2i, n = 10 per group, P < 0.05 by log-rank analysis).

Highlighting the need for effective combination therapies, DIPG cells that survive chronic panobinostat exposure demonstrated resistance when re-challenged with the drug (Supplementary Fig. 10a). As a recent elegant study demonstrated utility of the histone demethylase inhibitor GSKJ4[18], we tested GSKJ4 alone or in combination with panobinostat. Consistent with the findings of Hashizume and colleagues[18], GSKJ4 decreased cell viability in H3K27M DIPG cell cultures (Supplementary Fig. 10b). Panobinostat synergizes with GSKJ4 in H3.3K27M mutant DIPG cells (Supplementary Fig. 10 c – d). In contrast, panobinostat does not improve the efficacy of the commonly used but ineffective drug temozolomide (Supplementary Fig. 10e).

Discussion

The present study represents a collaborative effort to identify functionally relevant targets and develop an effective DIPG treatment strategy. Integration of genomic data and chemical screening data, validated by in vivo preclinical testing, converge to identify the potent HDAC inhibitor panobinostat as a promising therapy for DIPG. That the leading therapeutic candidate is an epigenetic modifying agent is congruent with the discovery that the majority (∼80%) of DIPG tumors harbor K27M mutations in the histone H3 genes HIST1H3B or H3F3A[4-6,19], which results in broad epigenetic dysregulation[7,8]. Indeed, epigenetic modifying therapies are emerging as the most promising class of agents for DIPG[18]. The FDA-approved drug panobinostat could be rapidly translated to the clinic. Other classes of epigenetic modifying agents, such as histone demethylase inhibitors[18], may prove clinically useful for DIPG, alone or in combination, when such agents become available for clinical use.

Key caveats of our study should be noted: the drug screen used was not exhaustive, and future screening efforts may elucidate additional useful drugs. Further, in vitro systems do not model the cellular complexity of the tumor’s microenvironment and thus do not screen for drugs that target key microenvironmental interactions. Nevertheless, by leveraging newly generated experimental model systems in an international collaborative effort, the data presented here define a promising therapy for a presently untreatable disease.

Materials availability

All cell lines described here can be obtained through an MTA with the originating institution (Stanford University, Johns Hopkins University, Institut Gustave Roussy, VU Medical Center, or Texas Children’s Cancer Center).

Online Methods

Primary cell culture establishment & propagation

For all human tissue studies, informed consent was obtained and Institutional Review Board (IRB) approval was granted at each participating institution. Cell culture models (SU, NEM, JHH. Li and VU cultures) were validated by DNA fingerprinting. SU, NEM and JHH lines were authenticated using short tandem repeat (STR, see Supplementary Table 2). Li cultures were authenticated using single nucleotide polymorphism (SNP) analysis. VU lines was authenticated by routine Sanger sequencing for histone gene mutation. All cell cultures were routinely tested for mycoplasma. Please see Supplementary Table 1 for age, sex and other clinical characteristics associated with each cell culture. Cultures derived from the tumors of both male and female individuals aged 2 – 9 years.

Stanford cell lines (SU) and primary tumor cell culture

As previously described[2,9,10], tumor tissue was harvested under sterile conditions in the autopsy suite, collected in serum-free DMEM/12 (Cellgro) and transported on ice to the laboratory. The tissue was then minced with a sterile No.10 scalpel and washed twice in HEPES-HBSS (HeHBSS) to remove debris. Minced tissue was then added to dissociation buffer (HeHBSS) with DNaseI (250U/ml) and collagenase type IV (1mg/ml). Minimal mechanical force was utilized and the cells were allowed to dissociate at 37C on a Nutator. The cells were then sequentially strained through a 100, 70 and 40 µM cell strainer and centrifuged and processed through a sucrose gradient. Dissociated cells were treated with ACK/RBC lysis buffer (0.15M NH4Cl, 1.0mM KHCO3 and 0.1mM Na2-EDTA) and plated at subclonal density of 100 cells/ml, for neurosphere formation in Tumor Stem Media (TSM) consisting of Neurobasal(-A) (Invitrogen, Carlsbad, CA), B27(-A) (Invitrogen), human-bFGF (20ng/ml) (Shenandoah Biotech, Warwick, PA), human-EGF (20ng/ml) (Shenandoah), human PDGF-AB (20ng/ml) (Shenandoah) and heparin (10ng/ml). Cells were subsequently cultured as tumor neurospheres and passaged every one to two weeks.

Johns Hopkins cell line (JHH-DIPG-1)

The JHH-DIPG-1 cell line was generated using the protocol described above for the Stanford cell lines.

Gustave Roussy cell lines (NEM) and primary cell cultures

Tumor pieces were collected into DMEM (Dulbecco’s Modified Eagles Medium, PAA Laboratories GmbH, Pasching, Austria). Biopsies were cut into 1 mm3 pieces and placed into either DMEM for immediate processing or into freezing medium (90% serum, 10% DMSO) prior to being progressively cooled to −80°C. Cell dissociation was performed mechanically by passage through increasingly finer needles (19G to 26G). Single cells were seeded in AmnioMAX™ C-100 complete medium containing gentamycine, L-glutamine and FBS (Gibco, Invitrogen, Paisley, UK) and maintained at 37°C in a 5% CO2 humidified atmosphere. The cells were further cultured until appearance of adherent cells and colony formation and then weekly passaged.

Texas Children’s Hospital xenograft models and recovered cultures (Li)

Tumor tissue was dissociated and immediately xengrafted in NOD SCID mice. Tumors were passaged by serial xenografting and cell culture experiments performed in recovered xenografts. The recovered xenograft cultures contain approximately 20% mouse cells and 80% human DIPG cells. Recovered xenograft cells were authenticated using single nucleotide polymorphism (SNP) analysis as described below.

SF7761 DIPG cell line

was cultured as previously described[3].

VU DIPG cell lines

VU-DIPG-A and VU-DIPG-B were cultured as previously described[11].

Chemical screens

A panel of 83 drugs (Fig. 1a) was developed for screening drug sensitivity of DIPG cells; the majority of these compounds are in clinical trials currently. Each compound was plated in four concentrations (10 µM, 1 µM, 100 nM, 10 nM) in triplicate in 384-well format. Early passage cells cultured from autopsy or biopsy tissue were plated at a density of 2,500 cells per well and incubated at 37C with humidified 5% CO2 for three days. Cell viability was determined by the CellTiterGlo luminescent assay (cat #G7572, Promega, WI) per manufacturer’s instructions and measured with the IVIS Lumina II imaging system. IC50 values were determined by a nonlinear best-fit method using Excel Solver.

DIPG orthotopic xenograft in vivo studies

All animal procedures were performed with approval from each institution’s animal care and use committee (the Stanford University Administrative Panel on Laboratory Animal Care [APLAC] and the Baylor Institutional Animal Care and Use Committee [IACUC]) and adhered to the NIH guide for the care and use of laboratory animals. Both male and female animals were used equally.

Patient-derived H3K27M DIPG orthotopic xenograft model (Stanford University)

DIPG pontine xenografts were generated as previously described [2]. Briefly, a single cell suspension was made of SU-DIPG-VI-luc neurospheres and 100,000 cells (50, 000 cells/µl, 2 µl) were stereotactically injected into the fourth ventricle/pons of NOD-SCID-IL2 gamma chain-deficient cold-anesthetized postnatal day 2 mouse pups by stereotactic injection through a 31 gauge (G) burr hole (stereotactic coordinates: 3 mm posterior to lambda suture and 3 mm deep.)

CED was performed two months post-xenograft using a 31 G Hamilton with a 27 G needle threaded on the outside in isoflurane-anesthetized animals[20]. (The inside needle is 1 mm longer than the outside needle in order to create the convection effects.) Infusion was performed using a digital pump set at a rate of 0.4 µl/min. Five µl of 2 µM panobinostat (SelleckChem) was delivered over a period of 12.5 min. Stereotactic coordinates (from surface of brain) used were from lambda AP −0.8mm, lateral 1 mm and 5 mm deep. Controls received vehicle (DMSO diluted in 5% dextrose). In vivo bioluminescent imaging was performed prior to CED of panobinostat or vehicle, and again 7 days later, using an IVIS imaging system (Xenogen) under isoflurane anesthesia.

Systemic administration of panobinostat was performed with intraperitoneal injection three days per week (M, W, F) for dose levels 1 mg/kg and 10 mg/kg, and administered once per week for the 20 mg/kg cohort. Controls were injected I.P with an identical volume of DMSO vehicle. Panobinostat was dissolved in 70 mg/ml DMSO, then serially diluted in water to a concentration of 0.1 to 2 mg/ml (depending on the target dose level) such that 10 µl/gram or ∼200 µl total was administered I.P each dose.

For IVIS imaging analyses of in vivo DIPG tumor growth, animals were imaged at baseline. Animals were excluded if no tumors were present, and then animals were randomized to control and treatment groups such that each group had equivalent distribution of initial tumor sizes.

Patient-derived orthotopic xenograft model IBs-W0128DIPG (Texas Children’s Hospital)

This is a serial xenograft model. SCID mice were implanted with 50,000 IBs-W0128DIPG cells (H3WT) into the brainstem (n=10 per group) and treatment (10 mg/kg, i.p.) was started 17 days post tumor implantation (5 day-on, 5 day off, followed by another 4 days on) for a total of 9 doses.

Blinding was not performed.

Liquid Chromatography Tandem Mass Spectrometry (LC-MS/MS) analysis of Panobinostat concentrations in tissues and serum

Sample preparation

A single 20 mg/kg dose was delivered intra-peritoneally in NOD-SCID-IL2 gamma chain-deficient mice, and tissue samples collected 30 minutes later for analysis using liquid chromatography/mass spectrometry (LC-MS/MS). Tissues samples were weighed and 1 volume of bullet blender beads (Next Advance) and 2 volume of Milli Q water were added. Tissues were homogenized by a bullet blender (Next Advance) at 4°C according to manufacturer’s instruction. Stock solution of neat panobinostat (Selleck Chemicals) and internal standard (IS) Panobinostat-d8 (Toronto Research Chemicals, Canada) were prepared in methanol and further diluted in 50% methanol. For spiked standards, 25 µl of neat panobinostat (5 ng/ml-20 µg/ml) and 10 µl of IS (100 ng/ml) were added to 25 µl of blank tissue homogenate or serum. For samples, 25 µl of 50% methanol was used instead of neat panobinostat. 150 µl of methanol/acetonitrile 50:50 (v/v) was added to the mixture and vortexed vigorously for 1 min followed by centrifugation at 10,000 g for 5 min. For brain standards and samples, the supernatant was diluted 3 times with 0.1% formic acid (v/v) in water. For serum and kidney standards and samples, the supernatant was diluted 10 times and 100 times respectively, in 30% methanol with 0.1% formic acid.

LC-MS/MS analysis

The LC-MS/MS system consists of a QTRAP 4000 mass spectrometer (AB SCIEX) coupled to a Shimadzu UFLC system. LC separation was carried out on an ACE C18 column (50 mm × 2.1 mm, 5 µm) with isocratic elution using a mobile phase composed of 65% water and 35% methanol. 0.1% of formic acid was added to both aqueous and organic phases. The flow rate was set at 0.25 ml/min. Column temperature was 25°C. The analysis time was 2.5 min. The injection volume was 10 µl. The mass spectrometer was operated in the positive mode with multiple-reaction monitoring (MRM). The m/z 350.2→158.1 and 358.3→164.2 transitions were used for Panobinostat and Panobinostat-d8 (IS). Data acquisition and analysis were performed using the Analyst 1.6.1 software (AB SCIEX).

Western blot analyses of protein levels and RT-qPCR analyses of RNA levels

For whole cell extract analyses, the cells were lysed by using RIPA buffer plus proteinase inhibitors. For chromatin extract analyses, chromatin extract was prepared as follows: Cells were lysed in Buffer A' (10 mM Tris-HCl, [pH 7.5], 10 mM KCl, 1.5 mM MgCl2, 0.34 M sucrose, 10% glycerol, 1 mM DTT and Protease inhibitors) with 0.2% Triton X-100 on ice for 10min. Nuclei were pelleted by low speed centrifugation (5 min, 1,300 × g, 4°C) and then lysed in buffer B' (3 mM EDTA, 0.2 mM EGTA, 1 mM DTT and protease inhibitors). Insoluble chromatin was collected by centrifugation (5min, 1,700 × g, 4°C). The chromatin extract was prepared by MNase digestion in Buffer A' with 2mM CaCl2 at 37°C for 5min. To stop the digestion, EDTA was added. NaCl and Triton X-100 were then added to concentrations of 150mM and 0.5% respectively. Protein concentrations were measured by using Pierce™ BCA™ Protein Assay (Thermo Scientific™) and equal amounts of cell extracts were loaded for western blot analysis. The primary antibodies used in this study are anti-mTOR (CST, #2972), anti-β-tubulin (Abcam, ab6046), anti-H3K27me3 (Abcam, ab6002), anti-H3Ac (H3-acetyl K9 + K14 + K18 + K23 + K27, Abcam, ab47915), anti-H3 (Abcam, ab1791), anti-HDAC1 (CST, #5356), anti-HDAC2 (CST, #5113), anti-EZH2 (Active motif, 39933), anti-HA (Santa Cruz, sc-7392).

For RNA analyses, RNA extraction was performed by using Qiazol Lysis Reagent (79306, Qiagen, Venlo, Netherland) according to the manufacturer’s instructions. Reverse transcription was performed by using the High Capacity cDNA Reverse Transcription Kit (4368813, Invitrogen) and real-time qPCR was performed by using 2x Maxima SYBR Green qPCR Master Mix (#K0251, Thermo Scientific). qPCR primers are listed as follows:

CCND1: TTCAAATGTGTGCAGAAGGA and GGGATGGTCTCCTTCATCTT

MKI67: TTGGAGAATGACTCGTGAGC and CGAAGCTTTCAATGACAGGA

HDAC1: ATCGGTTAGGTTGCTTCA and TCATTCGTGTTCTGGTTAGTC

HDAC2: GACAGTGGAGATGAAGATGGA and TTCTGATTTGGTTCCTTTGG

GAPDH: TGACTTCAACAGCGACACCCA and CACCCTGTTGCTGTAGCCAAA

Cell viability, cellular proliferation and apoptosis analyses

GSK-J4, BEZ235, AUY922, carfilzomib, panobonistat, vorinostat, flavopiridol, BKM120, AZD8055, RAD001 and temozolomide were purchased from Selleckchem (Houston, TX, USA). For generating drug dosage curves or time-course growth curves, cells were plated at a density of 5,000 cells per well in 96-well plates in at least triplicate and then subjected to drug treatment for 72 hrs as indicated. Cell viability was then measured by using CelltiterGlo assay according to the manufacturer’s instructions (G7571, Promega, WI, USA), and data were collected on a TECAN Infinite 200 plate reader. For testing combinatory effects of two drugs, the cells were treated with each drug individually or in combination before subjecting to CelltiterGlo assay. CalcuSyn software (Biosoft, Cambridge, UK) was used to calculate the combinatory index (CI). CI less than 1.0 was considered to be synergistic.

Cellular proliferation was measured by using the Click-iT EdU Alexa Fluor 488 Flow Cytometry Assay Kit according to the manufacturer’s instructions (C35002, Invitrogen, CA, USA). EdU+ population represented the proliferating cell population. DIPG cells were treated with 0.1% DMSO vehicle or panobinostat at 0.1 or 0.5 µM respectively for 24 hours and then incubated with 10 µM EdU for another 16 hours. Cells were then subjected to EdU FACS analysis to assess the proportion of cells in S-phase (EdU+). Cellular apoptosis was measured by using Annexin V-FITC Apoptosis Detection Kit II (556570, BD Biosciences, CA, USA) with some minor modifications. DAPI was used in the staining instead of PI. DIPG cells were treated with 0.1% DMSO or panobinostat at 0.1 or 0.5 µM respectively for 48 hours and then subjected to Annexin V/DAPI FACS analysis to assess the proportion of cells undergoing apoptosis. For FACS analyses, the stained cells were analyzed using BD Fortessa FACS machine (BD Biosciences, CA, USA). The data were analyzed using FlowJo software (FlowJo, LLC, OR, USA).

shRNA expressing lentivirus preparation and infection

shRNA expressing lentiviral constructs against human mTOR, HDAC1 and HDAC2 from the RNAi consortium human collection were purchased from Sigma. Lentiviral expressing constructs were co-transfected with packaging plasmids (pDelta 8.92 + VSV-G) into 293T cells to generate lentiviral particles. Lentiviral particles were then concentrated by polyethylene glycol (PEG) precipitation method. The precipitated lentiviruses were resuspended in PBS and aliquoted for −80 °C storage. For lentiviral infection, DIPG cells were incubated with shRNA expressing lentivirus for 16 h. 48 h post infection, puromycin (0.5 µg/ml) was added to select virally infected cells for further experiments. Cells were seeded in triplicate into 96-well plates and cell viability was assessed using the CelltiterGlo assay at 0, 3 and 5 days.

Testing for resistance to panobinostat in DIPG cell cultures

SU-DIPG-VI and SU-DIPG-XIII cells were exposed to panobinostat at ∼IC50 dose (0.05 µM for SU-DIPG-VI and 0.1 µM for SU-DIPG-XIII) chronically for 3 weeks. The old culturing medium was replaced with fresh drug-containing medium every 3–4 days.

DNA and RNA sequencing

For all human tissue studies, informed consent was obtained and Institutional Review Board approval was granted at each participating institution. Please see Supplementary Table 1 for age, sex and other clinical characteristics associated with each cell culture. Age range 2 – 21 years and approximately equal distribution of male and female individuals.

Generation of exome-capture libraries

Exome libraries of matched pairs of tumor / normal genomic DNAs were generated using the Agilent SureSelect XT kit and Agilent Automation Systems NGS system per manufacturer’s instructions. 1 ug of each genomic DNA was sheared using a Covaris E220 to a peak target size of 150 bp. Fragmented DNA was concentrated using AMPureXP beads (Beckman Coulter), and DNA ends were repaired using T4 DNA polymerase, Klenow polymerase, and T4 polynucleotide kinase. 3’ A-tailing with exo-minus Klenow polymerase was followed by ligation of Agilent paired-end oligo adapters to the genomic DNA fragment. Ligated DNA was PCR amplified for 8 cycles and purified using AMPure XP beads and quantitated using the Quant-It BR kit (Invitrogen). 500 ng of sample libraries were hybridized to the Agilent biotinylated SureSelect v4 Exome + UTR Capture Library at 65°C for 72 hr following the manuufacturer’s protocol. The targeted exon fragments were captured on Dynabeads MyOne Strepavidin T1 (Invitrogen), washed, eluted, and enriched by amplification with Agilent post-capture primer and an indexed reverse primer for multiplexing 12 additional cycles. After purification of the PCR products with AMPure XP beads, the quality and quantity of the resulting exome libraries were analyzed using an Agilent Bioanalyzer High Sensitivity chip.

Somatic point mutation identification by exome capture sequencing

The methodology described below is a minor modification of the approach described in [21]. All captured DNA libraries were sequenced with the Illumina HiSeq in paired end mode, yielding 80 base pairs from the final library fragments. The reads that passed the chastity filter of Illumina BaseCall software were used for subsequent analysis. Matepairs were pooled and mapped as single reads to the reference human genome (NCBI build 36.1, hg18), excluding unordered sequence and alternate haplotypes, using Bowtie [22] keeping unique best hits, and allowing up to two mismatched bases per read. Likely PCR duplicates, defined as reads with equal match intervals on the reference genome, were removed; and individual basecalls with Phred quality less than Q20 were excluded. A mismatched base (SNV) was identified as a variant when 1) it had at least six reads of support, 2) it constituted at least 10% of the coverage at that position, 3) it was observed on both strands, and 4) it fell within 50 bases of a region targeted for capture.

In matched samples (samples with matched normals), a variant was called somatic if 4) there was 8x coverage in the matched normal, 5) it did not occur in the matched normal sample in more than two reads and 2% of the coverage (or 4% of the coverage if the tumor variant fraction was at least 20%), and 6) if it had any support in the matched normal, it was not present in either dbSNP (v137) or the Exome Sequencing Project dataset (ESP6500).

In unmatched samples, a variant was called probable somatic if it affected the same codon as either a somatic point mutation observed in the matched samples or a somatic mutation from COSMIC (v66).

Finally, all somatic/probable somatic variants were screened against the full set of benign samples. Variants were removed from further consideration if they appeared in any benign sample with at least 3 reads and 10% of the coverage.

Identification of coding indels in exome capture data

The methodology for identifying indels in exome capture data was adapted from [23] with minor modifications. Reads for which Bowtie was unsuccessful in identifying an ungapped alignment were converted to fasta format and mapped to the target regions, padded by 200 bases on either side, with cross_match (v0.990329, http://www.phrap.org), using parameters –gap_ext −1 –bandwidth 10 –minmatch 20 –maxmatch 24. Output options were –tags –discrep_lists –alignments. Alignments with an indel were then filtered for those that: 1) had a score at least 40 more than the next best alignment, 2) mapped at least 75 bases of the read, and 3) had two or fewer substitutions in addition to the indel. Reads from filtered alignments that mapped to the negative strand were then reverse-complemented and, together with the rest of the filtered reads, remapped with cross-match using the same parameters (to reduce ambiguity in called indel positions due to different read orientations). After the second mapping, alignments were re-filtered using criteria 1–3. Reads that had redundant start sites were removed as likely PCR duplicates, after which the number of reads mapping to either the reference or the non-reference allele was counted for each. An indel was identified as a variant when 1) it had at least six reads of support, 2) it constituted at least 10% of the coverage at that position, 3) it was observed on both strands, and 4) it fell within 50 bases of a region targeted for capture. Indels were reported with respect to genomic coordinates. For insertions, the position reported is the last base before the insertion. For deletions, the position reported is the first deleted base.

In matched samples, an indel was called somatic if 4) there was 8x coverage in the matched normal, 5) if it had any support in the matched normal, it was not present in either dbSNP (v137) or the Exome Sequencing Project dataset (ESP6500). In both matched and unmatched samples, all remaining indels were screened against the full set of benign samples. Variants were removed from further consideration if they appeared in any benign sample with at least 3 reads. Supplemental Supplemental Table S2 provides the somatic point mutation and indels called in matched and unmatched samples in drug target genes, while Table S3 provides the somatic point mutations and indels called using matched samples.

Annotation

We annotated the resulting somatic mutations using CCDS transcripts wherever possible. If no CCDS transcript was available, we use the coding regions of RefSeq transcripts. HUGO gene names were used.

Exome copy number analysis

Copy number aberrations were quantified and reported for each gene as the segmented normalized log2-transformed exon coverage ratios between each tumor sample and its matched normal as described in [24]. For unmatched samples, we used the average coverage derived from the full set of benign samples as a “normal pool;” for chromosomes X and Y only, we formed separate normal pools from the XX and XY benign samples. We identified segments as focal if they contained ten or fewer CCDS gene annotations. To identify a gene as gained or lost we first considered the segmented data, requiring a copy number gain or loss of 30% (ratio >= 1.3 or <= 0.7) to call the gene as gained or lost. In addition, we considered the distribution of copy number estimates for each gene’s individual exons. If the mean exonic copy number showed a gain or loss of at least 30% and deviated from the null hypothesis by at least 2.5 standard deviations, we called the gene as gained or lost.

For the total copy-number call over all tumor tissue samples, we counted each gain or loss of at least 60% as +/−2 on its respective segment, and each gain or loss of 30–60% as +/−1, and summed these contributions across all tumor tissue samples. In order to omit germline CNVs from the total, we looked at segmented exon coverage ratios between the individual normal samples and the normal pool. Any segment that was gained or lost by 30% in at least two normal samples was flagged as a probable CNV, and was not included in the total copy-number sum.

Generation of RNAseq libraries

RNAseq libraries for transcriptome analysis were prepared using the TruSeq RNA Sample Preparation Kit (Illumina) and Agilent Automation NGS system per manufacturers’ instructions. Sample prep began with 1 µg of total RNA from each sample. Poly-A RNA was purified from the sample with oligo dT magnetic beads, and the poly(A) RNA was fragmented with divalent cations. Fragmented poly-A RNA was converted into cDNA through reverse transcription and were repaired using T4 DNA polymerase, Klenow polymerase, and T4 polynucleotide kinase. 3’ A-tailing with exo-minus Klenow polymerase was followed by ligation of Illumina paired-end oligo adapters to the cDNA fragment. Ligated DNA was PCR amplified for 15 cycles and purified using AMPure XP beads. After purification of the PCR products with AMPure XP beads, the quality and quantity of the resulting transcriptome libraries were analyzed using an Agilent Bioanalyzer High Sensitivity chip.

Gene expression in RPKM inferred from RNAseq data

All transcriptome libraries were sequenced with the Illumina HiSeq in paired end mode. The reads that passed the chastity filter of Illumina BaseCall software were used for subsequent analysis. We trimmed all reads to 85-mers and aligned them to the reference human genome (NCBI build 36.1, hg18), plus a splice junction set including 84 bases on either side of the 2008 Illumina splice junction set, using Bowtie [22] in single read mode keeping unique best hits and allowing up to two mismatched bases. Mate pairs from paired end runs were pooled and treated as single reads.

Next, all of the exons for a single gene were concatenated to form a single "transcript" for that gene. Reads that mapped to the exons in the gene, as well as reads that mapped to the splice junctions, were remapped to the "transcript". We then walked the transcript and summed the coverage at each position, then divided the result by the transcript length times the number of reads in the sample, and then multiplied the result by one million. This method is a modified version of the one described in [25].

DNA fingerprinting using common high frequency variants

We considered 147 genomic positions that tend to have good coverage in both whole-exome and RNAseq data and for which at least two alleles are widespread. For each sample, we constructed a "fingerprint" for those positions with at least six fragments of coverage; an allele is included in the fingerprint if it is seen in at least two fragments of coverage and at least 15% of the total fragments. Two samples are deemed a match if at least 85% of the positions with sufficient coverage in both samples have identical fingerprints. In practice, samples from the same patient match at greater than or equal to 90% of positions, while samples from different patients match at less than or equal to 75% of positions. For this study, all the included samples from the same patient met the greater than or equal to 85% match criteria, including all RNAseq and exome samples. In addition, no samples from different patients exceeded 74% similarity.

RNAseq of panobinostat-treated cultures

RNAseq was performed as above in SU-DIPG-IV and SU-DIPG-VI cell lines following exposure to 1 µM panobinostat (Selleck Chemicals) for 48 hours. RNAseq was also performed in SU-DIPG-IV, SU-DIPG-VI and JHH-DIPG1 cell lines following exposure to 0.1 µM panobinostat (Selleck Chemicals) for 24 hours in triplicate. Cells were then lysed and RNA extracted in TRIzol solution.

Post-panobinostat RNAseq data analysis: There was a good concordance for the ratio of panobinostat-treated to vehicle-treated cell gene expression and the distribution of RPKM values were similar in both cell lines and with the two conditions. In general, the variability in terms of ratio was low, usually below 1.

GSEA

RNA-seq profile was pre-processed by applying the following filters: minimum RPKM≥0.5 and maximum RPKM≥10. Gene set enrichment analysis (GSEA) was then performed as previously described (www.broadinstitute.org/software/gsea). MSigDB-C2 (Curated gene sets) database and MSigDB-C5 (GO gene sets) database were applied to gene set size filters (min=15, max=500). After that, 2557 gene sets from the MSigDB-C2 (curated gene sets) database and 586 gene sets from the MSigDB-C5 (Go gene sets) database were used in GSEA. Five additional H3.3K27M-related gene sets were also used. The five H3.3K27M-related gene sets were CHAN_NSC_ONLY (genes with H3K27me3 peaks only in NSC line but not SF7761 H3.3K27M DIPG line), CHAN_SF7761_ONLY (genes with H3K27me3 peaks only in SF7761 H3.3K27M DIPG line but not NSC line), CHAN_NSC SF7761_OVERLAP (genes with overlapping H3K27me3 peaks in NSC line and SF7761 H3.3K27M DIPG line), BENDER_K27M_UP (significantly upregulated genes in K27M pHGG tumors compared to WT pHGG tumors) and BENDER_K27_DOWN (significantly downregulated genes in K27M pHGG tumors compared to WT pHGG tumors). Enriched gene sets were selected using cut-off of FDR ≤ 0.1. Then Venn diagram analysis was used to identify the shared enriched gene sets of the four RNA-seq data sets (SU-DIPG-IV-DMSO or Pano 0.1 µM, SU-DIPG-VI-DMSO or Pano 0.1 µM, SU-DIPG-IV/IV-DMSO or Pano 1 µM and JHH-DIPG1-DMSO or Pano 0.1 µM).

Predictive analysis to identify relevant targets for therapeutic intervention

In this section, we discuss methodology and present results for a computational model of drug sensitivity designed to identify relevant functional targets in DIPG. We utilized the Target Inhibition Map (TIM) approach [12,13] to generate multivariate blocks of targets whose joint inhibition is predicted to increase drug sensitivity. This approach incorporates only functional drug screen sensitivity values with available RNAseq expression data to identify key functional targets in the disease cohort.

Data Pre-processing and Notations

For this analysis, we considered a total of 6 samples that had both drug sensitivity and RNA expression data. A subset consisting of 8 cell lines was tested on Drug Screen 1 (SU-DIPG-IV, SU-DIPG-VI, VU-DIPG.A, NEM-145, NEM157) and another subset of 7 cell lines was tested on Drug Screen 2 (NEM-163, NEM-165, NEM-168, NEM-175, JHH-DIPG1, NEM-157, NEM-215) with 2 cell lines tested on both screens.

To identify the mechanistic targets through which drugs derive their sensitivity, prior information on the kinase targets of the compounds are utilized. The Drug target inhibition information is contained in a matrix Δ which quantifies the ability of a drug to inhibit its kinase targets. Each row vector is relatively sparse, as most drugs inhibit a select few primary targets and a small set of secondary, but often relevant, side targets. Each entry of Δ is a real number between 0 and 1 denoting the inhibition of the kinase target. An entry of Δ is close to 1 when a small concentration of a drug can phosphorylate 50% of the kinase target. The matrix Δ is generated based on published drug target Kd’s or EC50’s from http://pubchem.ncbi.nlm.nih.gov/ and related publications. We used a hill curve interpolation with hill coefficient of 1 and drug concentration of 10,000nM to convert the EC50’s to normalized inhibitions between 0 and 1. For each drug tested over a cell line, the functional response is measured as cell viability, the percentage of living cancerous cells following 72 hr application of a drug. The cellular death rate for each drug, considered as the percentage of cancerous cells killed following 72 hr application of the drug, is simply (1 - cell viability). The cellular death rate values are tested with 3 replicates in 4 concentrations. The cellular death rate values are converted to IC50’s via hill curve interpolation. The IC50 values are then converted to sensitivities between 0 and 1 using a hill curve with hill coefficient of 1 and drug concentration of 10,000nM. A NaN or NA IC50 is converted to sensitivity 0. This set of sensitivities is a matrix denoted by Σ. Note that we have a matrix Σ of size 60×8 for drug screen 1 and a matrix Σ of size 60×7 for drug screen 2.

To incorporate the mutation characterization, we consider the set of kinase targets present in the drug screen (denoted by Ѳ1) and created a binary matrix M0 which denotes whether each of the drug target is mutated or not for each cell line. Considering the common 404 targets of the two drug screens, the size of M0 is 15×404 where 15 denotes the cell lines (two rows are repeated as there are unique 10 cell lines).

To incorporate RNA expression data generated from the RNAseq experiments, we consider the set of kinase targets present in the drug screen (denoted Ѳ1) and the targets with quantified expression present in the RNAseq data (denoted Ѳ2). The set of targets for which we have usable information is then Ѳ1∩Ѳ2, the intersection of the two target sets. The remaining targets in Ѳ1 are targets for which we do not have any RNA expression information and thus we are only able to gain information on these kinase targets from the drug screen data. Let the RNA expression be given by matrix G0. In this analysis, the set Ѳ1∩Ѳ2 has 403 targets and thus the size of G0 is 15 × 403. The RNA expression matrix is converted to a normalized matrix G by dividing the RNA expression of each gene by the mean expression of the matched normal samples or mean expression values. Thus a value ≫1 in G denotes that the gene is highly expressed for that sample compared to normal matched samples or the expression mean.

TIM predictive algorithm to identify multivariate target combinations for therapeutic intervention

In this section, we utilize the Target Inhibition Map approach for target selection and sensitivity prediction. We next provide a brief outline of the approach (for details, please refer to [13]). The primary goal of the TIM model is to generate hypotheses of set of targets that can potentially produce desirable treatment results based on the available drug screen and gene expression data. There are two main steps in the TIM method: step one is selection of the optimal drug targets based on functional data generated from drug screens. This is the model generation step in the TIM algorithm. The model generation is approached from a Boolean logic perspective, where we aim to generate logical relationships between targets that identify the mechanisms by which tumors persist. Due to the Boolean nature of the approach, the matrix Δ generated for this analysis consists of 0’s or 1’s. The drug inhibition profiles (in the form of target EC50’s) are used to generate the binarized inhibition vectors based on the drug’s IC50 value. The binarized vectors are generated based on the following equation: a log (IC50) ≤ log (EC50) ≤ b log (IC50), 0 ≤ a ≪ b This binarization approach uses the functional drug sensitivity data to generate vectors of relevant drug targets. We assume that any target sufficiently beyond the IC50 point is unlikely to have been the cause of the drug response, and thus we assume that we cannot gain information from it.

To generate drug sensitivity values from the IC50 response data, we use the following equation: where IC50, is the IC50,i of drug i, MaxDose is the maximum tested dosage of drug i, and y is the resulting sensitivity value of drug i.

To incorporate RNAseq data into the TIM analysis, we eliminate unlikely kinase targets from consideration based on G, the matrix of normalized RNA expressions over the set Ѳ1∩Ѳ2. A kinase target is kept in the dataset for consideration if G(i,j) > ρ, a threshold for minimum relative expression. Note that we are considering expressed targets only as the effect of a molecularly-targeted drug is to inhibit the target when it is expressed and thus, non-expressed drug targets will have limited effect on predicting the drug response. For the purpose of this project, the RNA expressions constraint was G(i,j) > ρ = 1 meaning the RNA expression of the cancerous line must be at least that of the matched normal expression. As G does not provide information about the remaining targets in Ѳ1 − Ѳ2, they are not eliminated from consideration as they could potentially be driver mechanisms for the tumor. We performed the primary analysis on the 16 sets of joint RNAseq-matched drug screen data and cell line drug screen data.

Subset-Superset Biological Constraint

We utilize the following biological constraints to identify consistent and inconsistent sets of protein targets [13].

For any drug d and any kinase target set S, Here, the ≼,≽ operators signify component-wise inequality in all components.

Constraint 1: If Δ(d1, S) ≼ Δ(d2, S), then Σ (i, d1) ≤ Σ (i, d2) (more inhibition of oncogenes should improve sensitivity)

Constraint 2: If Δ(d1, S) ≽ Δ(d2, S), then Σ (i, d1) ≥ Σ (i, d2) (less inhibition of oncogenes will reduce sensitivity)

The vectors and associated sensitivity values from matrix Σ are used to select the optimal targets for the TIM model. The optimal targets are selected by solving the following optimization problem that incorporates the inter-bin error (the error associated with multiple drugs having the same inhibition profile but differences in their sensitivity values) and the inconsistency error (the error associated with drug A being less effective than drug B but the targets of drug B are a subset of the drug A’s targets). In this context, a bin is any subset of a set of kinase targets selected to build the model. The equation for inter-bin error is given by ∑∑|P(s|T)–Y(s|T)|, where P(•) is expected sensitivity of s under target set T, and Y(•) is the experimental sensitivity of s under target set T. The inter-bin error amounts to ensuring that kinase inhibitors with similar inhibition profiles with respect to the model have similar experimental sensitivities, indicating the likely mechanisms have been identified. The inconsistency error incorporates the subset-superset biological constraint introduced previously and is given by the following equation: ∑∑χ(bin, drug)|P(s|T)|-Y(s|T)|, where P(•) and Y(•) are as above and X(•) is an indicator function which is 1 when the subset-superset constraint is violated. The total optimization problem we wish to solve then is ∑∑|P(s|T)-Y(s|T)| + ∑∑χ(bin, drug)|P(s|T)|-Y(s|T)|. We find the minimum solution to the optimization problem using a suboptimal search algorithm known as Sequential Floating Forward Search (SFFS). The TIM framework requires that the subset-superset biological constraints are satisfied during model generation. This allows for identification of unique target combinations that may not be identified as relevant through a standard linear modeling approach.

Step two in TIM analysis is sensitivity prediction based on selected targets from step one. This step is used to generate the TIM circuit, which provides a visual representation of effective treatments inferred from the functional data. To generate the target blocks, viable target combinations are explored systematically in a breadth-first search style approach. Here, desired sensitivity levels for each level of inhibition (in terms of targets inhibited) is specified by the user; based on this specification, combinations of targets that satisfy the desired sensitivity levels are selected and no further supersets of the combination are considered. The inferential steps are based on bounding the expected sensitivity by the closest subset of the target combination and the closest superset of the target combination. The equation for inferred sensitivity for a target set {t is given by , where y is the maximum sensitivity of the subsets of {t, y is the minimum sensitivity of the supersets of {t, α is weight for a kinase target, and β determines the order of the curve fit for y({t1,t2,…t}). For this analysis, for all targets meaning all targets have equal weight and α = 1 which provides a linear fit. The sum in the denominator is the sum of weights between the subset generating y and the superset generating y. In this way, the predicted sensitivity lies between y and y depending on the distance of {t1,t2,…t} to the subset and superset. Based on the number of additional targets needed to move from the subset to the desired target combination, the sensitivity of the drug combination is inferred. In the event, the predicted sensitivity y({t1,t2,…t}) is sufficient to satisfy the sensitivity level defined by the user, the target combination is designated a valid block combination, and the supersets are no longer considered.

In this supplemental material, Figure S3b is the TIM model generated for the set of available DIPG drug screen data. Note that TIM analysis generates potential targets from drug screen information and uses the gene expression to remove false positives. Since the drugs considered inhibit oncogenes, we keep the targets that are differentially expressed in tumor cell lines as compared to normal cell lines. If a TIM analysis target is differentially expressed in at least one tumor cell line, the target is considered relevant. The number of genes displayed in Figure S3b belongs to blocks whose scores[13] are above a threshold of 0.5.

Figure S4 contains the results of in vitro drug combination testing and synergy analysis for a target combination predicted by the TIM modeling approach to have in vitro synergy.

The Chou-Talalay Combination Index (CI) was used to detect synergy (CI<1). The CI is calculated based on the following equation where D1 and D2 denotes the combination dosages tested for drugs 1 and 2 to achieve δ% inhibition and D and D denote the corresponding dosages of individuals drugs 1 and 2 required to inhibit δ%.

To provide quantification for the predictive power of the TIM approach, 10-fold cross validation mean absolute errors (MAE) and correlation coefficients between experimental and predictive sensitivities were calculated. The first test used only the functional data generated from the Drug Screens for prediction. For this analysis the 10-fold CV error was 0.101 and a correlation coefficient of 0.848. The second test incorporated the RNA expression data to eliminate potential false positive before model development. Here the 10-fold CV error was 0.111 and correlation coefficient of 0.773.

Statistical analyses

Bioinformatics and computational methods are described above. IC50 values were determined by a nonlinear best-fit method using Excel Solver. For comparison of gene expression levels of the 429 genes associated with the targets of the drug screen, we generated RNAseq data and then compared tumor tissues to normal tissues to generate significance scores using a t-test with the Benjamini-Hochberg multiple comparison correction.

Sample size for in vivo experiments was based on variance estimated from pilot experiments with control pontine DIPG xenograft animals showing a standard deviation of approximately 30% of the mean. Sample size calculations to detect a 50% difference in tumor growth between two groups with an alpha of 0.05 and power 0.8 indicate a minimum group size of 3 animals.

Two-tailed t tests and two-way ANOVAs were used for in vitro experiments. A two-tailed t-test was used for comparison of tumor growth rate in treatment groups in orthotopic xenograft experiments; log-rank test was used for orthotopic xenograft survival analysis.

Outlying data points (>2 standard deviations from the mean in either direction) were excluded. Variance was similar between groups compared.