A model for the proliferation–quiescence transition in human cells

: The process of revitalising quiescent cells in order for them to proliferate plays a pivotal role in the repair of worn-out tissues as well as for tissue homeostasis. This process is also crucial in the growth, development and well-being of higher multi-cellular organisms such as mammals. Deregulation of proliferation-quiescence transition is related to many diseases, such as cancer. Recent studies have revealed that this proliferation–quiescence process is regulated tightly by the Rb − E 2 F bistable switch mechanism. Based on experimental observations, in this study, we formulate a mathematical model to examine the effect of the growth factor concentration on the proliferation– quiescence transition in human cells. Working with a non-dimensionalised model, we prove the positivity, boundedness and uniqueness of solutions. To understand model solution behaviour close to bifurcation points, we carry out bifurcation analysis, which is further illustrated by the use of numerical bifurcation analysis, sensitivity analysis and numerical simulations. Indeed, bifurcation and numerical analysis of the model predicted a transition between bistable and stable states, which are dependent on the growth factor concentration parameter ( GF ). The derived predictions conﬁrm experimental observations .


Introduction
The human body consists of approximately 10 13 -10 14 cells. A large proportion of these cells are quiescent, a biochemically distinct state of growth arrest from which cells can re-enter the cell cycle [1]. The proportion of quiescent cells consists of cells that can no longer be re-activated to re-enter the cell cycle, and cells that can be re-activated to re-enter the cell cycle in response to growth factor signals under normal physiological conditions [2]. The cell cycle comprises four phases, namely: G 1 (first gap phase), S (synthesis phase), G 2 (second gap phase) and M (mitosis), where gap phases refer to the time interval between the synthesis phase and mitosis [3,4]. Here, mitosis refers to the process by which a single cell divides into two identical daughter cells. The commitment of a cell to cell division is driven by growth factor signals. The availability of sufficient growth factors just beyond a certain point between G 1 and S, known as the restriction point, leads to the progression of the cell cycle process, otherwise the cell enters a reversible state of growth arrest (G 0 ) [5][6][7][8][9][10][11][12].
Mammalian cellular homeostasis (the state of steady internal, physical, and chemical conditions maintained by cells) depends on the ability of cells to reversibly switch between quiescence and proliferation [13]. This mechanism is crucial for tissue repair and regeneration and is fundamental to the growth, development and well-being of mammals. The decision by cells to exit or enter quiescence is dysregulated in cancer and degenerative diseases [14][15][16]. Therefore, understanding the molecular mechanisms that control the reversible transition between quiescence and proliferation is crucial and remains a challenging problem in biological and medical sciences. In the 1970s and 1980s, three classes of models were proposed to describe the transition between cellular quiescence and proliferation [17]. These include: the "transition probability" models [18][19][20][21], the "growth controlled" models [22,23] and the hybrid models [24,25]. Hybrid models were developed by integrating transition probability-and growth-controlled models [24,25]. Though these models are coherent with diverse experimental data, they lack depth, since they are descriptive. Modelling proliferation-quiescence transition has shifted from this approach since the discovery, in the contexts of molecular and cell biology, of certain genes that regulate proliferation-quiescence transition [17]. Continuum mathematical models describing the temporal dynamics to provide insights into proliferation-quiescence transition within the cell cycle using gene regulatory networks have been proposed in the last two decades [26][27][28][29][30][31][32][33][34][35]. Limit cycles [28,30,32], cell-mass-regulated bi-stable systems [25,34,36], bi-stable and excitable systems [33,37] and transient processes [26,29,38] are examples of the cell-cycle dynamics deduced through cell-cycle modelling. It has been demonstrated experimentally that proliferation-quiescence transition is controlled by the Rb − E2F signalling network, which acts as a bi-stable switch [17]. Here, Rb and E2F represent the Retinoblastoma and Transcription factors respectively. Previously, models for the regulation of the Rb − E2F have been put forward and simulated [29,32,38], with more details provided in the model proposed by Aguda and Tang [38]. At the heart of the regulation of the Rb − E2F signalling network are: cyclin D(CycD), cyclin E(CycE), cyclin A(CycA), cyclin-dependent kinases (CDKs), family of transcription factors E2F, Myc and the retinoblastoma (Rb) family of proteins [39,40]. The retinoblastoma (Rb) protein family is responsible for regulating proliferation in most cells. The E2F family of transcription factors is responsible for the regulation of genes involved in DN A replication and cell cycle progression [41,42]. Interactions among these regulators have been outlined and verified experimentally [17]. However, due to the complex nature of the Rb − E2F network, which consists of criss-crossing linkages, a clear description of the proliferation-quiescence transition is elusive. In this study, we simplify the Rb − E2F signalling network that has been theoretically and experimentally verified in [43][44][45][46], and formulate an activator-inhibitor model system following the same philosophy as that proposed by Tyson et al. in [25], to analyse the proliferation-quiescence transition. To investigate the dynamics of the Rb − E2F signalling network describing proliferation-quiescence transition, we have taken the approach of combining redundant and overlapping cellular activities and collapsing linear cascades, as presented by Yao et al. [42], to simplify the model to three nodes connected by activation and inhibition links, as experimentally verified in [45]. Although the terms activation and inhibition are generally used in the literature, in this study we will also refer to these molecular processes as the production/synthesis and degradation/removal of abundant species, respectively. In addition, we consider the conservation of the total concentration of the Rb family of proteins and show through bifurcation and numerical analysis that the resulting system generates a set of three dynamical behaviours shown in experiments [1,2,17]; namely, stability, bistability and stability when the growth factor concentration (GF) is varied. Bifurcation and numerical analysis also show that quiescence is achieved under low growth factor stimulation and proliferation requires strong growth-factor stimulation. This paper is organised as follows: in Section 2, we utilise the current understanding of activator-inhibitor systems and apply biologically reasonable assumptions to motivate a system of non-linear ordinary differential equations (ODEs) and describe the pertinent interactions. In the same section, we non-dimensionalise the system, outline positivity, boundedness and prove the existence and uniqueness of solutions. This section ends with a further study of the non-linear ODE system to gain insights into its stability behaviour. In Section 3, we conduct a bifurcation analysis of the model to analyse the effect of varying the growth factor concentration on the steady states of the system. We also conduct numerical simulations to illustrate, computationally, the theoretical results. The last part of Section 3 presents sensitivity analysis for parameters of the system. We discuss the interpretations and implications of our results in Section 4. Finally, in Section 5, we discuss the limitations of our study and present possible future research problems.

