Linear and structural stability of a cell division process model

The paper investigates the linear stability of mammalian physiology time-delayed flow for three distinct cases (normal cell cycle, a neoplasmic cell cycle, and multiple cell arrest states), for the Dirac, uniform, and exponential distributions. For the Dirac distribution case, it is shown that the model exhibits a Hopf bifurcation for certain values of the parameters involved in the system. As well, for these values, the structural stability of the SODE is studied, using the five KCC-invariants of the second-order canonical extension of the SODE, and all the cases prove to be Jacobi unstable.

2 Stability of a cell division process model (ii) E-the active complex formed by cyclin E and a cyclin-dependent kinase named cdk 2; (iii) R-the unphosphorilated retinoblastoma protein Rb. The time-evolution of these variables is governed by the following SODE: (1.1) The parameters involved are as follows: (i) a E is the rate of forming the cyclin E/ cdk 2 complex; (ii) a M and g are rates of activation of M proportional to E (g is an autocatalitic rate); (iii) p 1 and q 1 are parts of a rational term used to model the dephosphorilation of Rb; (iv) p 2 and q 2 are parts of a rational term used to model the phosphorilation of Rb; (v) d E is the rate of degradation of E promoted by M; (vi) d M is the rate of degradation of M; (vii) h is a constant related to the inhibition of E production by Rb; (viii) R T -total retinoblastoma protein Rb concentration. The parameter values for the mathematical model to obtain a normal cell cycle, a neoplasmic cell cycle, and multiple cell arrest states are given in Table 1.1 [15].
These parameter values were selected to obtain cell cycles for normal and neoplastic cells and to find multiple arrest states.

