Leveraging Mathematical Modeling to Quantify Pharmacokinetic and Pharmacodynamic Pathways: Equivalent Dose Metric

Treatment response assays are often summarized by sigmoidal functions comparing cell survival at a single timepoint to applied drug concentration. This approach has a limited biophysical basis, thereby reducing the biological insight gained from such analysis. In particular, drug pharmacokinetic and pharmacodynamic (PK/PD) properties are overlooked in developing treatment response assays, and the accompanying summary statistics conflate these processes. Here, we utilize mathematical modeling to decouple and quantify PK/PD pathways. We experimentally modulate specific pathways with small molecule inhibitors and filter the results with mechanistic mathematical models to obtain quantitative measures of those pathways. Specifically, we investigate the response of cells to time-varying doxorubicin treatments, modulating doxorubicin pharmacology with small molecules that inhibit doxorubicin efflux from cells and DNA repair pathways. We highlight the practical utility of this approach through proposal of the “equivalent dose metric.” This metric, derived from a mechanistic PK/PD model, provides a biophysically-based measure of drug effect. We define equivalent dose as the functional concentration of drug that is bound to the nucleus following therapy. This metric can be used to quantify drivers of treatment response and potentially guide dosing of combination therapies. We leverage the equivalent dose metric to quantify the specific intracellular effects of these small molecule inhibitors using population-scale measurements, and to compare treatment response in cell lines differing in expression of drug efflux pumps. More generally, this approach can be leveraged to quantify the effects of various pharmaceutical and biologic perturbations on treatment response.

Introduction

The parameterization of in vitro treatment response data is central to biomarker and drug discovery and the quantitative study of cancer therapies. With recent exceptions (Hafner et al., 2016; Harris et al., 2016), investigation of treatment response in vitro has been limited to cell survival assays that assess cell viability at a single, specified timepoint following treatment with a temporally constant concentration of drug. A range of drug concentrations are evaluated in these assays, and the results are conventionally summarized by Hill function parameters, which quantify cell survival with respect to applied drug concentration (Fallahi-Sichani et al., 2013). While this approach has yielded significant insights into cancer biology, it is fundamentally limited by the coarseness of parameters used to summarize treatment response. In particular, these parameters do not explicitly characterize the dynamics of treatment and subsequent response. Further, response metrics are reported relative to the extracellular concentration of drug in the assay, overlooking drug exposure times and variable cell line pharmacologic properties. This not only impairs analysis of in vitro treatment response data, but also presents a challenge in translating these therapies in vivo.

There are a host of biochemical processes that modulate a tumor cell's response to therapy. For example, the accumulation of drug within cells can be altered by drug metabolism or modification of surface proteins that regulate drug flux through the membrane (Larsen and Skladanowski, 1998; Larsen et al., 2000). Indeed, the multi-drug resistance protein 1 (MDR1) is a well-studied mechanism of resistance to cytotoxic therapies (Clarke et al., 2005). This ATP-dependent pump actively effluxes drug from cells, decreases drug accumulation within cells, and confers resistance to anthracyclines, taxanes, and several other agents (Mechetner et al., 1998). Similarly, pharmacodynamic response to therapies can be altered through modulation of signaling pathways downstream of the therapeutic target. With respect to DNA-damaging agents, changes in DNA repair pathways, which are activated in response to treatment, can alter sensitivity to those agents (Fink et al., 1998; Bouwman and Jonkers, 2012). For example, DNA-dependent protein kinase (DNA-PK) plays a major role in the repair of double strand DNA breaks via non-homologous end joining (Smith and Jackson, 1999). Increased expression of DNA-PK has been shown to confer resistance to doxorubicin, an anthracycline commonly used clinically (Shen et al., 1998). Fundamentally, cell line-specific pharmacokinetic and pharmacodynamic properties, such as those described above, drive observed treatment responses. Using conventional methods, these processes are conflated by the parameters used to summarize in vitro dose response data (Prentice, 1976; Fallahi-Sichani et al., 2013). The resulting parameters are imprecise measures of drug efficacy, which limits the biological insights to be gained from the data.

More precise technologies are required to advance systems approaches to studying cellular response to therapy (Anderson and Quaranta, 2008). We posit that a mechanistic, mathematical modeling framework is essential to maximize the knowledge gained through treatment response studies (Yankeelov et al., 2013, 2015). In this paradigm, biologically-motived mathematical models are constructed to describe observed behaviors of the system under investigation. The model is then fit to experimental data, yielding a set of parameter values that provide mechanistic insight into observed data. There exist several models in the literature that explicitly incorporate drug pharmacokinetics (PK) and pharmacodynamics (PD) to describe treatment response. In vitro, transit compartment models have been used to describe the temporal relationship between drug application and effects (Lobo and Balthasar, 2002). More biologically-motivated PK/PD models have been employed to study specific pharmacokinetic and pharmacodynamic parameters (Lankelma et al., 2003, 2013). PK/PD models have also been developed to investigate treatment response in vivo (Simeoni et al., 2004; Sanga et al., 2006; Wang et al., 2015). Recently, we proposed and validated a coupled PK/PD model of doxorubicin treatment response in vitro (McKenna et al., 2017). The model incorporates measured doxorubicin pharmacokinetics and pharmacodynamics and predicts response to a specified treatment timecourse on a cell line-specific basis. The model behaves consistently across a wide spectrum of treatment protocols and cell lines, thereby demonstrating that the response dynamics of cancer cell lines to doxorubicin is predictable within this framework. Specifically, the PK model-estimated concentration of doxorubicin bound to the cell nucleus is predictive of cell line pharmacodynamic rates. We further noted a mismatch of drug uptake and response among the investigated cell lines, suggesting that each cell line has an intrinsic sensitivity to stress induced by doxorubicin. By explicitly modeling both drug uptake and subsequent effect, these processes can be independently quantified to study each component of treatment response (McKenna et al., 2017).

It is the goal of the present effort to demonstrate the utility of a mechanistic, mathematical modeling framework in quantifying treatment response and PK/PD pathways. We leverage mathematical models to filter experimental data to yield quantitative measures of specific cellular processes. Specifically, we experimentally perturb doxorubicin pharmacokinetics and pharmacodynamics with chemical inhibitors of each process. We modulate pharmacokinetics in an MDR1 over-expressing cell line and modulate pharmacodynamics via DNA-PK in a BRCA1-mutated cell line. These data are analyzed with the proposed PK/PD model to yield quantitative measures of these pathways. We further illustrate the utility of our approach by proposing the equivalent dose metric, which we derived from the PK model. The equivalent dose is analogous to that in radiation therapy, which is used to compare radiation fractionation schedules (Fowler, 1992). In the context of chemotherapy, we define equivalent dose as a functional measure of drug exposure. We specify that for a given equivalent dose, treatment response dynamics are similar. As this approach accounts for variable pharmacologic properties, we posit that it allows for more precise comparisons among cell lines relative to metrics based on extracellular drug concentration. We demonstrate this experimentally through comparison of treatment response in cell lines differing only in MDR1 expression. The modeling-based framework proposed in this work can be leveraged to more precisely quantify the effects of various pharmaceutical and biologic perturbations on treatment response.

Materials and Methods

Mathematical Model of Doxorubicin Treatment Response

Doxorubicin is an anthracycline that remains standard-of-care therapy for several cancers (Tacar et al., 2013). Ultimately, doxorubicin induces a host of cellular stress responses which either inhibit further DNA synthesis allowing for cellular recovery, or initiate a cascade leading to cell death (Gewirtz, 1999). At high doxorubicin concentrations, extensive DNA damage often results in cell death via apoptosis. Low to moderate concentrations of doxorubicin induce cell senescence and cell death via mitotic catastrophe (Chang et al., 1999; Eom et al., 2005). Whereas, apoptosis is immediate (on the order of hours to days), mitotic catastrophe is a relatively protracted process (on the order of several days).

We previously developed and validated a parsimonious treatment response model to describe doxorubicin pharmacokinetics and pharmacodynamics (McKenna et al., 2017). Briefly, a three-compartment model was employed to describe the uptake and binding of doxorubicin in cancer cells. This process is modeled via mass conservation through Equations (1–3):

$\begin{array}{ll}\frac{d{C}_E(t)}{dt}=k_{FE}\frac{v_I}{v_E}{C}_F(t)-k_{EF}{C}_E(t)& (1)\end{array}$

$\begin{array}{ll}\frac{d{C}_F(t)}{dt}=k_{EF}\frac{v_E}{v_I}{C}_E(t)-k_{FE}{C}_F(t)-k_{FB}{C}_F(t)& (2)\end{array}$

$\begin{array}{ll}\frac{d{C}_B(t)}{dt}=k_{FB}{C}_F(t)& (3)\end{array}$

where C_E (t), C_F (t), and C_B (t) are the concentrations of doxorubicin in the extracellular, free, and bound compartments, respectively, at time t. The free compartment represents drug that has diffused into the cell, while the bound compartment represents drug that has bound to DNA. The k_ij parameters are rate constants that describe the movement of doxorubicin between the i^th and j^th compartments; for example, k_FE describes the rate of drug transfer from the free, intracellular compartment to the extracellular compartment. Similar definitions apply to k_EF and k_FB. The volumes of the extracellular and intracellular compartments are denoted by v_I and v_E, respectively (see Table 1 for a full list of model parameter definitions). We note that in this model, several intracellular processes, including doxorubicin metabolism and dissociation from DNA, are not explicitly considered. In previous work (McKenna et al., 2017), we evaluated the performance of several candidate models in describing our experimental data (described in section Doxorubicin Uptake Imaging and Image Processing) with the Akaike information criterion. Of the proposed models, we found that Equations (1–3) best balanced accuracy and the number of model parameters.

TABLE 1

Table 1. Model parameter definitions.

A logistic growth model, Equation (4), modified by either of two empirical time-dependent response functions, Equations (5, 6), reflecting distinct mechanisms of cell death, was proposed to describe population level response to doxorubicin therapy as follows:

$\begin{array}{ll}\frac{d{N}_{TC}(t)}{dt}=\left(k_p-k_d(t,D)\right)N_{TC}(t)\left(1-\frac{N_{TC}(t)}{\theta (D)}\right)& (4)\end{array}$

