## Abstract

Cell-to-cell variations contribute to drug resistance with consequent therapy failure in cancer. Experimental techniques have been developed to monitor tumor heterogeneity, but estimates of cell-to-cell variation typically fail to account for the expected spatiotemporal variations during the cell growth process. To fully capture the extent of such dynamic variations, we developed a mechanistic mathematical model supported by *in vitro* experiments with an ovarian cancer cell line. We introduce the notion of dynamic baseline cell-to-cell variation, showing how the emerging spatiotemporal heterogeneity of one cell population can be attributed to differences in local cell density and cell cycle. Manipulation of the geometric arrangement and spatial density of cancer cells revealed that given a fixed global cell density, significant differences in growth, proliferation, and paclitaxel-induced apoptosis rates were observed based solely on cell movement and local conditions. We conclude that any statistical estimate of changes in the level of heterogeneity should be integrated with the dynamics and spatial effects of the baseline system. This approach incorporates experimental and theoretical methods to systematically analyze biologic phenomena and merits consideration as an underlying reference model for cell biology studies that investigate dynamic processes affecting cancer cell behavior. *Cancer Res; 76(10); 2882–90. ©2016 AACR*.

### Major Findings

Cells are never isolated from their surroundings. Even in an experimental setting, cells always grow as part of a population. As such, they are constantly exposed to random fluctuations in external signals. One such fluctuating signal comes from local cell density (LCD). LCD and the ability of cells to move can dramatically and dynamically change cell-to-cell variation. Tracking the intrinsic spatiotemporal process of cell growth reveals heterogeneity in crucial features, such as cell-cycle duration, growth, proliferation, and apoptosis rates. These constitute intrinsic dynamic baseline variations that are characteristics of all cell populations. These baseline variations are needed to explain and predict cell heterogeneity and to determine the significance of any alteration resulting from cell-to-cell variability, especially in cancer.

### Quick Guide to Equations and Assumptions

We developed a spatial mathematical model of cancer cell growth in a two-dimensional spatial environment. The model is formulated as an agent-based model that incorporates spatial aggregation/repulsion forces, as well as transitions through three cellular states: proliferation (P), apoptosis (A), and quiescence (Q). Cell movement is governed by a stochastic differential equation (SDE) with no constraints to a specific lattice, and transition rates between cellular compartments depend on a cell's microenvironment via a local cell density (LCD) measurement. A more thorough description of the model is provided in the Supplementary Materials. An important feature of this work is that a cell interacts only with other cells in its local neighborhood, defined here as a sphere of radius ε_{a} centered at the particle's position [*X*_{k}(*t*) for the *k*^{th} particle at time *t*]. We model cells as hard spheres with centers *X*_{k}(t) and fixed radius ε_{r}. Thus, the LCD can be defined as
where is the number of distinct cells in the local neighborhood. The k^{th} cell's motility is then governed by a SDE combining both attractive and repulsive forces:
Here, describes the population LCD at time t, and () models the attractive (repulsive) forces experienced by a cell due to the presence of neighboring cells. Both force terms are modeled deterministically as a gradient of symmetric functions, with distinct supports. The attractive force is chosen to account for the observed phenomenon of cellular aggregation and generally has a large support relative to , whereas the repulsive force models volume exclusion. Finally, random effects are included in the standard Wiener term , where denotes the mean free path. Randomness is included to account for the experimentally observed random walk-like motion by cells in low-density conditions. Note also that the movement is not constrained to a lattice and that this framework is adapted from Morale and colleagues (1–3).