Proliferation-Quiescence Dynamics Model
Our mathematical model is formulated based on first principles by first simplifying the Rb − E2F model studied by Yao et al. [17], collapsing all the several reaction networks into three main nodes denoted by: R, M and E and thus reducing the reaction links from 10 to 8, as shown in Figure 1. The R node consists of the Rb family of proteins (Rbp107 and Rbp108), the E node consists of the family of transcription factors E2F, Cyclin E and Cyclin A, and the M node consists of cyclin D and Myc. Experimentally, it is observed that R proteins are conserved throughout, while those of M and E are not [40]. Due to lack of experimental justifications for some reactions, and following other published works [12,40], in formulating our model we will employ Michaelis-Menten kinetics, mass action as well as the Hill function [47]. Michaelis-Menten kinetics describe the rate of enzymatic reactions by relating reaction rates to the concentration of a substrate, while the law of mass action states that the rate of a chemical reaction is directly proportional to the product of the activities or concentrations of the reactants [48,49].
Yao et al. [17] represented the three respective nodes by symbols: RP, EE and MD. These nodes are connected by 10 regulatory linkages. In their model, they assumed activation of the three species RP, EE and MD by Hill functions of the form: A n K n +[A n ] , where A represents any of the three species RP, MD or EE, K is a Michaelis-Menten constant and 1 ≤ n ≤ 10. In addition, they assumed that the inhibition of RP, EE and MD was through mass action kinetics. Of interest was their exclusion of the constitutive synthesis of EE together with its self-activation and the conservation of mass of the RP family of proteins, which we will consider in our model. Moreover, they considered self-inhibition of EE and inhibition of EE by MD, which we removed in our model. The simplified network self-reorganises to form an activator-inhibitor network, as shown in Figure 1.
In our model formulation, we ignore the spatial localisation of the proteins; such an extension forms part of our current studies and is beyond the scope of this work. We note, however, that such an approach leads naturally to partial differential equations. The resulting model, when spatial effects are neglected, is given in terms of a system of non-linear ordinary differential equations (ODEs) describing the rate of change in the concentrations of M(t), R(t) and E(t). The main novelty of the model is the use of the Hill function with n = 2 to describe kinetics with saturation. The model is formulated based on the following assumptions: 1 The production of M(t) is through mass action by extra-cellular growth signals GF and by E(t), which is modelled using Michaelis-Menten kinetics, whereas its decay is modelled by mass action. Yao et al. [2] considered Michaelis-Menten kinetics for the production of M(t) through extra-cellular growth signals. 2 The activation of R(t) is enhanced by E(t) and its inhibition is intensified by M(t) and E(t), which are both modelled using mass action. It is pertinent to note that, in this study, we assume conservation of mass for the Rb family of proteins, which was not considered in [2]. In addition, we assume self-activation and inhibition of R(t) using the Hill function with n = 2. On the contrary, Yao et. al. [2] modelled the production and depletion of R(t) using Michaelis-Menten kinetics and mass action, respectively, and self-activation and de-activation were not considered. 3 E(t) is synthesised, cf. [40], with the use of Michaelis-Menten kinetics and its synthesis is enhanced by M(t) with the use of a Hill function with n = 2, while its decay is enhanced by R(t) with the use of mass action. On the contrary, the authors in [2] did not consider constitutive synthesis of E(t), which has been observed experimentally as indicated in [40]. The main novelty of our proposed model is that we consider mass conservation of the Rb family of proteins represented by R(t) existing in both hypo-phosphorylated form denoted by R u (t) and hyper-phosphorylated form denoted by R p (t); whereas in [42] it was assumed that the concentration of the Rb family of proteins was abundant and hence there was no mass conservation. This assumption makes our model more realistic, as the transition from cellular quiescence to proliferation is dependent on the cycling between the hyper-phosphorylated and hypo-phosphorylated forms of Rb family of proteins mediated by cyclin-dependent kinases. Hyperphosphorylation of Rb enhances proliferation, whereas hypo-phosphorylated Rb enhances quiescence [50]. In addition, some reaction pathways were not considered in [2], including constitutive synthesis of E and the simple mass action for the activation of R(t). Therefore, employing mass conservation, the total concentration of R(t) is such that: The state variable R(t) denotes the concentration of R p (t) at time t in the model. R u (t) corresponds to the inactive term R T − R used in Equation (1c). Therefore, the time evolution of M(t), E(t) and R(t) are, respectively, described by the following system of non-linear ordinary differential equations (ODEs): The system is closed with non-negative initial conditions M(0) = M 0 , E(0) = E 0 , and R(0) = R 0 , 0 ≤ M(t), E(t) and 0 ≤ R(t) ≤ R T . In our model, we consider M and E families of proteins to exist in abundance; hence, there is no conservation of mass, unlike the Rb family of proteins, which is conserved as outlined in the model formulation. This implies, that for some appropriate model parameter values, M and E might be unbounded. Model parameters are chosen based on a previous modelling study on the Rb − E2F pathway [17]. All model parameters are strictly positive and their values, units and physical meaning are described in Table 1.

