Rapid identification of optimal drug combinations for personalized cancer therapy using microfluidics

Functional screening of live patient cancer cells holds great potential for personalized medicine and allows to overcome the limited translatability of results from existing in-vitro and ex-vivo screening models. Here we present a plug-based microfluidics approach enabling the testing of drug combinations directly on cancer cells from patient biopsies. The entire procedure takes less than 48 hours after surgery and does not require ex vivo cultivation. We screened more than 1100 samples for different primary human tumors (each with 56 conditions and at least 20 replicates), and obtained highly specific sensitivity profiles. This approach allowed us to derive optimal treatment options which we further validated in two different pancreatic cancer cell lines. This workflow should pave the way for rapid determination of optimal personalized cancer therapies at assay costs of less than US$ 150 per patient.


Introduction
Increasing interest has been devoted to personalized (precision) medicine, which aims to match the best treatment to each patient. This is especially important for cancer, where there is often a high variability in the response of patients to treatments. A critical example is pancreatic cancer, which has a 5--year survival rate of only 6% and is projected to become the second leading cause of cancer death by 2030 1 . One of the reasons for the low efficacy of current treatments is that pancreatic cancer is genetically highly heterogeneous with most mutations occurring at a prevalence of less than 10% 2,3 . Therefore, pancreatic cancer patients would benefit immensely from the possibility of a personalized or stratified therapeutic approach.
Most efforts in personalized medicine have been focusing on tailoring the treatment to the specific patient based on genomic data, which are increasingly available due to advances in sequencing technologies. While there have been some impressively successful examples 4,5 , cancer genomics is generally very complex and, despite the increasing knowledge on occurring mutations, there is still limited understanding on how they affect drug response 6 . Multiple efforts have been devoted to the large--scale in vitro screening of drugs across cell lines 7-9 that have proven useful to identify some genomic markers associated with drug response.
However, molecular data alone has not proven sufficient to predict the efficacy 10 or toxicity 11 of a drug on an individual cell line in a reliable way. This predictability is likely to be even lower in patients, given the additional complexities when compared to cell lines.
Due to these limitations, genomics data has to be supplemented with other information in order to optimally guide the treatment for each patient, and systems for phenotypic stratification are urgently needed 6,12 . The need for new approaches is even more acute for the application of drug combinations.
Combinatorial targeted therapy has been shown to be a powerful tool to overcome drug resistance mechanisms, which can be due to tumor heterogeneity, clonal selection or adaptive feedback loops 13 , and seems to be a particularly promising 3 approach for treatment of pancreatic cancer 6 . However, strategies to identify effective combinations are still in their infancy 14 .
A powerful means to overcome these difficulties would be the ability to test the drug compound directly on patient samples. Despite recent progress towards functional testing of live patient tumor cells 12,15 , the application of standard drug screening technologies is currently limited by the need of large numbers of cells 15 .
Therefore, large scale drug screening of patient tumors has been so far limited to blood tumors 16 (where a much larger amount of malignant cells is easily accessible) or requires some ex vivo culturing steps 6 (e.g. patient derived cell lines, PDX models and organoids 17,18 ) that require long times to grow the cells and can cause changes in the phenotype of the cells. Hence, there is, hence, no technology to perform drug screenings with the number of cells that can be obtained with biopsies from a solid tumor.
Microfluidic technology can in principle overcome screening restrictions based on limited starting material. Making use of tiny assay volumes, microfluidic systems have recently been applied successfully to the testing of a few individual drugs on primary tumour cells, tumour spheroids and tissue slices [19][20][21] . However, these studies were based on single aqueous phase microfluidic systems which can process only small sample numbers (max ~96 including replicates, typically much less). A possible solution for further scale--up is the use of droplet microfluidics 22 .
In these systems aqueous droplets surrounded by oil serve as independent reaction vessels. Such systems have already been used for genetic assays of cancer cells 23,24 , but they have not yet been applied to personalized phenotypic drug screens. This is probably due to the fact that encapsulation of different soluble drugs (rather than just different cells) into droplets or plugs (sequential aqueous segments in a microfluidic channel, capillary or piece of tubing; spaced out by oil or air) is technically still very challenging. Switching between multiple fluid sources requires robotic systems (e.g. sequentially aspirating samples from microtiter plates), which are rather slow, or complex microvalve technology (e.g. to switch between fluids injected in parallel) 22 . While a combination of both approaches has been described 25 ; to date none of these plug--based systems has 4 been used for phenotypic assays with human cells. This may be owing to incompatibilities, e.g. due to the fact that mammalian cells sediment and secrete proteins and metabolites causing the plugs to stick and break at the channel walls (so--called "wetting"; also promoted by growth factors in the media). In addition, a robust method for sample tracking is required. In previous studies [25][26][27] plugs were simply stored in a sequential fashion (within a microfluidic channel or a piece of tubing), wherein the common fusion or splitting of plugs causes "frameshifts" resulting in loss of information on sample composition.
Here, we present a platform that can overcome these limitations, enabling the screening of drug combinations on mammalian cells and patient biopsies in a plug format. Our approach requires significantly less cells compared to conventional, non--microfluidic formats and provides one to two orders of magnitude higher throughput (in terms of samples per experiment) than existing systems. Our platform, based on microfluidic Braille valves 28,29 and an external autosampler, can rapidly generate plugs containing cells, reagents of an apoptosis assay and systematic combinatorial drug cocktails. All workflows have been optimized to guarantee compatibility with live cells, and we furthermore implemented a fully scalable sample barcoding system. We demonstrate application of this platform for screening cell lines and, more importantly, patient samples ex vivo. Furthermore, pathway modeling of this data allows us to define apoptotic pathways in a patient-specific way. This should open the way for further translational efforts and clinical applications in the near future.

