Local generation and efficient evaluation of numerous drug combinations in a single sample

We develop a method that allows one to test a large number of drug combinations in a single-cell culture sample. We rely on the randomness of drug uptake in individual cells as a tool to create and encode drug treatment regimens. A single sample containing thousands of cells is treated with a combination of fluorescently barcoded drugs. We create independent transient drug gradients across the cell culture sample to produce heterogeneous local drug combinations. After the incubation period, the ensuing phenotype and corresponding drug barcodes for each cell are recorded. We use these data for statistical prediction of the treatment response to the drugs in a macroscopic population of cells. To further application of this technology, we developed a fluorescent barcoding method that does not require any chemical drug(s) modifications. We also developed segmentation-free image analysis capable of handling large optical fields containing thousands of cells in the sample, even in confluent growth condition. The technology necessary to execute our method is readily available in most biological laboratories, does not require robotic or microfluidic devices, and dramatically reduces resource needs and resulting costs of the traditional high-throughput studies.


Introduction
The optimization of targeted therapies is based on the determination of drug dose-response and the selection of a regimen that is best suited for an individual patient. A critical barrier to optimization screening is the availability and amount of relevant tissue sample (e.g. biopsy), especially in the case of combination therapy. This limitation is a reflection of the fact that the number of possible drug combinations grows exponentially with the number of drugs and dosages (cf., Supporting Information). For example, the number of all possible combinations of three drugs tested at ten different concentrations (including zero) is 10 3 . This number does not include the biological replicates and controls that render the careful testing of even three drugs extremely challenging.
There exists another challenging aspect to the identification of optimal drug dosage and regimen: accurate reproducibility of in vivo and in vitro experiments is notoriously difficult to achieve in many biological applications, such as drug dose-response quantification. Many confounding factors contribute to the variability in sample readout Balázsi et al., 2011;Niepel et al., 2019; these factors can generally be grouped into two categories, which are usually described as intrinsic and extrinsic noise.
Originally, this noise nomenclature was used to describe temporal fluctuations in a single protein within a single cell (Elowitz et al., 2002;Raj and van Oudenaarden, 2008). Here, we adopted this classification to an autonomous system which is an ensemble of cells isolated from other samples, but still subject to fluctuations in the external environment. Instead of measuring temporal fluctuations in individual cells, one can record snapshot(s) of cellular phenotype(s) within a single sample. The source of noise in these measurements is variability between individual cells within the sample; we refer to this type of noise as intrinsic. Noise that is driven by external factors to this autonomous system (a particular sample) is referred to as extrinsic. Variability between different samples due to temperature, humidity, etc., are examples of extrinsic noise in this context.
Sources of intrinsic noise in drug uptake by and action within cells include fundamental biochemical processes, such as passive uptake, active transport, and degradation. By contrast, external factors that contribute to noise are principally manifest through convection and advection. Namely, these external factors are related to movement and gradient formation of fluid or gas. The initial cell seeding that eventually leads to unique microenvironment formation can also be thought as an example of extrinsic noise using this terminology.
While intrinsic cell variability can be significant, we believe that it is the extrinsic factor(s) that drive sample variability in most experimental cellular systems. This conclusion derives from the Law of Large Numbers (Feller, 1968), since the typical number of cells in a culture sample is at least of the order of 10 4 , while the number of biological/technical replicates is ∼3.
If the observed sample variability is due to extrinsic noise then it is very tempting to conduct all drug response experiments in the same tissue or cell culture sample. Indeed, under these conditions, it is guaranteed that external fluctuations in the environment are exactly the same for all testing conditions. Hence, quantification of the drug's response does not require averaging over many different samples to achieve an accurate, relative comparison of the drug's effectiveness.
In order to test different drug combinations in the same cell culture or tissue specimen, one must create these combinations locally. One of the smallest possible quantifiable locales in both cell culture and tissue samples is a single cell. Indeed, we can treat each cell as an individual chemical "reactor" in which testing of different drug combinations is performed. Thus, if we can (i) create different drug combinations in individual cells, and (ii) subsequently quantify these combinations and the ensuing outcome (phenotype) for each cell, the logistics of testing combinations in a single macroscopic sample would be effectively addressed.
There exists a straightforward way to resolve the quantification issue (ii) above. Drugs can be uniquely labeled (e.g. using fluorescent tags), in which case, the amount of drug taken up by each cell would scale as the fluorescence of its corresponding tag.
Let us now address a more delicate problem (i): how to create a broad range of drug combinations in different individual cells. At first glance, an herculean effort would be required to accomplish this task. Fortunately, cells themselves solve this problem effortlessly for us. Even cells with identical DNA typically exhibit deviations in drug uptake. Therefore, the random nature of drug(s) uptake by individual cells can be utilized to create different drug combinations. Here, we further enhance heterogeneity and independence in stochastic drugs mixing by creating individual drug gradients across the sample.
Rather than delivering drugs homogeneously, we rely on simple forced drug advection (by local injection) to create transient local drug gradients independently for each drug. After delivery and incubation of random local drug combinations to individual cells, we analyze the ensuing cellular phenotypes and corresponding drug optical barcodes by fluorescence microscopy or flow cytometry.
Cell 'barcoding' terminology in biology usually implies DNA identification, but here we use fluorescent signal(s) instead to encode drug regimen. This method was used by Krutzik and Nolan to distinguish between individual samples exposed to different drug regimens for high-throughput analysis (Krutzik and Nolan, 2006). The number of samples in their experiments was equal to the number of drug conditions, since the authors were barcoding samples, not individual cells.
The statistical analysis of the imaging data is used to predict the effect of ratiometric dosing on the outcome (cell phenotype). In this fashion, we require as little as a single cell culture sample to test all feasible drug combinations, bypassing the need for an otherwise enormous number of controls and technical replicates.
In this work, fluorescently labeled drugs are clustered into three 'barcoding' classes distinguished by the physical association of the label with the drug: either irreversibly linked or cleavable/separable chemistry. Additionally, we developed a labeling method that does not require any chemical modification of the drug(s), such as covalent fluorescent tagging. This particular approach involves local co-injection of an unconjugated pair of pre-mixed drug and fluorescent dye that serves as a correlative barcode. We utilize the fact that on short incubation time scales (less than an hour), the initial injection producing forced advection is mainly responsible for transient gradients in the sample. Immediately after local injection, the drug and its tag concentration profiles across the sample correlate. The subsequent diversion of these profiles due to differences in diffusion and transport rates between the drug and its unconjugated tag occurs slowly (di Cagno and Stein, 2019) over this comparatively short incubation period, ensuring persisting correlation between drug and tag uptakes throughout the course of the typical experiment. For all drug barcoding classes we designed and tested, we rely on the fact that the drug barcode is retained by the cells. Hence, total drug exposure over time (area under the curve) can be assessed by a single time point measurement.
We also developed segmentation-free microscopy-based image analysis and statistical data processing capable of handling large optical fields containing thousands of cells in the sample, even in confluent growth conditions. These algorithms allow us to bypass the standard segmentation process, typically requiring additional fluorescent marker(s), significant computational time, and resources.
To demonstrate a range of possible applications of our method using different drug classes, we include implementation of our stochastic drug mixing strategy covering different biological topics. We first consider a benchmark system in which 'drugs' are diester dyes with distinguishable fluorescent properties. To go beyond this 'toy' system, we assessed the effect of a combination of siRNAs on a direct and quantifiable target, viz., GFP fluorescence, as modified by siRNA-dependent expression of GFP protein. This system is characterized by a simple phenotype readout, such as an unique fluorescent signal, and the modulation of the phenotype is a result of the direct drug-target interaction, and not a secondary, downstream effect.
Lastly, we conducted multiple cell culture experiments treating cells with combinations of drugs commonly used in multi-drug chemotherapy regimens, such as cisplatin, docetaxel, cyclophosphamide, and 5-fluorouracil. We adopted our technology to accommodate quantification of cell death as a phenotype common to all of these agents.
Our approach has distinct advantages over conventional drug combination screening methods, including far greater efficiency, a smooth (continuous) variation of drug concentrations in the resulting combinations, and concomitant phenotype readout. We believe this study demonstrates the applicability of this method to practical challenges in combination drug development.