Mathematical Analysis of the Model
In this section, we non-dimensionalise system (1) so that it is unit free for purposes of mathematical analysis, including proving positivity, boundedness, existence and uniqueness of solutions as well as performing the sensitivity analysis. However, numerical bifurcation and numerical simulations were performed on the dimensional model. This is meant to preserve the original concentration of species GF, which appears in Equation (1a).

Non-Dimensionalisation
There are nineteen parameters in system (1). We can reduce the number of parameters to sixteen by defining new coordinates M = Xx * , E = Yy * , R = Zz * and t = Tτ, where X, Y, Z and t are constant (dimension carrying) scales, to be chosen and x * ; y * ; z * ; τ are the numerical (dimensionless) variables. Substituting these into (1) and letting and a 16 = , and, dropping * so that x * = x, y * = y and z * = z, we obtain the non-dimensional system: with positive initial conditions x(0) = x 0 , y(0) = y 0 and z(0) = z 0 and positive constant parameters a 1 , . . . , a 16 .

Positivity, Boundedness and Existence and Uniqueness of Solutions
We first show thatsystem (2), associated with initial conditions, has a unique local solution. (2), associated with the initial condition (x(0), y(0), andz(0)), has a unique local solution in the interval (0, T) for some T > 0.

Lemma 1 (Local existence and uniqueness). System
Proof. System (2) can be written in matrix form, as follows: It follows that f is Lipschitz continuous, i.e., there exists a constant L ≥ 0 such that for any bounded region D in R 3 + . Then, local existence and uniqueness of solutions is established by the classic theory, cf. [51].
Next, the positivity of solutions of (2) is established; hence, the feasibility of solutions for the study of the underlying biological problem is guaranteed.
Proof. We must prove that (x(t), y(t), z(t)) will remain in R 3 . Thus, recalling that all parameters used in system (2) are positive, we derive the following inequalities in the interval (0, t 1 ) Now, integrating the above differential equation inequalities in the interval [t 1 , T), we obtain for some positive constants A 0 , B 0 , C 0 . Thus, for all t ∈ [0, t 1 ], x(t), y(t), and z(t) will remain in R 3 + . Repeating this argument, we can prove that (x(t), y(t), z(t)) ∈ R 3 + for any t ∈ (0, T), and this completes the proof.