Microfluidic platform for drug screening
We developed a plug--based microfluidic platform allowing us to rapidly screen a large number of combinations of chemical compounds on a very limited number of cells. Specifically, with a starting material of less than 1 million viable cells more than 1100 samples could be screened (10 chemical compounds in all the pairwise combinations, with about 20 plugs (i.e. replicates) per condition), each containing about 100 cells (Fig. 1a).
Our platform is based on Braille valves 28,29 controlling individual fluid inlets of the microfluidic chip ( Fig. 1b--c). All reagents (including the cell suspension, drugs and assay components) are permanently injected into the device and, depending on the valve configuration, sent either to the waste or to a droplet maker with a T-junction geometry. This approach allows rapid switching (~300 ms) between 16 liquid streams and, upon injection of fluorinated oil at the T--junction, the generation of combinatorial plugs at high throughput. If necessary, the chemical diversity can be increased further by connecting an autosampler to one of the inlets (sequentially loading compounds from microtiter plates; Supplementary   Fig. 4), but this procedure also requires a higher number of cells.
To avoid wetting issues we integrated an additional mineral oil inlet into our chip design, enabling the insertion of mineral oil droplets in between all (aqueous) sample plugs (Supplementary Movie 3). Such three--phase systems are particularly efficient in keeping samples separated 30 , even under conditions that normally cause wetting. A further crucial factor for reducing wetting was the use of special, protein--free media. However, while a combination of these measures significantly improved reliability, plug integrity could still not be ensured for all samples. Hence a sample identification system that would still work if individual plugs break or fuse had to be implemented. Therefore, we introduced sequential barcodes to separate and identify samples of different composition, based on sequences of plugs with binary (high/low) concentrations of the blue fluorescent dye cascade blue. This way sample numbers can be encoded (e.g. high--low = binary number "10" = decimal number "2"). To demonstrate the power of this barcoding approach we converted a simplified EMBL logo (Supplementary Fig. 3) into a binary black and white image and translated all 2808 pixels into a sequence of plugs with two different fluorescence intensities. These barcodes were then detected using our laser spectroscopy setup (Supplementary Fig. 1) and converted back into the initial image, which did not show a single mistake. This clearly demonstrates the scalability and reliability of our strategy.
In each sample plug, cells are encapsulated together with one or two compounds and a rhodamine 110 (green--fluorescent dye) conjugated substrate of Caspase--3, which is an early marker of apoptosis 31 . Furthermore Alexa fluor 594 (orange-fluorescent dye) was added to the cell suspension for verifying dilution by all assays reagents and monitoring correct operation of all valves. For each drug treatment multiple replicates were generated (see screening specific details below), followed by a fluorescent barcode, and stored sequentially in gas permeable tubing with a total length of 5m (having a capacity for about 3500 plugs in total). After overnight incubation, the readout was performed by flushing the samples through the detection module. In our case, we used an optical setup (Supplementary Fig. 1) with three different excitation lasers (375nm, 488nm and 561nm) and performed a multiplexed readout at three different wavelengths (450nm, fluorescence barcodes; 521nm, Caspase--3 activity; >580nm, orange marker dye to monitor mixing of reagents) using photomultiplier tubes (PMTs).