Stochastic drug mixing
Stochastic mixing occurs naturally due to randomness in the drug(s) uptake process by individual cells. We enhance the random process of drugs' uptake by creating transient drug gradients across the sample.
In the simplest possible implementation of deterministic gradient mixing, one can create a drug A gradient along the sample's x-axis and a drug B gradient along its y-axis, generating (locally) drug A and drug B combinations in a continuous fashion. This approach, however, requires microfluidic devices that would maintain stable drug gradients, which can be technically challenging (Wei et al., 2019;Kuo et al., 2019). Moreover, the implementation of multiple (i.e. >2) drug gradients ex vivo is probably not feasible on a small spatial scale. It is also not obvious as to how to prepare efficient combinations of three or more drugs using this approach.
Here, we used a much simpler approach involving local point injection to create spatio-temporal drug gradients in (mammalian) cell culture. We do not homogeneously distribute drugs in the media but, rather, inject them locally into cell culture in separate, discrete locations. Point injection of multiple drugs ensures chaotic/stochastic mixing patterns across sample (Box 1, Figure 1).
One barcoding option is to label the drug directly with a fluorescent tag either by chemical modification of the drug or other irreversible binding methods. Here, we demonstrate the application of our method for irreversibly tagged drugs using cell-retained dyes as a benchmark drug model. The cell-permeable diester dye, carboxy fluorescein succinimidyl ester, and its derivatives are an example of this type of drug.
We first consider a benchmark system in which 'drugs' are diester dyes with distinguishable fluorescent properties. These dyes are well retained by the cell for an extended period and, hence, can be used to characterize total uptake by individual cells. Typical targets of these 'drugs' are intracellular thiol or amine groups (Materials and methods).
In Figure 2a, we demonstrate the typical results of a random drug mixing experiment using a homogeneous delivery method. HeLa cells were exposed to four different dyes mixed together simultaneously in suspension cell culture according to the manufacturer's protocol and imaged using confocal microscopy. The color of each dye (violet, green, yellow, and red) corresponds to its specific fluorescence emission maximum. The color and intensity of each cell in the image correspond to a unique 'drug' combination.
One reason as to why the majority of cells exhibit somewhat similar color is the high degree of similarity in chemical and, hence, transport properties of these compounds. There is a propensity for cells to take up similar molecules proportionally, with only atypical cells (outliers) exhibiting different behavior in this case.
We experimented with multiple strategies for increasing variability, such as decreasing incubation time and cell culture volume, factors that we showed previously drove heterogeneity in drug uptake (Elgart et al., 2018). We found that the most robust way to control variability in uptake is to generate transient drug gradients in culture. A typical outcome of these experiments with HeLa cells is shown Box 1. Barcoding individual cells using stochastic drug mixing One can think of transient drug gradients as local, short incubation processes that facilitate uptake heterogeneity across the sample.
The duration(s) of these local uptake processes are governed by diffusion, advection, and the drug-target reaction kinetics in cell culture.
Box 1-figure 1. The aesthetic visualization of the stochastic mixing method implemented on a muchlarger scale.
Ink In Motion by Macro Room, https://www.youtube.com/embed/ICxC5ekWnUc. We use the fluorescent intensity of each uniquely labeled drug within a cell as a proxy tocorresponding total drug exposure (area under the curve). This allows us to barcode each cellwith information about the drug regimen to which it was exposed. Elgart  in Figure 2b as a raw image. Here, we used a rather crude and simple method to establish transient spatial gradients in suspension cell culture by point injection (Methods). The quantitative comparison of homogeneous and point injection-based dye delivery in terms of the fluorescence intensity distribution and pairwise correlation relationship is shown in Figure 2c and Figure 2d, respectively. Additional heterogeneity introduced by point injections leads to a widening of the fluorescence intensity distributions and, thus, a significantly broader concentration range compared to the homogeneous mixing experiment.
We next conducted parallel experiments with adherent HeLa cells using the gradient mixing method by point injection of different drugs in different locations in the culture well (Methods). Figure 1. The workflow of the temporal local gradient method: Random drug mixing is achieved by creation of the local temporal gradients across a cell culture sample. Cells exposed to different drug combinations are fluorescently barcoded. The drug combination and corresponding outcome (phenotype) are analysed using coarse graining approach to reduce biological and technical noise.
Due to the dispersion process, the initial drugs' distribution is spread ('averaged') locally across finite regions of space (spread size depends on drug diffusion and uptake rates, and also incubation time). This observation allows one to simplify the image processing by averaging fluorescence signal intensities not over individual cells, but, rather, over the mesoscopic regions of space containing multiple cells (Materials and methods). This simplification is critical since confluent cell growth conditions (necessary for a large statistical ensemble) usually present a problem for image segmentation. The statistics obtained using this approach describes typical local cell drug uptake within the given mesoscopic region (tile).
The results of a typical experiment utilizing point injections are compared to homogeneous dye delivery in Figure 3. Here, we either co-delivered all dyes homogeneously or locally by point injection. Point injections were performed by local delivery of two dye pairs simultaneously: One pair of premixed dyes was codelivered in one point location and another pair of pre-mixed dyes was injected at a different point location on the microscopy cell culture slide. As in the case of cells grown in suspension culture, point injection greatly increased the range of drug uptake across the population compared to homogeneous delivery (Figure 3a). We observed correlation patterns between dye uptakes that were dependent on whether dyes were paired by point delivery or not. Locally co-delivered dyes (point-injection) exhibited a far greater degree of correlation compared to those in unmatched pairs of drugs delivered to different slide locations ( Figure 3b). This behavior is apparent from the scatter plots shown in Figure 3c where color represents local cell count (higher count corresponds to a redder color in the color spectrum scheme). Dye pairs delivered independently by the point-injection method (two different culture slide locations) led to the generation of different mixing ranges of dyes (first two columns, Figure 3c). The pre-mixed dyes co-delivered locally exhibit a high degree of correlation with coefficients of correlation greater than 0.8 (cf., right column in Figure 3c), even for structurally and chemically different agents such as CellTrace and nanoparticle QD dyes. The underlying cause of this surprising behavior is identical initial dispersion of pre-mixed co-delivered dyes throughout a sample by forced advection. The subsequent diversion of concentration profiles due to differences in dye diffusion, transport, and reaction rates occurs slowly over this comparatively short incubation period, ensuring persisting correlation between dye uptakes throughout the course of the typical experiment. In these two experiments, the locations of point injections were altered to demonstrate the independence of co-delivered dyes parameterized from the geometric setting. The linear model fits of coarse-grain data are R 2 0.99 and 0.97 for these pairs. (e, f) Point co-delivery of CellTracker dye pairs (CMTPX and Deep Red) in fixed cells: (e) coarse-grained intensity data and its linear regression fit (f) the combined scatter plot of single-cell intensity data for two separate experiments in which dyes were delivered at homogeneous concentrations. Here, we used the same molar ratio of dyes (1 : 1) as in the point injection delivery, but delivered dye combinations homogeneously at concentrations of either 5 µM or 10 µM . The bi-modality is due to the separation of individual distributions of single-cell intensity data for each of these experiments.

Conjugation-free drug barcoding
One can think of one of the dyes in co-delivered pairs as a drug and another as a label/tag. Since we can predict drug uptake utilizing the high degree of correlation with the corresponding tag, no physical conjugation of drug and tag is necessary.
Here, we rely on the similarity of temporal concentration profiles of drug and dye within cell culture. Drug (and dye) uptake is governed by transport and diffusive processes, as well as molecular interactions with its targets, both specific and non-specific (Elgart et al., 2018).
In many cases, strong, specific drug binding to the target leads to fast reaction rates which, in turn, result in the uptake process being diffusion/transport-limited. In this situation, media concentration profiles created by the initial dispersion drive intracellular drug and dye uptake. We pre-mixed reagents used for co-delivery in a small volume that is injected into a much larger sample volume. Hence, point co-delivery results in similar local concentration profiles of drug and dye, at least for short incubation periods.
Unlike chemically labeled drugs (where the one-to-one correspondence between tag fluorescence and drug uptake is assured by design), the conjugation-free barcoding method requires statistical processing. This step is necessary to filter out the noise due to intrinsic fluctuation in both drug and tag uptakes and improve correspondence between tag fluorescence and drug uptake.
A demonstration of the predictive accuracy in estimating drug uptake using this conjugation-free approach is shown in Figure 4. Here, coarse-grained uptakes of co-delivered dye pairs are shown in the scatter plots, Figure 4a and Figure 4b, with linear model fit and confidence bands depicted by blue and yellow lines, respectively. Average intensity values for the homogeneous dye delivery method are depicted by the large, red-filled circle on both graphs, demonstrating the greatly expanded range of combined dye concentrations by the point injection method. These results are also independent of the specifics of injection point location and geometry, as demonstrated in Figure 4c and Figure 4d, where we compare raw intensity data from two dyes in two different experiments involving two different pairs of injection site locations. and drug action proxy HyPer7 sensor ( H 2 O 2 fluorescent detector) were imaged by confocal microscopy. Auranofin was pre-mixed with CellTracker Orange dye ( 10 µM concentration each) and co-delivered locally using point injection. Quantification of the relationship between drug and dye uptake was performed using coarse grained analysis. The linear model fit for the coarse-grained HyPer ratio and 95% confidence intervals are shown by black and orange lines, respectively (adjusted R 2 = 0.89 ). The same scatter data are shown on log-log scale in the inset to demonstrate the broad log-linear range ( ∼ 100 fold) of tag and drug correlation.
We also observed a high degree of linear correlation in the uptake of co-delivered dyes in fixed cells ( Figure 4e). As in the case of live cells, homogeneous dye delivery to fixed cells can be predicted well by linear regression obtained from the point co-delivery data ( Figure 4e). The correlation in uptake ( ρ = 0.79 ) for point co-delivery in fixed (suspension) cells is comparable to the correlation observed in experiments with live cells ( ρ ≥ 0.8 ).
To go beyond the 'toy' system of conjugation-free co-delivery in which we treat one of the diester dyes or quantum dot (QD) nanoparticles as a drug, we co-delivered an authentic drug, auranofin, and a conjugation-free tag (Orange CellTracker dye) to human pulmonary artery endothelial cells (HPAEC). Auranofin is a thioredoxin reductase (TRR) inhibitor that has been shown to increase substantially intracellular H 2 O 2 levels. We utilized HyPer7, a genetically encoded fluorescent probe for H 2 O 2 detection, as a reporter for the drug's effect on cell phenotype (Pak et al., 2020). The co-delivery of auranofin and conjugation-free dye was performed by point injection. The results of the experiment are shown in Figure 5, where we quantified the fluorescence of the barcode and the H 2 O 2 sensor's ratiometric reading from the imaging data using a coarse-grained approach.
We observed a significant degree of correlation ( ρ ∼ 0.5) between the barcode tag for the drug and the drug effect as assessed by corresponding fluorescence signals ( Figure 5). The observed relationship was largely linear with saturation in the ratiometric signal for H 2 O 2 detection in cells with high auranofin uptake levels.

