A multi-state model of chemoresistance to characterize phenotypic dynamics in breast cancer

The development of resistance to chemotherapy is a major cause of treatment failure in breast cancer. While mathematical models describing the dynamics of resistant cancer cell subpopulations have been proposed, experimental validation has been difficult due to the complex nature of resistance that limits the ability of a single phenotypic marker to sufficiently identify the drug resistant subpopulations. We address this problem with a coupled experimental/modeling approach to reveal the composition of drug resistant subpopulations changing in time following drug exposure. We calibrate time-resolved drug sensitivity assays to three mathematical models to interrogate the models’ ability to capture drug response dynamics. The Akaike information criterion was employed to evaluate the three models, and it identified a multi-state model incorporating the role of population heterogeneity and cellular plasticity as the optimal model. To validate the model’s ability to identify subpopulation composition, we mixed different proportions of wild-type MCF-7 and MCF-7/ADR resistant cells and evaluated the corresponding model output. Our blinded two-state model was able to estimate the proportions of cell types with an R-squared value of 0.857. To the best of our knowledge, this is the first work to combine experimental time-resolved drug sensitivity data with a mathematical model of resistance development.

sensitive and resistant parameters for the simulated data sets in order to demonstrate the relative error in parameter identifiability of the LD50 values at the experimental noise observed. b. We display the ninety-five percent confidence intervals around 100 fitted sensitive and resistant slope parameters for the simulated data sets in order to demonstrate the relative error in parameter identifiability of the slope values at the experimental noise observed. We observe a higher relative error in the sensitive slope than in the resistant slope. c. We display the ninetyfive percent confidence intervals around 100 fitted fraction estimate parameters for the simulated data sets in order to demonstrate the relative error in parameter identifiability of the fractions at the experimental noise observed. We observe the error to be fairly consistent across all of the mixtures.

Figure S3: Parameter distributions around clustered input fraction parameters are
identifiable within 10% with current experimental noise. a. Model parameters of clustered low resistant fractions were input at intervals of 10% from 0-30% resistant. One hundred simulated data sets with experimental noise were fit to obtain 100 model parameter estimates for each fractional parameter. We compared each pairwise set of low resistant fraction parameter distributions with a multiple comparison t-test in Matlab and found significant differences with a p-value of 9.92×10 -8 . b. The same procedure was performed for high resistant fractions at intervals of 10% increase from 70-100% resistant. Again, the distributions of high resistant fraction parameter estimates were significantly different with a p-value of 9.92 × 10 -8 . \ Figure S4: Data from the model validation experiment is used to estimate sources of error.
The standard deviation among each set of replicate measurements is graphed against its average value, with each set of replicates represented as a single point.

Discussion of Intrinsic Stochastic Biological Variability and Sources of Error
While we often look at variation in measured data and mentally summarize it as "measurement error", the controls present in our model validation experiment allow us to decompose this variation into several sources of error. Error in the viability measurement performed using the Nexcelom Cellometer is approximately 0.5%, estimated using replicate measurements drawn from a single pool of cells. The variation in viability is larger, at an average of 2.2% across all samples, but is also non-uniformly distributed, as shown in Figure   S4.
Part of this variation can be explained as the result of differences in the proportion of cells of each subline plated into a given well; this variation has a standard deviation of approximately 1.7%. In Figure S4, however, the standard deviation ranges as high as 9%, indicating that an additional source of variation is present, and appears to be dependent on the population viability for that set of replicate measurements. We believe that this is best described as an actual variability in the response that replicate populations will display when given the same stimulus, due to downstream effects of stochastic variation in the early response. Among other things, we are aware that confluence has a large effect on cell survival -as a result, random fluctuations in early survival could snowball into these substantial differences over the duration of the drug perturbation. As each cell dies, it increases the probability that other nearby cells will die. (We speculate that this is mediated by the loss of pro-survival signals which the cells exchange.) Any fluctuation of survival in the early stages of the response is then propagated and amplified, magnifying the fluctuations into the pattern of variation seen here in supplemental Figure S4. This theory is consistent with the non-uniform distribution of variability; doses that are high or low enough to have more deterministic effects show minimal variation, more comparable to the error known to be present from variation in the initial seeding proportions and technical error in the assay while in doses where the probability of death is closer to 50%, a marginal change in the probability of survival is more likely to influence the outcome.