Simon
Gordonov
ab,
Mun Kyung
Hwang
b,
Alan
Wells
c,
Frank B.
Gertler
bd,
Douglas A.
Lauffenburger‡
*abd and
Mark
Bathe§
*a
aDepartment of Biological Engineering, Massachusetts Institute of Technology, Cambridge, MA, USA
bThe David H. Koch Institute for Integrative Cancer Research, Cambridge, MA, USA
cDepartment of Pathology, University of Pittsburgh, and Pittsburgh VA Health System, Pittsburgh, PA, USA
dDepartment of Biology, Massachusetts Institute of Technology, Cambridge, MA, USA
First published on 26th November 2015
Live-cell imaging can be used to capture spatio-temporal aspects of cellular responses that are not accessible to fixed-cell imaging. As the use of live-cell imaging continues to increase, new computational procedures are needed to characterize and classify the temporal dynamics of individual cells. For this purpose, here we present the general experimental–computational framework SAPHIRE (Stochastic Annotation of Phenotypic Individual-cell Responses) to characterize phenotypic cellular responses from time series imaging datasets. Hidden Markov modeling is used to infer and annotate morphological state and state-switching properties from image-derived cell shape measurements. Time series modeling is performed on each cell individually, making the approach broadly useful for analyzing asynchronous cell populations. Two-color fluorescent cells simultaneously expressing actin and nuclear reporters enabled us to profile temporal changes in cell shape following pharmacological inhibition of cytoskeleton-regulatory signaling pathways. Results are compared with existing approaches conventionally applied to fixed-cell imaging datasets, and indicate that time series modeling captures heterogeneous dynamic cellular responses that can improve drug classification and offer additional important insight into mechanisms of drug action. The software is available at http://saphire-hcs.org.
Insight, innovation, integrationShape is a quantifiable property reflecting the phenotype of a cell that can change over time in response to extracellular cues and pharmacological perturbations. Characterizing the temporal changes of individual cells is essential to understanding cellular response in asynchronous, heterogeneous populations. This work presents an integrated experimental–computational framework to profile live-cell shape from high-content fluorescence imaging data. Probabilistic time series modeling is applied to classify single-cell shape dynamics as discrete morphological states explored in time. The framework reveals important similarities and identifies key differences in shape dynamics following inhibition of molecular regulators of cytoskeletal organization in breast cancer cells. Accounting for temporal dynamics of individual cells within heterogeneous populations helps to distinguish and interpret cellular response to diverse drug perturbations. |
However, fixed-cell assays, while relatively simple to perform through fluorescent staining and imaging, suffer from several important limitations. Principal among these is the loss of information regarding cellular dynamics in response to long-term or transient drug treatments. In addition, imaging artifacts may occur due to cell fixation and permeabilization, which may distort spatially resolved protein distributions.16 For these reasons, live-cell imaging is increasingly being used to characterize cellular phenotypes, particularly in the subcellular analysis of cell shape dynamics and polarization. For example, computational tools for cell boundary tracking,17–19 morphodynamics profiling,20–23 measurement of fluorescent reporters,24,25 and quantitative morphology and subcellular protein distribution analyses26 in live cells have become an integral component of high-resolution analyses of cell shape and its regulation, particularly in the context of cell migration. In cell migration studies, live-cell shape and signaling analyses have been complemented by direct quantification of motility properties such as cell speed and persistence of motion to establish links between molecular mechanisms and migratory phenotypes.27–32
In these applications, the relative strengths of high-resolution, live-cell imaging versus fixed-cell HCI assays are apparent: the former captures rich, dynamic properties of single-cell behavior while the latter enables large-scale screening of hundreds to thousands of cells. In an effort to bridge this gap, several mathematical approaches have been developed to infer dynamic properties of cell populations from fixed-cell measurements in HCI studies. For example, ergodic rate analysis based on differential equation modeling has been used to infer transition rates through cell cycle stages from images of molecular reporters that define various mitotic phases in individual fixed cells.33 Additionally, Bayesian network modeling of shape parameters coupled with RNAi knockdown of cytoskeleton-regulatory proteins has been used to infer shape state transitions of migratory cells and reveal underlying regulatory signaling modules.34,35 However, these approaches assume quasi-steady-state of the cell population, assign cells into pre-defined phenotypic categories, and, in the case of Bayesian networks, face difficulties in modeling repetitive processes such as motility cycle stages in migrating cells. Moreover, they are not directly applicable to the analysis of live cells over time to monitor individual cellular responses to drug perturbations.
To address these limitations, here we present a live-cell HCI framework that captures the dynamics of a large number of cells on the scale of a phenotypic screen. The approach combines high-content live imaging, image processing, multivariate data analysis, and probabilistic modeling to characterize cell shape dynamics in a drug-screening context. Inspired by existing methods for modeling cell cycle stages from time series images,36,37 our framework employs hidden Markov modeling to describe shape dynamics as temporal sequences of morphological states that are observed as noisy, multivariate image data. Describing temporal trajectories of multivariate shape measurements as a finite set of states that is limited in number by Bayesian model selection that penalizes model complexity, results in an efficient description of time-dependent shape categories explored by a cell. This approach provides a means of modeling the shape dynamics of hundreds of individual cells, capturing temporal evolution in cell morphology directly from live images without assumptions of steady-state cell populations or predefined shape categories assumed by fixed-cell analyses. We show here in a proof-of-principle study that drug-response profiles derived from these models can be used for phenotypic drug comparisons and can reveal spatially distinct and pathway-specific roles that drug-targeted species play in modulating shape dynamics. Our computational framework is available as open source for the computational cell biology community to apply using MATLAB at http://saphire-hcs.org.
In one experiment, we included two controls in order to set a baseline for cell behavior under no drug treatment: one of growth media only and the other with 0.1% v/v DMSO in growth media to assess effects of DMSO alone since it was used to dissolve drug stock. We also included two different doses for the ROCK and MLCK inhibitors to assess drug concentration effects on cell shape dynamics characterized in our framework. In an additional experiment, we expanded the panel of drug treatments and used a different microscope, culture media, and cell passage number to assess these factors on reproducibility of experiments and subsequent analyses. Epifluorescence 2-D microscopy using a 10×/0.3NA air objective offered a suitable tradeoff between field-of-view and resolution, enabling us to maximize the number of cells captured with sufficient detail to characterize individual cell morphologies over approximately eighteen hours for each treatment condition.
Visual inspection of the compiled time series movies revealed asynchronous populations of dividing, apoptotic, and migratory cell phenotypes (Videos S1–S3, ESI†). For cells not undergoing division or death, we observed a broad range of shape changes, with some cells changing shape rapidly with pronounced cell body protrusions and retractions, while others exhibited less noticeable morphodynamic activity. Although the cells generally exhibited autonomous behavior that was largely independent of one another, we observed frequent contact and spatial overlapping of some cells in addition to entry and exit of cells into and out of the imaging fields of view. Collectively, the large quantity of imaging data and phenotypic cell behaviors rendered fully manual data parsing and analysis intractable, thereby motivating the need for an automated computational analysis pipeline for processing the time-series images.
To visualize how morphologies vary in different directions and radial distances away from the data mean in shape-space, we created a polar representation of the two PCs. The polar shape-space was subdivided into twelve equal angular bins, with each bin subdivided further into quartiles. Fig. 3C shows the representative cell shapes in each radial quartile of the angular bins as well as the contributions of the original eighteen shape features to the values of the two PCs. The highest variance in shape, along the first PC, predominantly captures cell elongation and branching, while the second PC captures cell spreading area and roundness. Features such as major axis length, geodesic diameter, and maximum Feret length capture similar length properties of cells and therefore have similar PC coefficients and directions in shape-space. These features are anti-correlated with solidity, circularity, and extent, which point in the opposite direction for both PCs, showing that branched and elongated morphologies have smaller solidities and circularities, as expected. As captured mainly by the second PC, cells with longer minor axes and larger areas have generally lower eccentricities and smaller ratios of major to minor axis lengths.
Furthermore, we found that cells are densely packed around the data mean at the origin in shape-space, with the point density dissipating radially (Fig. 3A). This suggests that although MDA-MB-231 morphologies are visually distinct, shape properties of the population as a whole vary on a continuum, with no clearly distinguishable groups. Interestingly, however, when we plotted cell shape trajectories one at a time we found increased grouping and clustering of points. We applied k-means clustering to individual cells and to the cell population as a whole for k = 2 to k = 10 clusters, and computed average silhouette values in order to quantify the degree to which points form well separated, or cohesive, groups (Fig. 3D). As may visually be discerned in Fig. 3A, example cell trajectory i forms 2 groups while cell trajectory j forms 3 groups, which is confirmed by the maximum silhouette value that occurs for k = 2 and k = 3 clusters, respectively. On average, points from individual cell trajectories had higher cluster cohesion than did random samples with the same number of points as the trajectories or the entire cell population as a whole.
To compare our approach with GMMs that ignore dependencies in time series data, we simulated shape-space trajectories from underlying hidden Gaussian states to explore differences with HMM-based annotation using SAPHIRE. Specifically, we sought to compare state identification capability of the two methods, with SAPHIRE using Bayesian model selection and GMM using the Bayesian information criterion to penalize model complexity. Our two-state model simulations revealed that SAPHIRE is more likely to choose the model with the correct number of states and has better accuracy in inferring locations of state means compared with GMM due to SAPHIRE's incorporation of temporal dependencies in the data during the state inference steps (Fig. S1, ESI†).
Next, we applied SAPHIRE to experimental observations of cells in the two-dimensional PC shape-space to infer a time series model for each cell individually. The inferred model specifies the most probable number of hidden shape states, which is unknown a priori, as well as the state parameters. Each state is a symmetric, bivariate Gaussian distribution in shape-space with two variables corresponding to the two PCs. The inferred parameters for each Gaussian distribution are the standard deviation and the mean, which capture temporal shape variation and average morphology of the cell within the state, respectively. The model also specifies a state transition probability matrix that describes the probability of transitions between shape states. Each cell is annotated in time with the most probable, or maximum likelihood, shape state sequence, enabling us to determine the most likely state that the cell exists in at each point in time. We first compared shape state inference of our approach with that of commonly-used GMM under different model constraints. Allowing for elliptically-shaped Gaussian states in the GMM led to undesirable grouping of cell shapes with visually-dissimilar morphologies into the same state (Fig. S2, ESI†), and, similar to the results of the numerical simulations (Fig. S1, ESI†), SAPHIRE was better able to capture distinct morphological phases of cells that gradually move between distinct regions of shape-space (Fig. S3, ESI†). In certain cases of cells moving between two visually distinguishable underlying states, the state annotations of cell trajectories were similar for SAPHIRE and GMM (Fig. S4, ESI†). Moreover, from all modeled cell trajectories we found that a larger fraction of cells exists in two to three underlying states, while fewer cells explore either a single state or four or more states (Fig. S5, ESI†).
Despite some similarities in the numbers of states explored by cells, the inferred parameters of transition dynamics and state annotations varied considerably between cells. These variations not only existed for cells across different drugs, but also across cells within a given treatment and for cells with the same number of shape states. Fig. 4 shows the inferred shape states for two example cells treated with the same drug having the same number of inferred states, but with notable heterogeneity in state parameters and transition dynamics. Some cells exhibited rapid back-and-forth switching between states corresponding to periodicity in actin protrusions and retractions, resulting in shapes that are elongated, rounded, or those in between (Fig. 4A). Other cells exhibited individual instances of state transitions, progressively changing shape in a given direction in shape-space, such as going from larger and more spread to smaller and rounder morphologies (Fig. 4B). The maximum likelihood state annotations capture different phases in the morphodynamic history of a cell, with state transition parameters derived from the state sequence providing useful insights into its dynamic behavior. For instance, the model for the cell in Fig. 4A reveals that the cell is more likely to stay in state 1 (p = 0.19) or state 2 (p = 0.29) than transition to other states, and that direct transitions between elongated (state 1) and a rounded (state 3) morphology without going through an intermediate shape (state 2) is not likely (p = 0).
We next assessed how the distribution of cellular shape states in polar PC space was affected by the experimental treatment conditions (Fig. 5B). Cells in the DMSO and growth media controls were positioned closely and were fairly evenly distributed around the mean of the shape-space data without exhibiting biased morphologies in any particular direction compared with the morphologies induced by drugs. The ROCK and myosin II inhibitors pushed cells predominantly towards highly elongated and branched morphologies with longer dwell times in these states. On the other hand, MLCK inhibition had the opposite effect, biasing cells toward smaller and rounder morphologies along the negative first PC axis. In the experiment with the expanded panel of drugs, both MEK inhibitors tested led to a broader variety of states, predominantly either rounder, or more elongated, shapes, although no strong biases of shape state location or state dwell time were observed upon MEK inhibition (Fig. S6A, ESI†). Similarly, EGFR inhibition led to a broader and relatively more uniform distribution of morphological states around shape-space, albeit with a noticeable shift towards less elongated morphologies, whereas calpain inhibition biased cells into more elongated morphologies. Further, doubling the dose of the myosin II inhibitor shifted cells slightly towards more branched morphologies while diminishing elongated states, whereas an increase in dose of the MLCK inhibitor appeared to accentuate the bias towards smaller and more rounded shapes (Fig. 5B). For either drug, however, the overall distributions of cellular shape states and state dwell times were similar across both doses tested, as well as between the two imaging experiments.
Live-cell analysis also enables the determination of whether drugs differentially affect the trajectories that cells take through shape-space when they transition between distinct morphological states. To our surprise, state transition directions were similar across the drugs tested, with cells moving between two relatively narrow angular ranges, 120 to 180 degrees or 300 to 360 degrees in shape-space (Fig. 5C and Fig. S6B, ESI†) that correspond to decreased elongation with increased roundness versus increased elongation and decreased roundness, respectively (Fig. 3C). This finding suggests that most of the dynamics in shape that MDA-MB-231 cells undergo are along this morphological axis of increasing or decreasing elongation, regardless of drug treatment, which may be indicative of cytoskeletal protrusion and retraction cycles, while transitions toward larger, spread morphologies, for example, are quite rare for this cell type. Despite the similarity in the directions of state transitions taken by cells across drugs, ROCK and myosin II inhibition led to larger magnitudes and slightly broader distributions in state transitions, suggesting that these drugs induce more pronounced variations along the roundness-elongation axis of cell shape, likely by reducing transcellular actomyosin tension.
Moreover, readouts of shape state locations and state transitions can be used to group treatments based on similarities in induced phenotype dynamics. We derived a dynamic shape state “profile” for each treatment from the shape state and state transition histograms (Fig. 5B, C and Fig. S6A, B, ESI†) and hierarchically clustered the treatments using profile similarities for each imaging experiment separately (Fig. 6A and B). As anticipated based on the preceding results, the two MEK inhibitors clustered closest together, as did the DMSO and growth media controls. These results serve as internal controls to validate the HMM annotation and phenotypic drug comparison procedure proposed here. The MEK inhibitors were also found to have dynamic shape state profiles more similar to those of the EGFR and calpain inhibitors than to the MLCK, ROCK, and myosin II inhibitors. Myosin II inhibition with Blebbistatin produced phenotypes most similar to those with ROCK inhibition. Perhaps surprisingly, MLCK inhibition, which is known to alter myosin II activity, induced shape dynamics more similar to those under MEK and EGFR inhibition than under myosin II and ROCK inhibition. Overall, the computational analyses from these imaging experiments demonstrate that temporal dynamics of individual cells can be combined into quantitative profiles that serve as useful readouts for phenotypic drug comparison and for inferring shape-regulatory roles of targeted signaling molecules.
For the three treatments tested that were amenable to the classification analysis approach as in ref. 15, of the four fixed-cell methods the “Factor Analysis + Means” profiles performed best, correctly classifying 5 of 6 treatments (Fig. 7A and B). This approach has also previously been shown to have good classification performance in a large-scale screen of dozens of drugs from different mechanistic classes.15,44 Profiles generated using Gaussian mixtures43 and the K–S statistic10 correctly classified 4 of 6 treatments, whereas the simplest “Means”42 approach only classified 2 of 6 treatments correctly. Treatment profiles generated using SAPHIRE resulted in correct classification of all 6 treatments, demonstrating improvement over the fixed-cell profiling methods tested. Moreover, the HMM annotation of morphological states using SAPHIRE appears to be critical for improving correct treatment classification, as profiles derived from single-cell temporal dynamics without the HMM misclassified 2 of the 6 treatments.
To further understand which properties of the state-space treatment profiles generated using SAPHIRE improve treatment classification, we separately considered the use of features that capture state transitions versus features that only consider properties of the states themselves (Fig. 7C). Treatment classification confusion matrices using these state features separately revealed that state transitions correctly classify all six treatments, whereas exclusion of temporal transition information leads to misclassification of a treatment, yielding a classification performance similar to the fixed-cell “Factor Analysis + Means” profiling approach.
Collectively, these results demonstrate that all of the methods implemented, with the exception of the “Means” approach, correctly classify the MLCK and myosin II inhibitors, but fail to differentiate between cells treated with MEK inhibitor versus control, or those treated with ROCK versus myosin II inhibitors. Thus, as shown in Fig. 7D, when the shape distributions of two treatments differ, both state-space modeling of individual cells over time and existing fixed-cell profiling methods that measure features of different cells at multiple time points, can correctly classify and resolve treatment effects on cell shape. When shape distributions of two treatments are similar, however, state-space modeling of single-cell temporal transitions within the population distributions can improve discriminability and classification of treatments compared with existing methods that only capture shape properties of distinct cells within a population. Therefore, annotation of cell morphologies with a state-space representation using HMM, and in particular capturing state transition dynamics on a per-cell basis, can improve the accuracy of classifying treatment conditions in an HCI experiment.
It has been proposed that cell populations can assume either discrete or continuous morphological landscapes, principally determined by cell type, extracellular environment, types of perturbations, and the particular phenotype captured.45 For example, in cases when particular morphologies are reached under effect of genetic perturbations, cells may exhibit more stable steady-state morphologies leading to discrete shape states.13 Our results suggest, however, that shorter-term responses to pharmacological perturbation of highly asynchronous cell populations fall on a more continuous, finely graded morphological spectrum of shapes that conceals more cohesive and well separated shape state clusters within temporal trajectories of individual cells (Fig. 3). Consequently, two different approaches to cell profiling, one using commonly-applied fixed-cell analyses of populations of distinct cells and the other using temporal modeling of individual live cells presented here, may collectively serve as a powerful combined approach to reveal idiosyncrasies in phenotypic effects of different molecular perturbations (Fig. 8).
Fig. 8 Schematic of the conceptual difference between fixed-cell analyses versus SAPHIRE, and their combined insights into molecular mechanisms of cell shape regulation. The distribution of morphologies of drug-treated MDA-MB-231 breast cancer cells falls on a continuum, illustrated in gray. The extent of cell elongation and branching is the dominantly varying morphological property of the breast cancer cell line analyzed in our drug screen (Fig. 3). Fixed-cell profiling captures the morphological landscape of distinct cells in the population at a given time point. In contrast, SAPHIRE captures the morphological evolution of individual cells over time. SAPHIRE uses an HMM framework to detect the presence of morphological states (e.g. orange and green states for one cell, blue and pink states for another) used to model the temporal dynamics of each cell independently within the population. Application of both of these approaches to live-cell shape profiling reveals that MLCK inhibition decreases, while ROCK and myosin II inhibition increases, cellular elongation and branching. |
In particular, within the context of image-based phenotypic profiling, our results illustrate that state-space temporal modeling of cell shape using HMM can improve the resolution of cellular classification in response to distinct drug treatments. This result, shown in Fig. 7D, highlights an important practical outcome of this work in the context of HCI. Existing phenotypic profiling methods correctly classified some treatments and misclassified others in the actomyosin cytoskeleton-focused drug screen. Our analyses demonstrate that modeling morphological transitions over time on an individual cell basis using HMM can provide useful live-cell phenotypic information for discriminating between treatment effects. Some of these effects are not captured by simply modeling cell shape distribution properties of different cells in a population, even if the same population is profiled over time. This suggests that the framework presented offers additional information beyond fixed-cell profiling in image-based classification using morphological data. An interesting application of temporal modeling would be to complement fixed-cell assays that utilize immunofluorescence staining in order to augment phenotypic profiling and better resolve underlying molecular differences driving distinct cellular drug responses assayed using live imaging.
In addition to investigating the utility of temporal shape modeling in HCI applications, we characterized differential cell shape dynamics induced by small-molecule inhibitors in order to infer relationships between actomyosin signaling mechanisms and shape phenotypes. Similarities in shape dynamics between DMSO and non-DMSO controls, between the two MEK inhibitors, and under different doses of the same drugs, validates the morphological states and state transitions as fundamentally reflecting target-specific effects. Similar dynamic responses in cell shape induced by EGFR, MEK, and calpain inhibition suggest that these species function along a common signaling axis of cell shape regulation in MDA-MB-231 breast cancer cells (Fig. 6B). This result is in concordance with previous molecular mechanistic studies performed in fibroblasts, which showed that m-calpain activity is directly regulated by Erk by altering its spatial localization and association with PI(4,5)P2 at the plasma membrane.46 This process can be mediated by EGFR signaling through PLC-γ, which depletes PI(4,5)P2 at the leading edge resulting in localization and activation of m-calpain at the trailing end that plays an important role in migratory cell polarization. The increased elongation phenotype induced by calpain inhibition revealed from our analysis (Fig. S6A, ESI†) reflects the role of m-calpain in adhesion remodeling, the inhibition of which leads to trailing end retraction defects.47,48
Moreover, our observation that ROCK and myosin II inhibitors biased cells toward branched and elongated morphologies is consistent with reported effects of increased actomyosin tension promoting cell polarization and limiting the formation of new protrusions.49–52 The increased elongation and branching dynamics we observed reflect the impairment in actomyosin contractility and loss of cell polarity upon ROCK and myosin II inhibition, leading to multiple competing protrusions extending outwards around the cell periphery. Given that MLCK and ROCK are both known to induce myosin II phosphorylation and activation, we may have expected their inhibition to promote similar shape dynamics. Instead, the surprisingly dissimilar and near-opposite morphologies that their inhibition produced (Fig. 5B) suggests that MLCK and ROCK differentially affect cytoskeletal dynamics; this notion is supported by previous studies which showed that MLCK and ROCK have spatially distinct roles in regulating myosin II activity, and whose activities have opposite effects on the number of lamellipodial protrusions, cell elongation, and polarization.50,53 The MDA-MB-231 cell shapes induced by MLCK inhibition in our work are reflective of single-front motile cell morphologies for which formation of additional lamellipodial protrusions is limited (Video S3, ESI†), as has been suggested to also occur in keratocytes at different stages of development.53 The differential effects of ROCK and MLCK inhibition on shape dynamics that we identified through our analyses (Fig. 8) may have consequential effects on overall cell polarization and migratory behavior in cancer that warrants further investigation in future work.
More generally, the plasticity and variability in shape revealed in our analyses reflects the highly dynamic cytoskeletal changes characteristic of protrusion and retraction events, such as those of migratory cells undergoing directional changes or motility cycles. In light of this, further exploration of interest would be to characterize and understand what cellular functions the morphological states captured by HMM annotation may represent. Qualitative visual inspection suggests that cells in rounder states have higher amounts of cortical actin, whereas those in more elongated states are polarized and may therefore be migratory (Fig. 2B and Videos S1–S3, ESI†). Future studies will focus on characterizing morphological state properties in more depth by quantifying cytoskeletal organization in addition to whole-cell shape, incorporating signaling status from live or fixed-cell fluorescent reporters in post-capture analysis,24,54 and combining model-annotated state-space dynamics with motility measurements to establish temporal relationships between morphology and migratory behavior.27,55,56
Notwithstanding, the current work demonstrates that phenotypic information derived from imaging-based models of temporal shape dynamics reflects underlying signaling pathway activities that regulate cell shape, which enables hypothesis generation for more in-depth, mechanistic follow-up studies. In vitro profiling of phenotype dynamics is also relevant in the context of early-phase drug discovery. In that context, time series modeling of live-cell phenotypes can complement existing fixed-cell imaging methods as well as genetic and biochemical approaches for understanding cellular drug responses, particularly with the advent of genome editing approaches such as CRISPR/Cas9.57–61 The collective implementation and application of these approaches will become an increasingly important component of integrative drug profiling and holistic understanding of cell function.
Two imaging experiments were performed: for the experiment with the expanded panel of drugs cells were imaged on an IncuCyte ZOOM system incubated at 37 °C and 5% CO2 under a standard scanning protocol (Essen Bioscience) in drug-containing DMEM growth media (see above). Cells were imaged over approximately 18 hours at 20 minute intervals, with three wells per drug condition and four fields of view per well, producing an image data set of 84 fields (4452 time series image frames), from which 293 individual cell trajectories, with an average of 53 time frames, were obtained for further analyses. For the other experiment, which contained the experimental controls, cells were imaged in Leibovitz's L-15 media (Life Technologies) supplemented with the same additives as the DMEM culture media above, on a Nikon Eclipse Ti microscope equipped with an Andor Zyla sCMOS camera. Cells were imaged over approximately 18 hours at 8 minute intervals, with three wells per drug condition and four fields of view per well, producing an image data set of 96 fields (16128 time series image frames), from which 435 individual cell trajectories, with an average of 100 time frames, were obtained for further analyses. All inhibitor stocks were dissolved in dimethyl sulfoxide (DMSO) (Sigma), other than Y-27632, which was dissolved in high-purity water. All final concentrations of inhibitors in culture media used for imaging were at 0.1% v/v DMSO or lower. In both experiments cells were imaged using a 0.30 NA Nikon Plan Fluor 10× air objective. The inhibitors used in this work and their vendor sources were as follows: AZD6244 (MEK inhibitor; Selleck Chem.), ML-7 (MLCK inhibitor; Enzo Life Sciences), Blebbistatin(+/−) (non-muscle myosin II inhibitor; Enzo Life Sciences); Gefitinib (EGFR inhibitor; LC Labs), PD0325901 (MEK inhibitor; LC Labs), Y-27632 (ROCK inhibitor; Enzo Life Sciences), PD150606 (Calpain inhibitor; EMD Millipore).
θk = [{μx,i,μy,i,σi}Ki=1, {πi}Ki=1, {ϕij}Ki,j=1] |
For each inferred Gaussian state si in shape-space from a temporally annotated sequence, we convert the coordinates of the state mean value in the Cartesian PCA axes (μx,si, μy,si) to polar coordinates (μr,si, μθ,si). As a result, rsi becomes the location vector of the state in shape-space and |rsi| is the radial distance of the state from the origin, where |·| denotes vector magnitude. θ is the angle rsi makes in the counterclockwise direction with the positive PC1 axis for state si. For each state si the state radial distance is described by |rsi|, or equivalently . The normalized state dwell time for si is calculated as , where δt,si = 1 if the state is si in the model-annotated state sequence at time t, and 0 otherwise, with T equal to the number of elements in the sequence (i.e. number of time frames in the cell trajectory). For a given cell trajectory, we derive an angular histogram with twelve bins on the range of 0 to 360 degrees, with the first bin 0 to 30 degrees, second bin 30 to 60 degrees, and so on until the twelfth bin that ranges from 330 to 360 degrees. For all states si ∈ S in the cell sequence, the state radial distance value of each bin is computed as 〈|rsi|〉 where 〈·〉 denotes the mean, for all states si with μθ,si falling within the angular range of the bin, to generate a histogram of state radial distances. A state radial distance profile for a treatment condition is then computed as the average within each bin in the twelve-bin state radial distance histograms for all cells in that treatment condition. The state dwell time profile of the treatment condition is similarly computed except averaging normalized state dwell time values for states of each cell instead of the state radial distances for twelve angular bins.
We additionally derive state transition profiles for each treatment condition. Consider for a given cell a state transition from si to sj. The state transition direction θsji for sj from state si is the angle that the vector rsj − rsi makes in the counterclockwise direction with the positive PC1 axis in shape-space. The state transition magnitude for sj from state si is computed as |rsj − rsi| scaled by the normalized state transition dwell time. The normalized state transition dwell time for sj from state si is calculated as the number of time points a cell spends in state sj after transitioning to it from state si, divided by the sequence length, T. The state transition dwell time value of each of the twelve angular bins is the mean of the normalized state transition dwell times of all states whose means (μr,si, μθ,si) are in the bin. For transition directions θsji for all i ≠ j and a given cell, the state transition magnitude value of each of the twelve angular bins is computed as 〈|rsj − rsi|〉 for all θsji falling within the angular range of each bin, to produce a histogram of state transition magnitudes. A state transition magnitude profile for a treatment condition is then calculated as the average within each bin in the twelve-bin state transition histograms for all cells in that treatment condition. A state transition dwell time profile is similarly computed except averaging normalized state transition dwell times across cells for the twelve angular bins instead of the transition magnitudes.
To quantify similarities in effects of treatment conditions on cell shape dynamics, clustering using average linkage of Euclidean distances between phenotypic profiles was performed for treatment pairs to generate cluster dendrograms. Internal nodes were re-sorted into optimal leaf order without dividing the clusters or changing overall tree connectivity so that treatments with similar profiles are next to each other by maximizing the sum of similarities (, where i is a leaf adjacent to leaf j) between adjacent tree leaves.
Permutation testing was used to assess significance of Euclidean distances (similarities) between all pairs of treatments as follows. Dynamic phenotype profiles (48-element vectors) based on state locations and state transitions and corresponding dwell times (see above) of individual cells within any two treatment conditions being compared were randomly reassigned into the two treatments while preserving the number of cells and their signatures. Signature values were averaged between the randomly-assigned cells for a given treatment to generate a treatment profile. A permutation test p-value for similarity between two treatment profiles was quantified by comparing the actual Euclidean distance between the treatment profiles relative to the null distribution of distances between profiles that were obtained by repeating the cell assignment randomizations and distance calculations 10000 times for each pair of treatments.
Briefly, the “Means” approach averages each of the 18 cell shape features for each treatment condition for the five time points to generate a phenotypic profile for a treatment. For the “K–S Statistic”, the cumulative distribution function (cdf) is compared between the treatment and DMSO control cells for each shape feature separately. A signed Kolmogorov–Smirnov (K–S) statistic is then computed for each treatment, which is equal to the maximum distance between the two cdfs and set to positive if the treatment cdf is above the control cdf and negative otherwise. A “K–S Statistic” profile for a treatment is a vector of concatenated signed K–S statistics for each shape feature. For the “Factor Analysis + Means” method we first performed factor analysis on cells from all treatments at the five time points, selecting the number of factors using the Kaiser criterion, and then computed the average value of the scores of each factor for the cells in each treatment condition to generate a profile for that treatment. Finally, for the “Gaussian Mixture” method, a GMM with no covariance matrix constraints was fit to cell observations from all treatments at each time point separately using all shape features as variables, testing models with 2 to 30 mixtures, and selecting the best-fitting model using BIC. For each treatment, the posterior probability of each cell under that treatment belonging to each of the mixtures was computed and averaged across cells for each mixture to generate a profile at a given time point. The averaged posterior probabilities of all mixtures and time points were combined to generate a phenotypic profile for each treatment.
Treatment classification accuracy was assessed for each of the profiling methods described above, as well as using the profiles derived from the single-cell temporal modeling using SAPHIRE. To assess the value of the HMM in treatment classification, treatment profiles extracted from shape-space analyses of cellular dynamics without the HMM annotations were generated identically as for SAPHIRE but by treating each time point in a cell trajectory as a separate “state” (i.e. no modeling of PCA coordinates from underlying hidden Gaussian states). Treatment condition classification accuracy was assessed by computing pair-wise Euclidean distances between phenotypic profiles of treatments using each of the profiling approaches separately. A treatment was designated as being correctly classified if its phenotypic profile was closest, in terms of distance, to that of the same treatment at a different dose, compared to the profiles from other treatments.15
BIC | Bayesian information criterion |
EGFR | Epidermal growth factor receptor |
GMM | Gaussian mixture modeling |
HCI | High-content imaging |
HMM | Hidden Markov modeling |
PC | Principal component |
PCA | Principal component analysis |
MEK | Mitogen-activated protein kinase kinase |
MLCK | Myosin light chain kinase |
ROCK | Rho-associated protein kinase |
SAPHIRE | Stochastic annotation of phenotypic individual-cell responses |
Footnotes |
† Electronic supplementary information (ESI) available: Fig. S1–S6, Table S1, Videos S1–S3. See DOI: 10.1039/c5ib00283d |
‡ 77 Massachusetts Avenue, Building 16, Room 343, Cambridge, MA 02139 USA. E-mail: E-mail: lauffen@mit.edu; Fax: +1-617-258-0204; Tel: +1-617-252-1629 |
§ 77 Massachusetts Avenue, Building 16, Room 255, Cambridge, MA 02139 USA. E-mail: E-mail: mbathe@mit.edu; Fax: +1-617-324-7554; Tel: +1-617-324-5685 |
This journal is © The Royal Society of Chemistry 2016 |