Next Article in Journal
State Machines and Hypergroups
Previous Article in Journal
A Fault Localization Method Based on Metrics Combination
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

A Model for the Proliferation–Quiescence Transition in Human Cells

by
Kudzanayi Z. Mapfumo
1,†,
Jane C. Pagan’a
2,†,
Victor Ogesa Juma
3,†,
Nikos I. Kavallaris
4,† and
Anotida Madzvamuse
5,6,7,*,‡
1
Department of Mathematics and Computational Sciences, University of Zimbabwe, Harare P.O. Box MP167, Zimbabwe
2
Department of Statistics and Mathematics, Bindura University of Science Education, Bindura P.O. Box 1020, Zimbabwe
3
Mechanical Engineering Department, University of Zaragoza, Edificio Betancourt, Campus Rio Ebro, E-50018 Zaragoza, Spain
4
Department of Mathematics and Computer Science, Faculty of Health, Science and Technology, Karlstad University, 651 88 Karlstad, Sweden
5
Department of Mathematics, University of Sussex, Pevensey III, Brighton BN1 9QH, UK
6
Department of Mathematics, University of Johannesburg, Johannesburg 2006, South Africa
7
Department of Mathematics, University of British Columbia, 1984 Mathematics Road, Vancouver, BC V6T 1Z2, Canada
*
Author to whom correspondence should be addressed.
These authors contributed equally to this work.
Conceived the research study.
Mathematics 2022, 10(14), 2426; https://doi.org/10.3390/math10142426
Submission received: 23 May 2022 / Revised: 4 July 2022 / Accepted: 8 July 2022 / Published: 12 July 2022
(This article belongs to the Section Difference and Differential Equations)

Abstract

:
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 R b 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 ( G F ). The derived predictions confirm experimental observations.

1. 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 R b E 2 F signalling network, which acts as a bi-stable switch [17]. Here, R b and E 2 F represent the Retinoblastoma and Transcription factors respectively. Previously, models for the regulation of the R b E 2 F 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 R b E 2 F signalling network are: cyclin D ( C y c D ) , cyclin E ( C y c E ) , cyclin A ( C y c A ) , cyclin-dependent kinases ( C D K s ) , family of transcription factors E 2 F , M y c and the retinoblastoma ( R b ) family of proteins [39,40]. The retinoblastoma ( R b ) protein family is responsible for regulating proliferation in most cells. The E 2 F family of transcription factors is responsible for the regulation of genes involved in D N 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 R b E 2 F network, which consists of criss-crossing linkages, a clear description of the proliferation–quiescence transition is elusive. In this study, we simplify the R b E 2 F 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 R b E 2 F 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 R b 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 ( G F ) 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.

2. Materials and Methods

2.1. Proliferation–Quiescence Dynamics Model

Our mathematical model is formulated based on first principles by first simplifying the R b E 2 F 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 R b family of proteins ( R b p 107   and   R b p 108 ) , the E node consists of the family of transcription factors E 2 F , Cyclin E and Cyclin A , and the M node consists of cyclin D and M y c . 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: R P , E E and M D . These nodes are connected by 10 regulatory linkages. In their model, they assumed activation of the three species R P , E E and M D by Hill functions of the form: A n K n + [ A n ] , where A represents any of the three species R P , M D or E E , K is a Michaelis–Menten constant and 1 n 10 . In addition, they assumed that the inhibition of R P , E E and M D was through mass action kinetics. Of interest was their exclusion of the constitutive synthesis of E E together with its self-activation and the conservation of mass of the R P family of proteins, which we will consider in our model. Moreover, they considered self-inhibition of E E and inhibition of E E by M D , 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 G F 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 R b 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 R b 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 R b 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 R b family of proteins mediated by cyclin-dependent kinases. Hyperphosphorylation of R b enhances proliferation, whereas hypo-phosphorylated R b 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:
R u ( t ) + R p ( t ) = R T = C o n s t a n t .
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):
d M d t = α 1 G F p r o d u c t i o n   o f   M   b y   G F + β 1 M E k 1 + E p r o d u c t i o n   o f   M   b y   E δ M d e c a y   o f M ,
d E d t = α 6 E   b a s e l i n e   v a l u e + α 4 E α 8 + E E   c o n s t i t u t i v e s y n t h e s i s + β 4 E M 2 k r 3 2 + M 2 P r o d u c t i o n   o f   E   b y   M α 5 R E r e m o v a l   o f   E   b y   R ,
d R d t = β 2 ( R T R ) E A c t i v a t i o n   o f   R   b y   E + β 3 ( R T R ) 2 k r 1 2 + ( R T R ) 2 R   b a s e l i n e   a c t i v a t i o n α 2 R M I n h i b i t i o n   o f   R   b y   M α 3 R E i n h i b i t i o n   o f   R   b y   E β 5 R 2 k r 2 2 + R 2 R   b a s e l i n e   i n h i b i t i o n .
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 R b 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 R b E 2 F pathway [17]. All model parameters are strictly positive and their values, units and physical meaning are described in Table 1.

2.2. 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 G F , which appears in Equation (1a).

2.3. Non-Dimensionalisation

There are nineteen parameters in system (1). We can reduce the number of parameters to sixteen by defining new coordinates M = X x * , E = Y y * , R = Z z * 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 a 1 = α 1 G F α 4 , a 2 = β 1 α 8 α 4 , a 3 = k 1 α 8 , a 4 = δ α 8 α 4 , a 5 = β 2 α 8 2 α 4 , a 6 = R T α 8 , a 7 = β 3 α 4 , a 8 = k r 1 2 α 8 2 , a 9 = α 2 α 8 2 α 4 , a 10 = α 3 α 8 2 α 4 , a 11 = β 5 α 4 , a 12 = k r 2 2 α 8 2 , a 13 = β 4 α 8 α 4 , a 14 = k r 3 2 α 8 2 , a 15 = α 6 α 8 α 4 and a 16 = α 5 α 8 2 α 4 , and, dropping * so that x * = x , y * = y and z * = z , we obtain the non-dimensional system:
d x d t = a 1 + a 2 x y a 3 + y a 4 x ,
d y d t = a 16 + y 1 + y + a 13 x 2 y a 14 + x 2 a 15 y z ,
d z d t = a 5 ( a 6 z ) y + a 7 ( a 6 z ) 2 a 8 + ( a 6 z ) 2 a 9 x z a 10 y z a 11 z 2 a 12 + z 2 ,
with positive initial conditions x ( 0 ) = x 0 ,   y ( 0 ) = y 0   a n d   z ( 0 ) = z 0 and positive constant parameters a 1 , , a 16 .