The time-delayed evolution flow
We will study the biological flow when one variable coordinate is subject to time-delay. In our case, we assume this to be R-which denotes the concentration of unphosphorilated V. Balan and I. R. Nicola 3 retinoblastoma protein. It is most common within physiological processes which take place in the living cells that certain time lags might interfere. This happens because the chemical reaction which succeeds in a chain reaction can be delayed by a lot of molecular mechanisms. Our motivation for choosing the variable R as a delayed one is not quite random. Namely, it is known that the change of a cell from normal to cancerous can be caused by a simple mutation in the sequence of a certain gene-called "oncogene," which can be activated by a lot of biological reasons. The Rb protein is the product of a gene from the class of the so-called "tumor supressor genes," which helps on preventing the development of a cancer. The Rb protein is involved in the cellular replication during the S-phase of the cell cycle and is active only when it is not phosphorilated. We are interested to see if a time lag in a process of phosphorilation or dephosphorilation of Rb protein could cause Hopf bifurcation, leaving still open the question of biological interpretations to their full extent.
To this aim, we determine first the equilibrium points of the system (1.1), obtained by setting to zero the right side of the differential system (1.1). The resulting nonlinear system has in general several solutions; only the positive ones can be accepted from physiological point of view. We will denote such a solution as (M * ,E * ,R * ).
Further, in order to obtain the dynamical system with delayed argument in the dependent variable R(t), it is known that for a given probability density g : (2.1) Applying the time-delay process to R, one changes the system (1.1) into the new SODE [5,12,17] dM with where the transform R(t) is defined by (2.1) and ϕ : [−τ,0] → R is a differentiable function which describes the behavior of the flow in the O direction. In other words, the initial SODE is replaced by a differential-functional system. Regarding the linearization of the SODE (2.2) we have the following statement [10].
Proposition 2.1. The linearized SODE of the differential autonomous system with delayed argument (2.2) at the equilibrium point (M * ,E * ,R * ) iṡ where the t-differentation is denoted by "dot," V (t) = t (M(t),E(t),R(t)), and , and ( f 1 , f 2 , f 3 ) are the components of the field which provides the SODE (1.1).
We will further study the subcases when the SODE describes a normal cell cycle, a neoplasmic cell cycle, and multiple cell arrest states. For the three cases of parameter sets and for the three types of distributions we further develop the basic results regarding the stability and bifurcation of the considered delayed SODE.
The values of the parameters which describe the three states are given in Table 1.1. We attach to the autonomous SODE (1.1) the initial condition In this case, we obtain the following.
Three notable distributions are worthy to consider: Dirac, uniform, and gamma. In these cases, the delayed R-component of the system has, respectively, the following forms.
(1) If g is the Dirac distribution of τ ≥ 0, that is, then the transform R(t) = R(t − τ) denotes the variable R with delayed argument.
V. Balan and I. R. Nicola 5   6 Stability of a cell division process model For m = 1, we obtain the exponential distribution and In the following, we examine the three distribution cases with the parameter sets corresponding to the subcases (a), (b), and (c).
(1) The Dirac distribution. In this case, R(t) = R(t − τ) and (2.6) becomes which has the form Let J(λ) be the characteristic quasipolynomial function in (2.13). The existence of τindependent solutions of (2.13) requires the condition a 5 a 1 a 2 − a 3 2 + a 2 2 a 5 − a 1 a 4 2 = 0. (2.14) A straightforward remark is that there exist no such solutions in this particular subcase. Looking for the critical values of the parameter τ at which there is an exchange of V. Balan and I. R. Nicola 7 structural stability, we note that the solutions of the characteristic equation (2.13) are of the form λ = λ(τ) = u(τ) ± iω(τ) ∈ C, and that (2.13) is equivalent to ReJ(λ) = Im J(λ) = 0.
A prerequisite for studying the Hopf bifurcation consists in finding the critical values of the parameter τ, by imposing u(τ) = 0 and ω(τ) = 0. Under these assumptions, we infer the nonlinear system in terms of ω and τ: We have the following proposition. (b) In all three cases, λ(τ) = ±iω(τ), and the complex values λ (τ 0 ) given in the table satisfy the transversality condition Re λ (τ 0 ) > 0.
Proof. The results are derived using the Maple 9.5 package.
Proposition 2.5. (a) For τ = 0, (2.13) has just three roots, which are the eigenvalues (see Table 2.1, (EV)) of the Jacobian matrix of the vector field at the considered equilibrium point.
(b) In the cases of normal cell cycle and two arrest cell cycles, (2.13) has at least one root with positive real part, so τ 0 is not a Hopf bifurcation parameter.
In the case of neoplasmic cell cycle, (2.13) has two imaginary conjugate roots λ = ±iω, an infinity of roots with negative real part, and no root with positive real part, so this is the only case for Dirac distribution when Hopf bifurcation occurs.
Proof. Maple 9.5 computations lead to the following results: (i) in the normal cell cycle case, for τ = τ 0 and ω ∈ (−0.1,0.1) for the first equilibrium point, we get Re(λ) ∼ 0.02, and for the second equilibrium point, we get Re(λ) ∼ 0.21; (ii) in the two arrest cell cycles case, for τ = τ 0 and ω ∈ (−0.1,0.1) for the first equilibrium point, we get Re(λ) ∼ 0.02, and for the second equilibrium point, we get Re(λ) ∼ 0.55; (iii) in the neoplasmic cell cycle, while τ passes through τ 0 , the function u(τ) changes from negative to positive values. It follows that the critical value of τ for which the Hopf bifurcation appears is exactly τ = τ 0 , which prove the claims in the statement.
Following the same steps like in the case of Dirac distribution, we obtain that there exists no solution for (2.18) for u = 0 in terms of τ 0 and ω 0 .
(3) The exponential distribution. In this case, we have where the coefficients are given by Table 2.2.
Following the same steps like in the case of Dirac distribution, we obtain that there exists no solution for (2.21) for u = 0 in terms of ω 0 and τ 0 .
In the case of the Dirac distribution-neoplasmic cell cycle, the initial dynamical SODE becomes subject to the following result, known as the Hopf bifurcation theorem [16].