Effect of siRNA combinations on gene expression
Our benchmark system (dye combinations) was designed to demonstrate random drug mixing capabilities. We next considered a biological system with a detectable drug mixing effect that could be monitored. To do so, we chose drugs that are direct binding partners of a common, measurable target. We also designed this system to be as portable as possible, that is, be adaptable to different cell types and culture conditions. We implemented this system as follows: The drugs' target is the exogenous green fluorescent protein (eGFP) gene (d2EGFP, a short-lived variant, cf., Materials and methods). The drugs are small interfering RNA (siRNA) duplexes targeting different regions of the GFP transcript.
We designed two classes of siRNA: (i) typical 21 nucleotide-long siRNA molecules targeting complementary binding sequences on eGFP mRNA; and (ii) 'msiRNA', that is, siRNA with mismatch(es) in its targeting binding sequence. The msiRNA guide strand may still bind target transcripts (and, hence, interfere with target gene expression) via partial sequence complementarity by a mechanism closely mirroring microRNA (miRNA) silencing.
By analogy with endogenous miRNA in eukaryotic cells, one can expect weaker target protein repression by msiRNA compared to siRNA (Carthew and Sontheimer, 2009). [Naturally occurring Figure 7. The effect of siRNA combinations on target gene (GFP) expression. Three different siRNAs were delivered to cells by optically distinguishable QD nanocarriers (Q 1 and Q 2 ) in pair-wise fashion. Two siRNAs, siRNA 1 (S 1 ) and siRNA 2 (S 2 ), were designed to bind non-overlapping sequences of the mRNA target in completely complementary fashion. The third siRNA, msiRNA 1 (M 1 ), has almost identical primary structure as siRNA 1 with a single nucleotide mutation close to the seed region. The GFP expression was normalized with respect to control samples. (a, b) Relief of GFP repression mediated by msiRNA 1 is shown for different concentrations of the repressors siRNA 1 and siRNA 2 , respectively. (c) Combination of non-overlapping repressors siRNA 1 and siRNA 2 results in additive/synergistic repression of target GFP expression. (d) Control sample of siRNA free QD combinations was used to account for QD toxicity for both carriers using a machine learning algorithm, and its predicative accuracy is shown here for different concentrations of Q 1 and Q 2 . miRNAs can typically achieve significant repression of the target gene only by utilizing multiple/ tandem binding sites on the target mRNA].
We created more diversity in the drugs' action by introducing two different classes of siRNAs targeting different sequence locations: (i) well separated and (ii) overlapping. The motivation behind this design followed theoretical considerations. For low drug concentrations, we expected additive drug interaction regardless of the regions' proximity since there is an abundance of target copy numbers to which siRNA molecules can bind. At high drug concentrations, we expected different drug interaction patterns that depend on target site locations. Different siRNA molecules targeting overlapping regions may, for example, antagonize each other due to binding competition.
For siRNAs targeting non-overlapping regions, one can expect additive or, perhaps, synergistic drug interactions owing to accelerated mRNA degradation and modulated translation rates in the presence of multiple RISC complexes bound to the mRNA. The siRNA design principles, specific sequences, and corresponding target region locations are shown in Figure 6a and b (Materials and methods for the details). Cell-penetrating peptide-conjugated, nanoparticles, Quantum Dots (QD), were used both for drug delivery and barcoding (Figure 6c and d). We delivered QD-conjugated siRNAs targeting GFP mRNA in a pairwise fashion. In these experiments, we point-injected all possible combinations of three different siRNAs: siRNA 1 , siRNA 2 , and msiRNA 1 . As expected, point injection of a drug conjugated and delivered in this way introduces a broad range of drug uptake across the sample (Appendix 1figure 10), as compared to the homogeneous random mixing approach (Appendix 1-figure 9b), and suppresses the proportionality bias in drug combinations.
The siRNA 1 and siRNA 2 species have perfect complementarity to non-overlapping regions of mRNA. The msiRNA 1 was designed to bind to the same region of mRNA as siRNA 1 but has a mismatched nucleotide pair within the targeting sequence. The siRNAs were barcoded using two optically different Local cell density is determined by the fraction of cell-covered area within a square of 10 4 µ 2 . High and low cell confluency corresponds to a fractional range from 1 (blue) and 0 (red), respectively. (c) The pairwise initial drug distributions.
QD nanocarriers, QD 605 and QD 655 . To demonstrate that the drug combination effect is independent of the tag, we exchanged the tags for each siRNA mixing experiment we performed.
Inherited randomness in cell behavior is not only responsible for variable drug uptake, but also all but guarantees variability in response to drug treatment by individual cells. Even under conditions of identical uptake of a drug and the same exposure time, one expects heterogeneity in individual cell response/phenotype. To overcome this source of heterogeneity, we perform averaging over many cells within each sub-populations characterized by a specific combination of tags (and hence drugs). The number of cells in each subset must be large enough to suppress the variability in drug response ensured by the Law of Large Numbers.
The results of siRNA combinations on target gene expression are shown in Figure 7. Here, we partitioned the QD tag intensities into the 2D bins spanning effective parametric space. GFP intensity was averaged across a subpopulation of cells in each bin. Since the QDs at high concentrations lead to cytotoxicity, GFP intensity for each condition was normalized to the control sample in which we The y-axis values represent cellular dye uptake for each stage of staining (panels 1.1, 1.2, and 2.1). The final dye combinations are shown in the bottom right panels for both delivery methods (the bar color and the amplitude reflect the molar ratio and total uptake of dyes, respectively). Each wash step is mimicked by reshuffling the cell position index (cf., panels 1.1 and 1.2). Variation in dye uptake due to intrinsic cellular heterogeneity is illustrated for each of the delivery methods by imposing multiplicative noise. For simplicity, the noise in uptake is assumed to be independent for each dye. (c) Graphic illustration of the point co-delivery method.
point-injected QD 605 and QD 655 nanoparticles not conjugated to siRNA. (This strategy is equivalent to normalization of repression levels using scrambled siRNA samples as controls.) We used a machine learning algorithm to predict the normalization factor for QD combinations using this control sample (see Supporting Information for details).
At high concentrations, the "mismatched" siRNA, msiRNA 1 (with incomplete complementarity in its target binding sequence), acts as an inhibitor for the strong repressors siRNA 1 and siRNA 2 , as shown in Figure 7a and b. This antagonistic effect is more pronounced in the siRNA 1 -msiRNA 1 combination as compared to the siRNA 2 -msiRNA 1 combination.
The outcome of pairwise combinations of the non-overlapping strong repressors siRNA 1 and siRNA 2 suggests additive or synergistic interaction for these drugs (Figure 7c). The bias in the titration range of the siRNA concentrations in this experiment is due to limitation of effective mixing conditions (see Supporting Information for details). The predictive accuracy of the machine learning algorithm for GFP normalization (due to toxicity of QDs) is shown in Figure 7d for a broad range of mixing conditions.
Here, we, once again, employed the tiled image renormalization approach, which is especially well suited for point injection data analysis because local gradients in drug concentrations persist over scales sized larger than that of a single cell. Since tiles in the image relate to cells in the same spatial proximity, tiled renormalization is equivalent to concentration(s) gating by the local microenvironment.
The results of the pairwise drug combinations in Figure 7 demonstrate both additive and antagonistic effects of the siRNA combinations. We note that direct binding competition is not the only source of possible antagonistic siRNAs interaction. The repression by siRNA is mediated by the RNAinduced silencing complex (RISC). At the high intracellular concentration of siRNA duplexes, one may expect competition for RISC integration by different siRNA species due to the limited pool of RISC in the cell. This type of competition may result in antagonistic siRNA interactions, regardless of (the abundance of) their target (Castanotto et al., 2007;Loinger et al., 2012).