Data extraction and quality assessment
In the fluorescence data, each plug corresponds to a peak in one or more channels (i.e. green, orange, blue), as shown in Fig.  1d. When processing the data, we first identified the peaks in the blue channel, representing the barcode, and we used them to separate and identify the peaks corresponding to each sample (Online Methods). Each sample is composed of multiple peaks (typically 12 replicates per run) with signals both in the green and in the orange channel. For each peak we considered two measures: height and width. The height of the peak is proportional to the measured fluorescence intensity: the intensity in the green channel represents the activation of Caspase--3 (thus apoptosis) while the intensity in the orange channel represents the concentration of the orange--fluorescent dye, which was added to the cell suspension. Since each plug was produced by mixing four components (cells, Caspase--3 substrate and two compounds/medium), the intensity in the orange channel allowed assessing the quality of this mixture and was used to discard samples/peaks with extreme values (i.e. outliers, see Online Methods). The width of the peak represents the length of the plug. Thin or wide peaks were discarded as they correspond to split or fused plugs respectively.
Additionally, the first peak for each sample was typically discarded to avoid the effect of cross--contamination between samples.

Combinatorial screening of pancreatic cancer cell lines
To optimize and validate our microfluidics platform we performed combinatorial screening of compounds on two pancreatic cancer cell lines with different genotype and phenotype 32 : AsPC1 and BxPC3. Ten compounds were screened alone and in pairwise combinations (Table  1) (Fig. 2b,c).

Efficacious drug combinations and validation for screened cell lines
Combinatorial screening using the microfluidics platform suggested potentially interesting drug combinations (Fig. 2b,c). Some combinations showed strong effects in both cell lines (e.g. GDC0941 and Cyt387, z--score = 0.46 for BxPC3 and 0.67 for AsPC1), while, more interestingly, some stronger effects were shown to be  (Fig. 3b,d), confirming the behavior observed in the microfluidic system.

Combinatorial screening of resected patient pancreatic tumors
The same set of ten compounds used on cell lines was then applied to screen tumor samples resected from 4 pancreatic cancer patients. For each patient sample, the solid tumor was dissociated to create a single cell suspension that was then used for the screening (protocol in Online Methods). Similar to the pipeline used to screen the cell lines, a total of 62 samples were produced, with 12 replicates each (each plug with about 100 live cells). The whole sample sequence was repeated two or three times (different runs) depending on the amount of viable cells that were obtained from the patient biopsy. As in the case of the cell lines, unreliable peaks/samples were discarded based on the intensity in the orange channel and on the length of the peaks. Only samples passing the quality assessment in at least two runs and with consistent values across runs were considered for further analysis (boxplot with data for each run shown in Supplementary Fig. 7, more details in Online Methods). Interestingly the best combination (highest z--score) was different for each patient:
Gemcitabine and Oxaliplatin, the current first--line treatments for pancreatic cancer showed increased apoptosis induction when administered in combination with each other (in patient #3, z--score equal 1.31) or with other drugs. In particular, Gemcitabine strongly induced apoptosis in combination with Gefitinib for patient #1 (z--score equal 1.58) and with MK--2206 for patient #4 (z--score equal 1.06), while Oxaliplatin is effective in combination with AZD6244 or GDC0941 for patient #3 (z--score equal 2.12 and 1.02 respectively) and with Cyt387 for patient #4 (z-score equal 1.05).

Comparison between patient samples and cell lines
Clustering the screening data by individuals (patient/cell line) and by samples (Fig   4e) showed how most of the single compounds and some combinations, e.g.
AZD6244 and TNFα, failed to show strong efficacy on any individuals. Other treatments displayed high efficacy both on patient samples and on cell lines (e.g. combination of PHT--427 and ACHP). Importantly, none of the treatments is effective across all individuals, suggesting that they are not widely toxic, but rather patient specific. This further demonstrates the need for individualized treatments.

Network--based interpretation
In order to better understand the mechanism of action of the screened compounds, we derived a logic model of pathways involved in apoptosis (Fig. 5a)  Cell line and patient specific models were built using a logic ordinary differential equations formalism 38 as implemented in CellNOpt 39 . The starting general model structure was fitted to the experimental data for each cell line and patient (see Online Methods). Estimated parameters are one life--time parameter for each species (e.g τEGFR for node representing EGFR) and a regulation parameter for each interaction (e.g. EGFR → JAK). For each cell line/patient a bootstrap distribution was obtained for each parameter (resampling experimental data with replacement 300 times).

11
In order to identify the pathway differences that might underlie the differential response to the drugs in our two tested cell lines, we compared the distributions of the parameters using Wilcoxon rank sum test for assessing the significance (p-values Bonferroni corrected for multiple hypothesis testing > 0.05) and r effect size (medium--high effect size > 0.3, see Online Methods); results are shown in (Fig.   5b). Similarly, distributions across all cell lines/patients were compared using Interestingly, it has been previously shown that alterations in the activity of the AKT pathway are quite common in pancreatic adenocarcinoma due to mutation and epigenetic alterations 40,41 .

Discussion
We describe here a microfluidic platform enabling the fast screening of many drug combinations in a personalized approach. The very small assay volumes (~100 cells per 0.5 µL plug) allowed us to perform comprehensive screens directly on patient biopsies. Importantly, these can be done without the need for any intermediate cultivation steps, which might introduce significant cellular changes and/or the selection of individual clones. Hence our approach opens the way for comprehensive screens that have so far been restricted to blood tumours, for which patient--derived cancer cells are available in large quantities 16 .
Furthermore, the low volume and the high level of automation allow such screens to be rapidly performed at low cost: Within 48h after tumour resection or biopsy, and with consumable costs of about US$ 150, which is considerably lower than routine procedures such as MRT scans or surgical interventions. Therefore, we envision rapid translation of this technology into clinical application.
Functional combinatorial screening has great potential in predicting personalized therapy as revealed by the data presented in this paper. Our cell line screening recapitulates the lower sensitivity of KRAS mutants (AsPC1, vs wild type BxPC1) to PHT--427 42 , (Fig. 1 b-- Perturbation experiments with targeted drugs (e.g. kinase--specific inhibitors) can also be used to infer the cell line or patient specific pathways 45 . Signal transduction models can also be linked to cell viability/apoptosis 46 in order to predict drug sensitivity. We used data from our perturbation experiments and prior knowledge on the signaling pathways impinging on apoptosis to build cell line and patient specific models which were then used to unveil potential biological reason for the differential response. Models could, in principle, be in the future also used to predict new potentially interesting targets. Compared to other personalized approaches our approach has specific advantages: While organoids and xenografts are particularly well--suited for mimicking three-dimensional tumour architecture and the in vivo microenvironment, the use of individualized tumour cells as shown here facilitates rapid, massively parallelized assays at low cost. This could be exploited further by implementing single--cell droplet assays 22,47,48 , which require fewer cells and also allow to take intratumoral heterogeneity into account. In addition to deriving an averaged readout over 50 single--cell replicates per treatment option, one could quantitatively determine the number of non--responding cells for each drug cocktail. Such data would probably be very valuable to overcome therapy resistances.
Apart from using it for personalized cancer therapies, our microfluidic platform should be of interest for further applications, such as the screening of cocktails for targeted stem cell differentiation 49,50 or combinatorial chemistry 51  For screening, the Braille valve chip (Fig. 1B) was mounted onto an SC--9 Braille display (KGS Corporation, Japan) using an in--house holder as shown in Fig. 1C.  Supplementary Fig. 1. The resulting sample data was analyzed using custom R--scripts (see paragraph "data processing" for details).

Setup of the fluidic system, choice of additives and oils
We used fluorinated oil without stabilizing surfactant as a carrier phase (only 0.5% of the anti--wetting agent PFO was added), which turned out to have two major advantages: upon reaching the outlet, small droplets generated at the T--junction fused and formed larger plugs that completely filled the collection tubing, thus allowing to incubate all samples in a sequential order. Furthermore the lack of surfactant prevented the formation of micelles, which can cause the exchange of reagents between droplets 52 . To increase stability of the arrays, plugs were furthermore spaced out using mineral oil (Sigma).

Integration of an autosampler
To allow for upscaling of the screens, we integrated an autosampler into our microfluidic platform 27 . One of the inlets of the microfluidics chip was connected to a Dionex 3000SL Autosampler, aspirating samples from 96--well plates and injecting them sequentially into a target tubing connected to an external Harvard Apparatus Syringe pump. While the resulting throughput for loading compounds from different wells is rather low (~90 sec per reagent), each of them can be mixed with all of the drugs injected directly into the Braille valve chip. Therefore, the maximal throughput in terms of pairwise combinations is much higher (e.g. 9 sec per combination when mixing with 10 further reagents directly connected to the Braille display chip). However, one effect had to be taken into account: the concentration of compounds coming from the autosampler varied due to dispersion of the samples in the miscible carrier phase (the buffer used in the robotic system) 53 . To overcome this effect, we implemented a feedback loop between the autosampler and the Braille valves: The relay signal of the autosampler was used to send the beginning and end of each sample from the microtiter plates to the waste, using the two Braille valves controlling the relevant inlet on the microfluidic chip. This process allowed to overcome sample dispersion and is illustrated in Supplementary  Fig.  4. This way only the centre part of each

Preparation of cells for microfluidic experiments
Cultures were trypsinized to detach cells, harvested and washed with PBS. The same procedure was used to prepare the single cell suspension from either Additionally, we also discarded peaks showing very high and very low width (based on empirical thresholds) in order to remove peaks corresponding to fused or split plugs respectively. After these steps, we considered the distribution of the intensity of the orange peaks across all samples and discarded the extreme values (i.e. the outliers). Where Q1 is the 25 th percentile, Q3 is the 75 th percentile and IQR is the interquantile range (Q3--Q1), outliers were defined as values lower than Q1 --1.5*IQR, or higher than Q3+1.5*IQR. These strict rules were applied to guarantee higher quality of the data used for the analysis described in the paper. Code used to process the data is available as an R package in GitHub (https://github.com/eduati/BraDiPluS).

Apoptosis pathways model
The logic model shown in Fig. 5 was derived by manual literature curation starting anti--apoptotic proteins, the latter binds to Apaf1 (Apoptotic protease activating factor--1) and pro--caspase9 which is converted to its active form of Caspase--9 (Cas9 in Fig. 5) and in turn activates Caspase--3 (Cas3 in Fig. 5). Both Akt and ERK have an anti--apoptotic effect by phosphorylating BAD 36 and thus unbinding it from BclX and this can be modelled as an OR gate 55 . We also included the pro-apoptotic effect of ERK as regulator of p53 37 . Additional cross--talks between the pathways (i.e. RAS → MEKK1 and RAS → PI3K) were included as described in 35 . The resulting logic model described in Fig. 5 was considered as a Prior Knowledge Network (PKN) and was then optimised based on the experimental data for each patient/cell line in order to generate patient/cell line specific models. We used a formalism based on logic ordinary differential equations 38 where ordinary differential equations are derived from logic models using a continuous update function for each species , which can assume continuous values : Where are the N regulators of and each regulation is defined by Hill functions with parameters and : Optimization was performed using CNRode add--on of the CellNOptR package 39 : parameters and were estimated, while parameters were fixed to 3. The notation used in the paper and in Fig. 5