Consider in a neighborhood of the stationary point x = x(τ) the linear approximation dx/dt = A(τ)x of the system dx/dt
For n > 2, additionally assume that Re(λ k (τ)) < 0, k = 3,...,n and that there exists an isolated value τ 0 ∈ I such that u(τ 0 ) = 0, ω(τ 0 ) = 0, and du/dτ > 0. Under these hypotheses, one of the following assertions holds true: (a) the stationary point x = x(τ 0 ) is a center; for τ = τ 0 neighbor to τ 0 , there exists no periodic orbit around x(τ); (b) there exists a number b > τ 0 such that for each τ ∈ (τ 0 ,b), there exists a unique induced orbit around the stationary point x(τ) in a neighborhood of this point. This 1parameter family of closed orbits split at the stationary point x(τ 0 ), that is, for τ → τ 0 , the diameter of the closed orbit varies with |τ − τ 0 | 1/2 . In this case, for τ ≤ τ 0 , τ ∈ I, there exists no closed orbit neighbor to x(τ); (c) there exists a number a < τ 0 such that for each τ ∈ (a,τ 0 ), there exists a unique closed orbit around the stationary point x(τ 0 ) in one of its neighborhoods. This 1-parameter family of closed orbits split at the stationary point x(τ 0 ), that is, for τ → τ 0 , the diameter of the closed orbit varies with |τ − τ 0 | 1/2 . In this case, for τ ≥ τ 0 , τ ∈ I, there exists no closed orbit neighbor to x(τ).
We note that though in the case τ = 0, in a certain neighborhood of the singular point, there exist no closed orbits, for the delayed system at a certain delay τ = 0, the SODE may exhibit one of the situations (b) or (c) in Theorem 2.6, hence the presence of closed orbits arbitrarily close to the singular point may occur, and then the concentrations of the proteins (M,E,R) may have a cyclic evolution, provided that exists a certain delay in the retinoblastoma protein Rb. The complete characterization of these alternatives, including a detailed study of stability of limit cycles is the subject of further research.
Still, one may take into account as well the alternative when there exists no periodic orbit around the isolated point x(τ). In this case, it is very interesting that the only cell in which the Hopf bifurcation occurs is the neoplasmic one, which is intuitively expected. This might mean that if someone induces a delay in the processes which involve the Rb protein, there exists a possibility for the cell to defeat the onset of cancer. The study of stability of limit cycles can provide more information in this matter as well.