Lemma 3 (Boundedness).
There exists x v , y v , z v > 0 such that lim sup t→T x(t) ≤ x v , lim sup t→T y(t) ≤ y v and lim sup t→T z(t) ≤ z v provided that min(θ, φ) > 0 where θ := a 4 − a 2 and φ := a 15 z v − a 13 . Here, T is the maximum existence time for system (2) given by Lemma 1.
Proof. Since all involved constants a i , i = 1, . . . , 16 are positive, thanks to Lemma 2 we can deduce that By making the substitution z = z v we also obtain that Setting θ := a 4 − a 2 and φ := a 15 z v − a 13 , provided that min(θ, φ) > 0, then it follows that where λ := a 1 + a 16 + 1. Solving this differential inequality using the method of integrating factor, we obtain: and thus lim sup entails that (x + y)(t) is bounded from above and, therefore, x(t), andy(t) are also bounded from above by
Proof. The existence of a local positive solution for system (2) is an immediate consequence of Lemmas 1 and 2. Thanks to Lemma 2, such a solution is also bounded in (0, T) and, thus, the classical theory, cf. [51] (Corollary 2.3.2), guarantees its existence as a (global) positive solution for all t > 0.
Next, we demonstrate the existence of at least three uniform states for the model system (1).

Steady States
The uniform steady state (M * , E * , R * ) of System (1) is a solution to the following system of nonlinear algebraic equations Analytically, solving this set of equations gives a cubic polynomial which, by the Fundamental theorem of algebra [52], yields at most three steady states. The illustration of this is shown by plotting numerical bifurcation analysis results shown in Figure 2. Therefore, for some parameter values, system (1) admits up to three steady states, if they exist. We remark that finding exact analytical solutions of the uniform steady states was not possible due to the intractability of the system. Instead, we used bifurcation analysis to unravel the existence of such states.

