Long-term p21 and p53 trends regulate the frequency of mitosis events and cell cycle arrest

Radiation exposure of healthy cells can halt cell cycle temporarily or permanently. In this work, two single cell datasets that monitored the time evolution of p21 and p53, one subjected to gamma irradiation and the other to x-ray irradiation, are analyzed to uncover the dynamics of this process. New insights into the biological mechanisms were found by decomposing the p53 and p21 signals into transient and oscillatory components. Through the use of dynamic time warping on the oscillatory components of the two signals, we found that p21 signaling lags behind its lead signal, p53, by about 3.5 hours with oscillation periods of around 6 hours. Additionally, through various quantification methods, we showed how p21 levels, and to a lesser extent p53 levels, dictate whether the cells are arrested in their cell cycle and how fast these cells divide depending on their long-term trend in these signals.


Introduction
The tumor suppressor p53 is often referred to as the "guardian of the genome" for its importance in regulating a large number of target genes. It has an essential role in regulating the cell cycle, inducing cell apoptosis or DNA repair, and senescence [18]. Because of the centrality of p53 in maintaining the integrity of the DNA and preventing mutated cells from proliferating, mutation to the TP53 gene is found in over 50% of cancers [28]. Inactivation of p53 can take various forms: viral infection, deletion of the p14 ARF gene, deletion of the carboxy-terminal domain, multiplication of the murine double minute 2 gene (Mdm2), and amino-acid-changing mutation in the DNA-binding domain [34].
The interaction between p53 and its negative regulator Mdm2 is of great importance in understanding its dynamic behavior. The activation of p53 as a result of cell stress or damage induces an increase in Mdm2 levels that in terms inactivates p53 through proteasomal-dependent degradation [16]. Because this feedback loop has a time delay of 30-40 min between the increase in p53 and the increase in Mdm2 transcription, this creates the oscillatory behavior that is often observed in the p53 dynamics [17]. This p53-Mdm2 interaction and the oscillatory dynamics of p53 have been the subject of extensive research [1,35,38,14,3,27]. Investigating the p53 pathway and how it relates to cell fate, such as senescence and apoptosis, is essential to developing cancer treatments [31]. In addition to Mdm2, the interactions of p53 with other elements of the pathway such as p14 ARF , Chk2, Wip1, and ATM are considered of importance [3,9,13].
Cellular senescence as well as temporary instances of cell cycle arrest is also of great interest as a potential therapeutic target through the prevention of cancerous cells from proliferating [37,11,20]. Notably, for cells to remain arrested, this signaling pathway needs to be maintained [4,12,23]. Among a large number of targets that p53 influences, the cyclin-dependent kinase (CDK) inhibitor p21 plays a crucial role in the ability of cells to undergo cell cycle arrest [15,36,7]. p21 interacts with CDK2 as a bistable switch with increased levels of p21 leading to cell cycle arrest, while low levels of p21 allow cells to escape G1 phase arrest [6,21,2,26,23].
Through the applications of various mathematical techniques to two single cell datasets that monitored p21 and p53 activity after exposure to various levels of gamma or x-ray radiation, the influence of p21 and p53 signaling on cell mitosis is analyzed.

Clustering of periodic signals using t-SNE and the Wasserstein distance
Various types of distance metrics are used when clustering data using t-Distributed Stochastic Neighbor Embedding (t-SNE) [30]. Common implementations include the application of the Euclidean distance, Chebychev distance, correlation, cosine, and Spearman's rank correlation. By computing pairwise distances between high-dimensional points in the data, the algorithm utilizes this information to compute 2D or 3D representations of the data in a way that aims to cluster the data in a meaningful, interpretable way. Thus, choosing the appropriate distance to employ in the t-SNE framework is of great importance for the quality of the clustering that one derives for the separation of the various subsets of a given dataset.
For time series data, each dimension represents the value of the measured variable at a given time t. In the case of oscillating signals such as p21 and p53, employing t-SNE with a distance that compares time series by considering differences across measurements taken at the same time t may not yield meaningful clusters. This is because periodic signals may not be synchronized in their phases. Thus, even if one was to compare almost identical cells, slight phase offsets may be enough to make this comparison challenging. This is not as much of a problem with p21 as it is for p53, since the trend in the data is much larger in magnitude than the oscillatory part.
To overcome the limitation in dealing with oscillatory signals, we use the Wasserstein distance (WD) from optimal mass transport (OMT). Through the WD, one considers time series as distributions, removing the dependence on time t. Thanks to its properties and natural ability to capture geometric information when comparing different types of data and also the development in computational tools, the WD has attracted widespread interest in many application fields during the last several years. We will just sketch the relevant definitions here. For full details and numerous references, we refer the interested reader to [22,32,33].
The Wasserstein-1 distance (Kantorovich formulation) between two d-dimensional probability distributions ρ 0 , ρ 1 defined on a subdomain Ω ⊂ R d is given as follows: where · , is the ground metric. Thus the infimum is taken over all joint distributions π : Ω × Ω → R whose marginals are ρ 0 , ρ 1 .
For the time series analysis of p21 and p53 in the present work, one has d = 1, and in this case, the Wasserstein-1 distance has the following very simple form: where F and G denote the cumulative distribution functions of the distributions ρ 0 and ρ 1 , respectively.

Analyzing delay between oscillatory signals
Signal detrending and amplitude normalization The preprocessing step for analyzing the relationship between p53 and p21 signals was to detrend the data. This can be done by fitting polynomials, preferably of low order, to the time series data. By finding an appropriate order that captures the trend of the data, one can then subtract the "trend" from the actual data to obtain the oscillatory component. In dealing with signals like p53 and p21 that can be subjected to abrupt changes when the cells are about or are undergoing division, the detrending should be applied to subsets of the overall time series data that avoids including these abrupt changes.
We were able to obtain cleaner results using a method introduced in our previous work [19]. The first part of this technique called the Detrended Autocorrelation Periodicity Scoring (DAPS) algorithm focuses on the detrending and amplitude normalization using a sliding window embedding, which can be briefly summarized in the following steps: 1. Given the sliding window of length M of an N -length time series x, the sliding windows are arranged into the columns of an M × (N − M + 1) matrix X, so that the i th column of X is In other words, the i th skew-diagonal of X is a constant value equal to the i th value of the time series x i . In this work, M is chosen to be 11 data points, or 5.5 hours, which is roughly the length of one period of p53 [14].
2. The mean of each column of X is subtracted from each element of that column. This normalizes for linear drift in the time series.
3. After point-centering, each column is normalized to be of unit norm. This controls for changes in amplitude in the signal.
4. Lastly, the following operation is performed on the matrix Finally, the i th value of the time series y is then simply the value of any of the elements in the i th skew diagonal of Y .
Quantifying delay using dynamic time warping and cross-correlation Dynamic time warping (DTW) [25] is a method that measures the similarity between two time series by finding an optimal time-ordered correspondence, known as a "warping path," that best maps every index in the first time series to every index in the second time series [5]. Once the signals have been properly detrended and normalized, we apply DTW to determine which of the proteins p53 or p21 leads or lags, as well as the actual delay between the two sequences. It is likely that, because of the periodic nature of these signals, the algorithm may pick up the real delay as well as this value adjusted by integer multiples of the period of these oscillations. We discuss this point in more detail in the Results section. The "warping path" of a pure delay is expected to be a straight line in this analysis with the true delay being a stronger signal than the "secondary" alignments that are offset by a certain amount of periods.
Similarly, cross-correlation may be used to determine a correspondence between two time series by shifting one relative to the other. This work introduces the concept of using DTW as an alternative to cross-correlation to capture signal delays, with the hope that other applications may benefit from it even though similar results were obtained using both techniques.

Long-term trend in the signals and occupancy density
The analysis of long-term trends in the signal is performed by using moving averages (m.a.) [8]. By removing the short-term fluctuations, this transformation acts as a low-pass filter that allows highlighting the trend from the overall signal. Using a relatively small window of 9 data points or 4.5 hours, all the time series shown in this work could be cleanly smoothed out. Since the amplitude of the oscillations can be quite large relative to the long-term trend, especially for p53, previous studies focused greatly on describing the behavior of the oscillations. A primary argument of this work is to show that focusing on the trend and not the oscillations for p21 and p53 provide important and novel insights about the current state of a cell as it pertains to cell cycle arrest.
Using moving average values of p21 and p53, 2-D maps of occupancy density were constructed to better understand where these values tend to lie with respect to their likelihood of undergoing mitosis or treatment conditions. By constructing a vector of these moving averages for all time points, each time point then takes up a (x, y) coordinate in space based on (p53 m.a., p21 m.a.). In this work, a grid was delimited by approximately the minimum and maximum values for each of the axis. Due to the large discrepancy of low and high values as well as the exponential nature of chemical reactions, a log 10 (p21 m.a.) vs. log 10 (p53 m.a.) was used and discretized by 0.1 increment on this scale. This discretized grid is used to count each time a coordinate lies within one of the cells. This is then normalized into the final occupancy density maps.

Experimental single cell datasets and general trends
The analysis in this work is based on two datasets of individual human cells that were published in [23,29]. The cells were monitored using fluorescence live-cell imaging for about a week after being subjected to irradiation. The single cells in the first dataset were subjected to gamma irradiation at 0 Gy, 2 Gy, 4 Gy, and 10 Gy. The culture media was replenished daily, except for a subset of the cells subjected to 10 Gy irradiation that had their culture media preserved throughout the experiment. Fluorescence markers for p21 and p53, as well as the mitosis events, were tracked throughout the experiment. The single cells in the second dataset were subjected to x-ray irradiation at 0 Gy, 0.5 Gy, 1 Gy, 2 Gy, 4 Gy, and 8 Gy. Fluorescence markers for p21, p53, and geminin, as well as mitosis events, were monitored throughout the experiment. The values used in this work are fluorescence intensities (a.u.) of p53-mNeonGreen reporter for p53, the mKate2 protein to form a p21-mKate2 complex, and the CFP-hGeminin(1-110) reporter for cell cycle progression [24].
The time series data is presented in Fig. 1 for the single cells exposed to various grades of gamma radiation and in Fig. 2 for cells exposed to x-ray radiation. As shown in Fig. 1(a) and Fig. 2 (a), between two-thirds and three-quarters of the unirradiated cells divide 3-5 times during the observation period of 5 days. As the radiation levels increase from 0 to 0.5 to 1 Gy in the x-ray dataset, the fraction of cells that divide often progressively diminishes and is mostly replaced by cells that do not divide or divide only once. The division results for gamma radiation levels of 4 Gy and x-ray radiation of 2 Gy are very similar, indicating the possible more harmful nature of x-ray radiation at equivalent exposure levels that translate into cells dividing less. Due to the two datasets being exposed to different radiation sources and the difference in intensity of the fluorescent light source at the time of the experiment, we reasoned that we would analyze these datasets independently and we focused on results that were conserved between the experimental conditions.  Looking at panels Fig. 1(b)-(c) and Fig. 2(c)-(d), p53 levels can be seen to be uniformly oscillating in irradiated cells with the exceptions of cells that divide only once during the 5 days and at least 3 days after radiation exposure. These cells tend to see dramatic increases in p53 levels observed after undergoing their sole mitosis event, consistent with previous reports [29]. On the other hand, p21 levels, exhibit much lower amounts of oscillations, and the trends in the data are of much greater importance to the local fluctuations. From Fig. 2(c) and (d), the relationship is clear that that the upticks of CFP-hGeminin(1-110) correlates with temporal drops in p21 levels. This is because this geminin reporter only accumulates during the S/G 2 /M parts of the cell cycle [30] and p21 is degraded in the G1/S transition [6,10]. Besides these "gaps" in p21 during the G1/S phases, p21 levels are extremely low in unirradiated cells. p21 levels also rapidly converge to very high values for cells that undergo cell cycle arrest and never divide in the course of the experiment.
4.2 p21 and p53 dynamics cluster by number of total observed divisions rather than radiation levels In the interest of classifying the dynamics of p21 and p53 time series, the t-SNE classification using the 1-Wasserstein distance was applied to the data to generate Fig. 3 For each dataset, we matched the t-SNE results on the dynamics with the appropriate labels into four subplots: p53 dynamics by levels of radiation exposure, p53 dynamics by the number of total observed divisions, p21 dynamics by levels of radiation exposure, and p21 dynamics by the number of total observed divisions. From these results, the relation between p21 dynamics and the number of total observed divisions creates the best separation between the groups in the data. The relationship between p21 dynamics and levels of radiation exposure is still significant, but there are a lot of more overlaps between these subcategories of cells. Similarly, the data clusters better relating p53 dynamics to total observed divisions than p53 dynamics to treatment conditions. However, it is clear that considering p53 and p21, the latter is the better indicator as to whether a cell decides to arrest, or divide less often than what one would expect from a healthy cell.
The most meaningful result from the t-SNE classification is that the total number of observed divisions is inversely correlated with the mean values of observed p21 and p53. The correspondence between a lower frequency of observed mitosis events and higher levels of mean p21 levels is particularly striking in the data. This correspondence is not as clear in p53, thus implying that high levels of p53 usually correspond to a lower number of divisions, but that there are also other factors a cell must consider in making the decision to undergo cell cycle arrest. Thus, p21 signaling, being a more downstream indicator, is a cleaner reflection of the cell's decision to undergo cell cycle arrest.  High levels of irradiation such as 8 Gy for the x-ray dataset and 10 Gy for the gamma-ray dataset do contain mostly cells that have high levels of p21 and p53. These cells tend to either not divide or divide only one time. For cells that are treated with "medium" levels of radiation, a fraction behaves similarly to healthy cells and the rest behaves like damaged cells that were subjected to "high" levels of radiation. Only the healthy cells that were not subjected to radiation damage separate cleanly from the rest of the data. The existence

of healthy behaving cells in irradiated conditions is most easily captured in Figs. 3(c) and 4(c).
There is almost always a number of cells that divide normally within 20-30 hours of irradiation exposure similar to unirradiated cells at lower levels of radiation. This attests to the heterogeneous nature of the cells when exposed to radiation and the reason why clustering cells by levels of radiation exposure does not cleanly separate the data into clusters. As expected, the ratio of damaged cells compared to healthy cells increase with increased irradiation exposure.

Quantifying the delay between p21 and p53 signaling
Given that p21 lies downstream of the p53 signaling pathway, part of the observed fluctuations in p21 is directly influenced by the oscillatory component of p53. By detrending and normalizing the two signals using part of the DAPS method as shown in Fig. 5(a) for the gamma radiation dataset, these transformed signals were used to calculate the delay between the two signals by warping or "aligning" them to one another using DTW. By performing this step for all time series in the gamma radiation dataset, the aggregated warping paths are shown in Fig. 5(b). There are mostly two main visible straight lines that cut across the graph, one indicating that the p21 signaling is lagging p53 signaling by about 3.5 hours and the other indicating that p21 leads p53 by 2.5 hours. Because the two signals are periodic, the later result is a natural artifact that arises from the fact that the patterns of transformed signals repeat themselves. Thus, the "correct alignment" offset by full periods of the signal will yield local maximum in terms of correlation between the compared signals.
Both the normalized cross-correlation and DTW results in Fig. 5(b)-(c) showed a stronger correlation when aligning with p53 leading p21 by about 3.5 hours. This holds whether looking at cells that did not divide, or divided 1, 2, or 3 times during the 5 days following radiation exposure. This is also identical to trends observed in similar results generated for the x-ray radiation dataset that are shown in Fig. S3. These observations do not hold when looking at cells that divide 3+ times over 5 days whether one looks at the gamma or x-ray radiation datasets. This can be partially explained by p21 signaling being suppressed during the G1/S phases, which represent a major fraction of the cell cycle, reducing the ability to compare the two signals. The oscillations are also noisier when cells are dividing at the pace of healthy cells. Thus, these single cell datasets containing a majority of cells for which the rate of mitosis was slowed down or for which cells are arrested made this analysis on the delay most appropriate.

Long-term p21 and p53 trends play a key role in determining cell behavior vis-à-vis of mitosis
In Fig. 6(a) for the gamma radiation dataset and Fig. S1 Videos of the datasets that show the evolution at every measured time point over 5 days are available at https://github.com/phongatran/p21p53/. There is a slight long-term drift in the two signals that makes the data points slowly disperse, which may potentially be due to a loss in fluorescence signal over time. By aggregating the m.a. values of these two markers, the occupancy densities are shown in Fig. 6(b) for the gamma radiation dataset and Fig. S1(b) for the x-ray radiation dataset. From these depictions, it is clear that p21 levels is the main determinant as to how often cells divide over the course of 5 days after irradiation. Note that the cluster of cells moves progressively upwards with the p21 m.a. as the y-axis, looking from > 3 down to 0 divisions. Additional plots are available in Fig. S2, illustrating the occupancy densities broken down by radiation levels for both datasets. To further understand the relationship between the trends in p21 and p53 with respect to cell cycle arrest, the m.a. values are centered around mitosis in Fig. 7 for cells that divide and are exposed to x-ray radiation. In Fig. 7(a), the cells that do not divide have a characteristic jump in p53 m.a. values that then subside within the first few hours after the cell receives the irradiation damage. The p21 m.a. values are constantly increasing for the next 20-30 hours until they reach a plateau of about ∼3000-5500. These non-dividing cells remain in this state with little variation in geminin, p21, or p53 m.a. values. Fig. 7(b)-(c), undergo the process of mitosis at significantly higher values of p21 m.a. values than cells that divide more often. These cells also seem to exhibit a substantial increase of both p21 and p53 m.a. values before mitosis. These elevated levels of p53 subsist long after these cells undergo their single mitosis event. Of note, the behavior between "early single-divider cells," defined as cells which undergo mitosis less than 3 days after exposure, exhibit distinct behavior from "late single-divider cells" (also called "escapers" in [23]). Early-single divider cells have noticeably lower levels of p21 and p53. There is also a significant jump in p53 levels in late single-divider cells near the mitosis event and afterward. Looking in detail at cells that divided exactly twice (Fig. 7(d)-(e)) or 3 times (Fig. 7(f)-(h)), the profiles are almost identical for the same number of divisions across all three markers. There does not appear to be a significant difference in the long-term trends of p21 levels before the mitosis process has taken place and afterward, once the values are stabilized. Cells that divide more than 3 times over five days in Fig. 7(i) exhibit in aggregate the lowest levels of p21, p53, and geminin. Note that this subplot shows the behavior across all observed mitosis events, unlike the other subplots.

Single-divider cells, shown in
Cells that do not divide despite high levels of geminin are shown in Fig. 8, split into 6 clusters using k-means with correlation as the distance. These subfigures show various instances in which cells that are in the process of dividing decide to forego division or have yet to divide. These cells had the characteristics of expressing geminin levels higher than 500 a.u. for a period of more than 12.5 hours, but for which no division event was recorded. There were 175 cells in the x-ray dataset. It is relatively clear that a high level of p21 is a necessary condition at the moment that the cell reverses this decision and the levels of geminin start dropping without mitosis occurring. All these cells also have an early spike in p53, which may be due to the radiation exposure at t = 0. When compared to cells that do divide in the previous panel (Fig. 7), it takes longer for geminin to reach a peak, > 20 hrs, from the moment geminin levels start noticeably rising as opposed to 15 − 20 hrs for dividing cells. At the end of this reversion, cells return to having high m.a. values of p53 and p21, similarly to non-dividing cells as shown in Fig. 7(a). Figure 8: The moving averages for the expression levels of geminin, p53, and p21 are plotted and centered around mitosis (at t = 0) for the x-ray irradiation dataset for cells that showed signs of division, but reversed course. These cells exhibited high levels of geminin for a period of more than 12.5 hours (>500 a.u.), but did not undergo mitosis. These cells were then split into 6 clusters using k-means clustering with the distance used being correlation. The middle line represents the 50 percentile value for a given relative time, while the top and bottom of the shaded region represent the 75 percentile and 25 percentile values.
In Fig. 9, the same results as Fig. 7 are shown for the gamma radiation dataset. The geminin levels were not measured in these single cells. The m.a. trends in p21 and p53 are nearly identical between the two datasets. A difference is that raw p21 levels in the gamma dataset are of lower values than the x-ray radiation dataset.
We should also note that the x-ray radiation profiles were very similar before and after mitosis. Looking at Fig. 9(d)-(e), one can notice elevated levels in the case of 2 divisions, after the first mitosis and before the second mitosis.   Fig. 10, there is a noticeable difference in distribution between cells that do not divide and for each additional division that these cells take over the course of the five days. This makes it clear that mean p21 n.d. levels are primarily responsible for the single cell decision making with respect to mitosis with lower values leading more divisions over the same period of time. Thus, cells do not simply undergo a binary decision of dividing or not dividing. Mean p21 n.d. levels in a given single cell strongly correlates with how frequently mitosis occurs in that cell. Plotting the mean of the non-dividing regions of p21 and p53 against one another also provides a cleaner correlation, though still noisy, between the two variables as shown in Fig. 6(a) or Fig. S1(a), where one observes higher levels of p21 in the presence of higher levels of p53.
Shown in Fig. 11 are swarm plots similar to Fig. 10, but limited to the ranges of 18-36 hours from the 5 days that the data was recorded after irradiation. The reason for picking this time frame is because it takes time for the long-term trends in p21 to settle in that may take up to a few days for some cells. Through these plots for both Fig. 11(a) and (b), we show that using these 15% of early data out of the full five days provide similar trends to those shown in Fig. 10. This provides the ability to predict the rate of division of a cell over a longer interval based on early observations because these levels of p21 tend to maintain before and after mitosis.

Discussion and Conclusion
Two sets of human single cells, treated with various levels of x-ray and gamma radiation, were monitored for cell division, p21, p53, and hGeminin(1-110) levels after radiation exposure [23,29]. From the data on mitosis occurrences, it is clear that the behavior of cells within a given treatment regimen, in terms of radiation levels and type, can be quite heterogeneous with cells dividing at different rates. Using t-SNE classification with the 1-Wasserstein distance, it was found that cells can be readily clustered using the number of cell divisions they undergo over five days, using the "histograms" of p21, and to a lesser extent, p53. This indicates that p21 dynamics plays a key role in the cell decision to undergo mitosis or not, and also affects the rate at which this occurs.
While it is known that biologically, p21 is located downstream of the p53 signaling pathway, this work shows how p53 signaling is an integral part of the p21 response by analyzing the oscillatory components of the two signals. Through detrending and normalizing the two signals using the DAPS technique, it was shown how it is expected for p53 to lead p21 by about 3.5 hours and that the period of the oscillations is about 6 hours. This analysis was made possible because of the rare occurrence of mitosis in heavily damaged cells, as p53 signaling is heavily suppressed in the G1/S phases.
From studying time-lapses of the moving averages of p21 and p53, it was found that p53 trends tend to settle within hours of cell damage, while p21 values may take up to a few days to fully settle into their long-term trends. Most cells that did not divide follow a pattern of expressing very high levels of p21 that builds up from the onset of the radiation damage. However, in rarer cases, cells may be attempting to divide through the high expression of geminin and the degradation of p21, but ultimately reversed course and these cells then express very high values of p21. It is also made clear through this work that the cell decision to divide is more complex than simply the idea that cells resume division after having repaired the radiation damage. In the studied datasets, cells that had high levels of p21 before mitosis oftentimes recovered similarly high p21 levels after mitosis. Finally, we found, looking at the mean values of p21 that exclude the parts where cells are actively dividing, that the pace at which cells divide is regulated by the long-term trends of p21. This would also suggest that cells do not necessarily decide to cease division when considered too damaged, but that very high levels of p21 make it unlikely for cells to divide within the observation period of five days.