Jacobi structural stability
From a physiological point of view, Jacobi structural stability adds a degree of accuracy to the classical linear analysis, by studying the robustness and fragility of a biological system. In the case of a cancerous cell, of main interest for medicine is to locate the region where the cell exhibits "robust arrest" (i.e., the region where one has both Lyapunov and Jacobi stability).
We recall first some basics of KCC-theory [1,2,14]. Let x = (x 1 ,...,x n ), and let t be the 2n + 1 coordinates of an open connected subset Ω ⊂ R n × R n × R 1 . We consider a second-order system of differential equations of the form where each function g i (x,ẋ,t) is smooth in a neighborhood of some initial conditions (x 0 ,ẋ 0 ,t 0 ) ∈ Ω. In order to find the basic differential invariants of the system (3.2) (see Kosambi [8], Cartan [6], and Chern [7]) under the nonsingular coordinate transformations we define the KCC-covariant differential of a contravariant vector field ξ i (x) on the open subset Ω via where ";" indicates partial differentiation with respect toẋ, and the Einstein summation convention is used throughout. Using (3.4), the system (3.2) becomes where ε i defined here is a contravariant vector field on Ω and is called the first KCCinvariant, which is interpreted as an external force [1]. For the normal cell cycle case, denoting (M,E,R) = (x, y,z), the first invariant is given in Table 3.1.
The functions g i = g i (x,ẋ,t) are 2 homogeneous ones inẋ if and only if ε i = 0. In other words, ε i = 0 is a necessary and sufficient condition for a semispray to be a spray. It is obvious that for the geodesic spray of a Riemannian or Finsler manifold, the first invariant vanishes.
It can be easily seen that, since the system is of the formẋ = X(x), the first invariant has the components ε i = (1/2)(∂X i /∂x r )ẋ r , and hence this vanishes for null velocities, that is, on the stationary points of the field X. The obtained strongly nonlinear equation does not depend on velocities, and admits solutions in the first octant of the position variables. In the general extended case, when (M,E,R) ∈ R 3 , there exists a region in TR 3 where the first invariant cancels.
On the other hand, it is known that if the trajectories x i (t) of (3.2) are varied into nearby ones with respect to x asx i (t) = x i (t) + ξ i (t)η with the parameter η small, one gets the variational equations where "," indicates partial differentiation with respect to x. Using now the KCC-covariant differential (3.4), one obtains (3.6) in the covariant form where P i is called the second KCC-invariant of the system (3.2), or deviation curvature tensor. Its eigenstructure is an alternative to the Floquet theory, with the eigenvalues of P i j replacing the characteristic multipliers (also called Floquet exponents, [3,13]). In our case, it has the generic form Note that (3.7) is the Jacobi field equation when the starting system (3.2) is geodesic equations in either Finsler or Riemannian geometry. This justifies the usage of the term Jacobi stability for KCC-theory.
On the other hand, the Jacobi (3.7) of the Finsler manifold (M,F) can be written in the scalar form where ξ i = v(s)η i is a Jacobi field along γ : x i = x i (s), η i is the unit normal vector field along γ, and K is the flag curvature of (M,F) [4].
It is also known that the sign of K influences the geodesic rays [4]. Indeed, if K > 0, then the geodesic rays bunch together (are Jacobi stable), and if K < 0, then they disperse (are Jacobi unstable).
Hence negative flag curvature is equivalent to positive eigenvalues of P i j , and positive flag curvature is equivalent to negative eigenvalues of P i j . The following is known. Theorem 3.1 [1,2]. The trajectories of (3.2) are Jacobi stable if and only if the real parts of the eigenvalues of the deviation tensor P j i are strict negative everywhere and Jacobi unstable, otherwise.
The notion of Jacobi stability presented until here can be extended to the general case of the SODE (3.2) using the theorem above as the definition for the Jacobi stability of the trajectories of a SODE. The third, fourth, and fifth invariants of the system (3.2) are, respectively, These invariants prove to cancel in our case. A basic result of the KCC-theory which points out the role of the five invariants is the following.  Based on Maple computations, we can infer straightforwardly that for our SODEsubject to the requirement of having real (positive) solutions (M,E,R), there exists no coordinate change such that the coefficients of the new second-order SODE semispray V. Balan and I. R. Nicola 13 do all vanish. According to [1], this implies that the trajectories of the second-order extended system (including the field lines of the initial SODE) can never be lines, whatever coordinate system one might choose.
The following results regarding the KCC stability of the SODE can be stated.
Proposition 3.3. In the "neoplasmic cell" subcase-that is, the parameter values in Table  1.1, the deviation curvature tensor P i j has (a) three positive eigenvalues: if the first equilibrium point is considered  ) is Lyapunov unstable; in conclusion the cell is in a state of "fragile arrest" in the first case and in a state of "fragile oscillations" in the second case.
Considering the parameter p 1 , which is the first parameter involved in the dephosphorilation of Rb protein, variable within the interval (0,1) and the other parameters taking the values as in the case of neoplasmic cell, we get the following results regarding the Jacobi stability.
Proposition 3.5. For p 1 ∈ (0,1), the discriminant of the characteristic polynomial of the Jacobi matrix attached to (1.1) is positive and depends continuously on p 1 . The Jacobi matrix has two complex conjugates and a positive real eigenvalue for p 1 ∈ (0,1) and, hence the field lines of the system are Jacobi unstable.
In the following, we consider the parameter q 1 , which is the second rate involved in the rational term modeling the dephosphorilation Rb protein, freely varying in the interval (0,1); the effect on this paths of the SODE (1.1) is described by the following result.
Proposition 3.6. For q 1 ∈ (0,1), the discriminant of the characteristic polynomial of the Jacobi matrix attached to (1.1) is positive and depends continuously on q 1 . In this case, the Jacobi matrix has two complex conjugates and a positive real eigenvalue; hence the field lines of the SODE are Jacobi unstable.
The study of Jacobi stability is complementary to the study of linear stability and is based on the study of Lyapunov stability of whole trajectories in a region. Therefore the perturbation yields trajectories closed by to the reference trajectory. It is well known that in the case of Lyapunov stability, the perturbations of a stable equilibrium point lead to trajectories which will be dumped out, provided that these are small enough so as not to escape from a basin of attraction (see [2,14]). By varying the parameters p 1 and q 1 which model the process of dephosphorilation of Rb protein, we obtain only cells in a state of "fragile oscillation."

Conclusions
An analysis of the delayed system was performed in Section 2, for three typical distributions: Dirac, uniform, and exponential; still, it was shown that significant results concern just the first one. In this case, for the neoplasmic cell cycle, a value τ 0 for the time-delay parameter was determined, such that for any τ < τ 0 , the positive equilibrium point is locally asymptotically stable, for τ = τ 0 the system exhibits Hopf bifurcation, and for τ > τ 0 , the equilibrium point is unstable. In Section 3, the Jacobi stability (which characterizes robustness/fragileness) of the nondelayed dynamical system was studied, with focus on the case of the neoplasmic cell. It was shown that for the parameters p 1 and q 1 varying in the interval (0,1), the field lines of the system are Jacobi unstable.