$\begin{array}{ll}{k}_d\left(t,D\right)=\{\begin{array}{cc}0& t<0\\ k_{d,a}\left(D\right)& t\ge 0\end{array}& (5)\end{array}$

$\begin{array}{ll}{k}_d\left(t,D\right)=\{\begin{array}{ll}\hfill 0\hfill & t<0\\ k_{d,b}\left(D\right)r\left(D\right)t{e}^{1-r\left(D\right)t}& t\ge 0\end{array}& (6)\end{array}$

where k_p and k_d are the proliferation and dose-specific death rates, respectively, D is the delivered dose [defined to be the bound concentration of drug, C_B, calculated with Equations (1–3)], r is a dose-specific constant describing the rate at which treatment induces an effect, θ is the dose-specific carrying capacity describing the maximum number of cells that can be observed in the experimental system, and N_TC(t) is the number of cells at time t. Logistic growth models have traditionally been used to describe growth of a variety of biological species whose total size is limited (Gerlee, 2013; Jarrett et al., 2018). This equation accurately describes our experimental system (described in section Doxorubicin Treatment Response Imaging), in which cell population is limited by the surface area of the experimental platform. Prior to treatment (i.e., t < 0), cells are modeled to have a constant proliferation rate, k_p. Following treatment at t = 0, Equation (5) assumes an immediate induction of a stable, post-treatment death rate (k_{d, a}). Equation (6) allows for a smooth induction of drug effect following treatment to a maximum death rate of k_{d, b}, while ultimately allowing for recovery of the cell population. The dynamics of this induction and decay is governed by r. A weighted averaging approach is used to incorporate both Equations (5, 6) in the treatment response model. This model was designed to describe cell death via apoptosis Equation (5) and mitotic catastrophe Equation (6) and fit experimental data well. Further details on the model can be found in McKenna et al. (2017).

Equivalent Dose

We define equivalent dose (D_eq) as a functional measure of therapy, a summary statistic connecting the amount of drug delivered with the biological effect of that drug. In the context of doxorubicin therapy, we define equivalent dose as the functional concentration of drug that is bound to the nucleus following therapy. To calculate D_eq for a specified treatment condition (i.e., extracellular drug concentration timecourse), Equations (1–3) are populated by cell-line- and treatment-specific k_EF, k_FE, and k_FB parameters derived from experimental data (described below). The model is then simulated using the experimentally-defined treatment condition. D_eq is the maximum concentration of bound drug (C_B) as predicted by the simulation. We hypothesize that the equivalent dose metric can account for variable cell line pharmacologic properties through its explicit consideration of k_EF, k_FE, and k_FB rates, and it can be leveraged to quantify the effect of agents that modulate those properties. Notably, in proposing the equivalent dose, we lump pharmacodynamic effects into the k_FB term. Specifically, k_FB is a mixed measure of doxorubicin binding and DNA repair and describes the functional net binding rate. The equivalent dose is illustrated in Figure 1.

FIGURE 1

Figure 1. Overview of equivalent dose metric. The response to therapy is determined by the applied drug concentration and cell-specific pharmacologic properties. Traditionally, therapeutic response is summarized relative to the applied extracellular concentration of drug. We propose the equivalent dose metric, D_eq, to summarize the contributions of various pharmacologic properties in shaping treatment response. We define equivalent dose as a measure of the functional drug concentration that enters the cell. The equivalent dose is calculated through a mechanistic biophysical model that considers several sources of variability in shaping treatment response. The metric consolidates variable drug uptake (quantified with k_EF), efflux (quantified with k_FE), and binding (quantified with k_FB) into a single descriptor of treatment. The equivalent dose summarizes pharmacologic properties to provide biological insight into treatment response and allows for more precise comparison of treatment response among cell lines.

Cell Lines

The MDA-MB-468 and SUM-149PT cell lines were obtained through American Type Culture Collection (ATCC, http://www.atcc.org) and maintained in culture according to ATCC recommendations. Cell lines were passaged no more than 30 times before being discarded. To facilitate automated image analysis for identifying and quantifying individual nuclei in time-lapsed microscopy experiments (described below), each cell line was modified to express a histone H2B conjugated to monomeric red fluorescent protein (H2BmRFP; Addgene Plasmid 18982) as previously described (Quaranta et al., 2009; Tyson et al., 2012).

To specifically modulate doxorubicin pharmacokinetics, the H2BmRFP-expressing MDA-MB-468 cell line (MDA-MB-468_H2B) was transduced to express a green fluorescent protein (GFP)-tagged MDR1 protein (ABCB1 gene, Origene Technologies, Rockville, MD). Following transduction, the cell line was cultured in 100 nM doxorubicin for 2 weeks to select a doxorubicin-resistant phenotype (MDA-MB-468_MDR1). These cells were serially imaged to ensure that all surviving cells stably expressed GFP.

The SUM-149PT cell line possesses a BRCA1 2288delT mutation (Elstrodt et al., 2006). BRCA1 is involved in maintaining genome stability through its role in repairing double strand DNA-breaks via homologous recombination (Gudmundsdottir and Ashworth, 2006). The BRCA1 mutation causes an increased reliance on alternate DNA damage repair pathways, such as non-homologous end joining (Farmer et al., 2005). The DNA damage repair pathway mediated by DNA-PK was targeted with a small molecule inhibitor to specifically modulate doxorubicin pharmacodynamics in the SUM-149PT cell line.

Chemicals

Doxorubicin was purchased from Sigma Aldrich (St. Louis, MO) and dissolved to a 1 mM stock concentration in sterile saline for subsequent experiments. Tariquidar (TQR) is a third-generation MDR1 inhibitor that non-competitively inhibits MDR1 function (Mistry et al., 2001). TQR is leveraged to modulate doxorubicin pharmacokinetics in the MDA-MB-468_MDR1 cell line. NU7441 is a DNA-PK inhibitor that has been investigated as a means to improve treatment response to DNA-damaging agents (Zhao et al., 2006; Ciszewski et al., 2014). NU7441 is used to modulate doxorubicin pharmacodynamics in the SUM-149PT cell line. TQR and NU7441 were both purchased from Selleckchem (Boston, MA). Each was dissolved to a 1 mM stock concentration in DMSO. We subsequently refer to these therapies (TQR and NU7441) as sensitizers. All solutions were stored in 250 μL aliquots at -80°C.

Doxorubicin Uptake Imaging and Image Processing

Time resolved fluorescent microscopy was employed to characterize the uptake of doxorubicin by each cell line (MDA-MB-468_H2B, MDA-MB-468_MDR1, and SUM-149PT) using a modification of the previously-published drug uptake assay (McKenna et al., 2017). The method leverages the intrinsic fluorescence of doxorubicin to quantify the movement of doxorubicin from the extracellular space into cells. Briefly, each cell line was introduced into 96-well microtiter plates at ~10,000 cells per well. Each well was imaged at 20–25 min intervals via fluorescent microscopy with a 20 × objective in 2 × 2 image montages on a BD Pathway 855 Bioimager (BD Biosciences, San Jose, CA). Imaging began 1 h prior to and continued for approximately 24 h following application of 1 μM of doxorubicin. After 8 h, doxorubicin was removed via media replacement. This timeframe allowed for an extended observation of drug uptake without inducing morphological changes and cell death that would limit the effect of the measurement. To measure the effect of TQR and NU7441 on drug uptake kinetics in the MDA-MB-468_MDR1 and the SUM-149PT cell lines, respectively, each sensitizer was applied over a range of concentrations (250–2 nM for TQR and 2 μM−16 nM for NU7441 both via a 2-fold dilution series) 1 h prior to doxorubicin application. At least three replicates of each treatment condition were collected.

The collected images were subsequently post-processed to correct for uneven background illumination and to isolate the contribution of each fluorophore in the experiment. First, the illumination function for each image was estimated (Jones et al., 2006). The image is defined:

$\begin{array}{l}I=L(C+b)\end{array}$

where I is the image, L is the illumination function, C is signal from cells, and b is the background. The signal from cells was removed from each image through use of a median disc filter with a radius of 50, isolating b. To estimate L, the background-only images in each well were averaged over all timepoints. A smooth surface was fit to this averaged image, and the surface was normalized to a maximum value of 1. Each image in the time series was divided by this surface (L) to correct for uneven illumination. Following illumination correction, a threshold-based approach was used to segment each cell.

To account for the various fluorophores in the experiment (H2BRFP, MDR1GFP, and doxorubicin), a linear unmixing approach was employed to isolate the signal from each fluorophore to more precisely quantify doxorubicin accumulation (Zimmermann, 2005). The approach leverages spectral imaging data collected at multiple excitation and emission wavelengths to isolate the signal from each fluorophore. This method can also be used for background subtraction by modeling the background (here, the signal from cell culture media) as an additional fluorophore. For these experiments, we define four fluorophores of interest: MDR1GFP, doxorubicin, H2BRFP, and background. The observed images are modeled as a linear combination of the signals from each of these fluorophores:

$\begin{array}{l}\left[\begin{array}{cccc}{S}_{H2B}& S_{MDR}& S_{Dox}& S_{background}\end{array}\right]T_{4\times n}=\left[\begin{array}{cccc}{I}_1& I_2& \dots & I_n\end{array}\right]\end{array}$

where S_H2B is the signal from the H2BRFP, S_MDR is the signal from the GFP-tagged MDR1, S_Dox is the signal from doxorubicin, and S_background is the background signal from cell culture media. T is the transformation matrix that estimates the contribution from each fluorophore in creating each image I. In this work, five images (n = 5) were collected at each timepoint. The excitation, dichroic, and emission filters for each image are listed in Table 2.

TABLE 2

Table 2. Filter settings.

To construct T, images of each fluorophore were collected from control samples. Specifically, control images of GFP, H2BRFP-positive cells, doxorubicin, and background were collected. For each fluorophore, the image with the highest intensity is assumed to be the true image; i.e., the corresponding entry in T is set to 1. The relative intensity of the other four images with respect to the true image are then estimated. This normalized spectrum is deposited into the row of T corresponding to the current fluorophore. T is estimated at each timepoint to compensate for any temporal changes in fluorophore intensity.

With an estimate of T and a spectral image set for each well at each timepoint, the underlying signals (i.e., S_H2B, S_MDR, S_Dox, S_background) can be estimated using QR decomposition [implemented in MATLAB (Mathworks, Natick, MA)]. This can be done on a per-pixel basis as shown in Supplementary Figure 1. However, as we are only interested in the intracellular and extracellular doxorubicin signals, the average value from each image in the intracellular and extracellular (I_{i, I}, I_{i, E}) space was calculated using a cell segmentation (as detailed above). Each signal can then be recovered:

$\begin{array}{l}\left[\begin{array}{cccc}{S}_{H2B,I}& S_{MDR,I}& S_{Dox,I}& S_{background,I}\\ S_{H2B,E}& S_{MDR,E}& S_{Dox,E}& S_{background,E}\end{array}\right]T_{4\times n}=\left[\begin{array}{ccc}{I}_{1,I}& \dots & I_{5,I}\\ I_{1,E}& \cdots & I_{5,E}\end{array}\right]\end{array}$

where S_{Dox, I} and S_{Dox, E} are the signals from doxorubicin in the intracellular and extracellular spaces, respectively. Similar definitions apply for the other signals S.

Finally, S_Dox is converted into doxorubicin concentration. We assume that doxorubicin signal is linearly proportional to its concentration, [Dox] (McKenna et al., 2017):

$\begin{array}{l}{S}_{Dox}=a\left[Dox\right]+b\end{array}$

To calibrate this model, images are collected on a series of wells containing a range of known doxorubicin concentrations. Estimates of a and b were obtained by fitting the doxorubicin signal equation to these control data. The image processing pipeline is illustrated in Supplementary Figure 1.

Doxorubicin Treatment Response Imaging

Using the previously-published dose-response assay, each cell line was treated with a range of doxorubicin concentrations (5,000–10 nM via a 2-fold dilution series) for 24 h as monotherapy. Additionally, the sensitizing effects of TQR and NU7441 in the MDA-MB-468_MDR1 and the SUM-149PT cell lines, respectively, were investigated by applying those therapies over a range of concentrations 1 h prior to application of doxorubicin. TQR concentrations in a 2-fold dilution series from 250 to 2 nM were used for the MDA-MB-468_MDR1 cell line, and NU7441 concentrations in a 2-fold dilution series from 2 μM to 15 nM were used for the SUM-149PT cell line. These combination studies were each performed at three doxorubicin concentrations. All drug (doxorubicin and sensitizer) was removed from each well via media replacement at 24 h. These cells were imaged daily via fluorescent microscopy for at least 15 days following treatment. For these studies, fluorescence microscopy images were collected using a Synentec Cellavista High End platform (SynenTec Bio Services, Münster, Germany) with a 20 × objective and tiling of 25 images. To generate images, the H2BmRFP fluorophores were excited with 529 nm light for 650 ms, and emissions were collected at 585 nm. Nuclei were segmented and counted in ImageJ (http://imagej.nih.gov/ij/) using a previously-described, threshold-based method (Frick et al., 2015) to quantify cell population. Six replicates of each treatment condition were collected. Media was refreshed every 3 days for the duration of each experiment to ensure sufficient growth conditions for surviving cells. Data were manually truncated when cell populations reached carrying capacity. At this point, signals from neighboring nuclei overlap, and the cell counting algorithm becomes unreliable.

Model Fits

The three-compartment model described in Equations (1–3) was fit to the uptake data under each treatment condition (doxorubicin monotherapy and doxorubicin combination with sensitizer) for each cell line using a non-linear least squares optimization implemented in MATLAB. Of note, each cell line is assumed to have a single set of compartment model parameters (k_EF, k_FE, and k_FB) for each sensitizer concentration; i.e., a parameter set for doxorubicin monotherapy and a set for each sensitizer concentration. The mean errors of the best-fit model across all timepoints and treatment conditions with respective standard deviations are reported. Similarly, the pharmacodynamic model described by Equations (4–6) was fit to the dose response data from all treatment conditions (i.e., doxorubicin monotherapy and doxorubicin combination with sensitizer) for each cell line. Each treatment condition in each cell line was fit independently, yielding cell line- and treatment condition-specific parameter values. This was also accomplished through a non-linear least squares optimization implemented in MATLAB, and we report the mean percent errors of the best-fit models across all timepoints. For additional details on the model fitting procedure see McKenna et al. (McKenna et al., 2017).

Measurement of Pharmacologic Properties With Equivalent Dose

We assume, by definition, that each unique treatment response timecourse corresponds to a specific equivalent dose. As the equivalent dose is perfectly known for doxorubicin monotherapy [i.e., the equivalent dose is simply C_B, which can be directly calculated with k_FE, k_EF, and k_FB values measured from drug uptake studies], the equivalent dose for co-treatment conditions can be estimated by comparing treatment response dynamics from co-treatment conditions to those from doxorubicin monotherapy treatments. With appropriate experimental design to isolate each equivalent dose parameter (i.e., k_FE, k_EF, and k_FB), this approach can quantify the effect of each sensitizing agent on their PK/PD pathway. Specifically, by assuming the effect of each sensitizing therapy is limited to a single equivalent dose parameter, the effect of TQR on k_EF and the effect of NU7441 on k_FB can be measured. As response under all treatment conditions (i.e., doxorubicin monotherapy and co-treatment with a sensitizer) can by summarized by the parameters in Equations (4-6) (i.e., p = [k_{d, a}, k_{d, b}, r]), we use model parameters to compare treatment response timecourses.

The response parameters (p) from doxorubicin monotherapy experiments are first interpolated with respect to equivalent dose via a local linear approach. This yields a continuous set of parameters (p_est) across all possible equivalent doses in the range from no treatment to maximal doxorubicin dose. The fit parameter values (p_fit) for each of the m co-treatment conditions are then matched to the interpolated parameters from doxorubicin-only treatment conditions (p_est) to estimate the equivalent dose (D_est) for each co-treatment condition. Specifically, D_est is the set of equivalent doses that correspond to the best matches between p_fit and p_est in the L₂ norm sense (i.e., min ||p_est − p_fit||₂). This process is illustrated in Supplementary Figure 2. The following constrained objective function, G (k_x), can then be used to estimate k_x (the equivalent dose parameter under investigation; e.g., k_EF and k_FB) for each of the n sensitizer concentrations:

$\begin{array}{l}G\left(k_x\right)={\displaystyle \underset{k_x}{min}}{\displaystyle \sum _{i=1}^m}{\left(D_{est,i}-D_i\left(k_x\right)\right)}^2\\ suchthat{k}_{x,q+1}-k_{x,q}\ge 0\forall q=\left[1,\dots ,n\right]& (7)\end{array}$

where the D_i is the equivalent dose calculated for each co-treatment condition as described below, and D_{est, i} is the estimated equivalent dose for the i^th co-treatment condition. Specifically, in calculating D_ifor the NU7441 experiments, we fix k_EF and k_FE values and optimize k_FB values corresponding to each sensitizer concentration in the co-treatment conditions. The constraints in the objective function ensure that k_FB increases monotonically with sensitizer concentration. Similarly, for the TQR experiments, we fix k_EF and k_FB and optimize k_FE for each sensitizer concentration. This objective function was minimized via a constrained optimization routine implemented in MATLAB. The non-parametric interpolation and optimization procedures were utilized as we did not assume any functional relationships between model parameters and equivalent dose. While this fitting procedure could have been made more robust by proposing such functional relationships, we implemented this non-parametric approach to allow for greater generalizability.

Comparison of Cell Lines With Equivalent Dose

As the MDA-MB-468_MDR1 line was engineered from the MDA-MB-468_H2B line, we hypothesize that the response of these cell lines to doxorubicin therapy is not significantly different when compared via equivalent dose. Specifically, the mechanism of action of MDR1 is to increase drug efflux, which effectively reduces the equivalent dose in the MDA-MB-468_MDR1 line for a given treatment timecourse. Indeed, the proposed equivalent dose metric was developed to account for the differing pharmacokinetic properties between these cell lines to more precisely compare their respective responses to therapy. To test this hypothesis, survival of the parental MDA-MB-468_H2B cell line is compared to that of the MDA-MB-468_MDR1 cell line. This comparison is made utilizing a conventional treatment response assay in which survival is assessed 72 h following treatment. Specifically, each cell line was treated with a range of doxorubicin concentrations (5,000–10 nM via a 2-fold dilution series) for 24 h as monotherapy, and survival was assessed via cell counting. Survival data for each cell line was fit with a pair of Hill functions. The first of these Hill functions assumed the dose to be the applied doxorubicin concentration. The second utilized the equivalent dose (D_eq) calculated with cell-line specific k_EF, k_FE, and k_FB values. We report the EC₅₀ (drug concentration at half-maximal effect) for each cell line as measured via extracellular doxorubicin concentration and equivalent dose.

Results

Treatment Response in MDA-MB-468_MDR1 Cell Line

The measured intracellular doxorubicin concentration timecourses for the MDA-MB-468_MDR1 cell line under doxorubicin monotherapy and combination therapy with TQR are shown in Figure 2A. Intracellular doxorubicin increases with TQR concentration. The average intracellular concentration at the end of each experiment, estimated with the last 10 timepoints, is significantly different among the treatment groups (one-way ANOVA, p < 1e-5). Equations (1–3) are fit to these data, and the best-fit models are overlaid on the timecourses. The mean error of the best-fit pharmacokinetic models was 45.6 (±47.4) nM across all timepoints and treatment conditions, and the corresponding model parameters are shown in Figures 2B–D. Increasing TQR concentrations decrease doxorubicin efflux in the MDA-MB-468_MDR1 cell line in a dose-dependent manner. For example, the efflux rate (k_FE) is decreased from 0.216 (±0.028) h⁻¹ to 0.046 (±0.008) h⁻¹ as TQR increases from 2 to 250 nM (the bounds here and below correspond to the 95% confidence interval of the parameter estimates). k_EF values varied with TQR concentration, all falling within [1.63, 3.46] × 10⁻⁶ h⁻¹.

FIGURE 2

Figure 2. Doxorubicin and TQR combination studies in the MDA-MB-468_MDR1 cell line. Timecourses of the mean intracellular concentration of doxorubicin with corresponding standard deviations are shown for each treatment condition in (A). Doxorubicin accumulation increases along with TQR concentrations. Equations (1–3) were fit to the data, and the best-fit models are overlaid on the data (smooth lines) in a. Model parameter fits corresponding to the best-fit models are shown in (B–D). Similar k_EF and k_FB vales are observed across all TQR concentrations. There is a trend of decreasing k_FE values with increasing TQR concentrations (C), consistent with MDR1 inhibition by TQR. Cell counts of MDA-MB-468_MDR1 following combination treatment with TQR and doxorubicin are show in panels (E-G). In each plot, a fixed concentration of doxorubicin is applied with variable TQR concentrations. These counts are fit with Equations (4–6) as described in section Model Fits, and the best-fit model is overlaid on the cell counts [smooth lines in panels (E–G)]. Error bars represent the 95% CI from six experimental replicates for each treatment condition. Model parameters with corresponding 95% CI are shown in (H–J) as a function of TQR concentration. For each doxorubicin concentration, the death rate (k_d,a and k_d,b) increased with TQR concentration (H,I).

Treatment response timecourses for the MDA-MB-468_MDR1 cell line under doxorubicin combination therapy with TQR are shown in Figures 2E–G. Equations (4–6) are fit to these data, and the best-fit models are overlaid on the observed cell counts. Model parameters are shown in Figures 2H–J. For a fixed concentration of doxorubicin, increasing concentrations of TQR incrementally sensitize cells to doxorubicin. For example, at a fixed dose of 156 nM doxorubicin, increasing the TQR concentration from 0 to 250 nM increased the death rate (k_{d, a}) from −0.16 (±0.23) × 10⁻² h⁻¹ to 2.21 (±0.1) × 10⁻² h⁻¹. TQR monotherapy did not affect the growth of these cells as shown in Supplementary Figure 3. Treatment response timecourses of the MDA-MB-468_MDR1 line to doxorubicin monotherapy are shown in Figure 3A. These data are fit with Equations (4–6), and the best-fit models are overlaid on the observed cell counts. The mean percent error of the best-fit model across all timepoints and treatment conditions is 10.3%. Prior to treatment, the MDA-MB-468_MDR1 line demonstrated a proliferation rate (k_p) of 2.12 (±0.03) × 10⁻² h⁻¹. Treatment response varied smoothly with doxorubicin concentration, and this response is quantified by the parameters in Figures 3B–D. Notably, high variance in parameter estimates is observed as values of r approach 0.05 h⁻¹ and values of k_{d, b} approach 0 h⁻¹. There exists intrinsic uncertainty at this limit as the rapid dynamics (r) coupled with small k_{d, b} effects cannot be resolved by the current data. This uncertainty in r for small k_{d, b} does not affect model predictions as demonstrated by a sensitivity analysis in previous work (McKenna et al., 2017).

FIGURE 3

Figure 3. Treatment response in MDA-MB-468_MDR1 (left column) and SUM-149PT (right column) cell lines under doxorubicin monotherapy. The top row [panels (A,E)] shows cell counts over time from treatment response studies for each cell line. For these studies, cells were treated with a fixed concentration of doxorubicin for 24 h. These counts are fit to Equations (4–6) as described in section Model Fits, and the best-fit model is overlaid on the cell counts [smooth lines in (A,E)]. Error bars represent the 95% CI from six experimental replicates for each treatment condition. Model parameters with corresponding 95% CI are shown in the bottom three rows as a function of doxorubicin concentration. Panels (B–D) show fits from the MDA-MB-468_MDR1 experiments, and panels (F–H) show fits from the SUM-149PT experiments. For each doxorubicin concentration for each cell line, the death rate (k_d,a and k_d,b) increased with increasing doxorubicin concentrations.

By leveraging the proposed mechanistic model and equivalent dose statistic, k_FE values for each TQR concentration can be estimated using the measured treatment response data and the optimization routine outlined in section Measurement of Pharmacologic Properties With Equivalent Dose. To make these measurements, the equivalent dose for each doxorubicin monotherapy condition was first calculated with the PK model parameters measured in the doxorubicin uptake studies. Specifically, k_FE, k_EF, and k_FB were measured to be 0.313 h⁻¹, 3.08 × 10⁻⁶ h⁻¹ and 0.0212 h⁻¹, respectively. The equivalent dose statistic was then estimated for each co-treatment condition. To perform this estimation, treatment response parameters from co-treatment conditions were matched to those from doxorubicin monotherapy conditions (Figures 4A–C). As the equivalent doses for all monotherapy conditions are perfectly known (i.e., C_B = D_eq for doxorubicin monotherapy), the equivalent dose for each co-treatment condition can be estimated with the matching process illustrated in Supplementary Figure 2. Briefly, parameter values are estimated across a range of equivalent doses utilizing parameters from the doxorubicin monotherapy experiments. The equivalent dose for each treatment condition can then be estimated by matching measured parameter values to those estimates. To demonstrate the efficacy of the parameter matching in comparing treatment response timecourses, a subset of responses from doxorubicin monotherapy and co-treatment conditions are color-coded to their estimated equivalent dose (Figures 4D–F). Note similar dynamics for similarly-colored data, indicating the efficacy of the parameter matching in comparing treatment response timecourses. With estimates of equivalent dose for all co-treatment conditions, the k_FE value for each TQR concentration was estimated with the optimization routine summarized by Equation (7). As we hypothesized that the effect of TQR is limited to k_FE (Figure 4G), k_EF and k_FB values were fixed to the values reported above in the optimization routine. The optimized k_FE values for all TQR concentrations are shown in Figure 4H. Decreasing k_FE values were observed with increasing TQR concentrations, matching the measurements from the uptake studies in Figure 2.

FIGURE 4

Figure 4. Leveraging equivalent dose to estimate the effect of TQR in the MDA-MB-468_MDR1 cell line. The equivalent dose for each doxorubicin monotherapy condition was first calculated with the PK model parameters measured in the doxorubicin uptake studies. The equivalent dose statistic was then estimated for each co-treatment condition by matching treatment response parameters from co-treatment conditions to those from doxorubicin monotherapy conditions. Parameter values from all doxorubicin monotherapy and co-treatment conditions are plotted as a function of equivalent dose (A-C). A subset of responses from doxorubicin monotherapy and co-treatment conditions are color-coded to their estimated equivalent dose (D–F). Similar dynamics are observed with similarly-colored data, demonstrating the efficacy of the parameter matching in comparing treatment response timecourses. As TQR impairs the function of the MDR1 pump, we hypothesized the effect of TQR is limited to the k_FE parameter (G). With estimates of equivalent dose for all treatment conditions, the k_FE value for each TQR concentration was estimated with the optimization routine summarized by Equation (7) (H). These values, calculated with treatment response data, agree well with direct measurements of k_FE reported in Figure 2. We note the large confidence intervals are a result of the optimization approach, in which the value (1/k_FE) was optimized.

The equivalent dose can summarize all treatment conditions in the MDA-MB-468_MDR1 cell line and is predictive of response. Further, this statistic can be leveraged to quantify the effect of TQR.

Treatment Response in SUM-149PT Cell Line

The measured intracellular doxorubicin concentration timecourses for the SUM-149PT cell line under doxorubicin monotherapy and combination therapy with NU7441 are shown in Figure 5A. NU7441 treatment did not affect intracellular doxorubicin accumulation following treatment. The average intracellular doxorubicin concentration at the end of each experiment, estimated with the last 10 timepoints, did not demonstrate significant differences at the p = 0.05 level (one-way ANOVA). Equations (1–3) are fit to the uptake data, and the best-fit model is overlaid on the timecourses. The corresponding model parameters are shown in Figures 5B–D. The mean error of the best-fit pharmacokinetic model was 77.9 (±71.4) nM across all treatment conditions and timepoints. Further, similar values of k_FE, k_EF, and k_FB are observed across all NU7441 concentrations (Figures 5B–D). Given its effect on DNA-PK, NU7441 is not expected to affect intracellular doxorubicin accumulation.

FIGURE 5

Figure 5. Doxorubicin and NU7441 combination studies in the SUM-149PT cell line. Timecourses of the mean intracellular concentration of doxorubicin with corresponding standard deviations are shown for each treatment condition in a. No significant difference in doxorubicin accumulation was observed as a function of NU7441 concentration. Equations (1–3) were fit to the data, and the best-fit models are overlaid on the data (smooth lines) in (A). Model parameter fits corresponding to the best-fit models are shown in (B–D). For each model parameter, similar vales were observed across all NU7441 concentrations, consistent with the similar intracellular doxorubicin timecourses in (A). Counts of SUM-149PT cells following combination treatment with NU7441 and doxorubicin are show in panels (E–G). In each plot, a fixed concentration of doxorubicin is applied with variable NU7441 concentrations. These counts are fit with Equations (4–6) as described in section Model Fits, and the best-fit models are overlaid on the cell counts [smooth lines in panels (E–G)]. Error bars represent the 95% CI from six experimental replicates for each treatment condition. Model parameters with corresponding 95% CI are shown in panels (H–J) as a function of NU7441 concentration. For each doxorubicin concentration, the death rate (k_d,a) increased with NU7441 concentration (H). The parameters shown in panels (I,J) are unable to be resolved with the current data as discussed in section Model Fits.

Treatment response timecourses for the SUM-149PT cell line under doxorubicin co-treatment with NU7441 are shown in Figures 5E–G. Equations (4–6) are fit to these data, and the best-fit models are overlaid on the observed cell counts. Model parameters are shown in Figures 5H–J. For a fixed concentration of doxorubicin, increasing concentrations of NU7441 incrementally sensitized cells to doxorubicin. For example, with a fixed dose of 156 nM doxorubicin, NU7441 concentrations increased the death rate (k_{d, a}) from 0.25 (±0.16) × 10⁻² h⁻¹ to 2.00 (±0.06) × 10⁻² h⁻¹. NU7441 monotherapy did not affect the growth of these cells as shown in Supplementary Figure 3. Treatment response timecourses of the SUM-149PT line to doxorubicin monotherapy are shown in Figure 3E. These data are fit with Equations (4–6), and the best-fit models are overlaid on the observed cell counts. The mean percent error of the best-fit model across all treatment conditions is 11.9%. Prior to treatment, the SUM-149PT line demonstrated a proliferation rate (k_p) of 2.58 (±0.03) × 10⁻² h⁻¹. Treatment response varied smoothly with doxorubicin concentration, and this response is quantified by the parameters in Figures 3F–H.

By leveraging the proposed mechanistic model and equivalent dose statistic, k_FB values for each NU7441 concentration can be estimated using the measured treatment response data and the optimization routine summarized by Equation (7). To make these measurements, the equivalent dose for each doxorubicin monotherapy condition was first calculated with the PK model parameters measured in the doxorubicin uptake studies. Specifically, k_EF, k_FE, and k_FB were measured to be 4.00 × 10⁻⁶ h⁻¹ and 0.165 h⁻¹, and 0.236 h⁻¹, respectively. These were calculated by fitting the SUM-149PT uptake studies assuming constant parameters for all NU7441 concentrations. The equivalent dose statistic was then estimated for each co-treatment condition. To perform this estimation, treatment response parameters from co-treatment conditions were matched to those from doxorubicin monotherapy conditions (Figures 6A–C). As the equivalent doses for all monotherapy conditions are perfectly known, the equivalent dose for each co-treatment condition can be estimated with this matching process. To demonstrate the efficacy of the parameter matching in comparing treatment response timecourses, a subset of responses from doxorubicin monotherapy and co-treatment conditions are color-coded to their estimated equivalent dose (Figures 6D–F). Note similar dynamics for similarly-colored data. With estimates of equivalent dose for all co-treatment conditions, the k_FB value for each NU7441 concentration was estimated with the optimization routine summarized by Equation (7). As we hypothesized that the effect of NU7441 is limited to k_FB(Figure 6G), k_FE and k_EF values were fixed to the values reported above in the optimization routine. The optimized k_FB values for all NU7441 concentrations are shown in Figure 6H. Increasing k_FB values were observed with increasing NU7441 concentrations, indicating the functional increase in drug with NU7441, mediated through its effect on DNA-PK. We note that the DNA repair pathway affected by NU7441 is not directly measured in the uptake studies. Recall from section Equivalent Dose that k_FB is a mixed measure of doxorubicin binding and DNA repair and describes the functional net binding rate. Thus, these values cannot be directly compared to the values extracted from the uptake study.

FIGURE 6

Figure 6. Leveraging equivalent dose to estimate the effect of NU7441 in the SUM-149PT cell line. The equivalent dose for each doxorubicin monotherapy condition was first calculated with the PK model parameters measured in the doxorubicin uptake studies. The equivalent dose statistic was then estimated for each co-treatment condition by matching treatment response parameters from co-treatment conditions to those from doxorubicin monotherapy conditions. Parameter values from all doxorubicin monotherapy and co-treatment conditions are plotted as a function of equivalent dose (A-C). A subset of responses from doxorubicin monotherapy and co-treatment conditions are color-coded to their estimated equivalent dose (D–F). Similar dynamics are observed with for similarly-colored data, demonstrating the efficacy of the parameter matching in comparing treatment response. As NU7441 impairs the function DNA-PK, we hypothesized the effect of NU7441 is limited to the k_FB parameter (G). With estimates of equivalent dose for all treatment conditions, the k_FB value for each NU7441 concentration was estimated with the optimization routine summarized by Equation (7). Increasing values of k_FB are observed with increasing NU7441 concentrations, indicating an increase in functional drug bound (H). These values cannot be directly observed with the uptake study, demonstrating the utility of the equivalent dose in estimating parameters that cannot be directly measured with current techniques.

The equivalent dose can summarize all treatment conditions in the SUM-149PT cell line and is predictive of response. Further, this statistic can be leveraged to quantify the specific effect of NU7441 with observed treatment response data.

Comparison of MDA-MB-468_MDR1 and MDA-MB-468_H2B

The measured intracellular doxorubicin concentration timecourses with accompanying best-fit models for the MDA-MB-468_H2B and MDA-MB-468_MDR1 cell lines are shown in Figure 7. Decreased doxorubicin accumulation was observed in the MDA-MB-468_MDR1 cell line relative to its parental line, MDA-MB-468_H2B. Notably, drug efflux was significantly elevated in the MDA-MB-468_MDR1 line relative to its parental line with k_FE values of 1.01 (±0.08) × 10⁻¹ h⁻¹ and 0.52 (±0.04) × 10⁻¹ h⁻¹, respectively (p < 0.05). The mean errors of the best-fit pharmacokinetic models across all timepoints were 44.7 and 58.7 nM for the MDA-MB-468_H2B and the MDA-MD-468_MDR1 lines, respectively.

FIGURE 7

Figure 7. Comparison of MDA-MB-468_H2B and MDA-MB-468_MDR1 cell lines using equivalent dose. The intracellular doxorubicin concentration with 95% CI for each cell line is shown in (A). The MDA-MB-468_H2B line demonstrates increased intracellular accumulation of doxorubicin relative to the MDA-MB-468_MDR1 line. Equations (1–3) are fit to the doxorubicin uptake data, and the best-fit models are overlaid on the data in a (smooth line). The corresponding parameters with 95% CI are shown in (B–D). The MDA-MB-468_H2B data are shown in red, and the MDA-MB-468_MDR1 data are shown in blue. Notably, the efflux of drug from the MDA-MB-468_MDR1 (k_FE) line is significantly greater than the corresponding rate in the MDA-MB-468_H2B line (p < 0.05). Treatment response is traditionally summarized by cell survival and plotted against applied drug concentration. The cell count relative to control for each cell line is shown as a function of extracellular doxorubicin concentration and equivalent dose in (E,F), respectively. While a significant difference is observed when comparing these cell lines via EC₅₀calculated with the extracellular doxorubicin concentration, no significant difference is observed when comparing the EC₅₀ statistic derived from the equivalent dose. The equivalent dose can account for the differing pharmacokinetic properties to reveal similar doxorubicin pharmacodynamics in these cell lines.

The survival of each cell line 72 h following treatment is compared as a function of extracellular doxorubicin concentration and equivalent dose in Figures 7E,F. The EC₅₀ as measured with the extracellular doxorubicin concentration for the MDA-MB-468_H2B and the MDA-MB-468_MDR1 are 101.6 (±28.9) and 350.6 (±109) nM, respectively. These measures indicate that there is a statistically significant difference between these cell lines (p < 0.05, t-test). The EC₅₀ as measured with the equivalent dose for the MDA-MB-468_H2B and the MDA-MB-468_MDR1 lines are 53.3 (±15.1) and 93.7 (±29.2) nM, respectively. These values are not different at p = 0.05 (t-test), indicating the similarity of these lines. Indeed, the only difference between cell lines is the overexpression of the MDR1 efflux pump. The intrinsic sensitivity of these cell lines to treatment should remain similar, and the equivalent dose reflects this similarity. The response of the MDA-MB-468_H2B and MDA-MB-468_MDR1 cell lines are not significantly different as measured by D_eq.

Discussion

We have proposed and demonstrated the utility of a mathematical modeling framework to quantify pharmacologic properties. We further proposed a new metric, the equivalent dose (D_eq), which provides a biochemically-based measure of treatment effect. With the data presented here, we show that a mechanistic mathematical model of treatment response can succinctly summarize a range of treatments to allow for more precise comparison of treatment response among cell lines. Further, we have shown how this model provides quantitative biological insight into the biochemical drivers of treatment response. We demonstrate that a mathematical modeling framework allows for quantification of pharmacologic processes through population-scale measurements.

Treatment response is driven by cell-line specific pharmacologic properties. Conventional summary statistics of treatment response data often conflate these pharmacologic properties, limiting their utility. To more effectively advance the study of treatment response, methods that explicitly consider this variability are needed to more precisely quantify biological drivers of treatment response. While previous treatment response assays provide insight in the relative sensitivity of a cell line to therapy (Fallahi-Sichani et al., 2013), the proposed approach quantifies specific drivers of treatment sensitivity. Through the approach proposed in this work, we demonstrate how intracellular pharmacologic properties can be quantified using limited data from population-level observations of treatment response.

This work is limited by its use of doxorubicin, which is intrinsically fluorescent, thereby allowing for the uptake model to be fit with experimental data. However, this approach need not be limited to fluorescent drugs. With appropriate experimental design, the approach summarized by Equation (7) can be leveraged to quantify any of the rates proposed in the model. Indeed, the optimized values of doxorubicin efflux in the MDA-MB-468_MDR1 line in Figure 4 are similar to those values measured by the uptake assay in Figure 2. Further, the effect of NU7441 in altering pharmacokinetics was quantified using only the treatment response data, as this effect cannot be directly measured in the uptake assay. Importantly, this work demonstrates that all treatment conditions collapse onto a single, smooth trajectory through parameter space as a function of equivalent dose, and this property can be leveraged to provide quantitative insight into the biological drivers of treatment response. While cell lines could not be compared without precise estimates of all model parameters, this approach can nevertheless be used to quantify therapeutic perturbations within a given cell line. It is straightforward to extend the proposed modeling approach as a means to more precisely quantify the effects of other parameters in the experimental microenvironment (e.g., how does pH or a specific nutrient concentration affect treatment response?). In this way, these variables can be mapped onto a unified treatment response framework to more efficiently advance precision medicine approaches. More generally, the approach outlined in this work demonstrates how mathematical modeling can be used as a “filter” to derive more specific measures from experimental data to advance systems biology.

Therapies that target PK/PD pathways offer the potential to sensitize cells to cytotoxic therapies, increasing the efficacy of therapy and allowing for lower doses of such therapeutics. The approach proposed in this work provides a means to quantify the respective contributions of PK/PD pathways, providing mechanistic insight into treatment response. This approach differs from current methods used to assess drug synergism and antagonism (Chou, 2006; Jones et al., 2014; Foucquier and Guedj, 2015; Chen and Lahav, 2016; Lederer et al., 2018). These methods have great utility in discovering and quantifying drug interactions; however, they cannot be leveraged to understand the mechanisms underlying the identified synergy/antagonism. While other methods have leveraged mechanistic data to identify synergy (Al-Lazikani et al., 2012; Gao et al., 2017; Yin et al., 2018), the proposed equivalent dose framework provides quantitative mechanistic insight into intracellular drug effects and allows for predictions of treatment response under a variety of treatment conditions. We posit that this mechanistic approach could facilitate clinical translation of combination therapies. Notably, therapeutic approaches intended to sensitize tumors to doxorubicin have demonstrated great preclinical activity; however, their efficacy has been limited in clinical trials. Specifically, negative results have been seen with TQR due to excess toxicities and inactivity (Pusztai et al., 2005; Fox and Bates, 2007). Similarly, DNA-PK inhibitors such as NU7441 have yet to demonstrate an effect clinically despite their preclinical promise (Zhao et al., 2006; Helleday et al., 2008; Davidson et al., 2013). We posit that the proposed modeling framework can be used to identify more effective strategies for dosing and assessing these therapeutics. In particular, the proposed modeling approach can provide precise guidance on the necessary dose adjustments to achieve a desired effect in the context of combination therapy. As we have demonstrated, a target equivalent dose can be achieved in a variety of ways. For example, the extracellular drug concentration timecourse can be tuned to reach a specified equivalent dose. Alternatively, the same equivalent dose can be achieved by altering cell line pharmacologic properties through sensitizers with concomitant changes in the extracellular doxorubicin timecourse. While realizing this goal in vivo will require a more complete model of treatment response (i.e., one that incorporates plasma pharmacokinetics and organ system toxicities), we have demonstrated the proposed model to be robust to various doxorubicin treatments and is general to sensitizing agents.

While the results of this study are promising, several limitations exist in the current approach. The first order pharmacokinetics model assumes static kinetic rates throughout the experiment, and the pharmacokinetic rates were investigated at only a single concentration. These rates are calculated as an average over all observed cells, not accounting for intercellular heterogeneity. Further, these kinetics may saturate as a function of doxorubicin concentration. This method remains to be validated in additional cell lines with other pharmacologic targets to address its generalizability. Additional properties of in vitro assays not explicitly considered in the current model have been shown to confound observed effects. For example, local cell densities have been found to affect treatment response (Greene et al., 2016). Finally, this model is deterministic and does not consider either population heterogeneity or cell evolution. Despite these limiting assumptions, we note the accuracy of the equivalent dose in summarizing population-level response to a range of doxorubicin treatment conditions.

In this work, we have demonstrated how mathematical modeling can be leveraged to quantify PK/PD pathways and more precisely compare treatment response among cell lines. It is the ultimate goal of precision cancer therapy to deliver the optimal therapy on the optimal schedule for the individual patient (McKenna et al., 2018). A necessary step toward this goal is to establish a robust functional relationship between applied treatment and subsequent response. The present study demonstrates the utility of the modeling framework and provides additional evidence that the response to therapy is predictable. In summary, analysis of treatment response data with mechanistic models can effectively quantify the effects of various biological and pharmaceutical perturbations on treatment response.

Author Contributions

MM and TY conceived the experiments. MM conducted the experiments under the guidance of TY. MM analyzed all results. JW aided MM in developing the numerical methods to fit the proposed models. All authors reviewed the manuscript.

Funding

We thank the National Institutes of Health for funding through: NCI R01 CA138599, NCI R01 CA186193, NCI U01 CA174706, NIGMS T32 GM007347, NCI F30 CA203220, NCI K25 CA204599, and NIBIB R21 EB022380. We thank the Cancer Prevention Research Institute of Texas (CPRIT) for RR160005; TY is a CPRIT Scholar in Cancer Research.

Conflict of Interest Statement

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Acknowledgments

We thank the funding agencies listed above for their support of this work. We would also like to thank the editors and reviewers for their critical evaluation of the manuscript.

Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fphys.2019.00616/full#supplementary-material

References

Al-Lazikani, B., Banerji, U., and Workman, P. (2012). Combinatorial drug therapy for cancer in the post-genomic era. Nat. Biotechnol. 2012:2284. doi: 10.1038/nbt.2284

CrossRef Full Text | Google Scholar

Anderson, A. R., and Quaranta, V. (2008). Integrative mathematical oncology. Nat. Rev. Cancer 8, 227–234. doi: 10.1038/nrc2329

PubMed Abstract | CrossRef Full Text | Google Scholar

Bouwman, P., and Jonkers, J. (2012). The effects of deregulated DNA damage signalling on cancer chemotherapy response and resistance. Nat. Rev. Cancer 12, 587–598. doi: 10.1038/nrc3342

PubMed Abstract | CrossRef Full Text | Google Scholar

Chang, B. D., Broude, E. V., Dokmanovic, M., Zhu, H., Ruth, A., Xuan, Y., et al. (1999). A senescence-like phenotype distinguishes tumor cells that undergo terminal proliferation arrest after exposure to anticancer agents. Cancer Res. 59, 3761–3767. doi: 10.1038/nrc2961

PubMed Abstract | CrossRef Full Text | Google Scholar

Chen, S. H., and Lahav, G. (2016). Two is better than one; toward a rational design of combinatorial therapy. Curr. Opin. Struct. Biol. 41, 145–150. doi: 10.1016/j.sbi.2016.07.020

PubMed Abstract | CrossRef Full Text | Google Scholar

Chou, T.-C. (2006). Theoretical basis, experimental design, and computerized simulation of synergism and antagonism in drug combination studies. Pharmacol. Rev. 58, 621–681. doi: 10.1124/pr.58.3.10

PubMed Abstract | CrossRef Full Text | Google Scholar

Ciszewski, W. M., Tavecchio, M., Dastych, J., and Curtin, N. J. (2014). DNA-PK inhibition by NU7441 sensitizes breast cancer cells to ionizing radiation and doxorubicin. Breast Cancer Res. Treat. 143, 47–55. doi: 10.1007/s10549-013-2785-6

PubMed Abstract | CrossRef Full Text | Google Scholar

Clarke, R., Leonessa, F., and Trock, B. (2005). Multidrug resistance/P-glycoprotein and breast cancer: review and meta-analysis. Semin. Oncol. 32, 9–15. doi: 10.1053/j.seminoncol.2005.09.009

PubMed Abstract | CrossRef Full Text | Google Scholar

Davidson, D., Amrein, L., Panasci, L., and Aloyz, R. (2013). Small molecules, inhibitors of DNA-PK, targeting DNA repair, and beyond. Front. Pharmacol. 4:5. doi: 10.3389/fphar.2013.00005

PubMed Abstract | CrossRef Full Text | Google Scholar

Elstrodt, F., Hollestelle, A., Nagel, J. H., Gorin, M., Wasielewski, M., van den Ouweland, A., et al. (2006). BRCA1 mutation analysis of 41 human breast cancer cell lines reveals three new deleterious mutants. Cancer Res. 66, 41–45. doi: 10.1158/0008-5472.CAN-05-2853

PubMed Abstract | CrossRef Full Text | Google Scholar

Eom, Y.-W., Kim, M. A., Park, S. S., Goo, M. J., Kwon, H. J., Sohn, S., et al. (2005). Two distinct modes of cell death induced by doxorubicin: apoptosis and cell death through mitotic catastrophe accompanied by senescence-like phenotype. Oncogene 24, 4765–4777. doi: 10.1038/sj.onc.1208627

PubMed Abstract | CrossRef Full Text | Google Scholar

Fallahi-Sichani, M., Honarnejad, S., Heiser, L. M., Gray, J. W., and Sorger, P. K. (2013). Metrics other than potency reveal systematic variation in responses to cancer drugs. Nat. Chem. Biol. 9, 708–714. doi: 10.1038/nchembio.1337

PubMed Abstract | CrossRef Full Text | Google Scholar

Farmer, H., McCabe, N., Lord, C. J., Tutt, A. N., Johnson, D. A., Richardson, T. B., et al. (2005). Targeting the DNA repair defect in BRCA mutant cells as a therapeutic strategy. Nature 434, 917–921. doi: 10.1038/nature03445

PubMed Abstract | CrossRef Full Text | Google Scholar

Fink, D., Aebi, S., and Howell, S. B. (1998). The role of DNA mismatch repair in drug resistance. Clin. Cancer Res. 4, 1–6.

PubMed Abstract | Google Scholar

Foucquier, J., and Guedj, M. (2015). Analysis of drug combinations: current methodological landscape. Pharmacol. Res. Perspect. 2015:149. doi: 10.1002/prp2.149

CrossRef Full Text | Google Scholar

Fowler, J. F. (1992). Brief summary of radiobiological principles in fractionated radiotherapy. Semin. Radiat. Oncol. 2, 16–21. doi: 10.1016/S1053-4296(05)80045-1

CrossRef Full Text | Google Scholar

Fox, E., and Bates, S. E. (2007). Tariquidar (XR9576): a P-glycoprotein drug efflux pump inhibitor. Expert Rev. Anticancer Ther. 7, 447–459. doi: 10.1586/14737140.7.4.447

PubMed Abstract | CrossRef Full Text | Google Scholar

Frick, P. L., Paudel, B. B., Tyson, D. R., and Quaranta, V. (2015). Quantifying heterogeneity and dynamics of clonal fitness in response to perturbation. J. Cell Physiol. 230, 1403–1412. doi: 10.1002/jcp.24888

PubMed Abstract | CrossRef Full Text | Google Scholar

Gao, H., Yin, Z., Cao, Z., and Zhang, L. (2017). Developing an agent-based drug model to investigate the synergistic effects of drug combinations. Molecules 22:2209. doi: 10.3390/molecules22122209

PubMed Abstract | CrossRef Full Text | Google Scholar

Gerlee, P. (2013). The model muddle: in search of tumor growth laws. Cancer Res. 73, 2407–2411. doi: 10.1158/0008-5472.CAN-12-4355

PubMed Abstract | CrossRef Full Text | Google Scholar

Gewirtz, D. A. (1999). A critical evaluation of the mechanisms of action proposed for the antitumor effects of the anthracycline antibiotics adriamycin and daunorubicin. Biochem. Pharmacol. 57, 727–741. doi: 10.1016/S0006-2952(98)00307-4

PubMed Abstract | CrossRef Full Text | Google Scholar

Greene, J. M., Levy, D., Herrada, S. P., Gottesman, M. M., and Lavi, O. (2016). Mathematical modeling reveals that changes to local cell density dynamically modulate baseline variations in cell growth and drug response. Cancer Res. 76, 2882–2890. doi: 10.1158/0008-5472.CAN-15-3232

PubMed Abstract | CrossRef Full Text | Google Scholar

Gudmundsdottir, K., and Ashworth, A. (2006). The roles of BRCA1 and BRCA2 and associated proteins in the maintenance of genomic stability. Oncogene 25, 5864–5874. doi: 10.1038/sj.onc.1209874

PubMed Abstract | CrossRef Full Text | Google Scholar

Hafner, M., Niepel, M., Chung, M., and Sorger, P. K. (2016). Growth rate inhibition metrics correct for confounders in measuring sensitivity to cancer drugs. Nat. Methods 13, 521–527. doi: 10.1038/nmeth.3853

PubMed Abstract | CrossRef Full Text | Google Scholar

Harris, L. A., Frick, P. L., Garbett, S. P., Hardeman, K. N., Paudel, B. B., Lopez, C. F., et al. (2016). An unbiased metric of antiproliferative drug effect in vitro. Nat. Methods 13, 497–500. doi: 10.1038/nmeth.3852

PubMed Abstract | CrossRef Full Text | Google Scholar

Helleday, T., Petermann, E., Lundin, C., Hodgson, B., and Sharma, R. A. (2008). DNA repair pathways as targets for cancer therapy. Nat. Rev. Cancer 8, 193–204. doi: 10.1038/nrc2342

PubMed Abstract | CrossRef Full Text | Google Scholar

Jarrett, A. M., Lima, E. A. B. F., Hormuth, D. A., McKenna, M. T., Feng, X., Ekrut, D. A., et al. (2018). Mathematical models of tumor cell proliferation: a review of the literature. Expert. Rev. Anticancer Ther. 18, 1271–1286. doi: 10.1080/14737140.2018.1527689

PubMed Abstract | CrossRef Full Text | Google Scholar

Jones, L. B., Secomb, T. W., Dewhirst, M. W., and El-Kareh, A. W. (2014). The additive damage model: a mathematical model for cellular responses to drug combinations. J. Theor. Biol. 2014:32. doi: 10.1016/j.jtbi.2014.04.032

CrossRef Full Text | Google Scholar

Jones, T. R., Carpenter, A. E., Sabatini, D. M., and Golland, P. (2006). “Methods for high-content, high-throughput image-based cell screening,” in Proceedings of the Workshop on Microscopic Image Analysis with Applications in Biology, 65–72. Available online at: https://pdfs.semanticscholar.org/d0a7/250de15140526b06ccd3ed8effceb77b04fe.pdf

Google Scholar

Lankelma, J., Fernández Luque, R., Dekker, H., and Pinedo, H. M. (2003). Simulation model of doxorubicin activity in islets of human breast cancer cells. Biochim. Biophys. Acta 1622, 169–178. doi: 10.1016/S0304-4165(03)00139-9

PubMed Abstract | CrossRef Full Text | Google Scholar

Lankelma, J., Fernández Luque, R., Dekker, H., van den Berg, J., and Kooi, B. (2013). A new mathematical pharmacodynamic model of clonogenic cancer cell death by doxorubicin. J. Pharmacokinet. Pharmacodyn. 40, 513–525. doi: 10.1007/s10928-013-9326-0

PubMed Abstract | CrossRef Full Text | Google Scholar

Larsen, A. K., Escargueil, A. E., and Skladanowski, A. (2000). Resistance mechanisms associated with altered intracellular distribution of anticancer agents. Pharmacol. Ther. 85, 217–229. doi: 10.1016/S0163-7258(99)00073-X

PubMed Abstract | CrossRef Full Text | Google Scholar

Larsen, A. K., and Skladanowski, A. (1998). Cellular resistance to topoisomerase-targeted drugs: from drug uptake to cell death. Biochim. Biophys. Acta 1400, 257–274. doi: 10.1016/S0167-4781(98)00140-7

PubMed Abstract | CrossRef Full Text | Google Scholar

Lederer, S., Dijkstra, T. M. H., and Heskes, T. (2018). Additive dose response models: explicit formulation and the loewe additivity consistency condition. Front. Pharmacol. 2018:31. doi: 10.3389/fphar.2018.00031

CrossRef Full Text | Google Scholar

Lobo, E. D., and Balthasar, J. P. (2002). Pharmacodynamic modeling of chemotherapeutic effects: application of a transit compartment model to characterize methotrexate effects in vitro. AAPS PharmSci. 4, 212–222. doi: 10.1208/ps040442

PubMed Abstract | CrossRef Full Text | Google Scholar

McKenna, M. T., Weis, J. A., Barnes, S. L., Tyson, D. R., Miga, M. I., Quaranta, V., et al. (2017). A predictive mathematical modeling approach for the study of doxorubicin treatment in triple negative breast cancer. Sci. Rep. 7:5725. doi: 10.1038/s41598-017-05902-z

PubMed Abstract | CrossRef Full Text | Google Scholar

McKenna, M. T., Weis, J. A., Brock, A., Quaranta, V., and Yankeelov, T. E. (2018). Precision medicine with imprecise therapy: computational modeling for chemotherapy in breast cancer. Transl. Oncol. 11, 732–742. doi: 10.1016/J.TRANON.2018.03.009

PubMed Abstract | CrossRef Full Text | Google Scholar

Mechetner, E., Kyshtoobayeva, A., Zonis, S., Kim, H., Stroup, R., Garcia, R., et al. (1998). Levels of multidrug resistance (MDR1) P-glycoprotein expression by human breast cancer correlate with in vitro resistance to taxol and doxorubicin. Clin. Cancer Res. 4, 389–398.

PubMed Abstract | Google Scholar

Mistry, P., Stewart, A. J., Dangerfield, W., Okiji, S., Liddle, C., Bootle, D., et al. (2001). In vitro and in vivo reversal of P-glycoprotein-mediated multidrug resistance by a novel potent modulator, XR9576. Cancer Res. 61, 749–758.

PubMed Abstract | Google Scholar

Prentice, R. L. (1976). A generalization of the probit and logit methods for dose response curves. Biometrics 32:761. doi: 10.2307/2529262

PubMed Abstract | CrossRef Full Text | Google Scholar

Pusztai, L., Wagner, P., Ibrahim, N., Rivera, E., Theriault, R., Booser, D., et al. (2005). Phase II study of tariquidar, a selective P-glycoprotein inhibitor, in patients with chemotherapy-resistant, advanced breast carcinoma. Cancer 104, 682–691. doi: 10.1002/cncr.21227

PubMed Abstract | CrossRef Full Text | Google Scholar

Quaranta, V., Tyson, D. R., Garbett, S. P., Weidow, B., Harris, M. P., and Georgescu, W. (2009). Trait variability of cancer cells quantified by high-content automated microscopy of single cells. Methods Enzymol. 467, 23–57. doi: 10.1016/S0076-6879(09)67002-6

PubMed Abstract | CrossRef Full Text | Google Scholar

Sanga, S., Sinek, J. P., Frieboes, H. B., Ferrari, M., Fruehauf, J. P., and Cristini, V. (2006). Mathematical modeling of cancer progression and response to chemotherapy. Expert Rev. Anticancer Ther. 6, 1361–1376. doi: 10.1586/14737140.6.10.1361

PubMed Abstract | CrossRef Full Text | Google Scholar

Shen, H., Schultz, M., Kruh, G. D., and Tew, K. D. (1998). Increased expression of DNA-dependent protein kinase confers resistance to adriamycin. Biochim. Biophys. Acta 1381, 131–138. doi: 10.1016/S0304-4165(98)00020-8

PubMed Abstract | CrossRef Full Text | Google Scholar

Simeoni, M., Magni, P., Cammia, C., De Nicolao, G., Croci, V., Pesenti, E., et al. (2004). Predictive pharmacokinetic-pharmacodynamic modeling of tumor growth kinetics in xenograft models after administration of anticancer agents. Cancer Res. 64, 1094–1101. doi: 10.1158/0008-5472.CAN-03-2524

PubMed Abstract | CrossRef Full Text | Google Scholar

Smith, G. C., and Jackson, S. P. (1999). The DNA-dependent protein kinase. Genes Dev. 13, 916–934.

Google Scholar

Tacar, O., Sriamornsak, P., and Dass, C. R. (2013). Doxorubicin: an update on anticancer molecular action, toxicity and novel drug delivery systems. J. Pharm. Pharmacol. 65, 157–170. doi: 10.1111/j.2042-7158.2012.01567.x

PubMed Abstract | CrossRef Full Text | Google Scholar

Tyson, D. R., Garbett, S. P., Frick, P. L., and Quaranta, V. (2012). Fractional proliferation: a method to deconvolve cell population dynamics from single-cell data. Nat. Methods 9, 923–928. doi: 10.1038/nmeth.2138

PubMed Abstract | CrossRef Full Text | Google Scholar

Wang, Z., Butner, J. D., Cristini, V., and Deisboeck, T. S. (2015). Integrated PK-PD and agent-based modeling in oncology. J. Pharmacokinet. Pharmacodyn. 42, 179–189. doi: 10.1007/s10928-015-9403-7

PubMed Abstract | CrossRef Full Text | Google Scholar

Yankeelov, T. E., Atuegwu, N., Hormuth, D., Weis, J. A., Barnes, S. L., Miga, M. I., et al. (2013). Clinically relevant modeling of tumor growth and treatment response. Sci. Transl. Med. 5:187ps9. doi: 10.1126/scitranslmed.3005686

PubMed Abstract | CrossRef Full Text | Google Scholar

Yankeelov, T. E., Quaranta, V., Evans, K. J., and Rericha, E. C. (2015). Toward a science of tumor forecasting for clinical oncology. Cancer Res. 75, 918–923. doi: 10.1158/0008-5472.CAN-14-2233

PubMed Abstract | CrossRef Full Text | Google Scholar

Yin, Z., Deng, Z., Zhao, W., and Cao, Z. (2018). Searching synergistic dose combinations for anticancer drugs. Front. Pharmacol. 9:535. doi: 10.3389/fphar.2018.00535

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhao, Y., Thomas, H. D., Batey, M. A., Cowell, I. G., Richardson, C. J., Griffin, R. J., et al. (2006). Preclinical evaluation of a potent novel DNA-dependent protein kinase inhibitor NU7441. Cancer Res. 66, 5354–5362. doi: 10.1158/0008-5472.CAN-05-4275

PubMed Abstract | CrossRef Full Text | Google Scholar

Zimmermann, T. (2005). “Spectral imaging and linear unmixing in light microscopy,” in Microscopy Techniques, ed J. Rietdorf (Berlin, Heidelberg: Springer), 245–265. doi: 10.1007/b102216

CrossRef Full Text | Google Scholar