Linear Stability Analysis
The following analysis establishes the stability of the positive steady states for ODE system (1). To investigate the stability of (1), we linearize it around the steady state (M * , E * , R * ) to obtain: with and, thus, recalling that, due to Theorem 1. it also holds that (M * , E * , R * ) ∈ R 3 + .
Proof. The Jacobian matrix for the system is given by Next, we compute the eigenvalues of J by solving the equation where I is the 3 × 3 identity matrix. The characteristic polynomial of J is then given by with coefficients, By the Routh-Hurwitz stability criterion, all the roots of (12) have negative real parts if p 1 > 0, p 3 > 0 and p 1 p 2 − p 3 > 0. Clearly p 1 > 0 by condition (A). In addition, and, thus, by conditions (A) and (B) it follows that p 3 > 0. It remains to show that p 1 p 2 − p 3 > 0. Now, After simplifying and grouping, we obtain Due to conditions (A)-(C) and in conjunction with the sign of the parameter θ 3 < 0, then p 1 p 2 − p 3 > 0. Hence, the system of ODEs (1) becomes asymptotically stable.

Bifurcation Analysis
Bifurcation analysis of the system of ODEs (1) was carried out using XPPAUT, a software program freely available from http://www.math.pitt.edu/$\sim$bard/xpp/xpp. html (accessed on 22 November 2020) [53]. Parameter values for our model are shown in Table 1, unless otherwise specified. The parameter values were chosen based on the literature whenever possible and some of them adjusted to illustrate qualitative dynamics [17,42]. As shown in Figure 2a, there exists a unique steady state corresponding to the quiescent state when GF < 0.3735 with M * < 1.225. For 0.3735 < GF < 0.4138, bistable dynamics result from a saddle node bifurcation marked by two black dots labelled as SN (saddle node). As GF increases further out of the bistable regime, M * jumps into its high steady state, where a cell undergoes proliferation. We also plot k 1 against GF and generate a two parameter bifurcation diagram shown in Figure 2c. In the latter figure, the bistable region is coloured green.

Numerical Simulations
Numerical simulations were carried out using the MATLAB ode45 solver [54]. This is a standard solver for ordinary differential equations (ODEs) that implements a Runge-Kutta method with a variable time-step for efficient computations, (see [55] for details). We simulated the system by selecting values of GF corresponding to different dynamic regimes deduced through the bifurcation analysis in the previous subsection. We solved the system of ODEs (1) for those three regimes and plot the solution for separate time intervals with a variable growth factor concentration value and for different initial conditions. First, we show the time evolution for M(t), R(t) and E(t) by solving numerically system (1) with parameters as outlined in Table 1 assuming the growth factor has concentration value GF = 0.3 with initial conditions M 0 = 0.2, E 0 = 1.2 and R 0 = 3, which bring the system to the lower stable region, as shown in Figure 3a. Second, we solved the system with the growth factor concentration GF value of 0.39 and two sets of initial conditions M 0 = 0.2, E 0 = 1.2, R 0 = 3 and M 0 = 20, E 0 = 5, R 0 = 5 while keeping all the other parameters fixed. The results obtained show that the system exhibits bistability, as shown in Figure 3b. Third, we solved the model after adjusting the growth factor concentration GF to 0.45 and, keeping all other parameters unchanged and initial conditions M 0 = 20, E 0 = 5 and R 0 = 5, the model evolves to a higher steady state, as shown in Figure 3c. These results confirm the conditions set in Theorem 2.
Next, we explored conditions given in Lemma 3 and Theorem 2, as demonstrated in Figure 4. Figure 4 illustrates numerically the dynamics of model system (1) when both Lemma 3 and Theorem 2 are violated. When both conditions of Lemma 3 and Theorem 2 are not satisfied, then there exist two steady states, one stable and another unstable, and, henceforth, the possibility of two different dynamics: either bounded (stable) or unbounded (unstable) solutions. The unstable steady state acts as a switch, determining the initial conditions under which the system may exhibit bounded or unbounded behaviour.  To illustrate such dynamics, we set α 5 = 0.0079, δ = 1 and GF = 0.3, while keeping the rest of the parameters fixed, as shown in Table 1. This choice of parameters ensures that conditions (B) and (C) outlined in Theorem 2 are violated while, at the same time, the min(θ, φ) = 0, which violates Lemma 3. With this choice of parameters, model system (1) exhibits two steady states, as predicted above, an unstable and a stable steady state, which leads to either unbounded (unstable) or bounded (stable) solutions, respectively. Figure 4d demonstrates the validity of our theoretical results; the system of ODEs (1) becomes unstable with solutions for M(t) and E(t) growing unboundedly while those of R(t) remain bounded. The biological justification for this type of model behaviour lies in our assumptions. Unlike the Rb family of proteins, which are assumed to be conserved for all time, M(t) and E(t) species are not conserved, because they are assumed to exist in abundance. Therefore, they are not constrained by laws of mass conservation. We also plot the bifurcation diagram as shown in Figure 4a, for which we have the lower stable steady state shown in red (M * , E * , R * ) = (0.7255, 1.4184, 2.0449) with α 5 = 0.0079) , followed by an upper unstable steady state also for the same parameter value α 5 = 0.0079 (shown in grey) (M * , E * , R * ) = (0.9755, 2.25168, 2.08855). For completeness, we have included the bifurcation diagrams for E(t) and R(t), respectively, as illustrated in Figure 4b,c. The eigenvalues of the above steady states are given by (−3.7119, −0.4209, −0.0015) and (−5.6216, −0.3165, 0.0013), respectively, confirming that the lower steady state is stable while the upper is unstable. The unstable steady state acts as a switch mechanism between the unbounded and bounded solutions, determining the initial conditions under which the system may either exhibit unbounded or bounded solutions when conditions in Lemma 3 and Theorem 2 are violated.
We complete this section by illustrating stable dynamics of model system (1), as shown in Figure 4f. To obtain bounded solutions close to the unstable steady state, we adjusted our initial conditions to M 0 = 0.4782, E 0 = 0.594 and R 0 = 1.8616, such that these are located below the unstable steady state. The model system solved with this set of initial conditions exhibits bounded solutions which converge to the stable steady state computed above.

