Assessing Combinational Drug Efficacy in Cancer Cells by Using Image-based Dynamic Response Analysis.

The landscape of translational research has been shifting toward drug combination therapies. Pairing of drugs allows for more types of drug interaction with cells. In order to accurately and comprehensively assess combinational drug efficacy, analytical methods capable of recognizing these alternative reactions will be required to prioritize those drug candidates having better chances of delivering appreciable therapeutic benefits. Traditional efficacy measures are primarily based on the "extent" of drug inhibition, which is the percentage of cells being killed after drug exposure. Here, we introduce a second dimension of evaluation criterion, speed of killing, based on a live cell imaging assay. This dynamic response trajectory approach takes advantage of both "extent" and "speed" information and uncovers synergisms that would otherwise be missed, while also generating hypotheses regarding important mechanistic modes of drug action.

Introduction

Numerous efforts have been dedicated to cancer therapeutics over the last few decades, primarily focusing on identifying driver mutations and other target types, whose activities could be severely disrupted via single-target drug inhibitors.1,2 The initiating breakthrough for this type of approach was the development of Gleevec® (Imatinib), a BCR–ABL oncogene inhibitor that achieves 98% complete hematologic response and 86% complete cytogenetic response in patients with early stage chronic myeloid leukemia (CML).3 Though there are some further prominent examples of how patients can benefit from these advances, successes are few and far between. It has become clear that cellular complexity produces at least two mechanisms that can lead to the failure of candidate drugs to deliver the expected therapeutic effects. These are redundant pathways and feedback circuits. These pathways can lead to cancer cells developing either de novo resistance or acquired resistance or both.4 For example, in melanoma, PLX4032 (vemurafenib), a BRAFV600E enzyme inhibitor, has been shown to stimulate RAS activities, even though it suppresses proliferation by inhibiting the MAPK (Ras-Raf-MEK-ERK) pathway.5–7 Melanoma cancer cells can also switch from BRAF-dependence to CRAF-dependence in response to an RAF kinase inhibitor, namely, AZ628.8 Similar acquired resistance is explained by ERBB3–PI3K–AKT activation through MET amplification in non-small cell lung cancer, when treated with gefitinib (an EGFR inhibitor).9,10 The finding of multiple, independently active targets naturally leads translational research to consider combination therapies.11

Drug efficacy is typically evaluated in vitro before testing in animals and humans. A reliable and informative evaluation method for evaluating efficacy in vitro can help produce drug candidates that have better chances of delivering appreciable therapeutic benefits in vivo. A traditional way of measuring a combination drug treatment is to see if it can exceed the capabilities of singleton treatments. To detect such synergistic effects, researchers often assume one of the two widely used models11: if the component drugs in the combination are assumed to exert their inhibition through independent mechanisms, one would use the Bliss independence model12 and expect synergism, if the percentage of survival when treated with combination is less than the product of percentages from individual drugs. On the other hand, if the drugs are assumed to operate in a similar mechanism, one would use the Loewe additivity model, which was first proposed by Loewe et al.13 and later developed into measures such as the combination index14 and its other variations. In this work, we will assume the Bliss independence model as we purposely selected drugs targeting independent molecules (details about the drugs we selected can be found in the “Materials and methods” section).

Drug efficacy can be measured based on the “extent” of their ability in inhibiting cells: the percentage of cells being killed after drug exposure. There is, however, one more critical type of measure that has largely been neglected, the speed of killing. Given two drugs with similar extent of efficacy, a faster working drug may be preferred as it may reduce the chance that cells will develop acquired resistance and minimize the likelihood that the drug is broken down and filtered out of the system before providing maximal killing. In addition, crucial cell behaviors, including cytostasis and cell death, manifest differently in cell response dynamics. But the traditional way of measuring efficacy, where the number of live treated cells is compared to the number of live untreated cells, specifies a way but is unable to attribute it to either inhibition of cell growth or induction of cell death. In fact, dynamics is such a vital element in the study of system biology that analyzing multiple component response trajectories over time after perturbation has long been an indispensable tool for gaining insight into complex cell functions and machinery.

Here, we will introduce how a live cell imaging-based assay can exploit information about temporal responses to evaluate the efficacy of drug combinations. The response dynamics between single and combination drug therapies can reveal important modes of mechanism and help uncover potent drug combinations that would otherwise be missed.