2.4. Positivity, Boundedness and Existence and Uniqueness of Solutions

We first show thatsystem (2), associated with initial conditions, has a unique local solution.
Lemma 1
(Local existence and uniqueness). System (2), associated with the initial condition ( x ( 0 ) , y ( 0 ) , a n d z ( 0 ) ) , has a unique local solution in the interval ( 0 , T ) for some T > 0 .
Proof. 
System (2) can be written in matrix form, as follows:
d x d t = d x d t d y d t d z d t = f ( x ) = f 1 ( x , y , z ) f 2 ( x , y , z ) f 3 ( x , y , z ) ,
for x = x ( t ) = x ( t ) , y ( t ) , z ( t ) T and
f 1 ( x , y , z ) : = a 1 + a 2 x y a 3 + y a 4 x , f 2 ( x , y , z ) : = a 16 + y 1 + y + a 13 x 2 y a 14 + x 2 a 15 y z , f 3 ( x , y , z ) : = a 5 ( a 6 z ) y + a 7 ( a 6 z ) 2 a 8 + ( a 6 z ) 2 a 9 x z a 10 y z a 11 z 2 a 12 + z 2 .
It follows that f is Lipschitz continuous, i.e., there exists a constant L 0 such that f ( x 1 ) f ( x 2 ) L x 1 x 2 , for all x 1 , x 2 D , 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.
Lemma 2
(Positivity). If the initial condition in system (2) satisfies x ( 0 ) > 0 , y ( 0 ) > 0 , z ( 0 ) > 0 , then the solution ( x ( t ) ,   y ( t ) ,   z ( t ) ) of (2) will remain in R + 3 for all t [ 0 , T ] , where T > 0 stands for the maximum existence time guaranteed by Lemma 1.
Proof. 
We must prove that ( x ( t ) , y ( t ) , z ( t ) ) will remain in R + 3 for any t ( 0 , T ) , where T > 0 is the maximum existence time. Since x ( 0 ) > 0 , y ( 0 ) > 0 , z ( 0 ) > 0 then, by a continuity argument, there exists t 1 > 0 such that ( x ( t ) , y ( t ) , z ( t ) ) R + 3 for all t [ 0 , t 1 ] . Thus, recalling that all parameters used in system (2) are positive, we derive the following inequalities in the interval ( 0 , t 1 )
d x d t = a 1 + a 2 x y a 3 + y a 4 x a 4 x ,
d y d t = a 16 + y 1 + y + a 13 x 2 y a 14 + x 2 a 15 y z a 15 y z ,
d z d t = a 5 ( a 6 z ) y + a 7 ( a 6 z ) 2 a 8 + ( a 6 z ) 2 a 9 x z a 10 y z a 11 z 2 a 12 + z 2 a 5 z a 9 x z a 10 y z .
Now, integrating the above differential equation inequalities in the interval [ t 1 , T ) , we obtain
x ( t ) A 0 e a 4 ( t t 1 ) > 0 , y ( t ) C 0 e a 15 t 1 t z ( t ) d t > 0 , z ( t ) B 0 e a 9 t 1 t x ( t ) d t ( a 5 + a 10 ) t 1 t y ( t ) d t > 0 ,
for some positive constants A 0 , B 0 , C 0 . Thus, for all t [ 0 , t 1 ] ,   x ( t ) ,   y ( t ) , a n d   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 m i n ( θ , ϕ ) > 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
lim z a 6 d z d t = a 6 a 9 x a 6 a 10 y a 11 a 6 2 a 12 + a 6 2 0 ,
which implies that there exists z v > 0 such that lim sup t T z ( t ) z v ; hence, z ( t ) is bounded in [ 0 , T ] . By making the substitution z = z v we also obtain that
d x d t + d y d t = a 1 + a 2 x y a 3 + y a 4 x + a 16 + y y + 1 + a 13 x 2 y a 14 + x 2 a 15 y z v . a 1 + a 2 x a 4 x + a 16 + 1 + a 13 y a 15 y z v , a 1 + a 16 + 1 x ( a 4 a 2 ) y ( a 15 z v a 13 ) .
Setting θ : = a 4 a 2 and ϕ : = a 15 z v a 13 , provided that m i n ( θ , ϕ ) > 0 , then it follows that
d ( x + y ) d t λ θ x ϕ y λ m i n ( θ , ϕ ) ( x + y ) ,
where λ : = a 1 + a 16 + 1 . Solving this differential inequality using the method of integrating factor, we obtain:
( x + y ) ( t ) λ m i n ( θ , ϕ ) + c 0 e m i n ( θ , ϕ ) t λ m i n ( θ , ϕ ) + c 0 ,
and thus
lim sup t T ( x + y ) ( t ) λ m i n ( θ , ϕ ) + c 0 .
Now, choosing
x v = y v = λ m i n ( θ , ϕ ) + c 0 ;
entails that ( x + y ) ( t ) is bounded from above and, therefore, x ( t ) , a n d y ( t ) are also bounded from above by x v = y v for all t [ 0 , T ] . □
Theorem 1.
Under the assumption of Lemma 3 and provided that x ( 0 ) ,   y ( 0 ) ,   z ( 0 ) > 0 , then system (2) has a global-in-time positive solution; that is, ( x ( t ) , y ( t ) ,   z ( t ) )   w i l l   e x i s t   i n   R + 3 and is bounded for any t > 0 .
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).

2.5. Steady States

The uniform steady state ( M * , E * , R * ) of System (1) is a solution to the following system of nonlinear algebraic equations
α 1 G F + β 1 M E k 1 + E δ M = 0 ,
α 6 + α 4 E α 8 + E + β 4 E M 2 k r 3 2 + M 2 α 5 R E = 0 ,
β 2 ( R T R ) E + β 3 ( R T R ) 2 k r 1 2 + ( R T R ) 2 α 2 R M α 3 R E β 5 R 2 k r 2 2 + R 2 = 0 .
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.