Cells also experience transitions between three compartments: P, A, and Q. These transition rates are determined by an individual LCD , defined as in Eq. A. Cells in the Q compartment are neither dividing nor dying, but instead are in a resting phase. Often, this state is referred to as G_{0} and is the compartment where all dividing cells originate. From Q, cells can either transition into compartments P or A. Cells in P are undergoing proliferation, that is, progressing through the cell cycle. Thus, each cell in P has a state variable *a*(*t*), the cell-cycle stage/age, which measures a cell's progress through P. *a*(*t*) increases between discrete time steps *t* and *t* + Δ*t* via the rate ω(, c):
Here, c is the spatially constant drug concentration, which acts to increase (and thus increase the cell cycle duration). Note that the rate of cell-cycle progression is also dependent upon the LCD . We assume that due to LCD constraints, proliferating cells do not move through the cell cycle at a constant speed, but instead decelerate as increases. More precisely, we assume a base length of μ (=15) hours, after which, if *a*(*t*) > μ, division is successful and both mother and daughter cell reenter compartment Q. If the amount of time in culture spent in P exceeds a specified limit, the cell transitions to apoptosis, A. Cells undergoing apoptosis are destined to complete cell death; that is, once cells enter compartment A there are no transitions back to P or Q. Once completed, the cell is removed from the simulation. As in our previous report (4), the amount of time spent in A is dictated by an *a priori* γ distribution, independent of any state variables. Transitions can be either explicit or implicit. Explicit transition rates are interpreted as probabilities per unit time, that is, continuous Markov chain transition rates, while implicit rates depend on state variables specific to individual cells.

## Introduction

Within an individual tumor, there are typically genotypic and phenotypic variations. This heterogeneity, due to both genetic and nongenetic alterations, can be either temporary or irreversible (5–10). Tumor heterogeneity has been identified as one of the causes of cancer therapy failure, contributing to drug resistance (11, 12). Great efforts have been made to identify and categorize the different subpopulations of cells within a tumor/patient, and to determine their importance in relation to treatment, with the hope of finding ways to efficiently target them. It is accepted that such an approach primarily aims to find genetically stable clones, and assumes that each clone consists mainly of a homogenous population of cells, with insignificant variations concerning the subject of study. Thus, the common goal is to focus on (and target) the identified genetic alterations. However, this approach does not take into account the importance of temporal changes that are not necessarily the result of genetic alterations (13). Determining, additionally, the extent of each clone's plasticity would result in a more pragmatic treatment protocol.

It has been long recognized that a single clone of cells may have significant phenotypic variations, even concerning drug sensitivity (14, 15). Perhaps the most easily observed evidence of intrinsic nongenetic heterogeneity regarding drug response occurs in virtually every survival curve for cancer cells exposed to drugs, as killing curves have two key features: (i) a continuous curve, that is, a gradual slope, (ii) distinct residual cells that survive even after administration of high doses of the drug (Fig. 1A; Supplementary Table S1). Different types of cell-to-cell variations have been experimentally observed for a single population in many complex cellular processes, such as duration of apoptosis (8, 16), cell size and age (17), and duration of cell cycle (18). These variations occur in many organisms, generated by a variety of mechanisms that are based on stochastic and/or deterministic (primarily external) signals in a given cell population. In cancer studies, predictions of the disease dynamics are highly dependent on the way those expected baseline variations are evaluated, prior to any additional new alterations, both experimentally and theoretically. So far, the baseline variations have been reported in a limited way, as short-term observations. However, normalization with the expected dynamic baseline cell-to-cell variation of the spatiotemporal growth process has not been included, and thus, the current statistical approach that determines the reference model may produce false conclusions.