Results

We tested four different drugs, MLN9708, SH-4-54, PIK-75, and KPT-330, and their combinations, pre-selected for their critical inhibition capabilities on distinct and independent molecules. (For better readability and more clarity, especially in figures and equations, wherever appropriate, these drugs are coded as drug A, B, C, and D, respectively.) We applied them on ten different canine cell lines and the images of cells were captured 46 times, once in an hour. We processed the images and extracted relevant morphological features for classifying dead/alive cells. Details of these procedures can be found in the “Materials and methods” section.

Limitation of traditional “end point only” synergism detection approach

In order to assess the effectiveness of a drug, we define the killing efficacy (KE) of a drug as the percentage of cells classified as dead after treatment. We use the notation to indicate KE for a certain drug X at a given time t. This is calculated by dividing the number of cells classified as dead by the total number of cells.

When two drugs X and Y are assumed to work through independent mechanisms, and therefore follow the Bliss independence model, we declare synergism, if the statistic is significantly positive (one-tailed Welch’s t-test with the null hypothesis that Equation (1) equals 0). We term this as a single time point synergism test.

We now show the limitation of the traditional “end-point only” approach for detecting synergism of drug combinations. KEs are calculated for drugs MLN9708 (A) and SH-4-54 (B) and their combination at time point t = 46 (45 hours after treatment) for cell line ABRAMS:

There is no synergism using the single time point synergism test, because , which is hardly >0, and the significance test fail to reject the null hypothesis at the 0.05 level using all duplicated wells and sites. This should not be surprising, since MLN9708 alone delivers almost the maximal killing efficacy and there seems no reason to add SH-4-54 to the mixture. However, we will shortly see that this is entirely due to the fact that we only tested for synergism at a single time point.

Comparing response trajectories can discover strong synergism not revealed by a non-dynamic approach

We will now calculate KEs for the full 46 hours and plot them over time in Figure 1. It is obvious that though SH-4-54 itself is not a potent inhibitor, it sensitized the cells, so that the combination works in an extremely efficient manner and accomplishes maximal killing much sooner. As discussed earlier, this speed is a very desirable property of drugs.

Figure 1

KE over time for drugs MLN9708 (black), SH-4-54 (blue), and a combination of MLN9708 + SH-4-54 (red).

By considering KEs over a period of time, we define combinational drug synergism, based on the single time point synergism test, relative to two drugs X and Y by where T is the last time point of a given experiment and I = 1, if the single time point synergism test is passed at time point t and I = 0 otherwise. This definition is not invariant to imaging frequency or experiment end point T, and so comparison for comboKEX,Y should only be carried out for experiments over the same time span. Furthermore, if one experiment is sampled with higher frequency, then one may wish to down-sample it by taking averages of for t in the lower frequency experiment’s each time interval.

Figure 2 shows the KE over time across different cell lines: ABRAMS, BKOS, and CML1. As comboKEX,Y is nonnegative by definition in Equation (3), it is evident from Figure 2B that coupling MLN9708 with PIK-75 generates zero synergism (it can be actually tested for antagonism, but that is not our interest here), which serves as strong evidence that the combination of MLN9708 and PIK-75 can be removed from consideration. Figure 2A, 2C, and 2D shows that more synergism between MLN9708 and SH-4-54 is observed for ABRAMS, followed by CML1, and much less on BKOS.

Figure 2

KE over time (hours on the x-axis) for drugs MLN9708 (A; black), SH-4-54 (B; blue), and MLN9708 + SH-4-54 combination (red), except in (B) which has drug PIK-75 (C; blue) and combination A + C (red). The green dashed line shows the “synergistic points” calculated as for each t ( replaced by in (B)).

Note: comboKEX,Y equals the area shadowed by yellow vertical lines.

Ranking drug efficacy using a two-dimensional drug efficacy measuring technique