Random mixing and response quantification of three and four drugs
To demonstrate that our approach can be applied to combinations of more than simply pairs of drugs, we considered the following application. We tested drugs commonly used in combination chemotherapy regimens to treat different forms of cancer, such as cyclophosphamide (CP), docetaxel (TXT), 5-fluorouracil (5-FU), and cisplatin (CIS). The desired phenotype in this application is cell death, which can be assessed by utilizing fluorescent biomarkers, such as Annexin-V stain. Alternatively, death phenotype can be assessed by monitoring cell detachment and/or changes in cell morphology, which may be resolved by bright field/phase contrast microscopy. Cell shrinkage due to dehydration and cell debris can also be detected by a decrease in intensity of the forward light scatter (FSC) signal in flow cytometry measurements (Wlodkowic et al., 2011). Unlike annexin-V staining, these cell death detection methods do not require, in principle, additional fluorescence biomarker(s) and, hence, do not reduce the necessary bandwidth required to resolve quantitation of drug combinations in individual cells.
We used the HeLa cell line as a benchmark system to evaluate the effect(s) of drug combinations in individual cells. We utilized either time resolved confocal microscopy, flow cytometry, or a hybrid of both detection methods in our experiments. In all experiments reported below, cells were treated with drugs for a short duration (1/2 hr). Hence, drug concentrations were chosen to be significantly higher than their corresponding typical IC 50 values. (In typical experiments, drug(s) incubation is 48 hr, roughly 100 times longer than in our setting.) In order to create a negative drug control, cisplatin was inactivated by DMSO (Hall et al., 2014).
We note that our goal was to demonstrate the technical feasibility of our stochastic mixing approach. A more practical, generic application (in terms of both an expanded list of drugs and of cell lines) and biological data interpretation are beyond the scope of this sudy.

Flow cytometry measurements of cell death phenotype, three-drug combinations
In order to assess drug combinations leading to cell death using flow cytometry, we harvested detached cells and debris in the sample media at different time points and analyzed the fluorescence intensities of the specific drug tags, which correspond to effective drug combinations. At the final time point of the experiment, we harvested the remaining attached cells and analyzed the fluorescence intensities of the specific drug tags similarly, which correspond to ineffective drug combinations.
We stained cells homogeneously with the same auxiliary dye (Deep Red CellTracker dye) not conjugated to any drug. We used an identical concentration of this dye in all samples, including controls. The CIS, CP, and TXT drugs were labeled with blue, green, and orange CellTracker dyes, respectively. All CellTracker fluorescence probes (except for Deep Red) contain a chloromethyl or bromomethyl group that reacts with thiol groups, utilizing a glutathione-S-transferase-mediated reaction.
In most cells, glutathione levels are high ( mM range) and glutathione transferase is ubiquitous. The CellTracker Deep Red probe contains a succinimidyl ester functionality, which reacts with amine groups present on proteins, and, hence, is also ubiquitous in cellular cytoplasm. The ratio R i of any CellTracker dye i signal to Deep Red intensity can be viewed as a proxy to cellular concentration of the CellTracker and its cognate drug, since the Deep Red signal is proportional to cell volume.
We assumed that for sufficiently large cell fragments, this normalization also holds, which allows an estimation of the drug concentrations in both intact cells and dead cells (debris). We note that due to stochastic uptake of any dye including an auxiliary dye, this normalization is only meaningful after averaging over a significant number of detection events required to suppress the signal noise.
In order to convert the ratio R i to absolute drug concentration, we calculated R h i for a homogeneously delivered (with mixing) drug sample where the drug concentrations are known exactly. The final concentrations of drugs in these homogeneously delivered controls were {1, 1, 0.1 µM} for NC, CP, and TXT, respectively.
The factor R i /R h i calculated for individual cell is a proxy for the volume concentration of the drug i compared to homogeneously delivered drug control. Given the ubiquitous nature of all CellTracker dye targets, one expects that this estimate will hold even after cell division. The absolute fluorescent signal from the parent cell will decrease due to dilution by cell division, but the ratio of homogeneously distributed targets should persist.
To demonstrate that this expectation holds, we considered dynamic changes in fluorescence signal intensity in the controls of homogeneously delivered drugs. The population average factors, R i for i={NC, CP, TXT}, are shown in Figure 8a at different time points.
Small time-dependent modulation of R h i values are noted and are a consequence of the stochastic nature of drug or dye uptake, a process that ensures the existence of cells with a distribution of drug concentrations, including importantly some that take up less drug(s) than the cell average. The corresponding R i values for these cells is smaller compared to cells that received higher drug dosage. Hence, this low R i subpopulation is more likely to survive longer. Since this dynamic drift is rather small, however, we neglect it in further analysis.
For samples created by gradient drug mixing, one naively expects qualitatively the following picture for the population average of drug concentrations: the concentrations of effective and ineffective drugs should decrease and increase, respectively, with time. This expectation derives from the fact that the subpopulation of cells that received ineffective drug concentrations or drug combinations becomes increasingly dominant over time.
This qualitative picture (evaluated on a whole cell population level) holds only if drugs act independently in individual cells. Drug pharmacodynamics may also be non-uniform over time. We observed these behaviors in the gradient-derived mixed sample data shown in Figure 8b.
The dynamics of inactivated cisplatin (labeled NC in Figure 8b), indeed, displays the time dependence expected from an ineffective drug (concentration). It is much more difficult to interpret the pharmacodynamics of CP and TXT because the time dependence is not monotonic. Drug-drug interactions, drug effect on cell cycle, and other time-dependent drug actions can all play a significant role in overall cell population behavior.
To visualize the interplay of drugs on individual cells, we first reduced data dimensionality by discarding information about the level of [NC] in cells (assuming DMSO inactivated cisplatin completely). Under these circumstances, one can plot the probability distribution, P(x, y) Figure 9b.
Our overall goal was to demonstrate the technical capability of the stochastic mixing approach. More practical application (in terms of both list of drugs and cell lines) and data interpretation are beyond the scope of this study.

Microscopy measurements of cell death phenotype, four-drug combinations
For microscopy-based measurements and image analysis, we used the point injection method to create local drug gradients across the sample of confluent HeLa cells on a glass slide. Followed by the initial imaging of the drugs' distributions across the sample, we monitored cell detachment (and changes in morphology, e.g. shrinkage) at different locations on the slide at different time points. The estimated intracellular concentrations of all chemotherapy drugs used were assessed using the nonconjugated fluorescent tag (dye) method described above.
The image analysis allows one to correlate the initial drug combination concentrations in each locale on the slide and subsequent death phenotype in the same position at different time points. In principle, cell migration between locales can complicate the analysis, but we used confluent seeding conditions to minimize this effect.
The outcome of a typical microscopy experiment is shown in Figure 10. Here, the initial distribution of each drug across the sample is shown in Figure 10a using four different colors (blue, green, red, deep red) corresponding to the fluorescent label for each drug. Local cell density at different time points is shown in Figure 10b, using an inverse rainbow color scheme (where blue and red correspond to fully confluent cells and no cells, respectively). All possible pairwise initial drug distributions are shown in Figure 10c, see also Appendix 1- figure 16 in Supporting Information.
We used a coarse-grain approach for image analysis, averaging fluorescence intensities over an area of ∼ 10 4 µm 2 at each locale. Cell density was assessed by cell segmentation using bright field imaging data and the calculation of the area fraction covered by cells in each coarse-grain locale. The absolute drug concentration can be determined using the homogeneous delivery method discussed above.

Discussion
Biologists often treat noise in samples as a nuisance that must be minimized, an important consideration for typical experimental measurements. Here, we rely on engineered local heterogeneity in drug uptake as a tool to encode drug treatment regimens and to predict the macroscopic response to drug perturbations. This approach was inspired by the celebrated fluctuation-dissipation theorem in statistical physics (Cardy, 1996) that predicts the response of a system to external perturbations based on the underlying properties of noise in the system at equilibrium.
Rather than testing thousands of samples, as little as a single-cell culture sample is required to perform the search for optimal drug combination using our approach. Applications of our method are not limited to the identification of optimal drug concentrations, as the approach is applicable to any quantifiable cell phenotype. The method is readily generalizable to a much broader set of problems involving complex biochemical or molecular perturbations applied in vitro and in vivo. Examples of these perturbations include combinations of signaling molecules, optimization of drug delivery strategies, lineage differentiation identification, stem cell manipulations and tracking, etc.
The best possible random drug mixing is achieved in cells in suspension. Suspension cells can be easily mixed (i.e. cell locations can be randomly re-distributed), and local drug gradients can be applied in a statistically independent serial fashion. By contrast, since the locations of adherent cells in the sample are fixed, we exclusively rely on the randomness of the drug distribution in the media by means of forced advection, which is challenging to control. Hence, further theoretical and experimental work is required to improve the stochastic drug mixing range and its smoothness for adherent cells. A possible way to improve the mixing range (and minimize reagent use even further) can be utilization of recently developed 'Dye Drop' method (Mills et al., 2022).
We would like to point out an important limitation of our stochastic drug mixing method (at least in its current state). Throughout this paper, we tacitly assumed that both drugs and their labels are retained by the cells and are not re-distributed from cell to cell after the washing step. This is certainly the case for siRNA drugs delivered via nanoparticles and CellTracker dyes in our applications, but this assumption does not necessarily hold in general.
Drugs can, in principle, leak from the cells into the culture media after the washing step and, eventually, be redistributed homogeneously across the entire cell population. This behavior will affect the 'area under the curve' for individual cells after prolonged incubation periods, especially cells that had initially very low (or no) drug exposure. We assumed that this effect is small and did not study the kinetics or the extent of this 'secondary' drug exposure here. (Multiple washing steps were performed to mitigate this effect.) More careful analysis and 'compensation' is required for drugs that are not retained by the cells.
Another possible drawback of our method is a limit on maximal drug concentration due to drug solubility and its solvent (e.g. DMSO) toxicity. Since we utilize short drug exposure times, high drug concentrations may be required to test a broad range of conditions. Typically, low concentrations of individual drugs are desired in combination therapies and, hence, our method can be implemented even for drugs with low solubility.
We expect possible adaptation of the point injection method to facilitate random mixing in vivo. The spatial distances between injection sites, geometry, needle gauge size, injection volume, and the specific measurement protocol (readout) will, of course, all need to be established and optimized for a particular application.
Throughout this paper, we focus exclusively on the measurement approach rather than on data interpretation, which is context-and phenotype-dependent. Further work is necessary to develop a methodology, interpretation, and dedicated numerical methods to identify optimal drug combinations from high dimensional experimental data.