Sensitivity Analysis
Following the works in [56][57][58], we present the local sensitivity analysis of system (1) with respect to baseline parameter values in Table 1. We used local sensitivity analysis to identify the set of model parameters whose percentage change from the baseline parameter set causes significant changes in the model output [59]. Hence, local sensitivity analysis provides a useful insight into identifying model components that contribute most to the variability in the model output [59,60]. In general, there are two types of sensitivity analyses: local and global. In this article, we focused on local sensitivity analysis, which will allow us to evaluate changes in the model output with respect to variations in a single parameter at a time [59,61]. Global sensitivity analysis will allow for simultaneous changes in model parameters and evaluation of the interaction between parameters, which is not in our interest at this stage. Our methodology involves increasing parameters one at a time by 5%, and 10% and decreasing them by 5% and 10%, respectively. Next, we calculated the local sensitivity indices (percentage changes) in the solutions for M(t), E(t) and R(t).
The results for 5% increase, 5% decrease, 10% increase and 10% decrease are shown in Tables 2-5, respectively. We proceed to illustrate comparatively the results corresponding to 5% increase and decrease in Figure 5 as well as 10% increase and decrease in Figure 6 for each variable. First, we observe that the steady state M * is mostly sensitive to parameters k 1 and GF, as shown in Figure 5a,b. Clearly, k 1 and GF have negative and positive effects to M * , respectively, corresponding to 5% increase and decrease in parameters. This means that increasing k 1 reduces M * , whereas increasing GF increases M * . Second, the steady state E * is mostly sensitive to α 2 and α 5 followed by k 1 , k r 1 and GF, as shown in Figure 5c,d. Third, k r 3 and α 2 have the greatest effect on R * followed by α 4 and α 6 , as shown in Figure 5c,d. It is also pertinent to note that the effect of 5% parameter increases on the steady states of M(t), E(t) and R(t) has a direct opposite effect to that of reducing the parameters by the same margin, as shown in Figure 5a-f. Furthermore, the response of M(t), E(t) and R(t) steady-state values to a 5% increase in parameters is almost similar, as shown in Figure 5a,c.
In addition, the response of the M(t), E(t) and R(t) steady-state values to a 5% decrease in parameters is also similar, as shown in Figure 5b,d, albeit with different magnitudes.
This can interpreted biologically as follows: increasing (i) the growth factor concentration GF, (ii) the activation rate of M(t), and (iii) the inhibition rate of E(t) by R(t) will result in an increase in the steady-state values of M(t) and E(t), leading to proliferation. Decreasing these parameters will reduce their respective steady states, leading cells to quiescence. In Figure 5e, all the parameters except α 2 , α 5 , β 2 , β 3 , k r 1 , and k r 2 have an opposite effect on the steady state of R(t) when compared to those of M(t)and E(t), as shown in Figure 5a,c. Thus, increasing these parameters would have the effect of suppressing both M * and E * , as well as increasing R * .   Next, we proceed to investigate the correlations between parameter sensitivities to parameter changes for M against E, M against R and E against R, respectively, by considering a 5% increase in parameters. The results are shown in Figure 7. From Figure  7a, it is clear that all parameters except k 1 , δ, β 1 and GF affect the steady states M * and E * in the same way, since they lie on the same straight line. The parameters δ and k 1 have a negative effect on M * , E * and a positive effect on R * , i.e., an opposite effect on R * . Parameter GF has a positive effect on both M * and E * and an opposing effect on R * . In Figure 7b,c, there is almost no correlation between the changes in the steady-state values M * against R * as well as E * against R * .
To proceed, we interpret the correlations shown in Figure 7. Generally, the order of the impact of parameters on M(t) and E(t) between a 5% and 10% increase in parameters is found to be similar, see Tables 2 and 4. This fact is also true for 5% and 10% decreases in parameters, as shown in Tables 3 and 5. This confirmed that none of the parameters of the system has a switching effect; that is, when they reach a particular threshold, their effect on the steady state is greatly increased, in relation to the other parameters. Having identified this, it is inferred that the dynamics of proliferation-quiescence transition remained unaltered with respect to small perturbations in parameters. In Figure 7a, most of the parameters are distributed linearly; thus, these parameters affect both M(t) and E(t) solutions in a similar fashion. Both Figure 7b,c show a similar pattern, in which parameter values are grouped in two well correlated sets. The first set contains parameters α 2 , k r 1 , k r 2 , β 3 , β 2 and α 5 , while the other contains k r 3 , α 6 , β 4 , α 3 , R T , α 4 and α 8 . The effects of changes in α 3 are absolutely negligible, as shown in Figures 5 and 6. Biologically, α 3 represents the inhibition rate of R by E which is a two-stage process involving the dephosphorylation of the Rb family of proteins [40] followed by the inhibition of these proteins. This effectively means that the inhibition process depends on the dephosphorylation of the Rb family of proteins, thus explaining the negligible sensitivity observed in this case. The parameter δ lies off the linear pattern, as shown in Figure 7a, which is a result of its role in the proliferation-quiescence transition. It is known that variation in cyclin D levels through the cell cycle determines the proliferating fate of a cell [40]. Parameters δ, k 1 , α 1 and GF lie off the well-defined linear pattern followed by the rest of the parameters. From this observation, it is inferred that a 10% increase in the inhibition rate of M (Cyclin D and Myc) causes a significant decrease in the activation of E(t) (E2F, Cyclin E, Cyclin A) and a small increase in the activation of the Rb family of proteins. Under these circumstances, cells would naturally enter the quiescence state. Additionally, an increase in the expression of the Rb family of proteins significantly reduces the expression levels of E (E2F, Cyclin E and Cyclin A), thereby preventing cells from entering the cell cycle.    . (a,b) show sensitivity of M * after 10 % increase and decrease in parameters, respectively, (c,d) show sensitivity of E * corresponding to 10% increases and decrease, respectively. (e,f) show sensitivity R * to 10% increases and decrease in parameters, respectively.

Discussion
Biologically, it was shown that the growth factor concentration plays a pivotal role in a wide range of cellular processes, including cellular growth and differentiation [40]. As such, growth factor (GF) stimulation has a profound effect on cancer development [18,19,21]. Previous studies have shown that the Rb − E2F bistable switch converts graded and transient serum growth signals into a binary and hysteric E2F activity in individual cells [2].
In this paper • Based on the concept of first principles, we investigated the dynamical potential of growth factors in the regulation of the cell-cycle entry. • A mathematical model for the simplified Rb − E2F network was constructed based on the model proposed in Yao et al. [2]. While previous studies modelled all links using Michaelis-Menten functions only [2], we used mass action, and Michaelis-Menten and Hill functions, resulting in a simpler model. In addition, we considered the R species to exist either in hyper-phosphorylated or hypo-phosphorylated form and that their total concentration is conserved. • By varying the growth factor signal values through bifurcation analysis, numerical simulations illustrated that the magnitude of the value of the growth factor plays a critical role in regulating cell-cycle entry. Through bifurcation analysis, we deduced the existence of three consecutive dynamical behaviours, namely, stability, bi-stability and stability.
• Numerical simulations performed with different growth factor values validated the results derived from the bifurcation analysis. In particular, the biological interpretation of the uniform steady state can be established as follows: 1.
For GF < 0.3735, System (1) is asymptotically stable, indicating the regime in which cells are in a quiescent state. In this state, cells feature low levels of Cyclin D, Myc and high levels of R species.