It is a common practice to rank drug efficacies, so that a candidate drug or drug combination can be prioritized for testing in the next phase. Having drug response dynamics at our disposal gives us two-dimensional drug efficacy parameters: extent and speed, denoted by (E, S). Inspired by the resemblance of the response trajectory to the engineering receiver operating characteristic curve, we define the efficacy of drug X as the total area under the curve (AUC): AUCX. It is straightforward to see that AUCX ∊ [0, 1], AUCX = 0 indicates that drug X manages absolutely no killing in any of the time points tested, and AUCX = 1 boasts about a perfect killer that accomplishes its task in no time.

As an example, AUCX is calculated for Figure 2A and 2B for ABRAMS and ranked high to low as follows:

This again shows the potency of the combination MLN9708 + SH-4-54 (A + B). In practice, this ranking can be limited to drugs that have achieved a predefined minimum KE at the end time point of the experiment (KET). It should be noted that we assume that cell death is irreversible. Hence, the response curves are monotonic: KE ≥ KE for any t < T.

There are scenarios when one may wish to decompose AUC values and take a look at (E, S) directly. To that end, we compute the first derivative of the trajectory curves and find the time point T when the response slows down or plateaus (or reaches the end; details can be found in the “Materials and methods” section). Some examples are found in Figure 3, where the red vertical dashed line marks T.

Figure 3

Response curves with T (when the response slows down or plateaus or reaches the end) marked as a vertical dashed line on ABRAMS for drugs (A) MLN9708, (b) SH-4-54, (C) PIK-75, (D) MLN9708 + SH-4-54, (E) MLN9708 + PIK-75, and (f) SH-4-54 + PIK-75. Hours are plotted on the x-axis.

We will estimate (E, S) using (KE, T). This is especially useful if we wish to prioritize drugs or combinations using data from multiple cell lines. Figure 4 shows (KE, T) pairs for all four canine malignant melanoma cell lines, CML1, CML6M, JONES, and PARKS; six canine osteosarcoma cell lines, ABRAMS, BKOS, BOOZA, MCKOS, SKOS, and UWOS2 (all cell lines for the same drug share the same marker in the graph for clarity); and all drugs tested. Figure 4B and 4C shows the zoomed-in top and right-hand side regions of the same plot in Figure 4A, respectively. Drugs in the left-hand side and top portions of the 2D scatter plot kill more cancer cells in a shorter period of time. Again, the drug combination MLN9708 + SH-4-54 resides in the region. It should be noted that SH-4-54, although it can have very limited killing effect on its own, acts very quickly. This partly explains the sensitizing effect when used in junction with drug MLN9708.

Figure 4

(T, ) pairs, with T on the x-axis and on the y-axis, for all four canine malignant melanoma cell lines, CML1, CML6M, JONES, and PARKS, and six canine osteosarcoma cell lines, ABRAMS, BKOS, BOOZA, MCKOS, SKOS, and UWOS2. Parts (b) and (C) show the zoomed-in top and right regions of the same plot in (A), respectively.

Discussion

We have shown that the efficacy of drugs and drug combinations can be more thoroughly evaluated using a live cell imaging-based dynamic response trajectory approach. As the landscape of translational research has been shifting to combination therapies,11 our aim here is to develop a method that prioritizes drug combinations that kill cancer cells in a more complete and efficient fashion. The addition of time-line information is crucial in uncovering synergism that will otherwise be missed by the traditional two-point system approach (taking measurements at the start and end of an experiment only).

It should be noted that even though we only tested a combination of two drugs in this work, our method could be directly applied (eg, in Figs. 1 and 3) or easily extended (eg, to Fig. 2) to combinations of more than two. For example, if we desired to evaluate any possible synergistic effects of a 3-drug (X, Y, and Z) combination over the singletons, Equations (1) and (3) become:

Furthermore, the therapeutic agents testable in combinations are not limited to chemical compounds. For example, immunotherapeutic cells can also be tested using the tools outlined in this work, prioritized combinations can be administered, and their effectiveness can be monitored, thanks to recent advances in real-time imaging in immunotherapy.15,16