Drug labeling and delivery
The workflow of the temporal local gradient method used in this work is depicted schematically in the Figure 1. The cell culture sample is placed on a horizontal vibration-free surface, and buffer containing labeled drug is injected in the periphery of the sample. The ratio of injection to sample volumes was typically 1 − 5 % and injection speed is slow to ensure local initial distribution of the drug in the sample. (In order to visualize and adjust injection volume, release speed, and dispersion dynamics, dyes visible by naked eye or food-dyes such as Brilliant Blue FCF can be used on training samples.) Care should be taken not to disturb the sample during the incubation period (to avoid sample homogenization). After the incubation period, culture media should be quickly replaced with fresh media in adherent cell culture samples.
In suspension cell experiments, cells are placed into narrow cuvettes using a buffer density lower than cell density. After cell sedimentation, their spatial distribution corresponds to a narrow rectangle on the bottom of the cuvette. Local point injection is delivered in the 'corner' of this rectangle sufficiently slowly so as to minimize the cells' movement in the cuvette. Mixing labeled drugs with iodixanol (OptiPrep), an inert liquid with a density higher than PBS/media, can be used to minimize dispersion of local injection in the vertical direction (Mills et al., 2022). Immediately after the drug incubation, large wash volume is used to minimize re-distribution of drug(s) within the sample prior to and during spin down step.
We used multiple different dyes to perform the multiplexing experiments: CellTracker dyes (blue, green, orange, CMTPX red, and deep red) and CellTracer dyes (violet, carboxy fluorescein succinimidyl ester [CFSE], yellow, red, and far-red). The CellTracker probes have been designed to pass freely through cell membranes, and once inside the cell, they are converted into impermeable fluorescent products. The probes are transferred to daughter cells, but are not transferred to adjacent cells in a population. The CellTracker fluorescence probes contain chloromethyl/bromomethyl or succinimidyl ester reactive groups that react with ubiquitous intracellular thiol or amine groups, respectively.
More flexible barcoding can be achieved by utilizing a tag bound to a drug by a cleavable linker. The tag can be designed to be retained by the cell, with the delivered cellular drug amount proportional to the tag's fluorescence signal. We used nanoparticle drug delivery to apply this labeling method, namely, cell-penetrating peptide-conjugated quantum dots (QD) were used both for drug delivery and barcoding. The fluorescent reporter (QD) is retained within the cell upon entry and the drug, siRNA in our case, is released locally into the cytosol via cleavage of a disulfide bond in the cytosolic reducing environment.
Finally, we utilized the point injection method to co-deliver locally both drug and conjugationfree fluorescent dye retained by the cells. We note that to achieve one-to-one (monotonic) correspondence between total cellular drug uptake and its conjugation-free tag, it is important to ensure non-saturation labeling and detection conditions. If the tag fluorescence intensity saturates beyond a certain exposure concentration (or time), the total drug uptake prediction will be inaccurate. Most of the dyes we used (CellTrace and CellTracker) have abundant intracellular target concentrations and saturate only at very high extra-cellular concentrations ( > 1 mM ). Therefore, the saturation in tag signal may be caused by microscopy detection limits (for either too low or too high fluorescent signals at the given laser power and sensitivity of the photomultiplier settings). Careful adjustment of these parameters is required to ensure a maximal dynamic detection range.
To avoid label and drug homogenization during the wash/collection steps using media (containing both fluorescent tags and drugs), the washing steps need be performed quickly using large dilution volumes. This homogenization, in turn, can affect local labeled drug gradients in the culture, reducing the dynamic range of drug concentrations within the sample. Cells that were initially exposed to very low local drug(s) concentration are most affected. We note that this short time scale 'leakage' effect occurs with both label(s) and drug(s) proportionally. Hence, label intensity should be proportional to the corresponding drug concentration, and can be still used as a proxy for drug concentration.
Long time scale leakage occurs if the drug within cells is released to neighboring cells either via media or through cell-to-cell contact. In general, this secondary drug exposure cannot be accounted for using fluorescent tags in our method because the tags/labels are designed to be irreversibly bound to the cells and do not leak (labels can, however, be diluted by cell division). Drug leakage from cell to cell mediated by media can be estimated using the volume ratio of cells and culture media. This ratio in most cell culture condition is negligibly small and leakage of this sort can be ignored.
Drug leakage from cell-to-cell mediated by direct contact can be much more significant. However, this effect is mitigated in adherent cell culture experiments by the coarse-graining procedure when the optical field containing multiple neighboring cells is averaged.

Effective drug combinations by stochastic mixing
Steps involved in the point injection method compared to homogeneous delivery are shown in Figure 11a and Figure 11b, respectively. Cells in suspension culture were placed in a cuvette and a dye A (Red Dye in the example shown in Figure 11b) was injected locally (just beneath the meniscus). The volume of the point injection must be small compared to the cell culture volume in the cuvette (we used ∼ 1 : 20 volume ratio). (N.B., the graphic illustration in Figure 11b is an idealization because the gradient created by point delivery is not spatially linear.) The cultured cells suspended in the cuvette were left undisturbed during the short incubation process (10-30 min). The transient gradient of dye A results in higher dye uptake by cells located closer to the point injection site. The subsequent wash step ensures removal of remaining dye A and also mixing of differentially stained cells. The process is then repeated with another dye B (Blue Dye in the example shown in Figure 11b) in exactly the same fashion as with dye A .
As a result, a broad range of drug combinations can be achieved. Consider, for example, cells in the proximity of the dye B point injection site. These cells exhibit a wide range of dye A uptake and high uptake of dye B . By contrast, cells remote from the dye B point injection site are also characterized by a broad range of dye A uptake but low uptake of dye B . The conjugation-free barcoding method to co-deliver locally both drug and fluorescent dye retained by the cells is illustrated graphically in Figure 11c.
The local injection method was used to create stochastic drug combinations in adherent cell culture as well. (Cell mixing, naturally occurring in suspension cell culture, is not feasible in this case unless cells are detached and re-seeded, actions that require time and may interfere with drug response.) The local drug gradients across the single sample created using this method are driven mostly by initial drug dispersion (forced advection) during the injection if the incubation time is short (cf., Supporting Information for technical details on diffusion rate of small molecules).
The results shown in Figures 3 and 4 correspond to experiments performed with 90% confluent HeLa cells cultured in four wells with a chambered coverslip Ph+ (Ibidi). Point injections were performed pairwise: pair (i) CellTrace Violet and Qtracker 655 Cell Labeling Kit (referred to as Blue and QD655 dyes in the main text); pair (ii) CellTrace CFSE and CellTracker Red CMTPX (referred to as Green and Red dyes in the main text).
The total culture volume of each chamber was 750 µl and 30 µl were replaced with a pre-mixed pair of dyes injected at one of the chamber corners. The process was repeated with the second pair of pre-mixed dyes injected into another chamber corner. The manual injection rate is difficult to control but the total release time was ∼ 1 s. The incubation was performed at room temperature for ∼ 20 min in the cell culture hood (with all possible precautions taken to avoid shaking the slide).
The graphic illustration in Figure 12 demonstrates this approach for a pair of drugs. Here, an idealized drug mixture distribution is depicted by color and intensity in Figure 12a. All fluctuations due to cell variability, advection, and non-independence in drug uptake are ignored in these illustrations for simplicity. Moreover, we assumed that transient drug uptake gradients are symmetric (radial) with respect to the point injection locations (left and right top corners of the rectangular slide). In reality, of course, all of these assumptions, including radial symmetry, are violated, resulting in much more random behavior, which is precisely what we wish to exploit to enhance the heterogeneity in drug uptake across the cell population.
In Figure 12b, the origin of the distribution as a linear superposition of individual drug uptakes is shown, which is designed to illustrate that unique mixing conditions, even in this idealized case, depend on the geometry of the matrix and point injection locations. Discretization (binning) of uptake for each drug has a simple geometric interpretation in this idealized case, shown as regions of annular intersections in which colors are uniquely mixed. For a real system, this symmetry is violated due to processes such as advection, and we, therefore, must rely on the fluorescence intensities of the tags rather than geometry to partition the ensemble into bins; however, as in the case of the idealized system in Figure 12b, one expects a finite-sized local region on the slide where drug concentrations/ uptakes can be assumed to be constant.

