Figures
Abstract
Apoptotic cell death can be initiated through the extrinsic and intrinsic signaling pathways. While cell cycle progression promotes the responsiveness to intrinsic apoptosis induced by genotoxic stress or spindle poisons, this has not yet been studied conclusively for extrinsic apoptosis. Here, we combined fluorescence-based time-lapse monitoring of cell cycle progression and cell death execution by long-term time-lapse microscopy with sampling-based mathematical modeling to study cell cycle dependency of TRAIL-induced extrinsic apoptosis in NCI-H460/geminin cells. In particular, we investigated the interaction of cell death timing and progression of cell cycle states. We not only found that TRAIL prolongs cycle progression, but in reverse also that cell cycle progression affects the kinetics of TRAIL-induced apoptosis: Cells exposed to TRAIL in G1 died significantly faster than cells stimulated in S/G2/M. The connection between cell cycle state and apoptosis progression was captured by developing a mathematical model, for which parameter estimation revealed that apoptosis progression decelerates in the second half of the cell cycle. Similar results were also obtained when studying HCT-116 cells. Our results therefore reject the null hypothesis of independence between cell cycle progression and extrinsic apoptosis and, supported by simulations and experiments of synchronized cell populations, suggest that unwanted escape from TRAIL-induced apoptosis can be reduced by enriching the fraction of cells in G1 phase. Besides novel insight into the interrelation of cell cycle progression and extrinsic apoptosis signaling kinetics, our findings are therefore also relevant for optimizing future TRAIL-based treatment strategies.
Author summary
TRAIL (TNF-related apoptosis-inducing ligand) induces cell death preferentially in cancer cells. Whether cell cycle progression notably affects extrinsic apoptosis has remained unclear, since systematic experimental and mathematical studies to quantitatively understand such interdependencies remained challenging. Here, we applied statistical and mechanistic modeling, linked with experimental analyses, to determine if cell cycle and apoptosis progression are interconnected. Using sample-based modeling, we demonstrate that times required to commit apoptotic cell death depend on the cell cycle position at the time of TRAIL exposure, with delays manifesting during S phase, around a time that can be defined as a point of apoptosis deceleration (PAD). Overall, cells receiving TRAIL in the G1 phase were more likely to die, providing scope to optimize TRAIL-based treatment strategies with respect to cell cycle dynamics.
Citation: Imig D, Pollak N, Allgöwer F, Rehm M (2020) Sample-based modeling reveals bidirectional interplay between cell cycle progression and extrinsic apoptosis. PLoS Comput Biol 16(6): e1007812. https://doi.org/10.1371/journal.pcbi.1007812
Editor: Attila Csikász-Nagy, King’s College London, UNITED KINGDOM
Received: June 7, 2019; Accepted: March 23, 2020; Published: June 4, 2020
Copyright: © 2020 Imig et al. This is an open access article distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited.
Data Availability: All relevant data are within the manuscript and the repository: https://github.com/DirkeImig/CellCycleApoptosis.
Funding: DI would like to thank the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) for financial support under Germany’s Excellence Strategy-EXC 2075-390740016 at the University of Stuttgart. MR additionally receives support from the German Research Foundation (FOR2036 (MO 3226/1-1)). The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript.
Competing interests: The authors have declared that no competing interests exist.
Introduction
The cell cycle regulates the expression of numerous proteins [1, 2], and consequently affects various cellular processes, including for example glycolysis and oxidative phosphorylation [3]. Fluorescent, ubiquitination-based cell cycle indicator (Fucci) reporter systems [4] allow distinction between G1 and S/G2/M phases by coupling fluorescent markers to proteins that accumulate in specific cell cycle phases. As a result, start and end points of Fucci phases can be tracked in parallel to other processes such as cell death of single cells in long-term time-lapse microscopy experiments [5]. A connection between cell cycle progression and apoptosis, for example induced by irradiation or chemotherapeutics, was proposed years ago due to similar morphology and involvement of mutual genes [6]. Long-term time-lapse experiments already provided insight into the cell cycle dependency of intrinsic apoptosis [7–12]. Concerning TRAIL-induced extrinsic apoptosis, the role of the cell cycle instead remains vague and underexplored. Canonical extrinsic apoptosis requires the activation of death receptors, such as TRAIL receptors 1 and 2, which is followed by a signaling cascade that leads to cleavage of executioner caspases (-3, -6 and -7) and subsequent cell death. Starting with binding of TRAIL and cluster formation of death receptors, cell-intrinsic signaling activates caspase-8. Apoptosis is either directly induced by caspase-8-dependent cleavage of pro-caspase-3 (type I) or relies on Bid cleavage and the permeabilization of the outer mitochondrial membrane (type II) [13]. Since de novo synthesis of proteins subsequent to TRAIL exposure is not required for apoptosis induction, independence between extrinsic apoptosis and cell cycle progression could be expected. On the other hand, expression, phosphorylation and localization of several proteins involved in signal transduction is controlled in a cell cycle-dependent manner [14–16]. To study if both dynamical processes, extrinsic apoptosis and cell cycle progression, are coupled, and due to substantial cell-to-cell heterogeneities even in isogenic cell populations [12, 17], the development and application of mathematical models and appropriate statistical tools is inevitable. Mathematical modeling of the cell cycle machinery has a long history (e.g. [18, 19]), including studies integrating time-lapse microscopy data of Fucci reporter cells (e.g. [20]), but modeling studies linking extrinsic apoptosis and cell cycle dynamics have not yet been conducted. Where initial work in this direction was attempted, modeling approaches linked proliferation and cell death to defined signaling networks [21, 22]). Although complex models are necessary for the understanding of signal transduction kinetics and the role of cellular heterogeneity and noise in cell populations [23], parametrization of high-resolution signaling models requires a significant amount of data and biological knowledge. A preceding step for describing and quantifying possible interconnections between cell cycle progression and extrinsic apoptosis signaling is defining parameters phenomenologically. Here, we therefore focused on statistical methods and phenomenological models to study the relationship of extrinsic apoptosis and cell cycle progression in NCI-H460/geminin cells [24] and HCT-116/geminin cells when these were exposed to a 2nd generation hexavalent TRAIL receptor agonist (IZI1551) [25].
Results
Cells in S/G2/M phase require longer to die than cells treated in G1 phase
To allow for an analysis of potential links between cell cycle phases and cell death timing (Fig 1), we first characterized cell cycle progression in NCI-H460 cells expressing mAG-hGeminin(1/110) as a fluorescent reporter of S/G2/M phases [24]. Durations for G1 (geminin negative) and S/G2/M (geminin positive) phases were recorded for approximately 400 cells and described by lognormal distributions (Fig 2A and 2B). Lognormal distributions outperformed gamma, weibull and normal distributions in describing these data, judged from comparison of Bayesian information criteria (BIC) [26, 27] (S1 Table). This criterion was chosen because it takes model fit and complexity into account. Previous studies showed that S/G2/M phases were relatively constant and mainly variability in G1 caused different cell cycle times [28]. Our data provide a different picture: magnitudes of mean and variance were comparable for both phases (Fig 2A and 2B) and we observed a strong linear correlation of both phases with cell cycle durations (Fig 2C and 2D). The Pearson correlation coefficient of 0.33 for the durations of G1 and S/G2/M (Fig 2E) indicated a low linear correlation of the two processes. These findings are in line with [29] and [20], except for the cell line specific phase proportions of the total cell cycle durations (45% and 39% for G1 in NCI-H460 and HCT-116 (see S4 Fig), respectively).
Arrows indicate information flow from experimental measurements to data normalization and analysis. The two key questions for exploring the apoptosis/cell cycle connection were: Does the cell cycle influence death timing and/or cellular fate? Is the cell cycle influenced by TRAIL? To answer these questions, the following measurements were required: A I) length of phases in untreated cells, A II) the history of cells that were stimulated with TRAIL (end of previous phase), A III) durations from TRAIL exposure to cell death or to the end of a cell cycle phase. On basis of this data, the influence of TRAIL on the phase length was analyzed (B). Information about the latter was incorporated to estimate the initial cell cycle position at TRAIL addition (C). Results were used to analyze the influence of the initial cell cycle position on cell fate and death timing (D).
(A,B) Length of G1 and S/G2/M phases in control cells. Parameters of the fitted lognormal distributions are μG1 = 1.96, σG1 = 0.27, μSG2M = 2.18, and σSG2M = 0.19 with mean values of 7.4 h and 9.0 h, respectively. Censored data (end of phase was not known or observable) and outliers are highlighted in red. (C-E) Correlations of phases in control cells. Similar to [20], the length of the cell cycle was plotted for individual cells against the respective length of the phases (C,D) and both phase lengths against each other (E). Pearson correlation coefficients (95% confidence interval) are as follows. (C): ρ = 0.83(0.82;0.84), (D): ρ = 0.8(0.76;0.83), (E) ρ = 0.33(0.24;0.42). (F,G) Tracking of cellular events in response to TRAIL. (F) Letters represent the following events in the exemplary video segments. a: G1 cell (geminin negative), b: S/G2/M cell (geminin positive), c: apoptosis, d: cell division. (G) Time of death (tdeath) box plots of cells treated with TRAIL in G1 (left panel) or S/G2/M (right panel). Boxes cover 25th to 75th percentiles. Median tdeath for cells exposed in G1 was 2.9 h, and for cells in S/G2/M 3.5 h, indicated by red bars. Red crosses indicate outliers.
Next, we studied cell fate and the time required for individual cells to execute cell death when treated with the EC50 dose of a hexavalent TRAIL receptor agonist [25] (Fig 1A). Approximately 550 G1 cells and 500 S/G2/M cells were monitored for at least 19 h prior to TRAIL addition to ensure that the mitotic history and thus age of each individual cell was known. The cell fate following TRAIL addition was recorded for another 19.25 h for randomly selected cells in the fields of view. Exemplary images are shown in Fig 2F. The probability of a cell to die when treated in G1 or S/G2/M was 0.76 and 0.72, respectively. At 95% confidence level the difference exceeded the maximal error calculated from the null hypothesis, indicating that the latter cells had a slightly increased likelihood to escape cell death induction. However, this difference was rather small and its functional significance is debatable. More interestingly, the time required until cells executed apoptosis (tdeath) differed substantially between cells exposed to TRAIL in G1 or S/G2/M phase, with S/G2/M cells requiring significantly longer to die (Δ median tdeath: 0.6 h, P-value: 4.0e-8, Fig 2G).
TRAIL prolongs the durations of G1 and S/G2/M cell cycle phases
Next, we studied if TRAIL exposure affects the duration of cell cycle phases, irrespective of whether cell death is induced or not (Fig 1B). Cell death events might distort measurements by which drug effects on the duration of phases are to be determined [30, 31]. To test if censoring due to cell death plays a role in accurately determining the durations of cell cycle phases in our scenario, we first described cell cycle progression by a mathematical model for a virtual cell population. As in [32], we expressed the cell cycle C on a scale from 0 to 1, so that the G1 phase is represented by 0–0.45 and S/G2/M phase by 0.45–1. Values >1 correspondingly represent cells after mitosis. To define the virtual cell population that served as reference, the initial positions of cells in the cycle were sampled from an idealized steady state distribution, adapted from [33] (1) The position Ci in the cell cycle progresses for each individual cell i, allowing for individual growth rates gi in G1 and S/G2/M phases so that (2) (3) (4) where represents the durations of the respective cell cycle phases. Next, we assigned a representative cell death distribution to this virtual cell population, based on our experimental data of death times. To this end, we described the times required to die after TRAIL exposure (tdeath, Fig 2G) by a lognormal distribution (5) taking into account that sibling cells dying after mitosis commit synchronous apoptosis shortly after mitosis [12, 17]. As such, cell death events after mitosis (f1 generation) were not inappropriately over-weighted. With help of Eqs 1–5, the cell cycle positions at the time of death were calculated with (6) We next compared the duration of G1 and S/G2/M phases between control cells and a cell population dying after TRAIL using this model. Cells that progressed through their cell cycle phase without dying first (C(tdeath) > f) do so faster than cells in the control population (approximately 0.5 h and 0.3 h for G1 and S/G2/M phases, respectively) (S1 Fig). This can be explained by removing cells that die before phase end in these types of analyses. Profiles of faster cell death distributions augment these effects further (S2 Fig). To accurately compare experimentally measured phase durations following TRAIL treatment to durations in untreated cells, we considered these time shifts. As seen in Fig 3A and 3B, TRAIL treatment prolonged both, G1 and S/G2/M phases, by 2.4 h and 1.3 h, respectively. For validation, we showed that the experimentally determined G1 durations in cells prior to TRAIL exposure was indistinguishable from the control distribution (S3 Fig). Next, we investigated if the prolongation of cell cycle phases is linked to the duration from phase start to TRAIL exposure (di) within G1 or S/G2/M phases (Fig 3C and 3D). Comparing the total phase duration against di revealed a negative slope during the first 4-5 h, indicating that the earlier a cell was stimulated in its current phase, the more pronounced was phase prolongation. As cells receiving TRAIL at later time points already exhibit longer phase durations, increasing values at later times were expected. Consequently, the most simplistic assumption is that TRAIL decelerates cell cycle progression in apoptosis competent cell populations. This can be expressed by introducing a prolongation term z for the respective cell cycle phases, resulting in a corrected growth rate (7) To determine z, we conducted a parameter estimation (for detailed description see method section) by comparing phase durations of experimental data to large sample populations that were constructed according to Eqs 1–7 (Fig 3E and 3F). Best values describing the experimental data were 3.6 h (3-4.3 h, 95% confidence interval) and 3.4 h (3-3.8 h, 95% confidence interval), respectively. Simulating G1 and S/G2/M durations based on these values provided distributions very similar to the experimental results (Fig 3G and 3H, compare to Fig 3C and 3D). This confirms that the mathematical model describes the influence of TRAIL on cell cycle progression sufficiently well. For the cell line HCT-116, we observed a similar behaviour, with slightly lower values (2.5 h (2.1-3 h), and 1.95 h (1.4-2.5 h) for G1 and S/G2/M, respectively) (S5 Fig).
(A,B) Lengths of G1 and S/G2/M phases of NCI-H460 cells that survived the phase of TRAIL exposure. The expected distributions were approximated under consideration of the fraction of surviving and uncensored, apoptotic cells. (C,D) Experimentally determined correlations of phase length and time of TRAIL addition. A linear gaussian process was assumed to highlight the 95% confidence interval. One outlier (S/G2/M > 25 h) was omitted from the display (D). (E,F) Negative log-likelihood for varying prolongations of phases (z-values). 6 experiments with 3e7 samples were conducted. Mean and standard deviations are shown and a polynomial of 5th grade was fitted to the resulting mean. The lowest negative log-likelihood values were obtained for (3 h-4.3 h), and (3 h-3.8 h), respectively. Confidence bounds are indicated with dashed lines. (G,H) Exemplary model simulations with -values are shown. The 95% confidence intervals were calculated with a linear gaussian process and highlighted in gray and green.
The time of TRAIL exposure in S/G2/M phase appears linked to the likelihood of apoptotic cell death
Next, we studied if the position of a cell within the G1 or S/G2/M phases at the time of exposure to TRAIL affects the probability and timing of cell death (Fig 1D I). To do so, we developed an approach to estimate the position of individual cells within the cell cycle reliably (Fig 1C). Therefore, we took into account that the durations of the cell cycle phases are heterogeneous across the population, that apoptosis events impede determining the end point of a phase in a fraction of the cell population, and that TRAIL exposure itself alters the progression of the cell cycle (see Fig 3). The concept of our approach was initially tested on a virtual cell population consisting of i cells with proliferation rates according to Eq 7 and a correction term z*. The initial cell cycle positions were sampled from a standard uniform distribution and death time () was sampled according to Eq 5. We extracted the following information from the virtual cell population: i) the duration from start of the cell cycle phase to TRAIL addition (di), ii) the cell cycle position at time of cell death C(tdeath), and iii) the duration of the cell cycle phase (pi) for cells that failed to die within this phase. We then tested if the initial cell cycle positions of cells in this virtual population can be estimated by (8) (9) where f represents the respective value of G1 and S/G2/M phases (compare to Eqs 2 and 3 and associated Fig 2A and 2B) and pdf(pu) represents the distribution of phase lengths in unstimulated cells (see Eq 2). For cells that died after the end of the cell cycle phase in which they were treated, a sample of a truncated distribution was used for normalization of the cell cycle position. This allowed taking into account that such cells already spent a defined time within the respective cell cycle phase, prior to TRAIL treatment and induction of cell death. In cells that continued to proliferate, the cell cycle position was calculated as a function of , the estimated TRAIL-induced cell cycle prolongation, the duration of the respective phase pi, and the time of TRAIL addition di (Eq 8). Our estimations of the cell cycle positions agreed very well with the real position of the virtual population (Fig 4A(I)). In contrast, estimations based solely on normalization with samples of the truncated population (Fig 4A(II)) or the unstimulated population (Fig 4A(III)) were more error prone. Reliable estimations of cell cycle positions were also possible for alternative correction terms z* for the virtual cell population (Fig 4B). Our approach is therefore suitable for arranging cells according to their cell cycle position. Consequently, we next applied this to experimental data to asses if initial cell cycle positions influence cell fate. When NCI-H460 cells were stimulated with TRAIL in G1, cell fate was independent of the initial cell cycle position (Fig 4C). With respect to S/G2/M, cells exposed to TRAIL earlier were more likely to survive (P-value: 0.043, Fig 4C).
(A) C0 indicates the real initial cell cycle position whereas represents estimated values. For calculating , the approximation of the initial cell cycle position, a data set of a virtual cell population was normalized with a sample of (III) the original distribution of G1 or S/G2/M phases or (II) the respective truncated distributions. In (I), either the truncated distribution was used for normalization or, for cells that completed the phase, was calculated as a function of and the length of the phase. (B) Indicated values of z* were assumed for constructing the virtual cell population. Sum of least squares are indicated for the different scenarios. (C) Impact of initial cell cycle positions on cell fate in experimental data. Gray bars represent median values. Apoptotic: cell death before division or both daughter cells died. Survivor: both daughter cells survived. Mixed: 2 in G1 and 29 in S/G2/M. P-values according to Wilcoxon rank sum test [35] are as follows. b/c: 0.95 (not significant, ns), e/f: 0.043 (significant, *).
Model analysis describes apoptosis deceleration in S phase
Next, we analyzed if the cell cycle position at the time of TRAIL addition influences the time between TRAIL exposure and cell death (Fig 1D II). Interestingly, when plotting the estimated cell cycle position at TRAIL exposure against the time to death (tdeath) for each individual cell (Fig 5A), the population divides. Cells stimulated in late G1 or early S phase, either died in f0 with death times comparable to stimulation in early G1, or alternatively with a significantly prolonged death time after the subsequent mitosis. Prolonged death times disappeared in cells exposed to TRAIL close to mitosis (Fig 5A). To identify the cell cycle position at which apoptosis suppression manifests, we developed a phenomenological model. The model consists of Eqs 2–4 and an additional, coupled ordinary differential equation describing apoptosis progression (10) (11) (12) (13) (14) Apoptosis progresses over time with a constant slope m and p represents inhibition of apoptosis progression. PAD (Point of Apoptosis Deceleration) refers to the unknown point in the cell cycle after which apoptosis progression will be delayed. and emphasize model and experimental death times, respectively. Since molecularly the delay in apoptosis execution might be related to cell cycle-regulated changes in protein amounts, the delay or inhibition of signal propagation is removed after cell division. First, we approximated tdeath of early G1 cells with a lognormal distribution (μtd = 1, σtd = 0.39), resulting in an average death time of 2.9 h in these cells. mi (Eq 10) is the inverse of a sample of this distribution ( 1/h). Next, the negative log-likelihood of death times in dependence on the initial cell cycle position C0 (Eq 13) for different values of p and PAD was used to calculate corresponding significance levels and approximate best parameters and (Eq 14) (for detailed description see method section). We conducted these computations for six independent sampling experiments, all of which provided similar best points and confidence regions of and (Fig 5B). Based on 95% confident intervals of PAD and p, the time point of apoptosis deceleration was reached at around 49% to 55% of the cell cycle. We therefore conclude that apoptosis progression was delayed approximately in mid S phase. Importantly, apoptosis progression was not fully blocked, since slope (Fig 5B). An example simulation based on the estimated parameters indeed reproduced the experimental data very well (Fig 5C). In Fig 5D, the idealized courses of apoptosis progression in relation to cell age are visualized, where young cells die faster than older cells upon TRAIL exposure. When comparing these results to different “null” models in which cell cycle-independent distributions were fitted to the death kinetics, substantially larger values for the Bayesian information criterion were obtained (S2 Table). These data further support cell cycle dependence of extrinsic apoptosis progression. Although cell death was slower in the case of the cell line HCT-116, parameter estimation suggests a similar point of apoptosis deceleration at the end of G1 or beginning of S-phase (S6 Fig).
(A) Death times in experimental data (n = 921) were plotted against estimated initial cell cycle positions. Symbols indicate cell cycle phases at time of cell death. Green shading indicates S/G2/M phase and dashed line represents division of an averaged, untreated NCI-H460/geminin cell. (B) Parameter estimation was conducted six times with independent samples of initial model conditions. The sample size for each simulation experiment was approximately 8e4. Significance levels are shown by color coding. Indicated values of and are mean and extreme values of 6 experiments. (C) An exemplary model-based simulation with best parameter values and is shown. (D) Average apoptotic progression in dependence on the ages of cells, representing the time from birth to TRAIL addition, is illustrated. The area of decelerated apoptosis progression is highlighted in gray. (E) Different synchronization strategies (C0 ∈ (0, 0.1), C0 ∈ (0.1, 0.2),…) impact death timing and number of divisions.
Since cells that proceed through mitosis as a consequence of apoptosis suppression ultimately might escape cell death, we conceptually studied if and which cell cycle synchronization strategies would maximize rapid apoptosis execution in cell populations. In a virtual cell population constructed from the presented experimental results, we analyzed how synchronization within G1 or S/G2/M affects tdeath and the number of cells that go through mitosis (Fig 5E). Cells died most rapidly when treated in the first 30% of the cell cycle. In contrast, TRAIL exposure between 50% to 70% of the cell cycle resulted in the longest and most variable death times (Fig 5E). Correspondingly, pre-treatment with the FDA-approved CDK4/6 inhibitor Abemaciclib synchronized NCI-H460/geminin cells in G1 phase and significantly reduced the death time in response to TRAIL (S7 Fig). Overall, these data therefore demonstrate that apoptosis progression is suppressed in favor of mitosis by mid S phase and that optimal treatment responsiveness can be expected for cells exposed to TRAIL early after division.
Discussion
In this study we combined microscopic fluorescence-based monitoring of cell cycle progression and cell death execution with sampling-based mathematical modeling to study if extrinsic apoptosis signaling in NCI-H460/geminin cells is influenced by cell cycle progression. Interestingly, we found that cell cycle progression, in particular transitioning through S phase, significantly prolongs the time required for cells to die in response to TRAIL. In a companion study, we found that these delays suffice to let sublethally damaged cells escape from apoptosis [36]. In reverse, TRAIL addition prolonged G1 and S/G2/M phases, a previously unreported phenomenon. Our study is the first to systematically and quantitatively analyze the interdependence of cell cycle progression and TRAIL-induced apoptosis without exteriorly interfering with either of these processes during the course of the analyses. This could be achieved by integrating non-invasive time-lapse imaging with mathematical modeling and parameter estimation. As experimental data sources, our study strongly relied on determining accurate kinetics for cell cycle progression and the times required for TRAIL to induce apoptosis execution. Prolonged time-lapse imaging, coupled with information on the mitotic history of the analyzed cells and fluorescence-based cell cycle phase monitoring, allowed us to confidently estimate cell cycle positions at minutes resolution. In addition, population heterogeneity in cell cycle synchrony could be taken into account within the modeling and parameter estimation framework. Our approach therefore circumvented the limitations imposed by conventional cell synchronization in experimental studies, such as arresting and releasing cells in distinct cell cycle phases, which would have prevented analyzing signaling interplay during regular cell cycle progression [37].
Indeed, conditions where cells are arrested in or released from distinct cell cycle phases can easily lead to the misinterpretation of results on TRAIL sensitivity, since cell cycle arrest in itself sensitizes to apoptosis [38]. Consequences on apoptosis susceptibility in such settings might therefore be independent of cell cycle progression. For example, tumor suppressor p53 plays a crucial role in cell cycle arrest but also induces the expression of various apoptosis proteins [6, 39], which in turn can synergize with TRAIL-induced caspase 8-dependent cell death. Despite the interesting insight such studies provide for TRAIL response potentiation, a comprehensive understanding of the interrelation between cell cycle progression and susceptibility to TRAIL-induced apoptosis can only be obtained by the approach chosen here, using unsynchronized cell populations. However, also studies in which cells were not artificially synchronized reported inconsistent findings for the relationship between extrinsic apoptosis susceptibility and cell cycle progression [12, 38, 40–42]. One problem in “observation only” studies arises from neglecting the possibility of bidirectional relationships between cell cycle progression on TRAIL responsiveness and TRAIL exposure on cell cycle progression, both of which we identified in our study. This aspect is non-trivial, since cells dying during the course of treatment pose a significant problem for the cell cycle analysis. We solved this problem by integrating mathematical modeling and parameter estimation into the workflows. To our knowledge, phenomenological mathematical models that define connections between the cell cycle and extrinsic apoptosis progression have not been reported yet. However, model-based studies that combine proliferation and cell death exist [43, 44]. In larger scale mechanistic models, extrinsically stimulated cellular signaling processes were implemented together with abstracted cell cycle signaling [21, 22], however, their complexity would impose parameter estimation to become computationally highly expensive and problems in parameter identifiability would arise. Focusing solely on a phenomenological description of the interconnection between cell cycle and cell death instead allowed us to appropriately estimate parameters and obtain novel biological insight. This might serve as a basis on which more complex mechanistic models could be implemented to dissect the molecular mechanisms that underlie the observed effects. The approach chosen here can likewise serve as a template to study other cell biological processes and their potentially bidirectional relationship to cell cycle progression.
Materials and methods
Time-lapse microscopy and cell tracking
Establishment of the NCI-H460/geminin cell line was described in [24]. To observe cell cycle history and cell death kinetics, cells were grown on 35 mm glass-bottom dishes (CellView Cell Culture Dish, Greiner Bio One) in phenol red free RPMI 1640 containing 10% FBS. Time-lapse images were acquired every 15 minutes on a Zeiss Cell Observer microscope equipped with a Plan-Apochromat 20x/0.8 objective at 37°C, 5% CO2. Fluorescence of mAG-hGeminin(1/110) was acquired with a 470 nm LED module combined with a 62 HE filter set (Carl Zeiss). After 19 h, Fc-scTRAIL (0.06 nM, [25]) was added and cells were imaged for another 19.25 h. Images were processed using the ZEN2.1 blue software (Carl Zeiss). Randomly chosen cells were tracked manually and the time of birth, onset of geminin expression (representing entry into S phase), cell division and the time from TRAIL addition to death (tdeath) was recorded. A cell was assigned as geminin positive as soon as geminin expression was above the background. Geminin positive cells were assigned to the combined S/G2/M phases. A cell was assigned to be apoptotic as soon as cell blebbing was visible.
With respect to the experiments shown in S7 Fig, NCI-H460/geminin cells were pre-treated with DMSO or Abemaciclib (2 μM, Selleckchem) for 3 h. Then, medium containing Fc-scTRAIL was added and cells were imaged for another 10 h. Randomly chosen cells were tracked and tdeath of cells was recorded.
Modeling and statistical analysis
For calculations, MATLAB (MATLAB R2017a, The MathWorks, Natick, 2017) and the Statistics and Machine Learning Toolbox were used. For computing the Bayesian information criteria (BIC) [26, 27], the length of phases that were not completely measurable was included as censored information. In order to approximate expected distributions for apoptotic cells (S1 Fig), independent samples of C0 (Eq 1) and experimentally determined tdeath (Eq 5) distributions were drawn. was calculated (Eq 6) and cells with were discarded. A lognormal distribution was fitted to the remaining phase lengths.
For approximating expected distributions (Fig 3A and 3B), the experimentally determined ratio of surviving and apoptotic cells was taken into account. The probability that two distributions of interest are samples from underlying distributions with the same median was estimated with help of the Wilcoxon rank sum test [35]. Except for results shown in S6 Fig, ordinary differential equations were solved analytically. For S6 Fig, the Matlab toolbox IQM Tools [46] was used. Computational methods and unprocessed data are available at Github: https://github.com/DirkeImig/CellCycleApoptosis.
Parameter estimation.
For parameter estimation of correlations between length of phases and time of TRAIL addition di (Fig 3E and 3F), the following algorithm was used: First, representative samples of C0, tdeath and g were drawn to construct a virtual cell population. Surviving cells () were added according to the fraction of surviving cells that was observed experimentally. Cells that died before reaching the end of the respective phase (, Eq 6) were discarded. Next, di was calculated (Eqs 2–4) for the remaining cells. Then, the overall lengths of the phases were computed (Eqs 2–4 with gi = gi,*, Eq 7). A lognormal distribution was fitted to the vector of model phase lengths that correspond to a specific time of TRAIL addition and the respective likelihood of describing the experimental data was evaluated. Lastly, the sum of all negative log-likelihood values was calculated. Negative log-likelihood values for different z−values of 6 separate simulation experiments with 5e7 sampled cells were averaged and shown in Fig 3E and 3F. The minimum value of the negative log-likelihood revealed the best parameter estimate .
For parameter estimation of death kinetics (Fig 5B), first, representative samples of C0, g and m were drawn and was calculated for given values of p and PAD (Eqs 10–14). and represent model and experimental data, respectively. We approximated the bivariate model density of C0 and with help of a kernel density estimator. For each data tuple {,}, the likelihood of being described by the model was evaluated and the overall negative log-likelihood was calculated. The process was repeated for different values of PAD, and p, originated from a predefined grid in the parameter space. To avoid rounding errors, the lowest value within the density was set to 1e-9. Confidence intervals were calculated according to χ2 statistics [34].
Supporting information
S1 Table. Information criteria for different distributions.
The negative log-likelihood and respective BIC values are shown for the indicated model distributions. The number of parameters equals 2 for lognormal, weibull and normal distributions and 4 in case of the gamma distribution. A lower BIC indicates a better model fit, with a difference higher than 2 being positive in evidence [45].
https://doi.org/10.1371/journal.pcbi.1007812.s001
(TIF)
S2 Table. Δ BIC for different null models.
Differences of BIC values between null models with indicated distributions and best BIC for the model in Eqs 10–14 is shown. A Δ in BIC > 2 means positive in evidence [45]. Here, all differences are significant.
https://doi.org/10.1371/journal.pcbi.1007812.s002
(TIF)
S1 Fig. Expected distributions of G1 and S/G2/M in apoptotic cells.
A sample of the fitted lognormal distribution of untreated cells is shown in black (G1, A) and green (S/G2/M, B). In white, lengths in the modeled virtual population are shown, given C(tdeath) > f (Eq 3).
https://doi.org/10.1371/journal.pcbi.1007812.s003
(TIF)
S2 Fig. Influence of different tdeath distributions on expected distributions of uncensored cells.
Different values for μ and σ regarding the lognormal tdeath distributions are evaluated. The mean differences between original distributions and calculated values given that C(tdeath) > f (Eq 3) are shown. The lower μ, the greater the difference between original and approximated distribution. Reason for this is that with a faster death impulse, the likelihood of reaching the end of a phase is lower for cells with relatively long phase lenghts. The underlying phase lengths distribution corresponds to G1 phase lengths in control cells. This is similarly valid for S/G2/M phases.
https://doi.org/10.1371/journal.pcbi.1007812.s004
(TIF)
S3 Fig. Length of G1, f0 phase for cells treated in S/G2/M.
In order to exclude a possible influence of experimental conditions on phase length variability, the G1, f0 phase lengths of cells that were treated with TRAIL in S/G2/M (n = 446) was compared to the control case. With a P-value of 0.54 [35], the null hypothesis of the two distributions originating from a distribution with the same median could not be rejected.
https://doi.org/10.1371/journal.pcbi.1007812.s005
(TIF)
S4 Fig. Length of G1 and S/G2/M in control HCT-116 cells.
(A,B) Parameters of underlying lognormal distributions are μG1 = 1.87, σG1 = 0.34 and μSG2M = 2.34, σSG2M = 0.22 with mean values of 6.9 h and 10.7 h. Censored data (end of phase is not known) and outliers are highlighted in red.
https://doi.org/10.1371/journal.pcbi.1007812.s006
(TIF)
S5 Fig. Lengths of G1 and S/G2/M phases in HCT-116 in response to TRAIL.
(A,B) Experimentally determined correlations of phase length and time of TRAIL addition. A linear gaussian process was assumed to highlight the 95% confidence interval. (C,D) Negative log-likelihood for varying prolongations of phases (z-values). 3 experiments with 3e6 samples were conducted. Mean and standard deviations are shown and a polynomial of 5th grade was fitted to resulting mean values. The lowest negative log-likelihood values were obtained for (2.1 h–3 h), and (1.4 h–2.5 h), respectively. Confidence bounds are indicated with dashed lines. (E,F) Exemplary model simulations with -values are presented. The 95% confidence intervals were calculated with a linear gaussian process, highlighted in gray and green.
https://doi.org/10.1371/journal.pcbi.1007812.s007
(TIF)
S6 Fig. Timing of death depends on cell cycle positions in HCT-116.
(A) Death times of experimental data (n = 538) were plotted against estimated initial cell cycle positions. Symbols indicate cell cycle phases at time of cell death. Experiment time was 24.5 h. Green shading highlights S/G2/M phase and the dashed line represents division of an averaged, untreated HCT-116/geminin cell. (B) Parameter estimation with six independent samples of initial conditions. Each sample size was approximately 1e4. m was estimated from cells dying in early G1. For simulation, the Matlab toolbox IQM Tools [46] was used and PAD was extended for the second generation so that p = 0 if Ci ∈ [0, PAD] or [1, PAD + 1] (see Eq 10). Significance levels are shown by color coding. Values of and represent mean and extreme values of six simulation experiments. (C) An exemplary simulation with best parameter values and is shown. (D) Average apoptotic progression in dependence on the ages of cells, representing the time from birth to TRAIL addition, is illustrated. The area of decelerated apoptosis progression is highlighted in gray.
https://doi.org/10.1371/journal.pcbi.1007812.s008
(TIF)
S7 Fig. Cell death after inhibition of CDK4/6 in NCI-H460/geminin cells.
A. Representative time-lapse images of NCI-H460/geminin cells treated with Fc-scTRAIL (0.06 nM) or Abemaciclib (2 μM). Scale bar represents 50 μm. Dying cells are highlighted by white arrowheads. B. Times required to die (tdeath) of NCI-H460/geminin cells following treatment with Fc-scTRAIL alone or after pre-treatment with Abemaciclib (2 μM, 3 h). Shown are medians with interquartile ranges, plus min to max range from n = 50 cells observed in three independent experiments (**** p < 0.0001, unpaired t-test). Na = not applicable.
https://doi.org/10.1371/journal.pcbi.1007812.s009
(TIF)
References
- 1. Spellman PT, Sherlock G, Zhang MQ, Iyer VR, Anders K, Eisen MB, et al. Comprehensive identification of cell cycle-regulated genes of the yeast Saccharomyces cerevisiae by microarray hybridization. Mol Biol Cell. 1998;9(12): 3273–3297. pmid:9843569
- 2. Cho RJ, Campbell MJ, Winzeler EA, Steinmetz L, Conway A, Wodicka L, et al. A genome-wide transcriptional analysis of the mitotic cell cycle. Mol Cell. 1998;2(1): 65–73. pmid:9702192
- 3. Salazar-Roa M, Malumbres M. Fueling the Cell Division Cycle. Trends Cell Biol. 2017;27(1): 69–81. pmid:27746095
- 4. Sakaue-Sawano A, Kurokawa H, Morimura T, Hanyu A, Hama H, Osawa H, et al. Visualizing spatiotemporal dynamics of multicellular cell-cycle progression. Cell. 2008;132(3): 487–498. pmid:18267078
- 5. Sakaue-Sawano A, Kobayashi T, Ohtawa K, Miyawaki A. Drug-induced cell cycle modulation leading to cell-cycle arrest, nuclear mis-segregation, or endoreplication. BMC Cell Biol. 2011 Jan 13;12(2). pmid:21226962
- 6. Pucci B, Kasten M, Giordano A. Cell Cycle and Apoptosis. Neoplasia. 2000;2(4): 291–299. pmid:11005563
- 7. Miwa S, Yano S, Kimura H, Yamamoto M, Toneri M, Matsumoto Y, et al. Cell-cycle fate-monitoring distinguishes individual chemosensitive and chemoresistant cancer cells in drug-treated heterogeneous populations demonstrated by real-time FUCCI imaging. Cell Cycle. 2015;14(4): 621–629. pmid:25551170
- 8. Marcus JM, Burke RT, DeSisto JA, Landesman Y, Orth JD. Longitudinal tracking of single live cancer cells to understand cell cycle effects of the nuclear export inhibitor, selinexor. Sci Rep. 2015 Sep 24;5: 14391. pmid:26399741
- 9. Tsuchida E, Kaida A, Pratama E, Ikeda MA, Suzuki K, Harada K, et al. Effect of X-Irradiation at Different Stages in the Cell Cycle on Individual Cell-Based Kinetics in an Asynchronous Cell Population. PLoS One. 2015, Jun 18;10(6): e0128090. pmid:26086724
- 10. Maeda S, Wada H, Naito Y, Nagano H, Simmons S, Kagawa Y, et al. Interferon-α acts on the S/G2/M phases to induce apoptosis in the G1 phase of an IFNAR2-expressing hepatocellular carcinoma cell line. J Biol Chem. 2014;289(34): 23786–23795.
- 11. Jin Z, Dicker DT, El-Deiry WS. Enhanced sensitivity of G1 arrested human cancer cells suggests a novel therapeutic strategy using a combination of simvastatin and TRAIL. Cell Cycle. 2002 Jan;1(1): 82–89. pmid:12429913
- 12. Spencer SL, Gaudet S, Albeck JG, Burke JM, Sorger PK. Non-genetic origins of cell-to-cell variability in TRAIL-induced apoptosis. Nature. 2009;459(7245): 428–432. pmid:19363473
- 13. Johnstone RW, Frew AJ, Smyth MJ. The TRAIL apoptotic pathway in cancer onset, progression and therapy. Nat Rev Cancer. 2008;8(10): 782–798. pmid:18813321
- 14. Li F, Ambrosini G, Chu EY, Plescia J, Tognin S, Marchisio PC, et al. Control of apoptosis and mitotic spindle checkpoint by survivin. Nature. 1998;396(6711): 580–584. pmid:9859993
- 15. Scaffidi C, Volkland J, Blomberg I, Hoffmann I, Krammer PH, Peter ME. Phosphorylation of FADD/ MORT1 at serine 194 and association with a 70-kDa cell cycle-regulated protein kinase. J Immunol. 2000;164(3): 1236–1242. pmid:10640736
- 16. Harley ME, Allan LA, Sanderson HS, Clarke PR. Phosphorylation of Mcl-1 by CDK1-cyclin B1 initiates its Cdc20-dependent destruction during mitotic arrest. EMBO J. 2010;29(14): 2407–2420. pmid:20526282
- 17. Rehm M, Huber HJ, Hellwig CT, Anguissola S, Dussmann H, Prehn JHM. Dynamics of outer mitochondrial membrane permeabilization during apoptosis. Cell Death Differ. 2009;16(4): 613–623. pmid:19136937
- 18. Gérard C, Goldbeter A. A skeleton model for the network of cyclin-dependent kinases driving the mammalian cell cycle. Interface Focus. 2011;1(1): 24–35. pmid:22419972
- 19. Tyson JJ, Novak B. Regulation of the eukaryotic cell cycle: molecular antagonism, hysteresis, and irreversible transitions. J. Theor. Biol. 2001;210(2): 249–263. pmid:11371178
- 20. Dowling MR, Kan A, Heinzel S, Zhou JH, Marchingo JM, Wellard CJ, et al. Stretched cell cycle model for proliferating lymphocytes. Proc Natl Acad Sci U S A. 2014;111(17): 6377–6382. pmid:24733943
- 21. Imig D, Pollak N, Strecker T, Scheurich P, Allgöwer F, Waldherr S. An individual-based simulation framework for dynamic, heterogeneous cell populations during extrinsic stimulations. J. Coupled Syst. Multiscale Dyn. 2015;3(2): 143–155.
- 22. Bertaux F, Stoma S, Drasdo D, Batt G. Modeling dynamics of cell-to-cell variability in TRAIL-induced apoptosis explains fractional killing and predicts reversible resistance. PLoS Comput Biol. 2014 Oct 23;10(10): e1003893. pmid:25340343
- 23. Suderman R, Bachman JA, Smith A, Sorger PK, Deeds EJ. Fundamental trade-offs between information flow in single cells and cellular populations. PNAS. 2017;114(22): 5755–5760. pmid:28500273
- 24. Kuritz K, Stöhr D, Pollak N, Allgöwer F. On the relationship between cell cycle analysis with ergodic principles and age-structured cell population models. J. Theor. Biol. 2017;414: 91–102. pmid:27908704
- 25. Hutt M, Marquardt L, Seifert O, Siegemund M, Müller I, Kulms D, et al. Superior properties of Fc-comprising scTRAIL fusion proteins. Mol Cancer Ther. 2017;16(12): 2792–2802. pmid:28904131
- 26. Schwarz G. Estimating the dimension of a model. The Annals of Statistic. 1978;6(2): 461–464.
- 27. Wit E, van den Heuvel E, Romeijn JW. ‘All models are wrong…’: an introduction to model uncertainty. Statistica Neerlandica. 2012;66(3): 217–236.
- 28. Smith JA, Martin L. Do cells cycle? Proc Natl Acad Sci U S A. 1973;70(4): 1263–1267. pmid:4515625
- 29. Chao HX, Fakhreddin RI, Shimerov HK, Kedziora KM, Kumar RJ, Perez J et al. Evidence that the human cell cycle is a series of uncoupled, memoryless phases. Mol Syst Biol. 2019;15(3): e8604. pmid:30886052
- 30. Hawkins ED, Turner ML, Dowling MR, van Gend C, Hodgkin PD. A model of immune regulation as a consequence of randomized lymphocyte division and death times. Proc Natl Acad Sci U S A. 2007;104(12): 5032–5037. pmid:17360353
- 31. Hawkins ED, Markham JF, McGuinness LP, Hodgkin PD. A single-cell pedigree analysis of alternative stochastic lymphocyte fates. Proc Natl Acad Sci U S A. 2009;106(32): 13457–13462. pmid:19633185
- 32. Imig D, Kuritz K, Pollak N, Rehm M, Allgöwer F. Death patterns resulting from cell cycle-independent cell death. IFAC-PapersOnLine. 2018;51(19): 90–93.
- 33. Powell EO. Growth rate and regeneration time of bacteria, with special reference to continuous culture. J. Gen. Microbiol. 1956;15(3): 492–511.
- 34. Wilks SS. The Large-Sample Distribution of the Likelihood Ratio for Testing Composite Hypotheses. Ann. Math. Statist. 1938;9(1): 60–62.
- 35. Wilcoxon F. Individual comparisons of grouped data by ranking methods. Journal of Economic Entomology. 1946;39(2): 269–270. pmid:20983181
- 36. Pollak N, Lindner A, Imig D, Kuritz K, Heinrich I, Stadager J, et al. Mcl-1 governs transmitotic resistance to extrinsic apoptosis. Submitted.
- 37. Cooper S. Rethinking synchronization of mammalian cells for cell cycle analysis. Cell Mol Life Sci. 2003;60(6): 1099–1106. pmid:12861378
- 38. Ehrhardt H, Wachter F, Grunert M, Jeremias I. Cell cycle-arrested tumor cells exhibit increased sensitivity towards TRAIL-induced apoptosis. Cell Death Dis. 2013 Jun 6;4:e661. pmid:23744361
- 39. Kastenhuber ER, Lowe SW. Putting p53 in Context. Cell. 2017;170(6): 1062–1078.
- 40. Márquez-Jurado S, Díaz-Colunga J, das Neves RP, Martinez-Lorente A, Almazán F, Guantes R, et al. Mitochondrial levels determine variability in cell death by modulating apoptotic gene expression. Nat Commun. 2018 Jan 26;9(1): 389, pmid:29374163
- 41. Hashimoto T, Juso K, Nakano M, Nagano T, Kambayashi S, Nakashima A, et al. Preferential Fas-mediated apoptotic execution at G1 phase: the resistance of mitotic cells to the cell death. Cell Death Dis. 2012 May 24;3:e313. pmid:22622132
- 42. Li Y, Dida F, Iwao A, Deguchi T, Azuma E, Komada Y. Cell cycle dependency of caspase activation in Fas-induced apoptosis in leukemia cells. Cancer Sci. 2007;98(8): 1174–1183. pmid:17561974
- 43. Tyson DR, Garbett SP, Frick PL, Quaranta V. Fractional Proliferation: A method to deconvolve cell population dynamics from single-cell data. Nat Methods. 2012;9(9), 923–928.
- 44. Sizek H, Hamel A, Deritei D, Campbell S, Ravasz Regan E. Boolean model of growth signaling, cell cycle and apoptosis predicts the molecular mechanism of aberrant cell cycle progression driven by hyperactive PI3K. PLoS Comput Biol. 2019 Mar 15;15(3):e1006402. pmid:30875364
- 45. Kass RE, Raftery AE. Bayes Factors. Journal of the American Statistical Association. 1995;90(430): 773–795.
- 46. Schmidt H, Jirstrand M. Systems Biology Toolbox for MATLAB: a computational platform for research in systems biology. Bioinformatics. 2006;22(4): 514–515. pmid:16317076