2.6. 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:
d ( M M * ) d t d ( E E * ) d t d ( R R * ) d t = J M M * E E * R R * ,
where
J = g 1 M g 1 E g 1 R g 2 M g 2 E g 2 R g 3 M g 3 E g 3 R ( M * , E * , R * ) .
with
g 1 ( M , E , R ) = α 1 G F + β 1 M E k 1 + E δ M , g 2 ( M , E , R ) = α 6 + α 4 E α 8 + E + β 4 E M 2 k 3 2 + M 2 α 5 R E , g 3 ( M , E , R ) = β 2 ( R T R ) E + β 3 ( R T R ) 2 k r 1 2 + ( R T R ) 2 α 2 R M α 3 R E β 5 R 2 k r 2 2 + R 2 ,
and, thus,
g 1 M | ( M * , E * , R * ) = δ + β 1 E * k 1 + E * : = ω 1 , g 1 E | ( M * , E * , R * ) = k 1 β 1 M * ( k 1 + E * ) 2 : = ω 2 > 0 , g 1 R | ( M * , E * , R * ) = 0 , g 2 M | ( M * , E * , R * ) = 2 M * β 4 E * k 3 2 ( k 3 2 + M * 2 ) 2 : = γ 1 > 0 , g 2 E | ( M * , E * , R * ) = α 4 α 8 ( α 8 + E * ) 2 + β 4 M * 2 k 3 2 + M * 2 α 5 R * : = γ 2 , g 2 R | ( M * , E * , R * ) = α 5 E * : = γ 3 < 0 , g 3 M | ( M * , E * , R * ) = α 2 R * : = θ 1 < 0 , g 3 E | ( M * , E * , R * ) = β 2 ( R T R * ) α 2 R * : = θ 2 , g 3 R | ( M * , E * , R * ) = β 2 E * 2 ( R T R * ) 2 β 3 k r 1 2 [ k r 1 2 + ( R T R * ) 2 ] 2 : = θ 3 < 0 ,
recalling that, due to Theorem 1. it also holds that ( M * , E * , R * ) R + 3 .
Theorem 2.
Assume that ( M * , E * , R * ) is the steady-state solution of ODE system (1). Given that
(A). 
ω 1 , γ 2 < 0 and θ 2 > 0 ,
(B). 
ω 1 θ 2 < ω 2 θ 1 and ω 2 γ 1 < ω 1 γ 2 ,
(C). 
ω 2 θ 1 γ 3 + ω 1 ω 2 γ 1 > 0 and γ 2 γ 3 θ 2 + ω 2 γ 1 γ 2 > 0 , then (M*, E*, R*) is asymptotically stable.
Proof. 
The Jacobian matrix for the system is given by
J = ω 1 ω 2 0 γ 1 γ 2 γ 3 θ 1 θ 2 θ 3 .
Next, we compute the eigenvalues of J by solving the equation
| J λ I | = ω 1 λ ω 2 0 γ 1 γ 2 λ γ 3 θ 1 θ 2 θ 3 λ = 0 ,
where I is the 3 × 3 identity matrix. The characteristic polynomial of J is then given by
λ 3 + p 1 λ 2 + p 2 λ + p 3 = 0 ,
with coefficients,
p 1 = ( ω 1 + γ 2 + θ 3 ) , p 2 = ω 1 γ 2 + ω 1 θ 3 + γ 2 θ 3 γ 3 θ 2 ω 2 γ 1 , p 3 = ω 1 γ 3 θ 2 ω 1 γ 2 θ 3 + ω 2 γ 1 θ 3 ω 2 γ 3 θ 1 .
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,
p 3 = ω 1 γ 3 θ 2 ω 2 γ 3 θ 1 + ω 2 γ 1 θ 3 ω 1 γ 2 θ 3 = γ 3 ( ω 1 θ 2 ω 2 θ 1 ) + θ 3 ( ω 2 γ 1 ω 1 γ 2 )
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,
p 1 p 2 p 3 = ( ω 1 + γ 2 + θ 3 ) ( ω 1 γ 2 + ω 1 θ 3 + γ 2 θ 3 γ 3 θ 2 ω 2 γ 1 ) ( ω 2 γ 1 θ 3 + ω 1 γ 3 θ 2 ω 1 γ 2 θ 3 ω 2 θ 1 γ 3 ) .
After simplifying and grouping, we obtain
p 1 p 2 p 3 = ω 1 2 γ 2 ω 1 2 θ 3 + ω 1 ω 2 γ 1 ω 1 γ 2 2 ω 1 γ 2 θ 3 γ 2 2 θ 3 + γ 2 γ 3 θ 2 + ω 2 γ 1 γ 2 ω 1 γ 2 θ 3 ω 1 θ 3 2 γ 2 θ 3 2 + γ 3 θ 2 θ 3 + ω 2 θ 1 γ 3 = ω 1 2 γ 2 ω 1 2 θ 3 + ω 2 θ 1 γ 3 + ω 1 ω 2 γ 1 ω 1 γ 2 2 ω 1 γ 2 θ 3 γ 2 2 θ 3 + γ 2 γ 3 θ 2 + ω 2 γ 1 γ 2 ω 1 γ 2 θ 3 ω 1 θ 3 2 γ 2 θ 3 2 + γ 3 θ 2 θ 3 .
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. □

3. Results

3.1. Bifurcation Analysis

Bifurcation analysis of the system of ODEs (1) was carried out using XPPAUT, a software program freely available from https://asmedigitalcollection.asme.org/appliedmechanicsreviews/article/56/4/B53/463898/Simulating-Analyzing-and-Animating-Dynamical (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 G F < 0.3735 with M * < 1.225 . For 0.3735 < G F < 0.4138 , bistable dynamics result from a saddle node bifurcation marked by two black dots labelled as SN (saddle node). As G F 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 G F and generate a two parameter bifurcation diagram shown in Figure 2c. In the latter figure, the bistable region is coloured green.

3.2. 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 G F 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 G F = 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 G F 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 G F 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 G F = 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 m i n ( θ , ϕ ) = 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.