2.
On the other hand, in the range 0.3735 < GF < 0.4138, System (1) exhibits bistability, marking the position of the restriction point, as deduced in Yao et al. [42]. This point sets a high threshold separating quiescence from proliferation and acts as a barrier against unregulated and accidental cell growth. In addition, it provides a low-maintenance mechanism ensuring that the cell cycle proceeds, albeit later due to changes in the extracellular environment which is crucial for maintaining genome integrity.

3.
For values of GF > 0.4138, the system generates a stable dynamical behaviour where a cell is in the proliferation mode marking the higher steady state value. This state features high levels of Cyclin D, Myc and low levels of the Rb family of proteins.
However, it remains to investigate the conditions under which the system exhibits excitable and oscillatory dynamics as observed in a different model proposed in [62], but that would be investigated in a subsequent work. While Yao et al. [42] identified a basic gene circuit underlying resettable Rb − E2F bi-stable switch controlling cell-cycle entry, we obtained a range of values of the growth factor concentration for the three dynamical regimes.
In this study, we focused our attention on the quantitative aspects of bi-stability, including ascertaining some constraints and the region for bi-stability, whereas, in [42], the focus was on the qualitative aspects of bi-stability.
Our consideration of the conservation of the R species did not affect the bistable nature of the system but revealed that the system becomes unstable at very low levels of the concentration of the R species. Local sensitivity analysis revealed that increasing parameters that enhance the availability of Cyclin A, Cyclin E, Cyclin D, Myc and E2F in the model results in the hyper-phosphorylation of Rb proteins. In its hyper-phosphorylated state, the Rb loses its affinity to bind free E2F, which results in transcription, leading to proliferation [40]. On the contrary, reducing these parameters by a small margin results in Rb dephosphorylation, which increases its affinity for E2F. Naturally, in its hypo-phosphorylated state, Rb binds to free E2F, impairing transcription. This results in cells being unable to pass through the restriction point and, hence, remaining in quiescence.

Limitations
In this study, we constructed an activator-inhibitor model to describe the Rb − E2F signalling interaction network and analysed its dynamics and biological implications. However, in our study there are some limitations. Firstly, the model does not seem to exhibit periodicity and excitable dynamics. Secondly, spatial localisation of the proteins was completely ignored. This will give rise to partial differential equations. Thirdly, we did not fit the model to data to infer model parameters; this forms part of our future studies. Though simplification of the model produced the interesting bistable behaviour consistent with experimental observations, there were consequences on other dynamics of the Rb − E2F signalling network such as oscillatory dynamics, as observed in [62]. Moreover, we assumed that all proteins other than the Rb family of proteins exist in abundance which may not be the case and hence must be investigated in future studies.

Data Availability Statement:
The authors confirm that the data supporting the findings of this study are available within the article.