Testing the behavior of actin cytoskeletal elements has been undertaken as a result of the accelerated kill rate response of MLN9708 after being combined with SH-4-54 (Fig. 1). This acceleration suggests the possibility that SH-4-54 might be eroding the integrity of some fundamental cell structures. This was confirmed by comparing actin cytoskeletal filaments in cells treated with MLN9708 or SH-4-54 alone, or both. From Figure 5, it is clear that the cytoskeletal structure is compromised for cells treated with SH-4-54 or with the combination, but not for those treated with MLN9708 only, and this occurs in concurrence with pyknosis and death of the cells. MLN9708 alone has about the same killing efficacy at 13 hours as that of SH-4-54 (Fig. 1); however, cells have dramatically different cytoskeletal structural morphologies (Figs. 5B and 5C), indicating that two different modes of mechanism are involved when these two drugs exert their individual inhibitory powers. Imaging the dynamics of various behaviors of cellular components in many individual cells over time in response to differing drugs is a way to analyze the systems active in cells. Recent advances in the development of reporters for a wide variety of cellular components make it possible to more fully interrogate cellular activities from a systems viewpoint.

Figure 5

Actin cytoskeletal filament images for ABRAMS at 13 hours for (A) untreated, (b) treated with MLN9708 only, (C) treated with SH-4-54 only, and (D) treated with MLN9708 + SH-4-54 combination.

Materials and Methods

Selection of drugs

Two proteins that are highly involved in generating prosurvival activities in a wide variety of cancer cells are Mcl1 and Stat3. Mcl1 binds with the proapoptotic Bax and Bak1 proteins, keeping them from forming multimers that breach the outer membranes of mitochondria, initiating apoptotic cell death. A number of drugs suppressing Mcl1 production or activity have been produced. The drug PIK-75 targets the suppression of Mcl1 production by inhibiting the ability of the Cdk9 kinase to phosphorylate RNA polymerase II, which is required for transcription of Mcl1.17 The drug MLN9708 (Ixazomib, Takeda Oncology) is currently in clinical use to reduce prosurvival activity of Mcl1 in cancer cells. This drug stimulates the production of Pmaip1 (NoxA) protein by inhibiting the chymotrypsin-like proteolytic (β5) site of the 20S proteasome, which ultimately stimulates transcription and translation of the PMAIP1 gene. Pmaip1 protein binds avidly to Mcl1 protein and prevents it from binding with the proapoptotic Bax and Bak1 proteins, allowing apoptosis to proceed. Although useful, its limitations restrict its approval to be used in combination with other drugs.18,19

Stat3 protein has long been studied as a likely target for cancer drug development.20 Drugs designed to antagonize Stat3 activity have mostly been designed to target kinases that phosphorylate Stat3 at a position that allows the phosphorylated protein to form homodimers. These types of drugs have not been very successful, which is likely due to several different kinases known to be active on Stat3 and probably more that have not yet been identified. More recently, two Stat3 antagonists with differing modes of repression have been developed. An inhibitor of CREB-binding protein-mediated Stat3 acetylation, KPT-330, effectively blocks Stat3 transcription of Mcl1.21 More recently, a Stat3 antagonist, SH-4-54, was developed that binds with phosphorylated Stat3 in a way that prevents homodimerization.22 This type of inhibition is more specific in its action and less likely to be confounded by alternate interactions.

A summary of these drugs is provided in Table 1.

Table 1

Drugs selected and their code names and description. All drugs were purchased from Selleckchem.com.

CODE NAME	DRUG NAME	SELLECKCHEM CAT. NO.	BRIEF DESCRIPTION
A	MLN9708	S2181	Inhibitor of chymotrypsin-like proteolytic (β5) site of the 20S proteasome
B	SH-4-54	S7337	Inhibitor of STAT3/5
C	PIK-75	S1205	Inhibitor of p110α (PI3k) and P-TEFb (CDK9)
D	KPT-330	S7252	Inhibitor of XPO1 (CRM1)

Cell culture

We have imaged ten different canine cell lines kindly provided by Dr. H.M. Wilson-Roble. Canine osteosarcoma cells (ABRAMS, BKOS, BOOZA, MCKOS, SKOS, and UWOS2) and canine malignant melanoma cells (CML1, CML6M, JONES, and PARKS) were grown in imaging medium (IM) having low levels of autofluorescence and were maintained at 37 °C with 5% CO2. IM contains 70% Medium 199 (Thermo Fisher Scientific, 11825015), 30% Roswell Park Memorial Institute-1640 medium (11875085) supplemented with 10% fetal bovine serum (16000044), 20 mM Hepes (15630080), 20 mM Glutamax (35050061), and 11 mM D-glucose (A2494001).