Statistical processing of single-cell data
In all single-cell experiments, multidimensional intensity data were compiled using averaged fluorescence intensity measurements of each optical drug barcode for each cell. Additionally, the ensuing phenotype for each individual cell was recorded using the corresponding fluorescence signal in functional experiments (conjugation-free barcoding, auranofin, and siRNA combinations).
For experiments utilizing a single barcoded drug, we used coarse-grained modeling to recover the relationship between drug and tag as follows: all measured single-cell fluorescence intensity data were ordered by increasing value of the intensity of the tag. The ordered data were grouped into bins (subpopulations) of equal length r (number of cells) large enough to suppress intrinsic noise. After averaging the intensity data in each bin, the dimension of the resulting data array is reduced by a factor r . This coarse-grained process ensures reduction in intrinsic noise due to the Law of Large Numbers (cf., Supporting Information). Alternatively, noise filtering can be implemented by moving the average approach, preserving data dimensionality.
A similar strategy was employed for coarse-grained analysis of barcoded multi-drug experiments. Cell intensity data were clustered hierarchically based on the barcodes' signals, and renormalized data were produced by averaging all optical parameters within each cluster (cf., Supporting Information).
Specifically, for siRNA experiments, after image segmentation, fluorescence of QD drug carriers and of GFP was assessed for individual cells. Each drug's uptake was assumed to be proportional to the intensity of the corresponding QD carrier in each cell. The data were typically recorded as a 3-tuple {QD 1 , QD 2 , GFP} i for each cell i .
We grouped cells into bins characterized by different QD combinations (in terms of fluorescence intensities {QD 1 , QD 2 } ). For pairwise drug combination experiments, binning was achieved by partitioning the 2d space of the corresponding QD-tags' fluorescence values into a finite mesh. The bins were designed to contain 30 or more cells each. We averaged a phenotypic marker (e.g. GFP intensity) over each binned cell population.

Cell culture and microscopy imaging
The great majority of the experiments were conducted with the HeLa cell line. We used human pulmonary artery endothelial cells (HPAEC) for quantification of auranofin delivery. All experiments utilizing the point injection method were performed in confluent cell culture. The cells were grown on culture plates (microscopy slides) as follows: we used rectangular 4 and 8 chamber slides (Ibidi μ-Slides) for siRNA experiments, and either rectangular or circular slides for dye labeling experiments.
In experiments requiring single-cell imaging, cells were either 'shrunken' or completely detached using Accutase treatment to facilitate image segmentation. Homogeneous dye and QD labeling were performed according to the manufacturer's (Thermo-Fisher) suggested protocol(s), typically at room temperature.
For point injection experiments, we used 10 x dye concentration (compared to the homogeneous staining condition). Delivery was carried out by injection of ∼ 0.05 − 0.1 cell culture volume locally. Cell imaging was performed with an inverted fully motorized Zeiss LSM-800 confocal scanning microscope utilizing 405, 488, 561, and 640 nm diode lasers for excitation of dyes and markers. Most of the microscopy experiments were performed using the 5 x objective to capture larger optical fields necessary for statistical analysis.
Imaging data from single-cell experiments were processed as follows: cell segmentation was performed using bright-field and/or dedicated cytosolic markers. We employed Cellpose (Stringer et al., 2021) for accurate cellular segmentation. For the auranofin co-delivery experiments, the ratiometric signal of the HyPer7 sensor was calculated as the ratio of cellular fluorescence intensities of the GFP channel excited by 488 and 405 nm laser sources.
Excitation and emission settings were chosen according to the dye manufacturer's suggestions. In all QD experiments, we used a 405 nm laser for excitation with corresponding emission spectra filters. We tested fluorescent labels/markers individually to detect possible cross-channel leakage/overlap. This analysis was especially critical in dye co-delivery experiments where the degree of correlation may be strongly affected by channel leakage. To ensure channel 'isolation', we utilized dyes/fluorophores with well separated emission spectra for these experiments.
HeLa cells were maintained in Dulbecco's modified Eagle's medium (DMEM) containing 0.11 g/ liter sodium pyruvate, 2 mM L -glutamine, 4.5 g/liter glucose, and 10% fetal bovine serum (FBS). Human pulmonary artery endothelial cells (HPAEC) were cultured at 37°C in modified EGM2 medium (Clonetics). Cells were maintained in continuous culture in an air-5% CO 2 atmosphere at constant humidity and used within three to five passages. For fixation experiments, after aspiration of the culture media, cells were incubated with 10% formalin for 5 min at room temperature followed by a PBS wash. For cell suspension experiments, we detached HeLa cells with Accutase solution (Sigma-Aldrich), washed the cells with culture media, and kept them at either room temperature or 4°C in dPBS in non-treated plastic cuvettes for a short duration of dye(s) labeling (typically 1/2 hour).
In siRNA experiments, point-injection of two different QD nanoparticles carrying different siRNAs was performed at two different corners of the cell culture plate (chamber slide). Samples were incubated at 37° for 30 min. (Of note, no special precautions were taken to suppress the advection process, which facilitated the creation of concentration gradients).

siRNA design and delivery method
We used short-lived green fluorescent protein (d2EGFP) to ensure correlation between fluorescent GFP signal and mRNA target concentration in cells. The d2EGFP plasmid (Warren et al., 2010) was a gift from Derrick Rossi, Addgene #26821. We used siRNAs that were designed and experimentally validated by Dharmacon (cf., Accell eGFP Control siRNA) and the Sharp laboratory with binding strand sequences GCCA CAAC GTCT ATAT CAT (siRNA 1 ) and GCAC CATC TTCT TCAA GGA (siRNA 2 ), respectively. The binding sites of siRNA 1 and siRNA 2 are separated by distance d = 153 nucleotides (in primary sequence). The msiRNA1 construct has a single nucleotide mismatch in the binding strand, GCCA CAAC GGCT ATAT CAT, and the separation distances to siRNA 1 and siRNA 2 were d = 0 and d = 153 nucleotides, respectively.
We utilized a fluorescent tag for each siRNA that was retained by the cell as follows. All siRNAs were designed with a cleavable disulfide linker (CL) attached to their 3' end ( Figure 6). We used streptavidin-conjugated Quantum Dots (QDs) as fluorescent tags for each class of siRNAs, creating uniquely labeled siRNA-QD pairs. Quantum Dots (QDs) are semiconductor nanoparticles, typically excitable by UV light, with fluorescence emission spectra that are fairly narrow. The emission peaks of QDs depend on the nanoparticle size, a property that makes QDs very useful for high-throughput microscopy measurements (Wegner and Hildebrandt, 2015). To facilitate siRNA-QD cellular delivery, we further functionalized the QDs with a cell penetrating peptide (CPP) [a nanopeptide (Arg) 9 ] that can facilitate cell uptake through endocytosis or non-endocytotic mechanisms (Verma and Stellacci, 2010) (cf., Supporting Information for technical details).

Coarse-grained image analysis
In addition to course graining of single-cell data, we develop a similar approach with respect to raw images. We performed effective local averaging over multiple confluent cells by partitioning the entire image into tiles. The size of the tile (typically, ∼50x50 microns) was estimated to contain at least a few cells on average.
Each tile in the image was replaced with an effective 'pixel' characterized by the average fluorescence intensity of real pixels within the tile. The only segmentation we performed in this case was identification of large lacunae (cell-free regions) in the optical field, with exclusion of these regions from image analysis. The funders had no role in study design, data collection and interpretation, or the decision to submit the work for publication. Appendix 1

Combinatorics of drug combinations
Test this Appendix 1- figure 1a and Appendix 1-figure 1c.
Appendix 1 The number of drug combinations, Nc , for a given number of drug doses can be expressed as where n i is the number of dosages of drug i, i = 1, · · · , N d . The Equation 1 accounts for all possible set of drug combinations by including an additional concentration condition c i = 0 for each drug i, i = 1, · · · , N d . For example, the number of all possible combinations of two drugs A and B at a single dose, N d = 1 , is (1 + 1) · (1 + 1) = 4 , namely { ∅ + ∅, A + ∅, ∅ + B, A + B }. The number of possible drug combinations Equation 1 scales exponentially with the number of drugs and dosages.
In the case of ex vivo testing, the size of the obtained (tissue or blood) specimen might not be sufficient to probe an enormous dosing parameter space. Another challenging aspect of the experimental identification of optimal combination therapy is logistics. Often, hundreds of samples need to be prepared, cultured, processed, stored, and analyzed. The whole undertaking cannot be entirely robotic, typically requiring human curation to ensure quality control.