3.3. 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 Table 2, Table 3, Table 4 and Table 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 G F , as shown in Figure 5a,b. Clearly, k 1 and G F 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 G F increases M * . Second, the steady state E * is mostly sensitive to α 2 and α 5 followed by k 1 , k r 1 and G F , 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 G F , (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 G F 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 G F 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 Table 2 and Table 4. This fact is also true for 5% and 10% decreases in parameters, as shown in Table 3 and Table 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 Figure 5 and Figure 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 G F 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 R b family of proteins. Under these circumstances, cells would naturally enter the quiescence state. Additionally, an increase in the expression of the R b family of proteins significantly reduces the expression levels of E (E2F, Cyclin E and Cyclin A), thereby preventing cells from entering the cell cycle.

4. 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 R b E 2 F bistable switch converts graded and transient serum growth signals into a binary and hysteric E 2 F 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 R b E 2 F 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:
    • For G F < 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.
    • On the other hand, in the range 0.3735 < G F < 0.4138 , System (1) exhibits bi-stability, 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.
    • For values of G F > 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 R b 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 R b E 2 F 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.

5. Limitations

In this study, we constructed an activator-inhibitor model to describe the R b E 2 F 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 R b E 2 F signalling network such as oscillatory dynamics, as observed in [62]. Moreover, we assumed that all proteins other than the R b family of proteins exist in abundance which may not be the case and hence must be investigated in future studies.

Author Contributions

K.Z.M. typed the manuscript, J.C.P. assisted in mathematical analysis as well as helped with editing the article. V.O.J. supervised the model formulation, numerical, bifurcation and sensitivity analyses as well as providing the Matlab codes. N.I.K. supervised most of the mathematical analysis as well as helped in editing the manuscript. A.M. conceived and supervised all aspects of the study, including model formulation, analysis and simulations, as well as overseeing the manuscript write up. All authors have read and agreed to the published version of the manuscript.

Funding

All authors (K.Z.M., J.C.P., V.O.J., N.I.K., A.M.) acknowledge support from the EPSRC grant (EP/T00410X/1): UK–Africa Postgraduate Advanced Study Institute in Mathematical Sciences (UK-APASI) and from the African Diaspora Mathematicians Program—Developing international relations in Mathematics between the University of Zimbabwe (Zimbabwe) and Sussex (UK), Commission for Developing Countries (awarded to A. Madzvamuse and S. Mukwembi). A.M.’s work was partially funded by the Leverhulme Trust Research Project Grant (RPG-2014-149) and the European Union Horizon 2020 research and innovation programme under the Marie Skłodowska-Curie grant agreement (No 642866). This work (A.M.) was partly supported by Health Foundation (1902431), the NIHR (NIHR133761) and by an individual grant from the Dr Perry James (Jim) Browne Research Centre on Mathematics and its Applications (University of Sussex). A.M. is a Royal Society Wolfson Research Merit Award Holder funded generously by the Wolfson Foundation. A.M. is a Distinguished Visiting Scholar to the Department of Mathematics, University of Johannesburg, South Africa.

Data Availability Statement

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

Acknowledgments

The authors wish to thank the University of Zimbabwe for funding the preliminary studies through the Graduate Teaching Assistantship Program (GTA). The authors also wish to thank G. Muchatibaya and colleagues in the University of Zimbabwe for supporting the research. AM thanks the International Mathematics Union, the Commission for Developing Countries and the UK Global Challenges Fund (through the Engineering and Physical Sciences Research Council) for supporting human infrastructure development and capacity building in Africa.

Conflicts of Interest

The authors declare that there is no conflict of interest regarding the publication of this paper.

References

  1. Wang, X.; Fujimaki, K.; Mitchell, G.C.; Kwon, J.S.; Croce, K.D.; Langsdorf, C.; Zhang, H.H.; Yao, G. Exit from quiescence displays a memory of cell growth and division. Nat. Commun. 2017, 8, 321. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  2. Yao, G. Modelling mammalian cellular quiescence. Interface Focus 2014, 4, 20130074. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  3. Harashima, H.; Dissmeyer, N.; Schnittger, A. Cell cycle control across the eukaryotic kingdom. Trends Cell Biol. 2013, 23, 345–356. [Google Scholar] [CrossRef] [PubMed]
  4. Miller, A.K.; Munger, K.; Adler, F.R. A mathematical model of cell cycle dysregulation due to human papilloma virus infection. Bull. Math. Biol. 2017, 79, 1564–1585. [Google Scholar] [CrossRef]
  5. Hartwell, L.H.; Weinert, T.A. Checkpoints: Controls that ensure the order of cell cycle events. Science 1989, 246, 629–634. [Google Scholar] [CrossRef] [Green Version]
  6. Kastan, M.B.; Bartek, J. Cell-cycle checkpoints and cancer. Nature 2004, 432, 316–323. [Google Scholar] [CrossRef]
  7. Naetar, N.; Soundarapian, V.; Litovchick, L.; Goguen, K.L.; Sablina, A.A.; Bowman-Colin, C.; Sicinski, P.; Hahn, W.C.; DeCaprio, J.A.; Livingstone, D.M. PP2A-mediated regulation of Ras signaling in G2 is essential for stable quiescence and normal G1 length. Mol. Cell 2014, 54, 932–945. [Google Scholar] [CrossRef] [Green Version]
  8. Nasmyth, K. Viewpoint: Putting the cell cycle in order. Science 1996, 274, 1643–1645. [Google Scholar] [CrossRef]
  9. Pardee, A.B. A restriction point for control of normal cell proliferation. Proc. Natl. Acad. Sci. USA 1974, 71, 1286–1290. [Google Scholar] [CrossRef] [Green Version]
  10. Qu, Z.; Weiss, J.N.; Maclellan, R. Regulation of mammalian cell cycle: A model of the G1-to-S transition. J. Physiol. Cell Physiol. 2003, 284, C349–C364. [Google Scholar] [CrossRef] [Green Version]
  11. Weinberg, R. The Biology of Cancer. In Garland Science; WW Norton & Company: New York, NY, USA, 2013; ISBN 978-0815342205. [Google Scholar]
  12. Zetterberg, A.; Larsson, O. Kinetic analysis of regulatory events in G1 leading to proliferation or quiescence of Swiss 3T3 cells. Proc. Natl. Acad. Sci. USA 1985, 82, 5365–5369. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  13. Pandey, N.; Vinod, P.K. Mathematical modeling of reversible transition between quiescence and proliferation. PLoS ONE 2018, 13, e0198420. [Google Scholar] [CrossRef] [PubMed]
  14. Chen, H.Z.; Tsai, S.Y.; Leone, G. Emerging roles of E2Fs in cancer: An exit from cell cycle control. Nat. Rev. Cancer 2009, 9, 785–797. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  15. Nevins, J.R. The Rb/E2F and cancer. Hum. Mol. Genet 2001, 10, 699–703. [Google Scholar] [CrossRef]
  16. Nevins, J.R. E2F: A link between the Rb tumor suppressor protein and viral oncoproteins. Science 1992, 258, 424–429. [Google Scholar] [CrossRef]
  17. Yao, G.; Lee, T.J.; Mori, S.; Nevins, J.R.; You, L. A bistable Rb-E2F switch underlies the restriction point. Nat. Cell Biol. 2008, 10, 476–482. [Google Scholar] [CrossRef]
  18. Burns, F.J.; Tannock, I.F. On existence of a G0-phase in cell cycle. Cell Tissue Kinet. 1970, 3, 321. [Google Scholar]
  19. De Maertelaer, V.; Galand, P. Some properties of a “G0” -model of the cell cycle. II. Natural constraints on the theoretical model in exponential growth conditions. Cell Tissue Kinet. 1975, 8, 11–22. [Google Scholar] [CrossRef]
  20. Shields, R.; Smith, J.A. Cells regulate their proliferation through alterations in transition probability. J. Cell. Physiol. 1977, 91, 345–355. [Google Scholar] [CrossRef]
  21. Smith, J.A.; Martin, L. Do cells cycle? Proc. Natl Acad. Sci. USA 1973, 70, 1263–1267. [Google Scholar] [CrossRef] [Green Version]
  22. Castor, L.N. A G1 rate model accounts for cell cycle kinetics attributed to transition probability. Nature 1980, 287, 857–859. [Google Scholar] [CrossRef] [PubMed]
  23. Koch, A.L. Does the variability of cell cycle result from one or many chance events? Nature 1980, 286, 80–82. [Google Scholar] [CrossRef] [PubMed]
  24. Friend, S.H.; Bernards, R.; Rogelj, S.; Weinberg, R.A.; Rapaport, J.M.; Albert, D.M.; Drtja, T.P. A human DNA segment with properties of the gene that predisposes to retinobalstoma and osteosarcoma. Nature 2004, 323, 643–646. [Google Scholar] [CrossRef] [PubMed]
  25. Tyson, J.J.; Chen, K.C.; Novak, B. Sniffers, buzzers, toggles and blinkers: Dynamics of regulatory and signaling pathways in the cell. Curr. Opin. Cell. Biol. 2003, 15, 221–231. [Google Scholar] [CrossRef]
  26. Aguda, B.D.; Tang, Y. The kinetic origins of the restriction point in the mammalian cell cycle. Cell Prolif. 1999, 32, 321–335. [Google Scholar] [CrossRef] [PubMed]
  27. Gardner, T.S.; Dolnik, M.; Collins, J.J. A theory for controlling cell cycle dynamics using a reversibility binding inhibitor. Proc. Natl. Acad. Sci. USA 1998, 95, 14190–14195. [Google Scholar] [CrossRef] [Green Version]
  28. Goldbeter, A. A minimal cascade model for the mitotic oscillator involving cyclin and cdc kinase. Proc. Natl. Acad. Sci. USA 1991, 88, 9107–9111. [Google Scholar] [CrossRef] [Green Version]
  29. Kohn, K.W. Functional capabilities of molecular network components controlling the mammalian G1/S cell cycle phase transition. Oncogene 1998, 16, 1065–1075. [Google Scholar] [CrossRef] [Green Version]
  30. Norel, R.; Agur, Z. A model for the adjustment of the mitotic clock by cyclin and MPF levels. Science 1991, 251, 1076–1078. [Google Scholar] [CrossRef]
  31. Nurse, P. The incredible life and times of biological cells. Science 2000, 289, 1711–1716. [Google Scholar] [CrossRef] [Green Version]
  32. Obeyesekere, M.N.; Knudsen, E.S.; Wang, J.Y.J.; Zimmernan, S.O. A mathematical model of the regulation of the G1 phase of Rb+/+ and Rb−/− mouse embryonic fibroblasts and an osteosarcoma cell line. Cell Prolif. 1997, 30, 171–194. [Google Scholar] [CrossRef] [PubMed]
  33. Thron, C.D. Bistable biochemical switching and the control of the events of cell cycle. Oncogene 1997, 15, 317–325. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  34. Tyson, J.J.; Novak, B. Regulation of the eukaryotic cell: Molecular antagonism hysteresis, and irreversible transitions. J. Theor. Biol. 2001, 210, 249–263. [Google Scholar] [CrossRef] [Green Version]
  35. Tyson, J.J.; Novak, B. Checkpoints in the cell cycle from modeler’s perspective. Prog. Cell Cycle Res. 1995, 1, 1–8. [Google Scholar] [CrossRef] [PubMed]
  36. Novak, B.; Tyson, J.J. Modeling the control of DNA replication in fission yeast. Proc. Natl. Acad. Sci. USA 1997, 94, 9147–9152. [Google Scholar] [CrossRef] [Green Version]
  37. Thron, C.D. A model for a bistable biochemical trigger of mitosis. Biophys. Chem. 1996, 57, 239–251. [Google Scholar] [CrossRef]
  38. Aguda, B.D. A quantitative analysis of the kinetics of the G2 DNA damage checkpoint system. Proc. Natl. Acad. Sci. USA 1999, 96, 11352–11357. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  39. Frolov, M.V.; Dyson, N.J. Molecular mechanisms of E2F-dependent activation and pRB-mediated repression. J. Cell Sci. 2004, 117, 2173–2181. [Google Scholar] [CrossRef] [Green Version]
  40. Heldt, F.S.; Barr, A.R.; Cooper, S.; Bakal, C.; Novak, B. A comprehensive model for the proliferation-quiescence decision in response to endogenous DNA damage in human cells. Proc. Natl. Acad. Sci. USA 2017, 115, 2532–2537. [Google Scholar] [CrossRef] [Green Version]
  41. Attwooll, C.; Denchi, E.L.; Helin, K. The E2F family: Specific functions and overlapping interests. EMBO J. 2004, 23, 4709–4716. [Google Scholar] [CrossRef]
  42. Yao, G.; Tan, C.; West, N.; Nevins, J.T.; You, L. Origin of bistability underlying mammalian cell cycle entry. Mol. Syst. Biol. 2011, 7, 485. [Google Scholar] [CrossRef] [PubMed]
  43. Blagosklonny, M.V.; Pardee, A.B. The restriction point of the cell cycle. Cell Cycle 2002, 1, 103–110. [Google Scholar] [CrossRef] [Green Version]
  44. Gierer, A.; Meinhardt, H.A. A theory of biological pattern formation. Kybernetik 1972, 12, 30–39. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  45. Kamps, D.; Koch, J.; Juma, V.O.; Campillo-Funollet, E.; Graessl, M.; Banerjee, S.; Mazel, T.; Chen, X.; Wu, Y.; Portet, S.; et al. Optogenetic Tuning Reveals Rho Amplification-Dependent Dynamics of a Cell Contraction Signal Network. Cell Rep. 2020, 33, 108467. [Google Scholar] [CrossRef] [PubMed]
  46. Sears, R.C.; Nevins, J.R. Signaling network that link cell proliferation and cell fate. J. Biol. Chem. 2002, 277, 11617–11620. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  47. Alon, U. An Introduction to Systems Biology: Design Principles of Biological Circuits; Chapman and Hall/CRC: New York, NY, USA, 2006; ISBN 9781439837177. [Google Scholar]
  48. Michaelis, L.; Menten, M.L. Die Kinetik der Invertinwirkung. Biochem. Z. 1913, 49, 333–369. [Google Scholar]
  49. Murray, J.D. Mathematical Biology I: An Introduction of Interdisciplinary Applied Mathematics; Springer: New York, NY, USA, 2002; pp. 257–271. [Google Scholar]
  50. Henley, A.S.; Dick, F.A. The retinoblsatoma family of proteins and their regulatory functions in the mammalian cell division cycle. Cell Div. 2012, 7, 1–14. [Google Scholar] [CrossRef] [Green Version]
  51. Hsu, S.-B. Ordinary Differential Equations with Applications, 2nd ed.; National Hua University: Hsinchu, Taiwan, 2013; ISBN 978-981-4452-92-2. [Google Scholar]
  52. Fine, B.; Rosenberger, G. The Fundamental Theorem of Algebra; Springer: New York, NY, USA, 1997. [Google Scholar]
  53. Ermentrout, B. Simulating, Analyzing, and Animating Dynamical Systems: A Guide to XPPAUT for Researchers and Students; SIAM: Philadelphia, PA, USA, 2002; ISBN 978-0-89971-506-4. [Google Scholar]
  54. Shampine, L.F.; Reichdt, M.W. MATLAB ODE Suite. SIAM J. Sci. Comp. 1997, 18, 1–22. [Google Scholar] [CrossRef] [Green Version]
  55. Dormand, J.R.; Prince, P.J. A family of embedded Runge-Kutta formulae. J. Comp. Appl. Math. 1980, 6, 19–26. [Google Scholar] [CrossRef] [Green Version]
  56. Juma, V.O. Data-Driven Mathematical Modeling and Simulations of Rho-Myosin Dynamics. Ph.D. Thesis, University of Sussex, Sussex, UK, 2019. [Google Scholar]
  57. Juma, V.O.; Dehmelt, L.; Portet, S.; Madzvamuse, A. A mathematical analysis of an activator-inhibitor Rho GTPase model. J. Comput. Dyn. 2021, 9, 133. [Google Scholar] [CrossRef]
  58. Zagkos, L.; Mc Auley, M.; Roberts, J.; Kavallaris, N.I. Mathematical models of DNA methylation dynamics: Implications for health and ageing. J. Theor. Biol. 2019, 462, 184–193. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  59. Zhang, X.-J.; Trame, M.N.; Lesko, L.J.; Schmidt, S. Sobol sensitivity analysis: A tool to guide the development and evaluation of systems phamacology models. CPT Pharmacomet. Syst. Pharmacol. 2015, 4, 69–79. [Google Scholar] [CrossRef] [PubMed]
  60. Saltelli, A.; Chan, K.; Scott, E.M. Sensitivity Analysis: Wiley Series in Probability and Statistics; Jon Wiley: Chichester, UK, 2000; ISBN 978-0-470-74382-9. [Google Scholar]
  61. Hampy, D.M. A review of techniques for parameter sensitivity analysis of environmental models. Environ. Monit. Assess. 1994, 32, 135–154. [Google Scholar] [CrossRef] [PubMed]
  62. Yan, F.; Liu, H.; Hao, J.; Liu, Z. Dynamical behaviors of RBE2F pathway including negative feedback loops involving miR449. PLoS ONE 2012, 7, e43908. [Google Scholar] [CrossRef]
Figure 1. An activator-inhibitor network for the R b E 2 F signalling network [42]. The solid arrows represent activation mechanism while the broken lines represent inhibition. Here, the growth factors G F activate M (cyclin D and M y c ) and, in turn M activates E (cyclin A , cyclin E and E 2 F transcription factors) while E activates M, forming a positive feedback loop. R ( R b family of proteins) inhibits E, while E activates and inhibits R and M inhibits R.
Figure 1. An activator-inhibitor network for the R b E 2 F signalling network [42]. The solid arrows represent activation mechanism while the broken lines represent inhibition. Here, the growth factors G F activate M (cyclin D and M y c ) and, in turn M activates E (cyclin A , cyclin E and E 2 F transcription factors) while E activates M, forming a positive feedback loop. R ( R b family of proteins) inhibits E, while E activates and inhibits R and M inhibits R.
Mathematics 10 02426 g001
Figure 2. Numerical bifurcation diagrams for the system of ODEs (1), with parameter values listed in Table 1. (a) One parameter bifurcation diagram for M * , with the growth factor G F concentration value as a bifurcation parameter. The saddle-node (SN) bifurcation points for G F are 0.3735 and 0.4138 and the corresponding M * are 3.826 and 1.225, respectively. The red line represents the stable steady states whereas the thin black line represent saddle points. The black dots labelled SN are the saddle-node bifurcation points, which span the region characterised by bistable dynamics. (b) A one-parameter bifurcation diagram showing a zoom of (a) around the two saddle-node bifurcation points. (c) Two-parameter bifurcation diagram with k 1 and G F as bifurcation parameters. The green coloured region is characterised by bistable dynamics, whereas the colourless region is characterised by stable dynamics.
Figure 2. Numerical bifurcation diagrams for the system of ODEs (1), with parameter values listed in Table 1. (a) One parameter bifurcation diagram for M * , with the growth factor G F concentration value as a bifurcation parameter. The saddle-node (SN) bifurcation points for G F are 0.3735 and 0.4138 and the corresponding M * are 3.826 and 1.225, respectively. The red line represents the stable steady states whereas the thin black line represent saddle points. The black dots labelled SN are the saddle-node bifurcation points, which span the region characterised by bistable dynamics. (b) A one-parameter bifurcation diagram showing a zoom of (a) around the two saddle-node bifurcation points. (c) Two-parameter bifurcation diagram with k 1 and G F as bifurcation parameters. The green coloured region is characterised by bistable dynamics, whereas the colourless region is characterised by stable dynamics.
Mathematics 10 02426 g002
Figure 3. Numerical simulations illustrating the time evolution dynamics of ( M ( t ) , R ( t ) , E ( t ) ) for system (1) corresponding to different dynamic regimes. As the value of G F increases, the system transitions from stable ((a) with G F = 0.3 and initial conditions M 0 = 0.2 , E 0 = 1.2 , R 0 = 3 ) through bistable ((b) with G F = 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 and R 0 = 5 ) then back to stable ((c) with G F = 0.45 and initial conditions M 0 = 20 , E 0 = 5 and R 0 = 5 ).
Figure 3. Numerical simulations illustrating the time evolution dynamics of ( M ( t ) , R ( t ) , E ( t ) ) for system (1) corresponding to different dynamic regimes. As the value of G F increases, the system transitions from stable ((a) with G F = 0.3 and initial conditions M 0 = 0.2 , E 0 = 1.2 , R 0 = 3 ) through bistable ((b) with G F = 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 and R 0 = 5 ) then back to stable ((c) with G F = 0.45 and initial conditions M 0 = 20 , E 0 = 5 and R 0 = 5 ).
Mathematics 10 02426 g003
Figure 4. Numerical bifurcations and simulations illustrating the dynamics of System (1) where we explore the conditions of Theorem 2. (ac) bifurcation diagrams for M ( t ) , E ( t ) and R ( t ) , respectively, with respect to α 5 ; other parameter values are fixed as shown in Table 1 and G F = 0.3 . (d) shows unbounded solutions of E and M with α 5 = 0.0079 , δ = 1 and G F = 0.3 , while keeping all the other parameters fixed as in Table 1 and initial conditions M 0 = 1.2 , E 0 = 3 , R 0 = 1.2 so that the conditions set in Theorem 2 and 3 are not met. To be more specific, γ 2 > 0 , which is against conditions (A), (B) and (C) of Theorem 2. The M ( t ) curve was amplified by a factor of 10 to visualize its behaviour with time and show that it grows unbounded with time while the R curve remain bounded because of conservation. (e) illustrates a zoom of the unstable dynamics of (d) for a short time t [ 0 , 2000 ] without scaling. (f) shows the stable time evolution of M ( t ) , E ( t ) and R ( t ) with α 5 = 0.0079 , δ = 1 and G F = 0.3 with initial condition ( M 0 , E 0 , R 0 ) = (0.4782, 0.594, 1.8616). With these parameters, the model converges to the steady state ( M * , E * , R * ) = (0.7255, 1.4184, 2.0449).
Figure 4. Numerical bifurcations and simulations illustrating the dynamics of System (1) where we explore the conditions of Theorem 2. (ac) bifurcation diagrams for M ( t ) , E ( t ) and R ( t ) , respectively, with respect to α 5 ; other parameter values are fixed as shown in Table 1 and G F = 0.3 . (d) shows unbounded solutions of E and M with α 5 = 0.0079 , δ = 1 and G F = 0.3 , while keeping all the other parameters fixed as in Table 1 and initial conditions M 0 = 1.2 , E 0 = 3 , R 0 = 1.2 so that the conditions set in Theorem 2 and 3 are not met. To be more specific, γ 2 > 0 , which is against conditions (A), (B) and (C) of Theorem 2. The M ( t ) curve was amplified by a factor of 10 to visualize its behaviour with time and show that it grows unbounded with time while the R curve remain bounded because of conservation. (e) illustrates a zoom of the unstable dynamics of (d) for a short time t [ 0 , 2000 ] without scaling. (f) shows the stable time evolution of M ( t ) , E ( t ) and R ( t ) with α 5 = 0.0079 , δ = 1 and G F = 0.3 with initial condition ( M 0 , E 0 , R 0 ) = (0.4782, 0.594, 1.8616). With these parameters, the model converges to the steady state ( M * , E * , R * ) = (0.7255, 1.4184, 2.0449).
Mathematics 10 02426 g004
Figure 5. Parameter sensitivities for system (1) corresponding to 5% increase and decrease in parameters. (a,b) show sensitivity of M * after 5 % increase and decrease in parameters, respectively; (c,d) show sensitivity of E * corresponding to 5% increase and decrease, respectively. (e,f) show sensitivity R * to 5% increase and decrease in parameters, respectively.
Figure 5. Parameter sensitivities for system (1) corresponding to 5% increase and decrease in parameters. (a,b) show sensitivity of M * after 5 % increase and decrease in parameters, respectively; (c,d) show sensitivity of E * corresponding to 5% increase and decrease, respectively. (e,f) show sensitivity R * to 5% increase and decrease in parameters, respectively.
Mathematics 10 02426 g005
Figure 6. Parameter sensitivities for system (1) corresponding to 10% increases and decrease in parameters. (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.
Figure 6. Parameter sensitivities for system (1) corresponding to 10% increases and decrease in parameters. (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.
Mathematics 10 02426 g006
Figure 7. Correlation between parameter sensitivity of M ( t ) , E ( t ) and R ( t ) corresponding to 10 % increase in model parameters. (a) shows the correction between M ( t ) and E ( t ) , (b) shows the correlation between M ( t ) and R ( t ) and (c) shows the correlation between E ( t ) and R ( t ) .
Figure 7. Correlation between parameter sensitivity of M ( t ) , E ( t ) and R ( t ) corresponding to 10 % increase in model parameters. (a) shows the correction between M ( t ) and E ( t ) , (b) shows the correlation between M ( t ) and R ( t ) and (c) shows the correlation between E ( t ) and R ( t ) .
Mathematics 10 02426 g007
Table 1. Parameters used for simulations and bifurcation analysis. Parameter values are taken from [17,42] and some of them adjusted to illustrate the qualitative dynamics. Here, U represents the unit of concentration and s the unit of time t.
Table 1. Parameters used for simulations and bifurcation analysis. Parameter values are taken from [17,42] and some of them adjusted to illustrate the qualitative dynamics. Here, U represents the unit of concentration and s the unit of time t.
ParameterDescriptionValueUnitsReference
α 1 Growth factors activation rate1 s 1 [17]
G F Growth factors concentrationvariesU[42]
δ Inhibition rate of M1.001 s 1 [17]
k 1 Michaelis–Menten constant1U[42]
β 1 Activation rate of M1 s 1 [42]
β 2 Activation of R protein family1U 1 s 1 [42]
α 3 Inhibition rate of R by E1U 1 s 1 Estimate
α 2 Inhibition rate of R by M1U 1 s 1 Estimate
β 3 R baseline inhibition rate1U s 1 Estimate
R T Total concentration of R5U[40]
α 4 E self activation rate0.02U s 1 [17]
β 4 Activation rate E by M0.02 s 1 [17]
β 5 R baseline inhibition rate1U s 1 Estimate
α 6 E constitutive activation rate0.001 s 1 [40]
α 8 Michaelis–Menten constant0.92U[42]
α 5 Inhibition rate of E by R0.01U 1 s 1 [42]
k r 1 Michaelis–Menten constant0.05U[42]
k r 2 Michelis–Menten constant1U[42]
k r 3 Michaelis–Menten constant1U[42]
Table 2. Percentage changes in the steady-state values M * , E * and R * after a 5% increase in parameter values of the non-linear model (1).
Table 2. Percentage changes in the steady-state values M * , E * and R * after a 5% increase in parameter values of the non-linear model (1).
5% Increase in Parameter% Change in M * % Change in E * % Change in R *
k 1 −10.330−10.1860.7016
δ −2.5797−2.6158−0.0973
α 6 −2.3573−6.62791.7049
α 4 1.62134.5563−1.1907
α 8 1.32123.7135−0.96101
β 4 −1.4816−4.16521.0916
k r 3 −3.6855−10.3612.78
α 5 3.675410.3310.9486
β 2 0.96892.72220.1993
R T 1.15863.2559−0.8912
β 3 0.46821.31500.0671
k r 1 −2.1313−5.9929−0.73 16
α 2 −4.4144−12.411−1.5102
α 3 0.00080.00217−0.0044
β 5 −0.50358−1.41710.2834
k r 2 −1.3792−3.8783−0.5012
G F 8.21818.6153−0.4477
β 1 4.83555.0458−0.30391
Table 3. Percentage changes in the steady-state values M * , E * and R * after a 5% decrease in parameter values of the non-linear model (1).
Table 3. Percentage changes in the steady-state values M * , E * and R * after a 5% decrease in parameter values of the non-linear model (1).
5% Decrease in Parameter% Change in M * % Change in E * % Change in R *
k 1 15.37916.574−0.74247
δ 3.05493.1573−0.23624
α 6 2.79897.8667−1.9380
α 4 −1.5130−4.25391.0604
α 8 −1.3083−3.67850.91636
β 4 1.48564.1750−1.1118
k r 3 4.694113.194−3.0811
α 5 −3.3974−9.5522−1.1635
β 2 −0.90100−2.5340−0.34861
R T −1,1643−3.27470.77670
β 3 −0.47058−1.3243−0.22096
k r 1 2.27506.39430.56813
α 2 5.652815.8891.4395
α 3 −0.00076−2.07070.00434
β 5 0.505061.4196−0.37425
k r 2 1.72254.84070.41722
G F −7.5303−7.48860.46275
β 1 −3.9715−4.01110.25848
Table 4. Percentage changes in the steady-state values M * , E * and R * after a 10% increase in parameter values of the non-linear model (1).
Table 4. Percentage changes in the steady-state values M * , E * and R * after a 10% increase in parameter values of the non-linear model (1).
10% Increase in Parameter% Change in M * % Change in E * % Change in R *
k 1 −18.792−18.0051.5896
δ −4.9970−5.02380.25366
α 6 −4.5627−12.8283.5291
α 4 3.54509.9639−2.4070
α 8 2.79187.8469−1.9470
β 4 −3.1078−8.73682.3944
k r 3 −6.9382−19.5055.7536
α 5 8.060522.6591.9911
β 2 2.12565.97360.52836
R T 2.42696.8214−1.7297
β 3 0.98002.75390.2105
k r 1 −4.3112−12.122−1.4863
α 2 −8.2964−23.324−3.0332
α 3 0.00170.0047−0.0093
β 5 −1.0545−2.96580.7076
k r 2 −2.6078−7.3322−0.8934
G F 18.29619.8−0.7814
β 1 11.81512.677−0.6242
Table 5. Percentage changes in steady-state values M * , E * and R * after a 10% decrease in parameter values of the non-linear model (1).
Table 5. Percentage changes in steady-state values M * , E * and R * after a 10% decrease in parameter values of the non-linear model (1).
10% Decrease in Parameter% Change in M * % Change in E * % Change in R *
k 1 43.28451.402−0.90841
δ 6.3407646.6298−0.40674
α 6 5.821616.364−3.7086
α 4 −2.7912−7.84752.0483
α 8 −2.4764−6.96311.7939
β 4 2.82647.9439−1.9852
k r 3 10.23628.774−5.9032
α 5 −6.2295−17.514−2.2045
β 2 −1.6616−4.6717−0.56234
R T −2.2183−6.23771.5914
β 3 −0.89676−2.5224−0.34685
k r 1 4.4444412.491.1354
α 2 12.39434.8392.8958
α 3 −0.0014110−0.00384230.0081792
β 5 0.959362.6969−0.73757
k r 2 3.681210.3 470.94703
G F −13.836−13.4701.0361
β 1 −7.0325−7.02610.41178
Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Share and Cite

MDPI and ACS Style

Mapfumo, K.Z.; Pagan’a, J.C.; Juma, V.O.; Kavallaris, N.I.; Madzvamuse, A. A Model for the Proliferation–Quiescence Transition in Human Cells. Mathematics 2022, 10, 2426. https://doi.org/10.3390/math10142426

AMA Style

Mapfumo KZ, Pagan’a JC, Juma VO, Kavallaris NI, Madzvamuse A. A Model for the Proliferation–Quiescence Transition in Human Cells. Mathematics. 2022; 10(14):2426. https://doi.org/10.3390/math10142426

Chicago/Turabian Style

Mapfumo, Kudzanayi Z., Jane C. Pagan’a, Victor Ogesa Juma, Nikos I. Kavallaris, and Anotida Madzvamuse. 2022. "A Model for the Proliferation–Quiescence Transition in Human Cells" Mathematics 10, no. 14: 2426. https://doi.org/10.3390/math10142426

APA Style

Mapfumo, K. Z., Pagan’a, J. C., Juma, V. O., Kavallaris, N. I., & Madzvamuse, A. (2022). A Model for the Proliferation–Quiescence Transition in Human Cells. Mathematics, 10(14), 2426. https://doi.org/10.3390/math10142426

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop