Bilateral Feedback in Oscillator Model Is Required to Explain the Coupling Dynamics of Hes1 with the Cell Cycle

Abstract

Biological processes are governed by the expression of proteins, and for some proteins, their level of expression can fluctuate periodically over time (i.e., they oscillate). Many oscillatory proteins (e.g., cell cycle proteins and those from the HES family of transcription factors) are connected in complex ways, often within large networks. This complexity can be elucidated by developing intuitive mathematical models that describe the underlying critical aspects of the relationships between these processes. Here, we provide a mathematical explanation of a recently discovered biological phenomenon: the phasic position of the gene Hes1’s oscillatory expression at the beginning of the cell cycle of an individual human breast cancer stem cell can have a predictive value on how long that cell will take to complete a cell cycle. We use a two-component model of coupled oscillators to represent Hes1 and the cell cycle in the same cell with minimal assumptions. Inputting only the initial phase angles, we show that this model is capable of predicting the dynamic mitosis to mitosis behaviour of Hes1 and predicting cell cycle length patterns as found in real-world experimental data. Moreover, we discover that bidirectional coupling between Hes1 and the cell cycle is critical within the system for the data to be reproduced and that nonfixed asymmetry in the interactions between the oscillators is required. The phase dynamics we present here capture the complex interplay between Hes1 and the cell cycle, helping to explain nongenetic cell cycle variability, which has critical implications in cancer treatment contexts.

1. Introduction

The transcription, translation and, ultimately, protein assembly of a specific gene is known as its gene expression. The temporal dynamics of gene expression are known to be a key factor in determining cell fate in many contexts [1,2,3]. Biological and technical advances (e.g., in fluorescent protein reporters and microscopy [4]) have meant that dynamic data has become increasingly abundant (for example, reviewed in [5]), yet dynamic gene expression remains difficult to measure and interpret [6,7,8].

The cell cycle is a biological process during which a cell undergoes an ordered series of predetermined procedures and checks, culminating in its division into two cells [9]. It has been well understood through the lens of gene expression and post-transcriptional dynamics that there are many proteins (e.g., cyclins) whose abundance peaks once and at specific times in each cell cycle [10,11] (Figure 1A bottom panel). It is also known that the cell cycle is closely linked to cell fate [12,13,14]. Hes1 (hairy and enhancer of split-1) is a transcription factor (TF) protein that is encoded by the Hes1 gene in mammals; it is known to regulate cell fate [15,16,17,18,19] but also interferes with the cell cycle through repression of cell cycle-related proteins [20,21,22]. Hes1 is also known to be important in cancer [23,24]; for example, it has been shown that elevation of Hes1 can contribute to cancer metastasis, and Hes1 overexpression can even lead to drug resistance [25].

In [26], we used single-cell live imaging to track Hes1 expression in a human ER+ breast cancer cell (BCC) line, known as MCF-7 [27,28], which contains breast cancer stem cells (BCSCs) [29] alongside tumour bulk cells. We discovered that the Hes1 protein exhibits a large oscillatory pattern averaging a circadian-like periodicity (≈25 h), while an oscillatory pattern with a shorter periodicity (≈4 h) was found when this larger oscillatory trend was removed. These nested ultradian (shorter than one day but longer than one hour in periodicity) oscillations, which are caused by transcriptional self-repression, a well-researched characteristic of Hes/Her genes [30,31,32], showed no correlation with the high levels of heterogeneity in cell cycle durations found in the MCF7 population. Conversely, the circadian-like oscillatory Hes1 trend within the cells did provide patterns related to longer and shorter cell cycles. We determined that it is the momentary position within these oscillations at the beginning of a cell’s life that has the greatest correlation with the eventual length of time taken before that cell divides (demonstrated in Figure 1). More specifically, when clustered by their Hes1 dynamics into three distinct groups, cells whose beginning mitoses were at or just after a trough in oscillation (cluster 3 as illustrated in Figure 1D) predominantly possessed longer cell cycles (Figure 1F) than those that began their life in a different part of the oscillatory wave (clusters 1 and 2 in Figure 1B,C). However, it remains unclear how this feature of Hes1 protein oscillations at the beginning of mitosis in a breast cancer cell results in specific mitosis to mitosis dynamics and, further, why this corresponds to distinguishable differences in cell cycle durations. Therefore, we ask why is it that our cluster 3 cells exhibit longer cell cycles and why do clusters 1 and 2 have similar and shorter cell cycle durations despite having different initial Hes1 dynamic behaviour? This paper aims to shed light on these questions.

Mathematics, and particularly dynamical systems theory, has had an impact on understanding gene expression dynamics [33,34,35] (specifically Hes1 dynamics [36,37,38,39]), the cell cycle [40,41,42,43,44,45] and their links to cell fate [46,47,48]. Although the relationship between Hes1 and cell cycle has previously been investigated in large, multiple parameter models (in [49], for example), thus far, no two-component (Hes1/cell cycle) model exists. Typically, models involving gene expression are very complicated and comprise vast multi-parameter spaces that are difficult to explore and relate to observations from real data. These models also require extensive experimental validation in order to verify the necessity of each component. Furthermore, we note that there has not yet been a study that models Hes1 dynamics with oscillations at the circadian scale, although coupled oscillator models such as [50,51] have been used to study circadian rhythms with the human sleep-wake cycle.

Our goal in this paper is to implement a minimal and thus, intelligible mathematical model in order to understand and make predictions on specifically the interaction between Hes1 and the cell cycle, which ideally can then reproduce the results in [26] (which are illustrated in Figure 1) without the need for many components. Discovery and critical evaluation of an intuitive mathematical explanation behind the relationship between the oscillatory expression of Hes1 and the cell cycle may provide clearer insight into the systems at play. By achieving this, one may be able to systematically probe this biological interaction via specifically designed experimentation. This can allow for causal effects between the two phenomena to be uncovered, resulting in a better understanding of the study of cell cycle elongation/dormancy, which, in a cancer context, may prove informative for clinical treatment.

2. Mathematical Modelling Strategy

While mathematical modelling approaches to gene expression dynamics often involve simulations of protein and/or mRNA abundances over time [52,53,54,55], we take the route of modelling the relative temporal progression of two coupled oscillators. By this, we mean mathematical models that are defined by systems of differential equations whereby the dependent variables each represent the phase angle of an oscillator. By “phase”, we mean the position at a given time, relative to some point of origin, within the progression of one cycle of a periodic object. We may define this mathematically as a real-valued angle $\theta \in \left[0,2\pi \right)$ of a point on the circumference of a circle from some predefined original point on the circumference. Models of this form have previously been applied to the cell cycle and the circadian clock oscillators in differing contexts in [56,57,58]: these publications provided particular inspiration for this paper. For the model we implement here, our two coupled oscillators shall be Hes1 gene expression dynamics (as described in [26] and displayed in Figure 1) and the cell cycle. Therefore, as such, we label our model variables as ${\theta }_{\mathrm{Hes1}}$ and ${\theta }_{\mathrm{C}.\mathrm{C}.}$ and periodic functions determining their temporal behaviour as $f_{\mathrm{Hes1}},f_{\mathrm{C}.\mathrm{C}.}$.

The decision to model Hes1 as an oscillator is justified since it was found that gene expression of Hes1 is oscillatory (and corroborated via multiple methods) in [26] on a near circadian scale. It may be argued though, that the cell cycle ought not to be considered an oscillator (as is done in [59]) due to its series of discrete points in time: checkpoints, potential stoppage events, etc. However, we shall take a more simplified view and assume that processes such as cell differentiation, apoptosis and cell cycle arrest, etc., do not occur. This allows us to assume that each time a cell division occurs, a new cell cycle begins, and thus, cell division can be seen as a repeated, periodic event in time. For the purposes of a mathematical model, we may, therefore, consider the cell cycle process as a continuous oscillator (like in [60,61]) with an inherent periodicity: i.e., the time between divisions. Furthermore, here we shall take a more abstract view of the cell cycle checkpoint events themselves in that we do not make any assumptions on the specific durations of the biological cell cycle phases: G1, S, G2 [9,62], which are present within the cell cycle. The results from multiple control experiments in [26] also suggest that Hes1 and the cell cycle are coupled oscillators further justifies our approach.

2.1. Phase Representation of Bcc Data

When choosing the specifics of our intended model, i.e., the nature of functions $f_{\mathrm{Hes1}}$ and $f_{\mathrm{C}.\mathrm{C}.}$, we examined the data we intend to model. One key aspect of the results from [26] is important here: in each case (Figure 1B–D), irrespective of the phase at the beginning of cell cycles, Hes1 invariably goes through a trough followed by an upward trend in dynamics before the cell cycle ends. This hints at a pattern forming between the phases of both oscillators as time progresses. To study instantaneous phase over time, our pseudo-time-series signals are plotted in Figure 1B–D respectively; we used hilbert.m and the angle.m functions in MATLAB. The function hilbert.m develops the analytic signal $\left(\xi \left(t\right)\right)$ of a signal vector $\mathit{s}$, which is defined in [63] as

$$\xi \left(t\right)=s\left(t\right)+i{s}_{\mathrm{H}}=A\left(t\right)e^{i\theta \left(t\right)}$$

(1)

where here, the function

$${\displaystyle s_{\mathrm{H}}\left(t\right)={\pi }^{-1}{\int }_{-\infty }^{\infty }\frac{s\left(\tau \right)}{t-\tau }d\tau }$$

(2)

is the Hilbert Transform of the signal vector $\mathit{s}$. In (1), $A\left(t\right)$ is the signal’s instantaneous amplitude, and $\theta \left(t\right)$ is the instantaneous phase. The function angle.m used upon an analytic signal vector returns $\theta \left(t\right)$ as a vector (with the same length as $\mathit{s}$) of real values modulo $2\pi$. The vector $\mathit{\theta }$ is now the estimated phase reconstruction of $\mathit{s}$. As an example, using functions hilbert.m followed by angle.m on a perfect sine wave will return a sawtooth wave. The phase reconstructions of the three clustered pseudo-time Hes1 signal traces (Figure 1B–D) are now vectors of length $T\in {\mathbb{Z}}^{+}$, which correspond to the position in the phase of each trace at all points in pseudo-time; we shall label them as ${\mathbf{\Theta }}_{\mathrm{Hes1}}^{\left(\mathbf{i}\right)}$ where $i=1,2,3$ denotes the cluster mean traces from which they were derived. We may then plot these against a uniformly spaced cell-cycle phase vector (${\mathbf{\Theta }}_{\mathrm{C}.\mathrm{C}.}$) between 0 and $2\pi$ of the same length to produce trajectories of our data in two-dimensional phase space. This phase space can then be drawn onto a square torus diagram (Figure 2B), where the phase of both oscillators are unfolded along the x and y-axes.

One noticeable feature in Figure 2 is that torus diagram trajectories of clusters 2 (blue line) and 3 (green line), despite beginning quite close together, appear to initially be moving away from each other as though repelled. Additionally, both trajectories, along with the cluster 1 trajectory (red line), appear to collect around a similar area (upper right-hand quadrant) in the phase space as the cell cycle oscillator phase angle reaches $2\pi$. This is also illustrated in Figure 2C, where in the far right-hand side column, all dots are close together around the top of the circle. This implies that in all cases, as the cell cycle reaches an end, the Hes1 oscillator is attracted to a particular phase. This indicates that a form of phase locking between the Hes1 and the cell cycle oscillators may exist. To that end, we choose the following model, which is known to exhibit this property.

2.2. Specifying the Coupled Oscillator Model

We model the interactions between cell cycle and Hes1 using the following previously characterised system of coupled oscillators in [64]:

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill \frac{d{\theta }_1}{dt}=& {\omega }_1+K_{1}sin\left({\theta }_2-{\theta }_1\right),\hfill \\ \hfill \frac{d{\theta }_2}{dt}=& {\omega }_2+K_{2}sin\left({\theta }_1-{\theta }_2\right),\hfill \end{array}\end{array}$$

(3)

where ${\theta }_i$ is the phase of oscillator i, ${\omega }_i$ is the free-running intrinsic frequency of oscillator i and $K_i$ is the strength of effect from oscillator j onto oscillator i for $i,j=1,2$ and $i\ne j$. This particular system of ordinary differential equations has previously been used to model circadian rhythms in humans [50,51], which makes it well-suited to be applied to our data since we see circadian-like Hes1 periodicity in human breast cancer stem cells in [26]. The above model is also an asymmetric, two-component version of the Kuramoto coupled oscillator model [65]. Although, in some cases, the Kuramoto model is used for purely theoretical studies, it has been previously used to good effect when used in conjunction with cell biology [66,67]. Analysis of the system Equation (3) will show that if $K_1,K_2>0$, then by considering the difference in phases between the oscillators $\varphi ={\theta }_1-{\theta }_2$ we get

$$\frac{d\varphi }{dt}=\frac{d{\theta }_1}{dt}-\frac{d{\theta }_2}{dt}={\omega }_1-{\omega }_2-\left(K_1+K_2\right)sin\left(\varphi \right).$$

This ODE has two fixed points $\left(\mathrm{i}.\mathrm{e}.,\frac{d\varphi }{dt}=0\right)$ when $\left|{\omega }_1-{\omega }_2\right|<K_1+K_2$ (shown by means of vector field diagrams in Equation (3)). Therefore, two phase-locked solutions will exist, one stable and one unstable, and the phase space trajectories will asymptotically approach the stable phase-locked solution in which both oscillators will run at a constant compromise frequency defined as

$${\omega }^{*}=\frac{K_1{\omega }_2+K_2{\omega }_1}{K_1+K_2}$$

(4)

and be separated by a constant phase difference $\left({\varphi }^{*}\right)$ [64]. These two phase-locked solutions will appear as diagonal lines with slope 1 when plotted on a square torus diagram and have also been called “attractor” and “repeller” lines in [57]. We now consider this model with the two oscillators representing Hes1 and the cell cycle; we aim to use this model to reproduce the phase data shown in Figure 2 and thus be able to find a mathematical explanation with the assistance of the previously-stated analysis. Examining Equation (3) closely makes us aware that an unwanted scenario may occur: suppose, for example, that ${\theta }_2-{\theta }_1=\frac{3\pi }{2}$ then $sin\left({\theta }_2-{\theta }_1\right)=-1$ and constants ${\omega }_1<K_1$, then we get $\frac{d{\theta }_1}{dt}<0$. Now, we know that the cell cycle cannot backtrack in its progression [62]. Therefore, allowing this to occur in our model would be unrealistic; hence we shall ensure that rates $\frac{d{\theta }_1}{dt},\frac{d{\theta }_2}{dt}\nless 0$, i.e., our oscillators may not progress backwards. We can impose this by using the following substitutions $K_1={\omega }_1{\kappa }_1,K_2={\omega }_2{\kappa }_2,\mathrm{where}{\kappa }_1,{\kappa }_2\in [0,1]$. Applying these substitutions, factorising and replacing variable names to our specifications yields the model with which we shall continue:

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill \frac{d{\theta }_{\mathrm{Hes1}}}{dt}=& {\omega }_{\mathrm{Hes1}}\left(1+{\kappa }_{\mathrm{Hes1}}sin\left({\theta }_{\mathrm{C}.\mathrm{C}.}-{\theta }_{\mathrm{Hes1}}\right)\right),\hfill \\ \hfill \frac{d{\theta }_{\mathrm{C}.\mathrm{C}.}}{dt}=& {\omega }_{\mathrm{C}.\mathrm{C}.}\left(1+{\kappa }_{\mathrm{C}.\mathrm{C}.}sin\left({\theta }_{\mathrm{Hes1}}-{\theta }_{\mathrm{C}.\mathrm{C}.}\right)\right),\hfill \end{array}\end{array}$$

(5)

where ${\theta }_{\mathrm{Hes1}},{\theta }_{\mathrm{C}.\mathrm{C}.}\in [0,2\pi )$ are continuous functions of time that represent the phase progression of Hes1 expression oscillations and the cell cycle, respectively. Constant parameters ${\omega }_{\mathrm{Hes1}}>0$ and ${\omega }_{\mathrm{C}.\mathrm{C}.}>0$ are the Hes1 and cell cycle oscillators’ intrinsic frequencies, respectively, and ${\kappa }_{\mathrm{Hes1}}$ is now a scaled coupling constant that represents the strength of effect that the cell cycle oscillator phase has upon the Hes1 phase and vice-versa for the parameter ${\kappa }_{\mathrm{C}.\mathrm{C}.}$. These coupling strength parameters are currently unassigned, but we allow the possibility for either or both to be equal to zero, implying the scenario of no coupling or of unidirectional coupling, respectively. As per our earlier definition, we make the following arbitrary choices:

• ${\theta }_{\mathrm{Hes1}}=0$ will be considered the peak of an oscillatory Hes1 wave, and ${\theta }_{\mathrm{Hes1}}=\pi$ is the trough.
• ${\theta }_{\mathrm{C}.\mathrm{C}.}=0$ will denote the beginning and the end of a cell cycle, i.e., a cell division event.

2.3. Numerical Implementation

In Section 3, the results are shown from our numerical simulations: solutions to the coupled system of Equation (5) are approximated in MATLAB via a fourth-order Runge-Kutta (RK) numerical scheme with a step-size of $h=0.1$. This scheme is described as follows: we set the cell cycle oscillator initial condition as ${\theta }_{\mathrm{C}.\mathrm{C}.}\left(t_0\right)=0$ (where $t_0$ is our initial time step) and set the Hes1 oscillator initial condition as one of the three initial phase values found from the phase reconstruction of the clustered Hes1 dynamics in the data: that is, ${\theta }_{\mathrm{Hes1}}^{\left(1\right)}\left(t_0\right)={\mathbf{\Theta }}_{\mathrm{Hes1}}^{\left(1\right)}\left(1\right)\approx 6.0830$, ${\theta }_{\mathrm{Hes1}}^{\left(2\right)}\left(t_0\right)={\mathbf{\Theta }}_{\mathrm{Hes1}}^{\left(2\right)}\left(1\right)\approx 2.0164$ and ${\theta }_{\mathrm{Hes1}}^{\left(3\right)}\left(t_0\right)={\mathbf{\Theta }}_{\mathrm{Hes1}}^{\left(3\right)}\left(1\right)\approx 3.3613$ (these values correspond to the radial coloured lines in Figure 1E). We choose the constant intrinsic frequency parameter values ${\omega }_{\mathrm{Hes1}}$ and ${\omega }_{\mathrm{C}.\mathrm{C}.}$ to be the reciprocals of the median Hes1 periodicity (24 h) and the median cell cycle length estimates (27 h), respectively, from the dataset of 164 BCCs found in [26]. At the nth iteration of the numerical scheme, we have the functions

$$\begin{array}{cc}\hfill f_{\mathrm{Hes1}}(t_n,{\theta }_{\mathrm{Hes1}}\left(t_n\right),{\theta }_{\mathrm{C}.\mathrm{C}.}\left(t_n\right))=& {\omega }_{\mathrm{Hes1}}(1+{\kappa }_{\mathrm{Hes1}}sin\left({\theta }_{\mathrm{C}.\mathrm{C}.}\left(t_n\right)-{\theta }_{\mathrm{Hes1}}\left(t_n\right)\right)),\hfill \\ \hfill f_{\mathrm{C}.\mathrm{C}.}(t_n,{\theta }_{\mathrm{Hes1}}\left(t_n\right),{\theta }_{\mathrm{C}.\mathrm{C}.}\left(t_n\right))=& {\omega }_{\mathrm{C}.\mathrm{C}.}(1+{\kappa }_{\mathrm{C}.\mathrm{C}.}sin\left({\theta }_{\mathrm{Hes1}}\left(t_n\right)-{\theta }_{\mathrm{C}.\mathrm{C}.}\left(t_n\right)\right)).\hfill \end{array}$$

Each time step iteration is determined by

$$\begin{array}{cc}\hfill {\theta }_{\mathrm{Hes1}}\left(t_{n+1}\right)& ={\theta }_{\mathrm{Hes1}}\left(t_n\right)+\frac{1}{6}h\left(k_1^{\mathrm{Hes1}.}+2k_2^{\mathrm{Hes1}.}+2k_3^{\mathrm{Hes1}}+k_4^{\mathrm{Hes1}}\right),\hfill \\ \hfill {\theta }_{\mathrm{C}.\mathrm{C}.}\left(t_{n+1}\right)& ={\theta }_{\mathrm{C}.\mathrm{C}.}\left(t_n\right)+\frac{1}{6}h\left(k_1^{\mathrm{C}.\mathrm{C}.}+2k_2^{\mathrm{C}.\mathrm{C}.}+2k_3^{\mathrm{C}.\mathrm{C}.}+k_4^{\mathrm{C}.\mathrm{C}.}\right),\hfill \\ \hfill t_{n+1}& =t_n+h,\hfill \end{array}$$

where

$$\begin{array}{cc}\hfill k_1^i& =f_i\left(t_n,{\theta }_i\left(t_n\right),{\theta }_j\left(t_n\right)\right),\hfill \\ \hfill k_2^i& =f_i\left(t_n+\frac{h}{2},{\theta }_i\left(t_n\right)+\frac{h{k}_1^i}{2},{\theta }_j\left(t_n\right)+\frac{h{k}_1^i}{2}\right),\hfill \\ \hfill k_3^i& =f_i\left(t_n+\frac{h}{2},{\theta }_i\left(t_n\right)+\frac{h{k}_2^i}{2},{\theta }_j\left(t_n\right)+\frac{h{k}_2^i}{2}\right),\hfill \\ \hfill k_4^i& =f_i\left(t_n+h,{\theta }_i\left(t_n\right)+h{k}_3^i,{\theta }_j\left(t_n\right)+h{k}_3^i\right),\hfill \end{array}$$

for $n=0,1,2,\dots$, $i,j=\mathrm{Hes1},\mathrm{C}.\mathrm{C}.\mathrm{and}i\ne j$. We consider one simulation to be complete at the time point $t_m,m\in {\mathbb{Z}}^{+}$ at which ${\theta }_{\mathrm{C}.\mathrm{C}.}\left(t_m\right)<2\pi$ and ${\theta }_{\mathrm{C}.\mathrm{C}.}(t_m+1)>2\pi$. This value of $t_m$ is recorded and the solution vectors ${\mathbf{\theta }}_{\mathrm{Hes1}}^{\mathbf{\left(}\mathit{i}\mathbf{\right)}}$, $i=1,2,3$ are then interpolated into pseudo-time (a vector of length T) and compared against their corresponding data vector (${\mathbf{\Theta }}_{\mathrm{Hes1}}^{\left(\mathbf{i}\right)},i=1,2,3$) as described in Section 2.4. Since ${\kappa }_{\mathrm{Hes1}}$ and ${\kappa }_{\mathrm{C}.\mathrm{C}.}$ are unknown parameters, the simulations are run many times over a finely gridded $(100\times 100)$ parameter space of ${\kappa }_{\mathrm{Hes1}}\times {\kappa }_{\mathrm{C}.\mathrm{C}.}$. This process is illustrated in the schematic in Figure 3.

2.4. Model Optimisation Strategy

We are interested in determining at which points in the parameter space ${\kappa }_{\mathrm{Hes1}}\times {\kappa }_{\mathrm{C}.\mathrm{C}.}$ our model is most successful at fitting the data. In a further test of the system, we shall impose a number of different parameter constraints (unidirectional coupling, symmetric coupling and (for completeness) no coupling) to investigate whether the data is reproducible with a more restricted model. These constraints will give us additional insight into the directionality of interaction between Hes1 and the cell cycle.

The success of our model is judged two-fold: first, it is determined by residual error (Euclidean distance) in the model solution vectors (${\mathbf{\theta }}_{\mathrm{Hes1}}$) and data vectors (${\mathbf{\Theta }}_{\mathrm{Hes1}}$); the smaller the error, the better the fit. Secondly, we can record the number of time steps ($t_m$) required for our model to complete one cell cycle simulation. Of course, these outputs are measured in iterative time steps and thus cannot directly correspond back to real-time cell cycle duration estimates. However, we may gauge the accuracy of our model by comparing these simulated predictions between cluster cases, similar to the bar chart in Figure 1F. The latter of these success gauges is the most interesting since we are only parameterising our model to fit the time-independent ${\mathbf{\Theta }}_{\mathrm{Hes1}}$ dynamics.

We may write out our optimisation problem more formally in the following manner:

$$\begin{array}{cccc}\hfill \underset{{\kappa }_{\mathrm{Hes1}},{\kappa }_{\mathrm{C}.\mathrm{C}.}}{\mathrm{minimise}}& \hfill & & R({\kappa }_{\mathrm{Hes1}},{\kappa }_{\mathrm{C}.\mathrm{C}.})\hfill \\ \hfill \mathrm{subject}\mathrm{to}& \hfill & \hfill \mathrm{Strategy}2:& {\kappa }_{\mathrm{Hes1}}={\kappa }_{\mathrm{C}.\mathrm{C}.},\hfill \\ \hfill \mathrm{or}& \hfill & \hfill \mathrm{Strategy}3:& 0\le {\kappa }_{\mathrm{Hes1}}\le 1,{\kappa }_{\mathrm{C}.\mathrm{C}.}=0,\hfill \\ \hfill \mathrm{or}& \hfill & \hfill \mathrm{Strategy}4:& {\kappa }_{\mathrm{Hes1}}=0,0\le {\kappa }_{\mathrm{C}.\mathrm{C}.}\le 1,\hfill \\ \hfill \mathrm{or}& \hfill & \hfill \mathrm{Strategy}5:& 0\le {\kappa }_{\mathrm{Hes1}}\le 1,0\le {\kappa }_{\mathrm{C}.\mathrm{C}.}\le 1.\hfill \end{array}$$

where

$${\displaystyle R({\kappa }_{\mathrm{Hes1}},{\kappa }_{\mathrm{C}.\mathrm{C}.})=\sqrt{\sum _{j=1}^T{\left({\mathbf{\Theta }}_{\mathrm{Hes1}}^{\left(i\right)}\left(j\right)-{\theta }_{\mathrm{Hes1}}^{\left(i\right)}\left(j\right)\right)}^2},\mathrm{for}i=1,2,3.}$$

Our resulting figures in the next section will illustrate square torus diagrams with trajectories for the best fit of our model found to the data for each parameter constraint. The resulting figures will also display the rebuilt Hes1 dynamics from that simulation: these are produced by inputting the resulting ${\mathbf{\theta }}_{\mathrm{Hes1}}$ solution vector into a cosine function that acts as the inverse action to the use of the Hilbert transform discussed earlier. Finally, they will show a bar chart of the $t_{\mathrm{end}}$ values recorded for the optimally fitted models. We emphasise here that a good result from our model to the data will have:

• A low residual fit (R) value.
• A torus diagram trajectory that closely tracks the coloured lines seen in Figure 2B, and subsequently, a rebuilt Hes1 expression trace that follows the mean clustered Hes1 pseudo-time series Figure 1B–D.
• Cell cycle simulation times resembling the bar chart in Figure 1F.

3. Results of Model Simulations and Comparisons to Biological Data

The following results are simulations of the system of Equation (5) and certain constraints on nonnegative coupling strength constants ${\kappa }_{\mathrm{Hes1}}$ and ${\kappa }_{\mathrm{C}.\mathrm{C}.}$. Initial conditions, methods and how results are displayed are described above in Section 2.

3.1. The Uncoupled Scenario

First, we may consider the seemingly trivial case of no coupling, i.e., ${\kappa }_{\mathrm{Hes1}}={\kappa }_{\mathrm{C}.\mathrm{C}.}=0$.

This constraint reduces Equation (5) to the following linear system:

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill \frac{d{\theta }_{\mathrm{Hes1}}}{dt}=& {\omega }_{\mathrm{Hes1}},\hfill \\ \hfill \frac{d{\theta }_{\mathrm{C}.\mathrm{C}.}}{dt}=& {\omega }_{\mathrm{C}.\mathrm{C}.}.\hfill \end{array}\end{array}$$

(6)

As discussed in [64], analysis can show us that all trajectories on the square torus diagram will be straight lines with slope ${\displaystyle \frac{{\omega }_{\mathrm{Hes1}}}{{\omega }_{\mathrm{C}.\mathrm{C}.}}}$. These straight-line trajectories can be seen quite clearly in Figure 4A–C and, interestingly, follow the cluster 1 and 3 trajectories reasonably well. Figure 4B,E, however, shows that this version of the model does poorly at capturing the cluster 2 dynamics. This could be expected because the data trajectory is not linear on the square torus diagram. In addition to the Hes1 dynamics not matching the data, the predicted cell cycle durations are all essentially identical (Figure 4G), which does not match the experimentation. Having examined this constraint, we now move on with the allowance of nonzero coupling to take place.

3.2. The Symmetric Interaction Strength Scenario

Let us first assume that the strength of effect that each oscillator has on the other is identical, i.e., ${\kappa }_{\mathrm{Hes1}}={\kappa }_{\mathrm{C}.\mathrm{C}.}$, we shall call this the symmetric case since the sine function within Equation (5) will mean that any increase in speed by one oscillator will be mirrored by an equivalent decrease in speed by the other and vice-versa. An interesting finding would be if the real-world data can be recapitulated by a simple model with only three parameters (two frequency constants and one coupling constant) compared to four.

Notably, here the cluster 1 trajectory (Figure 5A) is very similar to the corresponding image in the uncoupled case (Figure 4A) and matches the data well. This time, however, the cluster 2 trajectory also performs well in tracking the data (Figure 5B). The Cluster 3 trajectory fits the data better than in the uncoupled case, and this can be seen in the larger upward trend at the end of the cell cycle in Figure 5F than in Figure 4F. The symmetric model is a clear improvement upon the uncoupled case, and this is quantified in the residual values at the top of Figure 5A–C. However, we can see in Figure 5G that this model is insufficient since there exists a discrepancy between cell cycle predictions made by the model and the duration data seen in Figure 1F; comparisons between clusters are not maintained. The symmetric model predicts shorter cell cycles for cluster 2 than cluster 1, a feature that is not present in Figure 1F. This is likely because, for the dynamics to be matched in Figure 5B, then ${\theta }_{\mathrm{Hes1}}$ must slow down, and the symmetry of the system means that ${\theta }_{\mathrm{C}.\mathrm{C}.}$ must then speed up in an equivalent manner, yielding a faster cell cycle length simulation. Therefore, we can conclude that for our symmetric model, accuracy in either the simulation time or the dynamics must be sacrificed for cluster 2. Thus, the full extent of our data cannot be replicated. This logically leads us to our next step, where we may break the symmetry of our model and consider the asymmetric case where the ${\kappa }_{\mathrm{Hes1}}={\kappa }_{\mathrm{C}.\mathrm{C}.}$ constraint is removed. The results of these simulations are shown below.

3.3. The Unconstrained Interaction Strength Scenario

Figure 6A–F shows that solution trajectories $\left({\theta }_{\mathrm{Hes1}}^{\left(i\right)}\right)$ in the asymmetric case also match ${\mathbf{\Theta }}_{\mathrm{Hes1}}^{\left(i\right)},i=1,2,3$ well, but it is in Figure 6G where a more notable result lies: breaking the symmetry of our system allows for both phase dynamics and cell cycle duration comparisons to mimic that of the data seen in Figure 1. This is an important prediction from our model since we do not originally impose cell cycle lengths into the model, yet the expected lengths are consistent with our observations in [26] when the dynamics are fitted in pseudo-time.

We have found optimal solutions that recapture our data well. However, it may be plausible that a system with even fewer degrees of freedom is sufficient to imitate what we see in Figure 1. Therefore, we also tested (and excluded) the possibility that a one-directional interaction model can explain the data, as displayed in Appendix A. Below, we summarise our results by tabulating residual errors as well as collating the cell cycle length simulations for each constraint case in order to confirm, qualitatively, which of our models works best.

3.4. Asymmetry in Interaction Strength Is Predictive of Elongation in Bcsc Data

Table 1. Table of residual comparisons for each modelling strategy 1 to 5. Values are the minimised R values found by fitting the pseudo-timed model solution ${\mathbf{\theta }}_{\mathrm{Hes1}}^{\mathbf{\left(}\mathit{i}\mathbf{\right)}}$ to ${\mathbf{\Theta }}_{\mathrm{Hes1}}^{\left(i\right)}$ for $i=1,2,3$ within the given parameter constraints. Highlighted in cyan is the strategy with the lowest summed residuals across clusters.

Modelling Strategy	Coupling Strength Parameter Constraint	Cluster 1	Cluster 2	Cluster 3	Sum
1	${\kappa }_{\mathrm{Hes1}}=0,{\kappa }_{\mathrm{C}.\mathrm{C}.}=0$	12.1199	23.5206	6.8118	42.4523
2	${\kappa }_{\mathrm{Hes1}}={\kappa }_{\mathrm{C}.\mathrm{C}.}$	8.9520	3.0552	2.9236	14.9308
3	${\kappa }_{\mathrm{C}.\mathrm{C}.}=0$	9.9229	1.7894	3.0774	14.7897
4	${\kappa }_{\mathrm{Hes1}}=0$	8.8340	7.3753	2.8757	19.085
5	Optimal (no constraints)	8.8352	1.7967	2.8702	13.5021

shows that, overall, the case with the lowest residual error when compared to the data is the one with no constraints upon the coupling strength constants. What is also noticeable is that in Figure 7, the case that matches the data closest is also that of no constraints (modelling strategy 5). We may note that cluster 3’s simulation lengths in the ${\kappa }_{\mathrm{Hes1}}=0$ case (modelling strategy 4) predict long cluster 3 cell cycles and are promising; however, the cluster 2 simulation lengths are shorter than cluster 1 here, which is not what we see in the data (Figure 1F). We find that by constraining the interaction to be unidirectional, in one case, we have a good fit to dynamics but not to the cell cycle simulations, and in the other case, the opposite is true. Therefore, we determine, overall, that the case where the coupling is allowed to be asymmetric yields the best results, thus we may conclude that bidirectionality is a required feature of our system in order to replicate the data we have. Next, we examine the parameter values ${\kappa }_{\mathrm{Hes1}}$ and ${\kappa }_{\mathrm{C}.\mathrm{C}.}$, which provided the optimal fits and discuss why it might be necessary that they are unconstrained.

3.5. Cluster-Dependent Coupling Strength Points to Gene Expression and Cell Cycle Duration Differences

First, we inspect the cluster 1 case where Hes1 is initially around a peak in oscillation: Figure 8A shows that the distance between the model solution and data is very similar for all pairs of ${\kappa }_{\mathrm{Hes1}}$ and ${\kappa }_{\mathrm{C}.\mathrm{C}.}$, i.e., the data can be explained by a wide parameter region. This suggests that there are no possible combinations of coupling strength parameters that will divert Hes1 from its peak-to-peak behaviour within a cell cycle simulation when beginning with these initial conditions.

Secondly, we note that the initial Hes1 conditions for clusters 2 and 3 are similar in phase value $\left(\left|{\theta }_{\mathrm{Hes1}}^{\left(2\right)}\left(t_0\right)-{\theta }_{\mathrm{Hes1}}^{\left(3\right)}\left(t_0\right)\right|<\frac{\pi }{2}\right)$ yet, fascinatingly, require vastly different coupling strength parameters in order to fit the data. The optimal cluster 2 parameters show that a strong effect from the cell cycle onto the Hes1 oscillator and a weak effect in the opposite direction is necessary for the data to be closely replicated. Contrastingly, these effects are required to be weak in both directions for the cluster 3 case, albeit with a slight preference for a greater influence of the Hes1 phase on the cell cycle.

The heatmaps in Figure 8B,C appear to be almost inverted versions of one another; as though they are mutually exclusive cases. This suggests that there is something fundamentally different happening within cells when Hes1 begins a cell cycle before vs. after a trough in an oscillatory wave. Another noteworthy point is that when ${\kappa }_{\mathrm{Hes1}}>{\kappa }_{\mathrm{C}.\mathrm{C}.}$ (e.g., the cluster 2 case), the trajectory on the torus diagram (Figure 6B) is initially horizontally slanted. This implies that the effect from ${\theta }_{\mathrm{C}.\mathrm{C}.}$ onto ${\theta }_{\mathrm{Hes1}}$ is a repressive one, i.e., the cell cycle oscillator is slowing down the Hes1 oscillator. Similarly, we can deduce that only a repressive effect is being felt by the cell cycle oscillator in its progression in the cluster 3 case since the cell cycle predictions are being elongated in Figure 6G; this suggests that ${\theta }_{\mathrm{C}.\mathrm{C}.}$ is slowing down in its progression. Additionally, it appears to be the lack of effect of the Hes1 oscillator onto the cell cycle oscillator, which explains why we see the cluster 2 cell cycle length predictions from the model being similar to that of cluster 1 in Figure 6G: no elongation/shortening has occurred. Since the equations that make up our model are not overly complicated, we are able to analyse the system Equation (5) with our optimal parameter values to gain mechanistic insight into why we see these results mathematically.

3.6. Mathematical Analysis and Long-Term Behaviour of Our Model

Analysis of the ODE system Equation (5) (see [64]) shows us that there exists a stable and unstable phase-locked trajectory on the torus if

$$\left|{\omega }_{\mathrm{Hes1}}-{\omega }_{\mathrm{C}.\mathrm{C}.}\right|<{\kappa }_{\mathrm{Hes1}}{\omega }_{\mathrm{Hes1}}+{\kappa }_{\mathrm{C}.\mathrm{C}.}{\omega }_{\mathrm{C}.\mathrm{C}.}.$$

(7)

We may now check to see whether this inequality holds with the parameters discovered in the modelling in this paper. Succinctly,

Table 2. Table summarizing the optimal ${\kappa }_{\mathrm{Hes1}}$ and ${\kappa }_{\mathrm{C}.\mathrm{C}.}$ values for the model Equation (5) for each cluster case.

Cluster	1	2	3
Initial Hes1 phase (1 d.p.)	6.1	2.0	3.4
${\kappa }_{\mathrm{Hes1}}$	0	0.86	0.01
${\kappa }_{\mathrm{C}.\mathrm{C}.}$	0.99	0	0.13
${\omega }_{\mathrm{Hes1}}$	0.0417	0.0417	0.0417
${\omega }_{\mathrm{C}.\mathrm{C}.}$	0.0370	0.0370	0.0370
Condition (7) met?	Yes	Yes	No

shows the parameter values that developed the results from the optimal strategy, as seen in Figure 6.

By plugging the values from Table 2 into the condition of Equation (7), we can analytically show that the trajectories with the coupling parameter values in the cluster 1 and 2 cases will asymptotically approach a stable, periodic phase-locked solution in which both oscillators (${\theta }_{\mathrm{Hes1}}$ and ${\theta }_{\mathrm{C}.\mathrm{C}.}$) run at a so-called compromised frequency $\left({\omega }^{*}\right)$ [64]. However, since the parameters in the cluster 3 case do not meet the condition of Equation(7), and, in fact, the inequality faces the other direction, this means that the trajectory will exhibit a phenomenon known as quasiperiodicity, in which phase-locking between the cell cycle and Hes1 can never occur. Below, we illustrate simulations of long-term solutions to Equation (5) in the square torus diagram with the parameter values from Table 2 to show this visually.

We see in Figure 9A,B that over many cell cycle simulations, the phase space trajectories in both cluster 1 and cluster 2 scenarios will become phase-locked to a solution that is close to the main diagonal, corresponding to cell divisions occurring at a peak indefinitely. This will not happen with the cluster 3 case; Figure 9C shows quasiperiodic behaviour for this solution in which the trajectory will become arbitrarily close to every point on the torus as time tends to infinity. This would mean that cell division are not fixed to any particular Hes1 phase value. The results of this analysis suggest a qualitative difference in the solutions between cluster 3 (the long cell cycle case) and clusters 1 and 2 (the short cell cycle cases), as seen in Figure 9. It is very clear from the data in [26] that cluster 3 cells are particularly interesting due to the long cell cycle durations they possess. Our model predicts these longer cell cycles; the fact that our model now also predicts that phase-locking can not occur in these cells may be a contributing factor to the cell cycle elongation of these cells.

To end this investigation, we are interested to see what happens to solution trajectories when the oscillators maintain the coupling strength interactions but when Hes1 begins at different initial phases.

Figure 10 illustrates that when Hes1 begins at a certain phase, the bidirectional interaction strengths must be specific for that case for the data to be reproduced. For example, inputting the cluster 3 initial conditions into a simulation with the cluster 2 coupling parameters will yield a trajectory that corresponds to a Hes1 wave that is in anti-phase with cluster 1. We know from [26] that Hes1 does not have these specific mitosis to mitosis dynamics. Therefore, we conclude that it is the initial conditions and the specific bidirectional interaction that is required in each case to reproduce the results seen in [26].

4. Discussion

Our results show that the experimentally observed relationship between the expression dynamics of transcription factor (TF) Hes1 and the cell cycle can be described mathematically by a simple, deterministic model of coupled oscillators. Briefly, the experimental results that this work was based on were that, in a Human Breast Cancer cell line, the oscillatory protein expression of the Hes1 TF and the cell cycle fall into three clusters of related dynamics. Namely, one cell cycle from each cluster consists of high-low-high Hes1 dynamics in cluster 1, medium-low-high dynamics in cluster 2 and low-high-low-high dynamics in cluster 3. Cluster 1 cells divide at the peak of the Hes1 wave, which has been described as the preferential in-phase relationship of the majority of the cells, with a concomitant regular 24 h cell cycle length. By contrast, cluster 3 cells divide at the trough (“out-of-phase”); these include cells with significantly longer cell cycles. Cluster 2 cells exhibit an intermediate behaviour, dividing at the downward arm of Hes1 expression and having a 24 h cell cycle length.

Here, we have shown that when the two processes, Hes1 and the cell cycle, are represented only by their phase variables, the parameterised system of nonlinear ODEs from Equation (5) succeeds at generating phase space trajectories. These, in turn, predict the dynamic Hes1 expression behaviour found in the experimental data mentioned above. Furthermore, we have shown here that when the model is parameterised to the time-independent Hes1 data only, the number of time steps required to complete cell cycle simulations mimics the cell cycle duration patterns that are seen in the three clusters in human BCCs.

Our mathematical model shows that both Hes1 dynamics and cell cycle length information taken from real-world data are replicated if, and only if, the following allowances are made: (1) bidirectional coupling is allowed, (2) the symmetry of the system is allowed to be broken such that the coupling strength between cell cycle and Hes1 can differ in each direction and, furthermore, (3) the coupling strength parameter values are not maintained across each cluster case (i.e., they are different between clusters). The first prediction of bidirectional coupling is consistent with our previous experimental data [26]. However, the idea that the interaction strength between the oscillators may differ in either direction and that this happens on a cluster-dependent basis is unexpected and novel. Thus, this work offers a mathematical explanation for the complex outcome of the relationship between Hes1 and the cell cycle [26] and, furthermore, offers some novel parameterisation that will inform future work.

An intriguing product of the simulations is that the parameterisation of the model suggests that neither of the oscillators increases in their speed of progression. In fact, in the cluster 2 and 3 cases, one oscillator will slow down in order for the other to catch up to it. This is evidenced in cluster 2 by a plateau in Hes1 dynamics, yet cell cycle lengths remain similar to that of cluster 1. We may make some thought-provoking hypotheses from these findings: first, Hes1 dynamics cannot induce cell division, but it can delay it. Secondly, the cell cycle machinery can cause Hes1 protein abundance to plateau but not produce quicker oscillations. Future laboratory experimentation will be needed to test these hypotheses. A specific hypothesis to be tested is that of Hes1 being an integral controller, which has been found to be the case in other cell cycle systems such as the Whi5-Cln-cdk1 system [68]. Indeed this would work well with the result found in [26] of a positive correlation between cell cycle duration times and the shape of the Hes1 protein peak and, specifically, the “area under the curve” of Hes1 protein before its first peak was reached.

The benefits of using a simple mathematical model, such as Equation (5), is that it is intuitive, analytically manageable, easy to simulate computationally and yet still offers up surprising and worthwhile results. A more complicated model in the literature, focusing on neural progenitors instead of BCCs, takes into consideration many more interactive molecular species. Interestingly, although Hes1 operates on a much shorter, ultradian periodicity scale [49], it offers a similar conclusion, in part, suggesting that the behaviour of Hes1 is involved in the elongation of the G1-phase and hence the cell cycle. However, the model of [49] differs from ours in that it is unidirectional, meaning that no component of the large model feeds onto Hes1 expression offering no opportunity for Hes1 to be affected by the cell cycle. The interaction terms ${\kappa }_{\mathrm{Hes1}}$ and ${\kappa }_{\mathrm{C}.\mathrm{C}.}$, which we predict to be important here, may relate to the strength by which Hes1 represses cell cycle molecules, such as p21 and cyclin D [49], both of which are known molecular targets of Hes1. In the opposite direction, the cell cycle may be coupled to Hes1 via cell cycle-dependent phosphorylation of the Hes1 protein. Since Hes1 is a Notch signalling target, this could be due to variable Notch signalling during the cell cycle [69].

This is a point at which we may concede that the model in [49] has an advantage over ours as it contains many molecular components and is thus more detailed, allowing for closer inspection of the cell cycle’s inner workings. The ease of parameterisation of our functional model, on the other hand, is an advantage since the dynamics are uncomplicated to study. Going forward, this model could be easily adapted or modified to include more biological realism, e.g., including the addition of some external force acting upon one or both of the oscillators or offering some uncertainty quantification. Should our current hypothesis of the importance of changing bidirectionality be experimentally refuted, then it may be that the data is reproducible with the addition of the concept of a phase delay as is discussed in [70,71]. This may also be true with the introduction of oscillatory concepts such as dead zones (see [72]); this may be particularly fruitful since it is apparent that Hes1 is not responding to the cell cycle progression at the beginning of the cluster 2 trajectory in Figure 6. Further experimentation would be required to justify the implementation of these ideas.

This paper has been focused solely on Hes1’s relationship with one full cell cycle and has ignored the existence of previous or subsequent generations of cells. In future work, it would be interesting to see whether the coupling strength between the oscillators is capable of changing between generations. Changes over time within a cell lineage may align well with the fact that BCSCs are known to be plastic in nature, and interconversion between stem cell subtypes is known to occur [29,73]. The allocation of cells into Hes1 dynamic clusters also shows conversion between generations. Multigenerational analysis of Hes1 dynamics/cell cycle length will also be needed to elucidate the hypothesis of a “phase-skipping“ in Hes1, which may occur when the cell divides [26]. The concept of a phase-skipping event (also described as a hlphase-slip in [63]) of an oscillator due to stochasticity has been studied in [57], where it is stated that the “phenomenon occurs when a fluctuation in the non-deterministic elements exceeds the stability domain of the attractor”. Such an event with the mammalian circadian oscillator has been found to be linked to DNA damage in [74,75]. Phase skipping may also take place here since the periodicity of Hes1 studied here is circadian-like [26], and it would be interesting to examine the connection to DNA damage, especially since the data have been collected in MCF-7 cells, which are a model for breast cancer. In a cancer model, phase skipping and the associated elongation of the cell cycle may have important ramifications for the development and progression of the disease and could be a topic for future experimental and theoretical work on multigenerational gene expression of Hes1.

In summary, we have demonstrated in this theory paper that complex biological data, such as gene expression dynamics, can be modelled effectively using a straight-forward concept and very few core assumptions. By considering phases rather than protein abundances over time, we have been able to view the complicated landscape of genetic cell biology in a different light and strip back potentially superfluous detail to uncover key factors in the relationship between Hes1 and the cell cycle. This will allow us to propose new hypotheses (e.g., plasticity in cell cycle/Hes1 coupling strength between generations, phase-skipping of the Hes1 oscillator at a cell division) and promote exciting future work endeavours that would have otherwise been hidden.