Noise and law of large numbers
Fluctuations in measured biological parameters around their typical values in cellular systems can be described as a consequence of intrinsic and extrinsic noise Elowitz et al., 2002;Raj and van Oudenaarden, 2008. Here, we will refer to intrinsic noise as phenotype(s) fluctuations across isogenic cell populations cultured under the same conditions Kaern et al., 2005. Measurement noise in some cases can also be thought as intrinsic noise. Fluctuations in cellular phenotype(s) driven by the global environment will be referred to as extrinsic noise.
In cell culture, there exist several extrinsic factors that drive heterogeneity in the cells' response to the drug treatment. Among these factors are included cell density, local cell clustering in the well, passage number, time in culture prior to drug exposure, and position on a multi-well plate (leading, for example, to temperature and evaporation rate gradients). All of these factors contributes to extrinsic noise and, paradoxically, the fate of drug response experiments can be (pre-)determined (in terms of fluctuations in the outcome) even before they begin as a consequence. Hence, quantification of drug effects requires multiple replicates, and very carefully designed protocols and handling of the biological specimens.
There also exist many microscopic parameters that drive intrinsic noise in drug uptake. The drug diffusion coefficient, cell membrane potential and capacitance, active transport properties, and intracellular metabolism/degradation of the drug are among parameters controlling drug uptake by individual cells. Due to differences in cell shape and size, cell cycle state, and the microenvironment, these parameters can vary from cell to cell, further driving intrinsic noise in drug uptake.
While intrinsic cell variability can be significant (with coefficients of variation of the order CV ∼ 1 ), we believe that it is extrinsic factors that drive overall replicate sample variability in typical cell biological experiments. This conclusion is a consequence of the Law of Large Numbers Feller, 1968, since the typical number of cells in a culture sample is at least of the order of 10 4 . According to the law, the average of the results obtained from a large number of trials should be close to the expected value and will tend to become closer to the expected value as more trials are performed. Deviation from the expected value scales inversely proportionally to the square root of the number of trials. Hence, even if intrinsic cell variability is of the order CV ∼ 1 , averaging over 100 randomly chosen cells will decrease this variability to CV ∼ 0.1 , corresponding to 90% noise suppression. By contrast, the typical number of replicates in most biological applications is 3, averaging over which will only suppress noise by roughly 40%.
Throughout this paper, we utilize the Law of Large Numbers in different contexts, and, hence, it is worth demonstrating its action in a concrete biological example. To illustrate the application of the law, we assessed a few commonly used cellular characteristics, such as a gene expression level and dye uptake, at the single cell level. We also measured uptake by individual cells of such different fluorescent species as the Hoechst bisbenzimidine nuclear stain or quantum dot fluorescent nanoparticles (QDs). Detectable fluorescence signals from these agents can be thought of as a correlate to drug delivery: Fluorescence intensity of transient GFP-expressing cells is a direct consequence of Lipofectamine transfection of the corresponding plasmid. Hoechst dye becomes fluorescent only upon binding to the minor groove of DNA. QDs (functionalized by cell penetrating peptides) can be designed to carry a real drug payload. A large degree of intracellular heterogeneity in uptake is shown in Appendix 1-figure 1a for these diverse fluorescent species with corresponding reagents delivered sequentially to a single HeLa cells sample (cultured on a microscopy slide).
Using single-cell data as a starting point, we randomized cell ordering and uniformly partitioned data into bins. Each bin contains an equal number of elements, r , corresponding to measured parameters (attributes) of r distinct cells. We will refer to r as a coarse graining or scaling factor, borrowing the terms from the statistical physics approach to model reduction (renormalization) developed in 1960-1970s by Kadanoff, Wilson, and Fisher (Cardy, 1996. We averaged cellular attributes (e.g. GFP fluorescence) over each bin to produce a renormalized distribution of the phenotype for different values of the parameter r (Appendix 1-figure 1b).
A reduction in noise for all three phenotypes was observed after distribution renormalization, as shown in Appendix 1-figure 1b where narrowing of the distribution of Hoechst dye fluorescence intensity is demonstrated as a function of the scaling factor r (similar results hold for the other reporters, GFP and QDs). As expected, the degree of variability behaves precisely according to the Law of Large Numbers, as shown in Appendix 1- figure 1c where the coefficient of variation ( CV ) is shown as a function of r .
While it is reassuring that the Law of Large Numbers works flawlessly in a single cell statistical ensemble, we next need to consider its effect on a range of drug combinations. It is rather obvious that the renormalization process applied uniformly (i.e. randomly choosing cells for each bin) to a cell population results in a reduction of possible combinations precisely as a consequence of this law. This effect is illustrated in Appendix 1-figure 2a , where we show a scatter plot of QD and GFP cellular fluorescence intensities before (red) and after (blue) a uniform renormalization procedure. We can also apply renormalization to sub-populations of cells sharing a common characteristic(s). A trivial example of this process is illustrated in Appendix 1-figure 3, where we choose as a common characteristic a color. For example, we can partition the entire population of cells into bins designed to span the whole range of a particular reporter such as Hoechst dye fluorescence intensity. The partition into N bins is given by a sequence of intervals [b i−1 , b i ] , i = 1, · · · N covering the entire range [b 0 , b N ] of Hoechst dye fluorescence intensity, where b 0 and b N correspond to the minimal and maximal fluorescence intensity values, respectively, for the entire cell population. The resulting partition bins do not contain cells typical for an entire population but, rather, are organized by the degree of closeness with respect to the reporter used for cell sorting. The partitioning for each bin in this way is analogous to gating used in flow cytometry. Hence, we will refer to this partition strategy and subsequent coarse graining as a gated renormalization process (cf., uniform renormalization by binning cells randomly).
Gated renormalization is designed to integrate out fluctuations in attribute(s) of the cell that are not involved in the partitioning criteria. As a trivial example shown in Appendix 1-figure 3, the 'cellular' attribute not involved in the clustering and partitioning is 'cell' size, which has large fluctuations among each partition bin in this example. By contrast, the gating parameter(s) for each bin is(are) expected to lie somewhere within the gating criteria [b i−1 , b i ] for each bin i after averaging. If cell attributes are all statistically independent, gated renormalization will result in coarse grained, but still statistically independent, attributes; however, if some attributes correlate with one another, the gated renormalization is expected to make this dependence more apparent by reducing biological and measurement noise.
The variance across data renormalized in this fashion is expected to be comparable to that of raw data for gating attributes by construction. The variance in the cellular attributes not involved in the partitioning criteria after gated renormalization will depend on the degree of correlation between these attributes and the gating parameters(s). As an example of this point, in Appendix 1- figure  2b we show a scatter plot of QD and GFP cellular intensities with renormalization performed after sorting in ascending order of fluorescence intensity. This strategy can be generalized to design mesoscopic partitioning of a cell population into multi-dimensional bins characterized by a range of multiple reporters (e.g. by choosing Hoechst and GFP fluorescence intensities as gating parameters).
The gated partition data analysis for identification of the optimal drug combinations is, then, rather straightforward. First, we rely on natural biophysical phenomena to provide a broad range of random mixing in drug uptake. Next, we use the fluorescence intensity of drugs (or their corresponding tags) to partition all cellular attributes into bins. These bins are designed to be sufficiently capacious to contain dozens or hundreds of cells. All cellular attributes are averaged within these bins, which map out the space of drug concentrations and corresponding cellular responses/phenotypes. Finally, we treat data from each bin as an individual culture dish with delivered drug concentrations proportional to fluorescence intensities of the corresponding tags. (According to the Law of Large Numbers, averaging drug response over many cells in each bin suppresses intrinsic noise by the factor √ n i , where n i is the number of cells in bin i ).

Point injection and conjugation-free method
We labeled each drug uniquely with fluorescent tags, and assumed that the amount of drug taken up by each cell scales as the fluorescence of its corresponding tag. The quantification of each cell response is matched to fluorescence data by utilizing tags that are retained in cells. The inherited randomness in cell behavior is not only responsible for variable drug uptake, but also drives response to drug treatment by individual cells. To overcome this source of heterogeneity, we perform averaging over many cells within the sub-population binned by fluorescence intensities of drug tags. There are potential pitfalls in the random mixing method. For example, the uptake of different drugs may be correlated due to the similarity in chemical structures or delivery methods. Cells that take up one drug may, therefore, be likely to take up a similar drug, as well, in proportional fashion. Hence, drug mixing becomes not very effective since we are not exploring the entire drug mixing range by introducing a bias in the mixing process.
Multiple approaches can potentially mitigate this "proportionality" bias. One possibility is to utilize different delivery mechanisms for different drugs. This approach, however, is too cumbersome and does not guarantee bias elimination. Randomness in drug uptake might be driven by similar factors (such as membrane geometry and physical properties), even for different delivery methods.
Here, we developed a method to facilitate independent drug mixing by the creation of local drug gradients across the sample using point injection of individual drugs in different sample locations. The modification of the point injection method can be used for conjugation-free drug labeling by local co-delivery of a pre-mixed drug and its corresponding label. We confirmed that the conjugation-free drug labeling method works well, not only with dyes or small molecule drugs, such as auranofin, but also with much larger nanoparticles, such as QDs.
To assess the predictive power of the conjugation-free method, we compared the CV for each cell subpopulation (bin) used to produce the coarse-grained data shown in Section 6 and Section 6 to the CV of the entire population in a uniform dye delivery experiment. The results of this comparison are shown in Appendix 1- figure 4a and Appendix 1-figure 4b for dyes we treated as 'drugs'. Here, the large red circles represent the CV of the entire cell population in the uniform delivery method for these drugs.
Appendix 1-figure 4. Predictive accuracy in drug uptake using the point co-delivery method. (a) and (b), Predictive accuracy in drug uptake using the point co-delivery method; the coefficient of variation in 'drug' uptake for each bin of co-delivered dyes is shown for pairs {Blue, QD 655 } and {Green, Red}, respectively. Here, the Appendix 1-figure 4 continued on next page CV of 'drug' uptake is plotted against the average dye uptake for each bin comprised of 100 cells characterized by similar dye intensities. Variability in the uniform 'drug'-dye co-delivery experiment across the entire cell population is shown as the large, filled, red circles for each condition above.

Coarse-grain image processing by tiling
In order to utilize stochasticity in drug uptake (or any other cellular phenotype) as a tool to study biological correlations, one has to have enough cells to make statistical analysis feasible. This requirement is essential for two primary reasons. First, a large statistical ensemble allows one to collect enough information about atypical (rare) cellular responses to drug treatment that may nevertheless be critically important (e.g. resistance to treatment). Second, renormalization requires partitioning and subsequent bin averaging to suppress intrinsic noise, which is best achieved using a large ensemble.
While flow cytometry-based experiments can easily be designed to handle large numbers of cells, this is usually not the case for fluorescence microscopy. The experimental protocols, data storage, and analysis are usually non-trivial technically. Confluent cell culture conditions can present a problem for segmentation analysis both in terms of accuracy and run-time. Hence, many tiled optical fields might be necessary to image a large number of cells. Cell segmentation may also require a dedicated optically active biomarker resulting in bandwidth reduction and possible leakage to other channels in high-throughput experiments designed to assess multiple cellular attributes.
To overcome these issues with confocal fluorescence microscopy methodologies, we employed yet another coarse graining method for image analysis. Similar to renormalization of the cell population, microscopy images are renormalized as follows. All channel images are identically partitioned into square tiles of an area large enough to contain at least one cell on average. Fluorescence intensities for each of these tiles are averaged, and the resulting mean intensities are treated as a single-cell attribute. Of course, the data collected in this fashion would include intensities from incomplete cells, fragmented by the tile limits, and also intensities from multiple cells averaged together.
To study these effects on tiled renormalization and assess its quality, we generated synthetic images of geometrically regularized but dense objects (circles), as shown in Appendix 1-figure 3. (Appendix 1-figure 5a) or varying radii (Appendix 1-figure 5b). We assigned to each circle a unique pixel intensity (uniformly across each circle). The images were partitioned into tiles (we will refer to its tile size as a scaling parameter) as illustrated in Appendix 1-figure 5b by the red lines. We calculated average pixel intensity by normalizing total intensity to the foreground area for each tile. To this end, we created a binary image (mask) by separating foreground (circles) and background in the image. Unlike traditional quantification approaches that rely on identification of individual objects and involve separation of touching objects in the image, we used only masks for tile normalization.
Appendix 1-figure 4 continued contribution to the overall single pixel distribution will be A identical values I k where I k is an intensity assigned to the object k . While I k varies from object to object, pixel size A is constant. Hence, a single pixel distribution should be identical to that for objects' intensities { I k }, k = 1, 2, · · · . (Histogram binning may slightly affect the perfect match with an overlap of 1.) Therefore, we expect that as long as fluctuation in object intensity dominates variations in object size, the dependence of the distribution overlap on scaling length should monotonically decrease. This intuition, however, only holds for uniform objects in the image. Here, we synthesized more realistic images by introducing noise across individual circles' intensities. We applied a multiplicative noise function, G ϕ (I) , for each pixel of value I , as follows where a random variable ζ is drawn from a uniform distribution (for example), and the strength of noise is controlled by parameter ϕ .
The results of these simulations are shown in Appendix 1-figure 3. The typical synthetic image and its traditional segmentation are shown in Appendix 1-figure 7a. Here, we used the following statistical parameters to create this image: The comparison of the measured and exact intensity distributions' overlap as a function of scaling parameter is shown in Appendix 1- figure 7b for different values of noise strength ϕ . The red flat line corresponds to the overlap value obtained by comparison of segmented and exact distributions' overlap for ϕ = 0.5 . We used local adaptive segmentation and a smoothing approach to obtain this value.
It is clear from Appendix 1-figure 7b that unlike the case of noise-less (uniform) objects, tile partitioning may have an optimal scaling parameter value. This conclusion is because a consequence of the fact that a single object in the image does not contribute A identical values to the overall distribution in this case. Averaging over finite sized regions suppresses multiplicative noise in the model we introduced. Naively, we expect similar behavior in real microscopy data.
We found that the tile renormalization (and segmentation-free) approach, indeed, results in fairly accurate estimates of the statistical properties of cellular attributes (see Appendix 1- figure  3 for comparison between segmented single cell and renormalized [tiled] image statistics for the phenotypes we tested). We must note that, in principle, accurately determined statistics of each individual random variable in the multidimensional ensemble does not necessarily guarantee a correlative relationship between these variables. For example, some cellular markers (such as QDs) tend to be localized in sub-cellular compartments, and too aggressive image partitioning (i.e. too small a scaling parameter) may result in loss of the correct relationship between these localized markers. By contrast, too large a scaling parameter results in loss of sensitivity and range in the determination of the quantitative relationship between the markers. We, therefore, relied on a comparison of segmented single cell and tiled data to determine the optimal trade-off in the partition strategy empirically (using sparsely seeded control samples where traditional segmentation is possible).
For confluent cell culture conditions (or tissue), the application of tiling renormalization is straightforward. However, when lacunae are present in the culture due either to seeding inhomogeneities or cell death, one has to use masking to prevent averaging over empty space. Masking does not require segmentation but, nevertheless, may require an additional fluorescence marker or phase-contrast imaging to define the approximate limits of the lacunae.

siRNA delivery by quantum dots nano-carriers
After conjugation of streptavidin-QDs with biotinylated siRNA, we further functionalized QDs with an arginine-rich cell-penetrating peptide (CPP) that facilitate cell entry Medintz et al., 2008. The peptide, biotin-labeled (Arg) 9 , was purchased from Anaspec (AS-64078).
As in the case of the diester dye experiments described above, we observed bias towards proportional drug uptake when the QD nanocarriers were delivered uniformly across samples (Appendix 1- figure 9a and Appendix 1-figure 9b), which is, almost certainly, a consequence of identical drug delivery and uptake mechanisms. All of our drugs were siRNAs, and all of them were delivered in exactly the same fashion via QD nanoparticles functionalized by the same CPP. Hence, one expects the uptake kinetics of each drug to be quite similar for each individual cell (the kinetics, of course, varies widely among different cells in the population). Point injection of QDconjugated siRNAs produced a concentration gradient that resulted in much broader range of drug combinations across the samples (Appendix 1-figure 10).

Bin partitioning
We first assign to each cell n a set of d measured attributes ⃗ x n = {x 1 , x 2 , · · · , x d } n where n is an index denoted by a superscript, and the subscript is used to distinguish attributes. The bin partitioning is performed with respect to a subset of the attributes. To be specific, let us define the partition and subsequent subpopulation average based on a first attribute, x 1 .
⟨⃗ x⟩x 1 = ⟨{x 1 , x 2 , · · · , x d }n⟩ x1⊂ [x1, x1+∆x1] (4) The generalization to a higher dimension is straightforward. For example, for pairwise drug combinations, one has to perform binning and averaging over attributes {x i , x j } simultaneously satisfying conditions To decide what mesh sizes ∆x i should be used for each attribute i, 1, · · · d , we applied the following strategy. We set as a goal to partition the underlying parameter space into the mesh in such way that each bin contains a similar number of cells. The average number of cells, L d , in each bin is pre-determined from considerations of noise suppression. The number L d depends on the width of the joint distribution of d attributes comprising the parameter space; we typically chose L d ∼ 30-100 for single-cell data, and L d ∼ 3-10 for tiled image data.
To find a space partition mesh satisfying these criteria, we developed the following approximation. First, we assumed that binning and averaging are performed over independent attributes, such as QD carriers delivered independently by point injections. This very strong assumption reduces the problem into finding a partition in 1d space for each attribute.
For each attribute, x , we quantile the intensity distribution of collected N values of attribute x to ensure equal number of events, L , for each bin. If our assumption holds that all d attributes x are statistically independent, the number of events L d in a d -dimensional partition mesh is given by: Therefore, an estimate for the value L we must use in 1d partitioning is: For example, for L d = 100 and N = 10 4 , the 1d partition should incorporate L ∼ 10 3 events per bin. We optimized the partition scheme for one particular specimen, namely, the control sample.
(Alternatively, all sample data can be pooled, and partitioning can be optimized for the pooled intensity data.) We used this strategy because all samples were normalized to control data and, at least intuitively, we believe that mesh design should perform best for this sample.
In Appendix 1- figure 11, we demonstrate the 1d partition outcome for the same partition scheme applied to seven collected samples simultaneously. Here, we display event numbers L in each bin. Optimal performance corresponds to a line with zero slope (equal event numbers in each partition bin, L = const ). It is entirely possible to find an optimal bin partitioning (mesh) in 2d directly, but we found fairly satisfying results in terms of the number of events per bin using the attribute independence assumption (Appendix 1-figure 12).
Appendix 1-figure 11. Independent application of the 1d partition mesh to QD 605 and QD 655 intensity data collected in 7 different samples shown by different symbols and line styles, as indicated in the left and right panels, respectively. The partitioning was optimized for a control sample (7) and applied uniformly to all seven samples. To eliminate bleeding between channels ('PB450', 'FITC', 'PE', and 'APC') we used data from control samples (single drug/dye). The fitting function for each detection channel was determined after preliminary data processing by moving average filter. The results of the pre-processed data fitting are shown in Appendix 1- figure 15 for compensation due to CellTracker Deep Red dye used as an homogeneous label dye in our experiment. The leakage between other fluorescent signals was less than 5% and was neglected. The compensation was implemented by utilizing a look-up table.