Sample preparation

Cells of each type were plated at a density of 2,000 cells/well in 30 µL/well of IM on a 384-well microtiter plate (Greiner Bio-One 781091) pre-coated with 10 µg/mL rat tail collagen type I (BD Biosciences 354249) and incubated overnight at 37 °C with 5% CO2. Prior to taking baseline time points, nuclei were stained with 0.00667% Vybrant® DyeCycle™ violet stain (ThermoFisher Scientific V35003), which is a live cell permeable stain and produces blue fluorescence (<437 nm) when bound to double-stranded DNA and stimulated with a violet excitation source (<369 nm). Each well also contained 0.0125% CellTox™ green (CTG) dye for measuring cell death (Promega Inc. G8742). This cyanine dye preferentially stains the dead cell DNA when cell membrane integrity is compromised, producing bright green fluorescent signal (<509 nm) when stimulated with a 488 nm excitation source.

Treatment and imaging

After taking baseline time points (typically the first one or two time points), cells were treated with different drugs or drug combinations, and each well was scanned every hour for 1 or 2 days using an Image Xpress Micro XLS Widefield High-Content Analysis System (Molecular Devices) sampling at three different imaging sites within each well.

An example composite image is given in Figure 6A for a baseline time point and Figure 6B for 44 hours after being treated with drugs; the full-time sequence images can be found in the Supplementary Files (S01.avi for the nucleus channel and S02.avi for the CTG channel). Phase-contrast cell images were also captured to provide information about overall cell-structure changes over time. A plug-in tool based on ImageJ (National Institutes of Health) was developed that allows users to more efficiently view and manage the large number of images collected. The interface of the tool is demonstrated in Supplementary Figure S03.pdf.

Figure 6

Sample composite images of nuclei (blue) and CTG (green), showing (A) early time point before being treated with drugs and (b) a time point about 45 hours after being treated with drugs.

Image processing

The images were sent into an in-house pipeline developed in Matlab (Mathworks) using the SDC morphological toolbox.23 After image segmentation on both nuclear channel and green fluorescent channel,24 we collected a large number of measurements for each individual cell. Details on image segmentation method were previously reported.24,25 As our primary goal here is to assess the ability of various drugs and drug combinations in killing cancer cells, we used the following morphological features on each cell in dead/alive labeling (the rationales are provided in the following sections):

Nucleus size (Snuk)

Nucleus mean intensity (µInuk)

CTG mean intensity in nucleus (µICTG/nuk)

In addition, nuclei boundaries are also saved and can be easily overlaid on the original images using the ImageJ plug-in tool.

Identifying apoptotic cells morphologically

Debates are still ongoing about the most reliable cell morphological manifestations for classifying dead cells versus live cells.26,27 However, for our practical purpose, we can safely classify a cell as dead if its membrane is compromised, at which time the CTG dye in the media can enter the cell and bind to the DNA and produce a significant signal when compared to intact cells. Therefore, we can classify a cell as dead when where µI0,CTG/nuk is the threshold for the CTG mean intensity in the nucleus. It should be noted that we used the intensity in the nucleus only (as opposed to the intensity in the cell), as this yields the most distinction between the compromised membrane and intact cells.

Pyknosis, or condensation of chromatin in nucleus, is typically a hallmark of apoptosis and is irreversible. As such, the start of the process can be easily detected on nuclear channel images when the nuclei become much smaller and brighter, as determined by another classification rule, which is as follows: where S0,nuk and µI0,nuk are the thresholds for Snuk and µInuk, respectively.

It has been observed that pyknosis measured by nucleus intensity and size well precedes the loss of cell membrane integrity as detected by CTG and the exact time depending on the cell lines and drugs used. Figure 7 shows a couple of examples (the complete set of comparisons, including multiple cell lines and drugs, can be found in the Supplementary File S04.pdf). In Figure 7, blue lines are the percentages of cell death based solely on µICTG/nuk and black lines indicate the percentages of cell death based solely on nucleus intensity µInuk. The peak activities measured by either method are separated by more than 24 hours. This makes CTG unsuitable for early death detection. On the other hand, as drug-responsive cells move deeper into apoptosis and experience structural breakdown, the brightness observed at the beginning of pyknosis starts fading away and the nuclear boundary becomes harder to detect. This intensity deficit makes measuring of nuclear size and intensity unreliable for apoptotic cells in late stages. To optimize live/dead classifying accuracy late in response, we combined Equations (4) and (5). Our final classification rule for dead cells is as follows:

Figure 7

Blue lines are the percentages of cell death (y-axis) based solely on µICTG/nuk and black lines solely on nucleus intensity µInuk. The peak activities measured by either method are separated by more than 24 hours. Hours are plotted on the x-axis.

Selection of morphological thresholds

As discussed earlier, the determination of µI0,CTG/nuk is rather trivial as the CTG signals are well separated in intensities for dead and live cells. To find S0,nuk and µI0,nuk, we generally select cells in untreated control wells, and cells treated with drugs at time points when pyknosis is expected to be at its peak. Histograms on µInuk or Snuk for these cells typically show a bimodal distribution, which can be fitted with a two-class Gaussian mixture model. Here, we assume nonpooled variances (the two mixtures have different variances), and the threshold is taken so that the Bayes error is minimized. Two histograms are shown in Figure 8A and B for cell lines PARKS and CML1, respectively.

Figure 8

Sample histogram plots of nucleus mean intensity for cell lines PARKS (A) and CML1 (b). Values are in log2 scale on the x-axis. Groups of normal cells and cells in pyknosis in (b) are less distinguished when compared to (A).

If a more precise threshold is desired, as maybe in the case of Figure 8B, the threshold is refined by a secondary step by using a semi-automatic selection tool, which we have developed in ImageJ (CellTagger):

Identify cells close to decision boundaries as determined by histograms and mixture Gaussian fit;

Run CellTagger. The above identified cells are marked with yellow outlines (an example is shown in Figure 9A; outlines are obtained from the saved nuclei boundaries);

Figure 9

(A) Example nucleus images. Cells near decision boundaries are marked with yellow outlines. Only a portion of the full-sized image that contains these marked cells is shown. (b) Same image as in (A) after the user has marked dead cells (red), live cells (green), or undecided (remaining yellow outlined).

As the example shows in Figure 9B, an experienced molecular biologist will go through all these cells and mark them in red (dead cells), green (live cells), or leave them as yellow (undecided/ignored);

Repeat this for all the images in the sequence (example in Supplementary File S05.avi);

Collect all new classification labels (dead or live) and refine the thresholds.

After all thresholds are established, we use a rule in Equation (6) and for each cell line and drug, we compute the percentage of dead cells (%Apop) as

The KE of a drug is calculated as average %Apop over all imaging sites for all replicated wells.

Validating classification

To evaluate the thresholds selected and classification accuracy, we use the same ImageJ tool CellTagger after selecting all cells in the image (that is, all cells except those touching image borders, since they will have imprecise size information) and manually assign dead or alive labels by experienced molecular biologists. Due to the large amount of data, this validation process has been carried out randomly, and Figure 10 shows a typical comparison between the percentages calculated from our computer approach and from our two biologists. The decisions are quite comparable.

Figure 10

A typical comparison between percentages calculated from our computer approach and by our two biologists.

As one of the most sought-after properties of candidate drug or drug combination in the ability to kill cancer cells completely, we are particularly critical to the accuracy of calling at or close to 100% KE. To that end, we carried out a validation experiment, where we removed the drugs at the end of experiment run and watched whether the remaining cells, if any, started to grow again in 2 days. For drug combination MLN9708 + SH-4-54 on ABRAMS, where all cells are calculated as dead (Fig. 1), the dead cells are washed away (Supplementary Fig. S06B and E) and none of the remaining cells grew back in 2 days (Supplementary Fig. S06C and F). As a comparison, there were plenty of live SH-4-54-only treated ABRAMS cells, 48 hours after SH-4-54 was removed (Supplementary Fig. S07).

Finding time point T

The time point T is defined as the time point when the response starts to slow down or plateau, or if neither happens, the end time point of the experiment. This is calculated as the time point when the first derivative of the response curve first falls below ε and stays below ε for a minimum of τ time points. In our experiment, we chose the following values: ε = 0.5% and τ = 5.