Cells are never isolated from their surroundings. Even in an experimental setting, cells always grow as part of a population. As such, they are constantly exposed to random fluctuations in external signals. One such fluctuating signal comes from local cell density (LCD). LCD and the ability of cells to move can dramatically and dynamically change cell-to-cell variation. We hypothesized that the intrinsic spatiotemporal process of cell growth should result in intrinsic considerable cell-to-cell variations. These constitute intrinsic dynamic baseline variations that are characteristics of all cell populations. We argue that the initial global cell density (GCD) is not sufficient to determine the baseline variations of the system. To quantify and mechanistically understand the impact of spatial LCD and the ability of cells to move, we developed a spatial mathematical model of cancer dynamics that can be used to theoretically investigate whether a specific set of local cellular principles can explain the observed distributions in cell-cycle duration (CCD), growth, and drug response. We show that this mechanism of spatial LCD can explain *in vitro* experimental cell growth data. We experimentally tested two representative geometric seeding configurations from a range of low to high LCD. Both seeding configurations shared the same GCD. We seeded OVCAR-8 cells in a simple 2D culture system and examined the spatiotemporal cell-to-cell differences in growth rate, proliferation, and paclitaxel-induced apoptosis percentages. Experimentally and theoretically, we demonstrate the extent of those cell-to-cell dynamic baseline variations, based solely on differences in geometrical seeding configurations, yet sharing the same global conditions. Understanding and quantifying these dynamic baseline variations can help to determine the significance of any alteration resulting from cell-to-cell variability, especially in cancer.

## Materials and Methods

### Growth curve

Cells (OVCAR-8 parental human ovarian carcinoma cells) were seeded with culture medium (RPMI1640 medium + 10% FBS + 100 U/mL penicillin streptomycin + 2 mmol/L glutamine) in 100 mm dishes in three different GCDs (about 7%, 20%, and 30%), with each density in two different geometrical configurations Figs 1E and 2A). The OVCAR-8 cell line was obtained from the NCI, NIH (Bethesda, MD). It is part of the NCI-60 panel of cancer cell lines, which were authenticated by the NCI using short tandem repeat markers to confirm cell identity. The cells were passaged for approximately two weeks in our laboratory after resuscitation and then used for experiments. Although the random seeding was mixed very well, the seeding pattern was 5 dots, spaced uniformly on the plate, containing 2/3 × 10^{5}, 2 × 10^{5}, and 10/3 × 10^{5} cells per dot, respectively, to reach the desired GCD. Two plates per condition were seeded. The plates were incubated at 37°C in humidified 5% CO_{2}. Every 24 hours, each plate was counted using a Cellometer Counter (Nexcelom Bioscience) four times (Fig. 2B; Supplementary Table S2).

### Time series imaging using confocal microscopy

A Zeiss 710 confocal microscope was used to image over time the entire well on an 8-well slide (ibidi) and form a 24-hour time series movie of cell movements and growth (several snapshots are provided in Supplementary Fig. S1). The tool employed on the Zeiss software was "Tiling," which allowed the microscope to take a series of spatial-sequence snapshots of the well and combine them together to create a single whole matrix image with a wide view of hundreds of cells (OVCAR-8-DsRed2 cell line; ref. 19). The plates were imaged using a Zeiss 710 confocal microscope with 20× objective, 10× *z*-stacks, and a 5 × 5 tile, to form two Time Series Imaging movies.

### Proliferation detection

#### Detecting Ki67 using FACS.

The percentage of Ki67-mediated cell proliferation as a function of LCD and GCD was measured by a staining method using the Ki67 protocol using an Anti-human Ki67 FITC, Flow Cytometry Staining Buffer, and Foxp3 Fixation/Permeabilization Solution (eBioscience) according to the manufacturer's instructions. OVCAR-8 were seeded in 150 mm dishes in drug-free medium and incubated at 37°C for 24 hours to allow cells to attach. Three different global cell densities (7%, 20%, and 33%) were seeded, each density in two different geometrical configurations. Although the random formation was mixed very well, the circular formation was 15 dots containing 2/3 × 10^{5}, 2 × 10^{5}, and 10/3 × 10^{5} cells per dot, respectively. The percentages of proliferation were then evaluated by FACS (BD LSRFortessa) every day for 4 days. Data from 10,000 gated events per sample were collected. Data were collected three times. For experiments with drug, after 24 hours with drug-free media, random and circular plates from each density were dosed with IC_{50} paclitaxel, 0.001 ± 6.20E−6 μmol/L and with IC_{50} × 4 paclitaxel. After 3 days of incubation with the drug, data were collected. See Supplementary Table S3 for the results.

#### Detecting Ki67 using confocal microscopy.

Cells were seeded at the desired density in an 8-well plate (ibidi). To form a circular structure for imaging, cloning cylinders (purchased from Sigma-Aldrich) were placed using forceps in the center of each well before seeding. Cells (4 × 10^{4}) per well were seeded, and after 24 hours, the cylinders were removed and the circular arrangement of the cells could be observed. After 48 hours, we followed the manufacturer's instructions to add the primary antibody (see details in Supplementary Materials). PBS (200 μL) was added to each well prior to imaging. The plates were imaged using a Zeiss 710 confocal microscope with 20× objective, 10× *z*-stacks, and a 3 × 10 tile to form Fig. 3B.

### Survival assay (MTT)

Cell survival was measured by 3-(4,5-dimethylthiazol-2-yl)- 2,5-diphenyltetrazolium bromide (MTT) viability assays. OVCAR-8 were randomly seeded at two different densities, 10^{4} and 10^{5} cells per well in 24 well plates. After 24 hours, all cells were dosed with 20 mmol/L paclitaxel in 16 different dilutions (Fig. 1A; Supplementary Table S1). Three hundred microliters of each dilution was added to each of the 24 wells. Once all the plates were dosed, they were incubated for another 72 hours. The MTT assay was performed according to the manufacturer's instructions (Molecular Probes). Absorbance values were determined at 570 nm on a Spectra Max 250 Spectrophotometer (Molecular Devices). All MTT assays were performed three times in triplicate.

### Apoptosis

Cells were seeded as in the Ki67 FACS experiments. The percentage of apoptosis-mediated cell death was measured by a double staining method using a Dead Cell Apoptosis Kit with Annexin V Alexa Fluor 488 & PI (Invitrogen) according to the manufacturer's instructions. The percentages of apoptosis and necrosis were then evaluated by FACS (BD LSRFortessa). Data from 10,000 gated events per sample were collected. Cells in early stages of apoptosis were positively stained with Annexin V, whereas cells in late apoptosis and necrosis were positively stained with both Annexin V and propidium iodide (PI). The number of cells positively stained with Annexin V was used to estimate the death rate. Data were collected three times; see Fig. 3C and Supplementary Table S4 for the results.

## Results

### Dynamic variation in cell-cycle duration

In this study, we focused on cell-to-cell dynamic baseline variations in relation to three key features of cancer cells: cell cycle, growth, and drug sensitivity. First, we explored the variation between cells in their CCD, as this affects growth and drug sensitivity. We sought to identify a set of cellular principles that would explain our experimental patterns of OVCAR-8 growth, reproduce the common normal-like distribution of CCD with its parameters, and analyze the implications of spatial variations in proliferation rates. To study these issues, we took a mathematical modeling approach to study the entire growth and treatment processes as a function of both time and space. We hypothesized that the local spatial component significantly affects the characteristics of this duration distribution.

Our spatial mathematical model is formulated as an agent-based model (ABM) that incorporates spatial aggregation/repulsion forces, as well as transitions through three cellular states (see Supplementary Materials; Fig. 1B). Cell movement is described via a stochastic differential equation, and transition rates are functions of LCD . Using our *in silico* model, we investigated the effect of imposing a deterministic cell-cycle structure based on the spatial location of the cell. In particular, we were interested in the distribution of the CCD as a function of time. The initial simulated GCD was *ρ*(0) = 10%, and the population of cells was allowed to grow until a time *t*, such that a plate was fully covered, that is, *ρ*(*t*) = 100%. The entire population of cells was initially assumed to be quiescent synchronized and randomly distributed in space. By calibrating the model with growth curves over 4 days (initial 10% and 80% GCD; Fig. 1C) and experimental time series imaging data (random, 10%) for a period of 24 hours, we imaged a wide spatial range of the cultured plate with hundreds of cells (Fig. 1D and Supplementary Fig. S1A and S2A). In this way, we were able to determine the parameter space of the model (see Supplementary Materials, Supplementary Fig. S1 and S2).

Our *in silico* results clearly showed that the basic assumptions of an intrinsic CCD as a deterministic process with growth limitation based merely on LCD can produce the normal-like distribution of CCD in randomly seeded cells. Moreover, as illustrated in Fig. 1E (where the CCD is modeled using equations S.9 and S.10 in the Supplementary Material), a full 3D histogram for the lengths of successful divisions for 10-hour windows in a single simulation demonstrated normal-like distributions that are functions of both time and space. For example, the first histogram labeled 15 on the “occurrence interval” axis denotes the distribution for all divisions occurring in simulation time *t* such that *t*
[15, 25] hours. We observed a traveling wave of normal-like distributions, starting with short cell cycles at low GCDs, which gradually increased as the population grew. Averaging over all time intervals yields a distribution with mean = 23 hours, similar to our previously observed experimental value (24.41 hours; ref. 4). Figure 1E illustrates the interesting temporal dynamics of the variance distribution (simulation), with a relative maximum variation occurring at approximately *t* = 80 hours, that is, at this point the plate reaches a 50% to 70% value of GCD.

### Local conditions create different global behavior

As LCD affects the distributional dynamics, we experimentally tracked the growth of cells over time with different initial geometric configurations, while the initial total seeding number of cells held constant. Throughout this article, we describe our work concerning two main geometrical configurations that represent the extreme cases in a range of commonly used seeding methods in any experimental cell culture laboratory (illustrated in Fig. 1E), which are random seeding (uniformly distributed across the plate), or circular seeding structure (can be achieved by direct application with a pipette without mixing the plate, illustrated in Fig. 2A). Common growth curves with shared initial GCDs of about 7% and 20% demonstrate the dramatic variations in growth over 10 days (Fig. 1F). The cells with an initial 20% GCD in random versus patterned circular seeding grew about 70% faster [i.e., (# cells_{Rand}/# cells_{circular} −1) × 100 = 70%] even during the first 4 days (Supplementary Table S2). These results show that not only the known global, but also LCD conditions can produce distinctive cell-to-cell baseline variations in growth.

Given that cell spacing had major effects on growth, we used our mathematical model to ask how the previous distribution of CCD would vary when the system grew in a specific geometrical configuration of circular (Fig. 2A, top), instead of random seeding (initially with 10% GCD, as before). The simulation predicts a wide and stable CCD distribution, with the variation values of most time points approximately the same. Interestingly, a binormal distribution with two distinct average values was obtained (Fig. 2A, bottom). This reflects the two subpopulations in the P state, the slow growing cells closer to the center and the faster, more proliferative cells growing at the periphery, when the area occupied by each subpopulation varies based on the GCD. In the case of low GCD, most cells are in a proliferative state, even ones closer to the center. However, they have different CCDs. On the other hand, once the GCD increases, most of the cells in the center move to a quiescent state, and the CCD for the rest of the P population changes based on the circular bands that the cells occupy. Thus, the two mean values of the binormal distribution approach one another and may be considered as a single distribution at later stages of the process. This process can be seen in Supplementary Movie S1, in which we simulated the spatiotemporal heterogeneity in CCD for a circular structure. Contact inhibition is well known to regulate the progression of the cell cycle. The contribution of this work is mainly to mechanistically and systematically predict the entire process of cancer cell dynamics as a population, based on local interactions, beginning with low cell density, growing with or without geometrical constraints. For cancer cells, cell-to-cell interactions slow cellular growth, but in many cases do not necessarily completely stop cell divisions, leading to heterogeneity. The unknowns are the extent and implications of such dynamic heterogeneity.

Using confocal microscopy imaging, for the circular pattern, after 48 hours, we observed more proliferative cells in the surrounding area, as determined by Ki67 staining (Fig. 2B). Note that it took 48 hours to reach this spatial distribution of proliferation. Prior to 48 hours, most cells showed random spatial proliferation (data not shown), as was predicted by the simulated movie. We further measured the cell-to-cell differences in the cell population distribution of Ki67 protein expression using flow cytometry, for two seeding patterns for 4 days at three different GCDs: 7%, 20%, and 30%. Three key FACS measurements were closely studied, the percentage of Ki67-positive cells, the Ki67 median intensity of all cells (positive and negative), and the shape of the Ki67 intensity distribution. These three measurements include information about the entire population, proliferative, nonproliferative state, and the transition between them over time. On average, the Ki67 distributions (example plotted in Fig. 2C) in populations with circular seeding at low GCDs were wider than when cells were seeded randomly, as was predicted by the simulation. In higher GCDs, we can observe higher fluctuations in proliferation over time at the random case (Fig. 2D). We found that the Ki67 percentages over 4 days for low (less than 30%) GCDs reached a maximum in random seeding, while a minimum was observed in circular seeding (Fig. 2E and F; Supplementary Table S3); altogether, the highest variation in proliferation percentages between the two geometrical configurations, can reach up to 36% over time [e.g., (#Ki67PositiveCells_{Rand_72h}/#Ki67PositiveCells_{circular_72h}−1) × 100 = 36%]. In other words, the baseline variations of both seeding methods are dynamic over time and the highest variations occur 72 hours after seeding. This is a growth feature related to the cells' adjusting process after initial seeding, which is also a function of local spatial constraints.

For cells that initially were seeded at low GCDs of less than 30%, there was a pattern that included the LCD and GCD. During the first 24 hours, cells are mainly in the process of attaching themselves to the plate. During the second 24 hours, that is, 24to 48 hours, cells mostly proliferate. During the third 24 hours, that is, 48 to 72 hours, cells increase their movement to contact other close neighboring cells. During the fourth 24 hours, that is, 72 to 96 hours, the plate was much more dense, changing the initial geometry. As a result, cells that were locally dense act differently from the random case over time. These variations can be seen in the changes of median fluorescence intensity as a function of seeding methods with varying GCDs (Fig. 2E) and variations as a function of time (Supplementary Table S3) or variations in Ki67-positive percentage as a function of seeding methods (Fig. 2F).

Taken together, this evidence implies that there are two main mechanisms that affect the population dynamic and thus its expected variations: the transition process between the P and Q states and the CCD (for cells in P state). We showed that the fluctuation in the transition process between states is as profound as the variations in CCD, especially for the geometrical cases, such as the circular rearrangement at low GCDs. However, in the case of randomly seeded cells, changes in the CCD are the main factors causing variations in spatial proliferation. These experimental results validate our understanding from the simulated predictions of a traveling wave versus stable distribution of CCD illustrated in Figs. 1E and 2A, which show the different mechanisms of the expected variations.

### Cell-to-cell baseline variation in drug-sensitive cells

Finally, we examined whether or not our finding concerning the role of LCD and cell movement as an intrinsic inducer of spatial proliferation variation can also explain the cell-to-cell baseline variation in drug sensitivity. We found that both killing curve features, the continuous curve and the distinct residual surviving cells after high drug doses, can be largely explained on the basis of the LCD mechanism, at least for chemotherapy drugs, such as paclitaxel. Paclitaxel is a mitotic inhibitor and hence acts on cells in the proliferative compartment P. More precisely, its primary chemotherapeutic effect is on the cytoskeleton, targeting tubulin (20). During cell division, chromosomal microtubules are necessarily destabilized and broken down as the mother cell divides. Paclitaxel acts as a stabilizer for these microtubules, thus preventing metaphase spindle configuration. At high paclitaxel doses, the cell is unable to complete mitosis and remains at the mitotic checkpoint of the cell cycle for a prolonged period of time. Eventually, prolonged activation of the mitotic checkpoint triggers apoptosis due to unsuccessful cell division. Using FACS, we quantified the fluorescence intensity of Ki67 staining for cells that were treated with an IC_{50} dose of paclitaxel for 72 hours. We observed that in cases of low GCD (less than 30%), there were differences in the median Ki67 values of the entire population (all cells included) between the two seeding methods: random versus circular (Fig. 3A; Supplementary Table S3). Randomly seeded cells expressed a much higher level of Ki67, consistent with increased proliferation (or at least staying in "cell cycle") at low GCDs. When high doses (IC_{50} × 4) were applied, much higher differences in the median occurred (Fig. 3A, first row). The range of median ratios for IC_{50} was 17%–47% (7% and 20% GCDs, respectively), whereas for IC_{50}×4, the range was 73%–88%. In addition, when we examined the percentage of only Ki67-positive cells, we found the same trend, approximately a 10% range in linear differences (Fig. 3A, second row), and about 13% in their ratio. Again, with randomly seeded cells, there are more positive cells than occur with the circular seeding method.

Furthermore, we hypothesize that different geometrical configurations sharing the same GCD can produce variations in actual apoptosis rates, and this intrinsic resistance can be traced back to cell-cycle dynamics. In general, increasing GCD is a well-documented cause of decreasing drug sensitivity. An example can be seen in Fig. 1A (Supplementary Table S1), where we performed MTT assays on initially randomly seeded cells with 10% versus 80% GCD of cells treated with paclitaxel. Killing curves showed the percentage of surviving cells after 72 hours of treatment as a function of the drug dosage, normalized by nontreated cells growing in parallel. The observed increase in IC_{50} at high GCDs is slightly more than one order of magnitude. To examine our theory about significant variations due to fluctuations in LCD, we utilized our ABM to mechanistically study the killing curves for three different spatial seeding configurations: *ρ*_{1}(0) = 10% (of GCD), and *ρ*_{2}(0) = 80% in random structure and *ρ*_{3}(0) = 10% in circular structure. The mathematical model qualitatively showed that the randomly seeded plate can reach a higher relative survival percentage for most drug concentrations and experience a more gradual decrease with increasing drug concentrations, that is, a more gradual slope (Fig. 3B). The simulations also predicted that the compact configuration (i.e., circular structure) would exhibit the same trend as high GCD, meaning that the LCD can produce significant higher survival with a more gradually sloped curve and higher number of surviving cells after exposure to high drug concentrations than randomly seeded cells that share the same GCD of 10%.

The model could not predict differences for low drug concentrations and was not sensitive enough to provide information about the observed IC_{50} variations. However, the model did provide direction based on the underlying dynamics to examine the general survival variation in a more detailed way, namely to focus on the relationship between processes related to CCD, the transition between the P and Q states, and the apoptosis rates as a function of the LCD. To experimentally quantify the cell-to-cell baseline variations in paclitaxel drug sensitivity, we further directly measured the apoptosis percentages under different local conditions: random and circular, using Annexin V staining, with the same GCD (20%). For the IC_{50} and IC_{50} × 4 settings, the variation between the random versus circular structures reached 12% to 14% linear differences in actual apoptosis percentage after normalizing each geometrical seeding method with their nontreated apoptosis rates (Fig. 3C; Supplementary Table S4). Therefore, our results suggest that features related to the changes in growth (e.g., CCD and transition between the P and Q states) and apoptosis rates are functions of the geometric cell configuration, that is, LCD, in addition to GCD.

## Discussion

Recently, Sandler and colleagues (21) addressed a central question concerning the origin of cell-to-cell variation in CCD by studying the differences between generations of cells derived from a parental cell. The main debate is whether observed variation originates from stochastic or deterministic processes. Although stochasticity, for instance, results from random bursts of gene expression or noisy partition of cellular components at division, deterministic characteristics may result from the robust design of the dynamic cellular system that can produce a stable oscillator. They concluded, on the basis of an experimental setting with very low cell densities, that the cell cycle is a deterministic and not a stochastic process. This leads to the conclusion that the internal noise, without high external noise, does not seem to add enough stochasticity to cause cell-to-cell variation in CCD.

Most experimental settings involve higher GCDs (30%–50%), yet have low LCDs (by randomly seeding the cells), especially in high-throughput screening experiments. In our current study, we suggest a method to integrate the conclusions from single-cell studies (such as those from Sandler and colleagues) with experimental conditions that include higher cell densities into a single unified model. This unified model can serve as a general reference model for cell biology studies that focus on dynamic features of single cells, population of cells, and the relationship between them. Theoretically and experimentally, we show that the LCD, in addition to the GCD, can lead to significant differences in spatial growth, proliferation, and paclitaxel-induced apoptosis rates. Another important aspect of this work is that it addresses a phenomenon that all cell biologists take for granted: killing curves are not discrete in nature, but rather are gradual, reflecting intrinsic heterogeneity in a cell population that can be attributed to differences in LCD and cell cycle. Furthermore, our model allows us to systematically explain why, and specifically how, one population can lose synchronization much faster than another.

This type of mechanistic model can be used to determine whether the observed system behavior or characteristics are the result of specific principles. For example, a future extension to our model could test the following questions: to what extent would additional internal noise (neglected here due to the Sandler and colleagues findings) change the population dynamics? Or what is the threshold level of collective internal noise that would begin to affect the entire system dynamics, in relation to spatial growth, proliferation, or drug response? Moreover, beyond a cell's response to a drug, our model has implications concerning other cellular processes, such as adult stem cell regeneration. Lei and colleagues have recently suggested the origin and consequence of heterogeneity in healthy tissues (22). They showed that heterogeneous proliferation, when not all cells have the same probability of proliferation, increases robustness in the process of tissue regeneration and replacement of defective cells and hence provides improved long-term performance.

In summary, a comprehensive reference model is essential in cell biology studies to determine significance. We stress the value of including in the reference model not only statistical tests, but also dynamic features of the process as an integral method.

## Disclosure of Potential Conflicts of Interest

No potential conflicts of interest were disclosed.

## Authors' Contributions

**Conception and design:** J.M. Greene, D. Levy, M.M. Gottesman, O. Lavi

**Development of methodology:** J.M. Greene, D. Levy, O. Lavi

**Acquisition of data (provided animals, acquired and managed patients, provided facilities, etc.):** S.P. Herrada, O. Lavi

**Analysis and interpretation of data (e.g., statistical analysis, biostatistics, computational analysis):** J.M. Greene, D. Levy, O. Lavi

**Writing, review, and/or revision of the manuscript:** J.M. Greene, D. Levy, S.P. Herrada, M.M. Gottesman, O. Lavi

**Administrative, technical, or material support (i.e., reporting or organizing data, constructing databases):** D. Levy, M.M. Gottesman, O. Lavi

**Study supervision:** D. Levy, M.M. Gottesman, O. Lavi

## Grant Support

This research was funded, in part, by the Intramural Research Program of the NIH and also supported, in part, by a seed grant from the UMD-NCI Partnership for Cancer Technology. The work of D. Levy was supported by the John Simon Guggenheim Memorial Foundation.

The costs of publication of this article were defrayed in part by the payment of page charges. This article must therefore be hereby marked *advertisement* in accordance with 18 U.S.C. Section 1734 solely to indicate this fact.

## Acknowledgments

The authors thank George Leiman, Carol Cardarelli, Susan Garfield, Langston Lim, Poonam Mannan, Karen Wolcott, Subhadra Banerjee, Mirit Aladjem, and Haiqing Fu for their assistance and for insightful discussions.

## Footnotes

**Note:**Supplementary data for this article are available at Cancer Research Online (http://cancerres.aacrjournals.org/).

- Received December 2, 2015.
- Revision received February 4, 2016.
- Accepted February 17, 2016.

- ©2016 American Association for Cancer Research.