Mathematical analysis of a generalised p53-Mdm2 protein gene expression model

Applied Mathematics and Computation 328 (2018) 26–44

Contents lists available at ScienceDirect

Applied Mathematics and Computation journal homepage: www.elsevier.com/locate/amc

Mathematical analysis of a generalised p53-Mdm2 protein gene expression model Monika J. Piotrowska a, Agnieszka Bartłomiejczyk b,∗, Marek Bodnar a a

Institute of Applied Mathematics and Mechanics, Faculty of Mathematics, Informatics and Mechanics, University of Warsaw, Banacha 2, Warsaw 02–097, Poland b ´ sk University of Technology, G. Narutowicza 11/12, Gdan ´ sk 80–233, Poland Faculty of Technical Physics and Applied Mathematics, Gdan

a r t i c l e

i n f o

MSC: Primary 34K11 34K13 34K18 34K20 49K25 Secondary 92B05 93C23 Keywords: Delay differential equations Stability analysis Hopf bifurcation p53 protein

a b s t r a c t We propose the generalisation of the p53-Mdm2 protein gene expression model introduced by Monk (2003). We investigate the stability of a unique positive steady state and formulate conditions which guarantee the occurrence of the Hopf bifurcation. We show that oscillatory behaviour can be caused not only by time lag in protein transcription process, but also can be present in the model without time delay. Moreover, we investigate the stability of new born periodic solutions. Theoretical results are illustrated by numerical simulations and interpreted from the biological point of view. © 2018 Elsevier Inc. All rights reserved.

1. Introduction The p53 protein was discovered in 1979 and its name is connected with its molecular weight determined by the analysis of SDS-PAGE (which in fact has been incorrectly calculated). The p53 protein is involved in many biological processes including cell cycle regulation, cell differentiation and apoptosis. As a tumour suppressor it is one of the subjects of interests of scientists working on the designing and improvement of cancer therapies. It also plays an important role in the treatment of neurodegenerative diseases such as cerebral ischemia, traumatic brain injury, epilepsy and Alzheimer’s disease. In [1,2] it is demonstrated that an inhibition of p53 prevents cell death in different models of neurodegenerative disorders. For the first few years it was thought that p53 is an oncogene, that is a gene that may cause the growth of cancer cells. Only later it turned out that its role is rather reverse and it is a tumour-suppressor that is protein which is necessary to keep the cell division under control. Just like a car brakes regulate its speed, the tumour-suppressor gene acts as brakes to the cell cycle, the DNA replication and the division of cells. If these proteins do not act properly, the uncontrolled growth defining the feature of cancer cells can appear. It is so in the case of p53 protein which was found to not act correctly in the most of human cancers, for more details see e.g., [3]. In 1992, Lane [4] called the p53 protein “guardian of the genome” because of its ability to guard the cells from the malignant transformations. When DNA is damaged (e.g., by ionizing radiation or chemicals), the appropriate p53-mediated pathways are activated. This yields the arrest of the cell in the cycle, which prevents the proliferation of the cell containing ∗

Corresponding author. E-mail addresses: [email protected] (M.J. Piotrowska), [email protected] (A. Bartłomiejczyk), [email protected] (M. Bodnar).

https://doi.org/10.1016/j.amc.2018.01.014 0 096-30 03/© 2018 Elsevier Inc. All rights reserved.

M.J. Piotrowska et al. / Applied Mathematics and Computation 328 (2018) 26–44

27

damaged DNA. Next, the biochemical processes involved in DNA repair are triggered. If they are successful, the cell resumes its progression and the cell division can take place. If the repair is unsuccessful due to excessive damage, the p53-mediated apoptotic pathway becomes functional, leading to apoptosis, i.e., the process of programmed cell death. To understand the action of the p53 protein it is crucial to take into account the interaction between p53 and Mdm2 proteins, see [5–10]. That interaction is the main negative regulator of the p53 protein in the cell because in the normal cells the p53 pathway is switched off, that is p53 activity is kept low, such that the cell cycle progression is not disrupted by a negative feedback loop consisting of the p53 and Mdm2 genes. Additionally the p53 protein acts as a transcription factor and regulates the expression of several target genes including Mdm2 gene. On the other hand, the Mdm2 protein inhibits the p53 activity, forming a negative feedback loop and this ability of Mdm2 to keep p53 under control is essential for the normal cell function. The repression of Mdm2 operates via three main mechanisms: (i) Mdm2 binds p53 at its DNA binding domain in such a way that the latter can not function as a transcription factor; (ii) p53 with binding Mdm2 is labelled for degradation; (iii) Mdm2 is responsible for an export of p53 from the nucleus to the cytoplasm changing its transcriptional activity. In the mammalian cells DNA damaged activates of a protein called ATM kinase [11,12] and phosphorylation of the p53 protein at a specific site, preventing the binding of the Mdm2 protein to p53. In the absence of Mdm2 mediated degradation of p53, quantity of p53 proteins stabilizes at a higher level than under the presence of Mdm2. In fact, the interactions between the feedback loops in p53 gene expression system are much more complex, see e.g., [13]. One of the first models of the p53 gene expression system was proposed by Bar-Or et al. in [3], where it is numerically shown how the oscillatory behaviour of both p53 and Mdm2 protein may vary depending on the system parameters. Later, experimental results showed that p53-Mdm2 interactions proposed in [3] lead, under particular experimental conditions, to undamped oscillations in the concentrations of the considered proteins [14–17]. More recently, Gordon et al. [18] and Sturrock et al. [19] studied models postulating that the p53-Mdm2 oscillations have their source in a nuclear transport of the proteins between the cytoplasm and cell nucleus. In that models transport was modelled by a passive diffusion and the impact of the diffusion rates on the behaviour of the concentration of the p53 protein was investigated. Many models that include the stochastic effect were also considered, see e.g., [15,20–23]. For example in [15], the negative and positive feedback loops with some random noise added into protein’s production terms were studied. On the other hand, models taking into account a transcriptional time delay were considered by a number of the research groups [24–26], but none of these examples includes the mathematical analysis. The exception is the article [26] in which the existence and the direction of the Hopf bifurcation for simple p53-Mdm2 system with time delay is investigated. It is worth pointing out that some much more complex models of the p53 gene regulatory network were also proposed. For example in [24] the feedback loop involving five components was considered. Some models consisting of several ODEs that include both negative and positive feedback loops were considered, [27–29]. Also large models that consist of more than ten differential equations [30] were studied. Clearly, in reality the dynamics of the p53 signalling pathway is more complex, and it involves some additional, experimentally confirmed, feedback loops (both positive and negative), see [13,31], which were neglected in the considered model. As a consequence in this model it is not possible to have a second positive steady state as in [31], where the positive feedback through PTEN is taken into account. Thus, the dynamics of the studied system clearly differs from that more complex models. However, we wish to point out that in this paper we consider p53 system right after the DNA damage when the longer positive feedback is not activated. The advantage of such complex models is that they describe the studied phenomenon more precisely, but their main disadvantage is that the mathematical analysis of theirs properties (even only quantitative) is hardly possible. Moreover, such large complex models contains a large number of parameters and estimation of them usually is a very challenging task. However, all these models include the p53-Mdm2 complex which consists of three main components: the p53 and Mdm2 proteins and Mdm2 mRNA. Thus, following [14] we consider a generalised mathematical model of p53-Mdm2 with two Hill functions, which describe the influence of Mdm2 protein concentration on the rate of degradation of p53 and the influence of p53 protein concentration on the rate of transcription of Mdm2 mRNA, respectively. The paper is organized as follows. In Section 1.1 we formulate a generalisation of the p53-Mdm2 model proposed by Monk. In Section 2 we prove analytically the existence of positive steady state and provide the conditions guaranteeing the occurrence of the Hopf bifurcation. Moreover, we investigate the stability of new born periodic solutions. We also analyse the direction of the Hopf bifurcation for the p53-Mdm2 system without time delay by computing the first Lyapunov coeffcient. Finally, in Section 3 we present the numerical simulations illustrating our analytical results and discuss our findings. 1.1. Model presentation In the presented paper we consider a generalisation of p53-Mdm2 model, which original version was proposed by Monk in [14]. We assume that: the p53 protein is produced with a constant rate α p and it is degraded with the rate μ p1 , the p53 degradation rate is increased by the Mdm2 presence and it is described by the term with function fp with the maximal Mdm2-induced degradation rate μ p2 . The Mdm2 mRNA production is upregulated by p53 and in the considered system it is described by the term αm1 fr ( p(s − τ˜ )) assuming the maximal initiation rate at the level of αm1 and taking into account the delay of the transcript elongation, splicing and processing, as well as the protein transport time between the cytoplasm

28

M.J. Piotrowska et al. / Applied Mathematics and Computation 328 (2018) 26–44

Fig. 1. A diagram of the p53-Mdm2 interactions modelled by system (1.1).

and the nucleus. Moreover, degradation rate for mRNA and the Mdm2 protein are denoted by μm1 and μm2 , respectively. Parameter αm2 describes the synthesis rate of the Mdm2 protein. In the considered model we focus only on the main interaction i.e., a negative feedback loop in the p53 compartment. We do not consider any post-translational modifications or other proteins that may have an effect on the p53-Mdm2 interaction system. For simplicity we also assume that a continuous stress is acting on the system (reflecting some physiological stresses such as hypoxia), and as a result of this we do not model a DNA damage in particular. The relationships described above are schematically illustrated in Fig. 1. The considered model has the following form

p˙ (s ) = α p − (μ p1 + μ p2 f p (m(s )) ) p(s ), r˙ (s ) = αm0 + αm1 fr ( p(s − τ˜ ) ) − μm1 r (s ),

(1.1)

m˙ (s ) = αm2 r (s ) − μm2 m(s ), where p, r and m denote the concentrations of the p53 protein, Mdm2 mRNA and Mdm2 protein, respectively. The switching functions fp and fr are assumed to be continuous and strictly increasing with value zero at zero and converging to 1 at +∞. In [14, supplementary data], the particular functions fp and fr were considered

f p (m ) =

k2m

m2 , + m2

f r ( p) =

knp

pn , + pn

(1.2)

where n ∈ N is the Hill coefficient which describes the character of the binding process. It should be noted that in [14] the value of n was not indicated. Additionally, it is assumed that all the model parameters are positive. To close the model we postulate the non-negative initial condition

φ = (φ1 (t ), φ2 (t ), φ3 (t )) ∈ C ([−τ˜ , 0], [0, +∞ )3 ), where C (I, E ) denotes the set of continuous functions defined on the interval I with values in E. 2. Model analysis To simplify the mathematical analysis we rescale model (1.1) in a proper manner which is, however, differs from the one proposed by Monk. It turns out that the rescaling proposed by us gives a simpler form of a positive steady state thus, it is more convenient. First we prove the existence and uniqueness of the positive steady state of system (1.1). 2.1. Positive steady state ¯ ) of system (1.1). Moreover Proposition 2.1. There exists exactly one non-negative steady state ( p¯ , r¯, m

αp αp < p¯ < . μ p1 + μ p2 μ p1

(2.1) μm

¯ . From the first two equations of (1.1) we obtain a non-linear system Proof. From the last equation of (1.1) we have r¯ = αm 2 m 2 ¯ of two equations for p¯ and m

αp p¯

¯ ), = μ p1 + μ p2 f p ( m

μm 1 μm 2 ¯ = αm0 + αm1 fr ( p¯ ). m αm2

(2.2)

M.J. Piotrowska et al. / Applied Mathematics and Computation 328 (2018) 26–44

29

From the first one we calculate

p¯ =

αp . μ p1 + μ p2 f p (m¯ )

(2.3)

Plugging (2.3) it into the second one we obtain

  αp μm 1 μm 2 ¯ − αm0 = αm1 fr . m αm2 μ p1 + μ p2 f p (m¯ )

(2.4)

¯ that varies from −αm0 (for m ¯ = 0) to +∞. Note that the left-hand side of (2.4) is an increasing continuous function of m Because the functions fp and fr are strictly increasing and continuous, the right-hand side of (2.4) is decreasing. Thus, there ¯ ) < 1 for any m ¯ > 0, exists exactly one non-negative solution to (2.4). Inequalities (2.1) follow from the fact that 0 < f p (m and equality (2.3). This completes the proof.  ¯ be the positive solutions to (2.2). We use the following rescaling which is very convenient for the analysis Let p¯ and m of the model dynamics:

p , p¯

x=

r

y=

,

ραm1

μx = ρμ p1 , η = ρμ p2 ,

z=

m , ¯ m

τ=

τ˜ , ρ

t=

s

ρ

,

αx =

ρα p p¯

α=

,

μy = ρμm1 , μz = ρμm2 , ρ =

αm0 , αm 1

¯ m

αm1 αm2

(2.5) ,

and

¯ z ), f x (z ) = f p ( m

fy (x ) = fr ( p¯ x ).

After this rescaling system (1.1) reads

x˙ (t ) = αx − (μx + η fx (z(t )) )x(t ), y˙ (t ) = α + fy (x(t − τ )) − μy y(t ),

(2.6)

z˙ (t ) = y(t ) − μz z(t ). Proposition 2.1 implies that system (2.6) has the unique positive steady state (1, μz , 1) and the following holds

αx = μx + η fx (1 ),

α + f y ( 1 ) = μy μz .

(2.7)

Note that after proposed scaling (2.5) functions considered in [14] read

f x (z ) =

z2

β 2 + z2

and

f y (x ) =

xn

γ n + xn

,

(2.8)

¯ and γ = k p / p¯ . where β = km /m 2.2. Existence, uniqueness and non-negativity of the solutions Let



 αx α α+1 α α+1 = (x, y, z ) ∈ R : ≤x≤ , ≤y≤ , ≤z≤ . μx + η μx μy μy μy μz μy μz 3

αx

Theorem 2.2. The solutions to (2.6) exist, are unique, and defined for all t ≥ 0. The set







C = (x, y, z ) ∈ C [−τ , 0], R3 : (x(t ), y(t ), z(t ) ) ∈ , t ∈ [−τ , 0]



is invariant under the evolution of system (2.6). Proof. The right hand side of (2.6) is Lipschitz continuous, which implies the local existence and uniqueness of the solutions. It is easy to see that if x(t ) = 0, then x˙ (t ) = αx > 0, thus x is non-negative for all t for which it is defined. Similar arguments yield the non-negativity of y(t) and z(t). Because the inequalities 0 ≤ fx (z) ≤ 1 hold for all z ≥ 0 and x ≥ 0, we have the following estimates for the first equation of (2.6)

αx − (μx + η )x(t ) ≤ x (t ) ≤ αx − μx x(t ), which gives the following bounds on the function x for all t ≥ 0



min x(0 ),

αx

μx + η





≤ x(t ) ≤ max x(0 ),

αx . μx

Similarly, the assumption made on the function fy gives estimations

α − μy y(t ) ≤ y (t ) ≤ α + 1 − μy y(t ),

(2.9)

30

M.J. Piotrowska et al. / Applied Mathematics and Computation 328 (2018) 26–44

which imply



min y(0 ),

α μy



 ≤ y(t ) ≤ max y(0 ),

 α+1 . μy

(2.10)

The upper bound for y yields the estimation for z. This proves that C is invariant. From the existence of the invariant set C we easily deduce the global existence of the solutions to system (2.6).



Remark 2.3. In fact, the upper bound for y and z can be improved if we use directly the upper bound for x. That is it can be proven that the set









˜ , t ∈ [−τ , 0] , C ˜ = (x, y, z ) ∈ C [−τ , 0], R3 : (x(t ), y(t ), z(t ) ) ∈ where

     α + fy μαx +x η α + fy μαxx α + fy μαx +x η α + fy μαxx α x ˜ = (x, y, z ) ∈ R : ≤x≤ , ≤y≤ , ≤z≤ μx + η μx μy μy μy μz μy μz 3

αx

is invariant under the evolution of system (2.6). 2.3. Stability of the positive steady state First we prove a technical lemma that we use in the future. Lemma 2.4. Assume that the characteristic function for a steady state reads

W (λ ) = P (λ ) + be−λτ = λ3 + a2 λ2 + a1 λ + a0 + be−λτ with positive constants ai , b, i = 0, 1, 2. Assume also that the auxiliary function defined by

F ( ω ) = P ( i

√

ω )2 − b2 = (a0 − a2 ω )2 + ω (a1 − ω )2 − b2

(2.11)

has exactly one simple positive root ω0 . If the steady state is stable for τ = 0 then the steady state loses its stability due to the Hopf bifurcation at

1

τcr = √ arccos ω0

a ω − a  2 0 0 b

.

Proof. Due to the assumptions of the lemma the stability change and the existence of the Hopf bifurcation is obvious in the context of Theorem 1 from [32]. Clearly, the critical τ cr fulfils the following system

cos(τcr

√

ω0 ) =

a2 ω0 − a0 , b

sin(τcr

√

ω0 ) =

a1 − ω0 √ ω0 . b

√ To complete the proof it is enough to show that sin(τcr ω0 ) > 0. Due to the fact that F has one simple positive root ω0 , we deduce that F(ω) < 0 for 0 < ω < ω0 and F(ω) > 0 for ω0 > ω. The assumption of the stability of the steady state for τ = 0 implies that a2 a1 > a0 + b holds. This is equivalent to the condition F (a1 ) = (a2 a1 − a0 )2 − b2 > 0 implying ω0 < a1 , which completes the proof.  Clearly, condition a0 < b is a necessary one for the function F to have exactly one simple positive root however, it is not a sufficient one. Of course one is able to pose a direct form of the sufficient condition, however the final formula would be complex and not manifesting any clear biological interpretation. Below we formulate a theorem that gives conditions for the stability of the steady state and the existence of the Hopf bifurcation. Let define

F (ω ) = ω 3 + ω 2



  μ2y + μ2z + αx2 + ω αx2 μ2y + μ2z + μ2y μ2z + αx2 μ2y μ2z − η2 ζ22 ,

(2.12)

where

ζ1 = (μy + μz )(μy μz + αx2 ) + αx (μ2y + μ2z ),

ζ2 = fx (1 ) fy (1 ).

(2.13)

Theorem 2.5. (i) If ηζ2 > 2αx μy μz + ζ1 , then the steady state (1, μz , 1) of system (2.6) is unstable for all τ ≥ 0. (ii) If ηζ 2 < α x μy μz , then the steady state (1, μz , 1) of system (2.6) is locally asymptotically stable for τ ≥ 0. (iii) If αx μy μz < ηζ2 < 2αx μy μz + ζ1 , then there exists a critical value τ cr > 0 such that the steady state (1, μz , 1) of system (2.6) is locally asymptotically stable for 0 ≤ τ < τ cr and it is unstable for τ > τ cr . Moreover, for τ = τcr , the Hopf bifurcation √ occurs and periodic solutions arise with period close to 2π / ω0 , where ω0 is a unique simple positive root of the function F defined by formula (2.12), and

 (μy + μz + αx )ω02 − αx μy μz τcr = arccos . ω0 ηζ2 1



(2.14)

M.J. Piotrowska et al. / Applied Mathematics and Computation 328 (2018) 26–44

31

Proof. Linearising system (2.6) around the steady state (1, μz , 1) we get

x˙ (t ) = −(μx + η fx (1 ) )x(t ) − η fx (1 )z(t ),

y˙ (t ) = −μy y(t ) + fy (1 )x(t − τ ),

(2.15)

z˙ (t ) = y(t ) − μz z(t ). Using definition (2.13) of ζ 2 we write the characteristic function for the steady state in the following form



W (λ ) = − (αx + λ )(μy + λ )(μz + λ ) + ηζ2 e−λτ





=−

λ3 + λ2 (αx + μy + μz ) + λ(μy μz + αx (μy + μz ) ) + αx μy μz + ηζ2 e−λτ ,

(2.16)

where (2.7) holds since (1, μz , 1) is the steady state of system (2.6). To determine the stability of the steady state for τ = 0 we check the Routh–Hurwitz conditions for characteristic function given by (2.16) obtaining

 (μy + μz ) μy μz + αx (μy + μz ) + αx2 > ηζ2 .

(2.17)

Inequality (2.17) is equivalent to 2αx μy μz + (μy + μz )(μy μz + αx2 ) + αx (μ2y + μ2z ) > ηζ2 or, according to definition (2.13), equal to 2αx μy μz + ζ1 > ηζ2 . Hence point (i) is proved. Calculating the auxiliary function F given by (2.11) we obtain the function F defined by (2.12) which is a polynomial of the third degree. Applying Descartes’ Rule of Signs, we immediately deduce that F has a real positive root if and only if

ηζ2 > αx μy μz

(2.18)

holds and in such a case the root is simple. This yields point (ii). Point (iii) follows directly from Lemma 2.4.



For the particular forms of fr and fp considered in [14] and given explicit by (1.2) the conditions from Theorem 2.5 can be rewritten as conditions for n. To do that first we use the following lemma. Lemma 2.6. Parameters of system (2.6) with fx and fy given by (2.8) fulfil

2β 2



β2 + 1

2 ·

nγ n

(γ n + 1 )2

= ζ2 ,

α < μy μz < α + 1 ,

μx < αx < μx + η,

(2.19)

where ζ 2 is given by (2.13). Proof. Because of the scaling and the fact that (1, μz , 1) is the steady state of system (2.6) the parameters of model (2.6) fulfil the following identities

αx − μx =

η , β2 + 1

μy μz − α =

1

γn +1

.

(2.20)

Due to identities (2.20) and the assumption that all parameters are positive we immediately obtain inequalities (2.19). Calculating β 2 and γ n from (2.20) and plugging them into (2.13) we arrive at equality (2.19), which completes the proof.  Remark 2.7. A direct calculation of fx (1 ) fy (1 ) leads to the following formula

ζ2 =

(γ

n

2 nβ 2 γ n = 2n(η − αx + μx )(αx − μx )(1 + α − μy μz )(μy μz − α )η−2 . + 1 )2 ( β 2 + 1 )2

(2.21)

For the particular forms of fx and fy functions given by (2.8), that is when we consider the scaled system, and the results from Lemma 2.6 we rewrite the conditions given in Theorem 2.5. Clearly, the value of ζ 2 depends on the Hill coefficient n of the function fy . It turns out, that there are two significant thresholds for this parameter, namely

n1 =

αx μy μz η

2(η − αx + μx )(αx − μx )(1 + α − μy μz )(μy μz − α )

,

n2 =

(2αx μy μz + ζ1 ) n1 , αx μy μz

(2.22)

where ζ 1 is given by (2.13). For n < n1 , the positive steady state is stable for all τ ≥ 0 (compare Theorem 2.5 (ii)). For n > n2 the positive steady state is unstable for all τ ≥ 0 (compare Theorem 2.5 (i)), while for n1 < n < n2 , it is stable for τ < τ cr , unstable for τ > τ cr , and at τ = τcr the Hopf bifurcation occurs (compare Theorem 2.5 (iii)). We illustrate that outcome graphically in Fig. 2. Theorem 2.5 determines stability condition for the dimensionless model. It turns out that for n large enough the steady state of system (2.6) is unstable for all τ ≥ 0. However, the nondimensionalization procedure causes that the parameter n of ¯ thus in the dimensionless parameters. Therefore, in order to determine the the original model (1.1) is hidden in p¯ and m stability of the steady state for the parameters considered by Monk [14] and for large n we write the conditions (i)–(iii) of Theorem 2.5 in terms of original parameters and for the particular forms of functions fp and fr given by (1.2). In order to emphasise that p¯ depends on the parameter n in the following we write here p¯ n . Thus, till the end of this subsection we consider a sequence { p¯ n } of the second coordinates of the steady state of system (1.1) for n = 0, 1, 2, . . . Clearly, considering

32

M.J. Piotrowska et al. / Applied Mathematics and Computation 328 (2018) 26–44

Fig. 2. Stability of the steady state (1, μz , 1) of system (2.6) for functions fx and fy given by (2.8), where n1 , n2 are given by (2.22).

n = 0 we neglect the nonlinear influence of p53 protein on the production of Mdm2 mRNA assuming that it is constant and protein concentration independent in result obtaining a simple linear equation. However, such a case alow us to derive coherent conditions used later on, thus from now on we also consider n = 0. Let us define

( p¯ n ) = C

α

p

p¯ n

 n 

2

− μ p1

n



p¯ n kp

1+

 n  p¯ n kp

αm0 + (αm0 + αm1 )

 n p¯ n kp

3 ,

C =

2αm1 μ3m1 μ3m2 k2m

μ p2 αm2 2

.

(2.23)

Proposition 2.8. For the non-scaled system (1.1) the condition (i) of Theorem 2.5 is equivalent to



( p¯ n ) > (μm1 + μm2 )

μm 1 μm 2 +

αp  αp p¯ n

p¯ n

+ μm 1 + μm 2



,

(2.24)

while the condition (ii) of Theorem 2.5 reads

α p μm 1 μm 2

( p¯ n ) <

p¯ n

.

(2.25)

Proof. The condition (i) of Theorem 2.5 in the original parameters and for the particular forms of functions fp and fr given by (1.2) reads

2αm1 αm2 μ p2

(

  nknp p¯ nn ¯ αp  αp k2m m > ( μm 1 + μm 2 ) μm 1 μm 2 + + μm 1 + μm 2 . n n 2 2 2 p¯ n p¯ n ¯ ) (k p + p¯ n ) +m

k2m

(2.26)

¯ and p¯ fulfil (2.2), which for the functions fp and fr given by (1.2) reads The coordinates of the steady state m

αp

μm 1 μm 2 p¯ n ¯ = αm0 + αm1 n n n , m αm2 k p + p¯ n

p¯ n

= μ p1 + μ p2

k2m

¯2 m . ¯2 +m

(2.27)

¯ 2 from the second equality of (2.27) and introduce it into the first equality of (2.27) obtaining We calculate m

Ch2 h2 ( p¯ n ) = where

 h( p¯ n ) =

 αm0 + αm1

p¯ nn n k p + p¯ nn

2

α p − μ p1 p¯ n , (μ p1 + μ p2 ) p¯ n − α p

,

(2.28)

Ch =

μm 1 μm 2 k m . αm2

(2.29)

¯ from (2.26) and we obtain that inequality (2.26) is equivalent to inequality (2.24). Using (2.27) we eliminate m The condition (ii) of Theorem 2.5 in terms of the original parameters reads

2αm1 αm2 μ p2

(

nknp p¯ nn ¯ k2m m < n ¯ 2 )2 (k p + p¯ nn )2 +m

α p μm 1 μm 2

k2m

p¯ n

.

(2.30)

Because the left-hand side of (2.30) is exactly the same as the left-hand side of (2.26), thus (2.30) is equvialent to (2.25).



For any ξ ∈ [0, 1] denote by g(ξ ) an unique solution (denoted by p¯ ξ ) of the equation

Ch2 h2 ( p¯ ξ ) = (αm0 + αm1 ξ ) , 2

that is

g( ξ ) =

αp

 μ2

μ2m2 k2m + (αm0 + αm1 ξ )2 αm2 2

m1

(μ p1 + μ p2 )(αm0 + αm1 ξ )2 + μ p1



μ2m1 μ2m2 k2m αm2 2

.

Proposition 2.9. If p¯ 1 < k p , then p¯ n < k p for all n ≥ 1. If p¯ 1 > k p , then p¯ n > k p for all n ≥ 1. If p¯ 1 = k p , then p¯ n = k p for all n ≥ 1. Moreover, the sequence p¯ n is convergent to p¯ ∞ which is given by

⎧ ⎪ ⎨max {g(1 ), k p } if p¯ 1 > k p , if p¯ 1 = p¯ 2 = · · · = p¯ n = · · · = k p , p¯ ∞ = k p = g(1/2 ) ⎪ ⎩ min {g(0 ), k p } if p¯ 1 < k p .

M.J. Piotrowska et al. / Applied Mathematics and Computation 328 (2018) 26–44

33

Fig. 3. The illustration of the proof of Proposition 2.9. The dashed curve is a graph of the left-hand side of (2.28), while the solid curves are the graphs of the right-hand side of (2.28) for different values of n. Solutions p¯ 0 , p¯ 1 and p¯ 20 (i.e., the terms of the sequence p¯ n ) to Eq. (2.28) are also indicated.

Proof. For fixed n the coordinate p¯ n of the steady state is a solution to (2.28). Here, we study the behaviour of a sequence { p¯ n }. In order to avoid confusion in this proof we consider quality (2.28) with p¯ instead of p¯ n . The left-hand side of this equation is an decreasing function of p¯ for α p /(μ p1 + μ p2 ) < p¯ < α p /μ p2 and the graph of this function (as a function of p¯ ) does not depend on n. On the other hand, the graph of the right-hand side of (2.28) changes with the change of n in such a way that for any fixed p¯ , the value of the right-hand side of (2.28) at the point p¯ is an increasing sequence if p¯ < k p , and decreeing one for p¯ > k p (see Fig. 3). Thus, for n → +∞ the right-hand side of (2.28) tends to (αm0 + αm1 )2 if p¯ > k p , tends 2 for p ¯ < k p and is equal to (αm0 + αm1 /2 )2 for p¯ = k p . Thus if for n = 1 p¯ < 1 holds, then p¯ n remains below kp for all to αm 0 2 or to k for the solution larger than n ≥ 0 and tends either to the solution of (2.28) with the right-hand side equal to αm p 0 kp . This completes the proof. 

Remark 2.10. Observe that calculating the term ( p¯ n /k p )n from (2.28)



p¯ n kp

n

=

Ch h( p¯ n ) − αm0 αm0 + αm1 − Ch h( p¯ n )

and substituting it into (2.23), we obtain the alternative formula of the function ( p¯ n ) namely

( p¯ n ) =

n C

2 Ch2 αm 1

α

p

p¯ n

− μ p1



μ p1 + μ p2 −

αp  p¯ n

 (αm0 + αm1 − Ch h( p¯ n ) ) 1 −

αm0  . Ch h( p¯ n )

(2.31)

Corollary 2.11. The following statement are true: (i) if p¯ ∞ < k p then for n large enough the steady state of (1.1) is stable for all τ ≥ 0, (ii) if p¯ ∞ ≥ k p then for n large enough the steady state of (1.1) is unstable for all τ ≥ 0. Proof. Note that if p¯ n → p¯ ∞ < k p , then ( p¯ n ) converges to 0 while the right-hand side of (2.24) remains positive. This yields that the condition (ii) of Theorem 2.5 is fulfilled for large n. On the other hand, if p¯ n → p¯ ∞ > k p , the situation is reverse. The sequence {( p¯ n )} tends to +∞ and the condition (i) of Theorem 2.5 is fulfilled for large n. If p¯ n → k p as n → +∞ we can not determine a priori the limit of ( p¯ n /k p )n . However, using Remark 2.10 we deduce that ( p¯ n ) ∼ C n as n → +∞, while the right-hand side of (2.24) is asymptotically constant. This implies the condition (i) of Theorem 2.5 is fulfilled for large n.  In our paper we consider only a single delay (τ ) in the term that describes production of Mdm2 mRNA upregulated by p53. However, if we consider time lag (σ ) in the protein production from mRNA i.e., the time lag in translation process described by term fp (m(s)) of system (1.1), then analogous calculus to that presented in the paper leads to the characteristic function determining the stability of steady state of the same form as in Theorem 2.5, but with τ + σ instead of τ . Thus, the results regarding the change of the stability of the steady state would be similar, however they would depend on the sum of the delays. On the other hand, whenever we would like to study the type of the existing Hopf bifurcation the calculations would differ and an additional parameter σ (appearing separately from τ ) would be present in a nonlinear part of the system. 2.4. Analysis of the type of the Hopf bifurcation with increasing delay In this section we derive analytical formulae allowing to investigate the stability of new born periodic solutions following Diekmann et al. approach, [33]. First we need to introduce notation. We define

dx = fx (1 ),

dy = fy (1 ).

For ω0 being a unique positive root of F defined by (2.12) and τ cr given by (2.14), we define

ϑx = αx2 + ω02 ,

ϑy = μ2y + ω02 ,

ϑz = μ2z + ω02 ,

(2.32)

34

M.J. Piotrowska et al. / Applied Mathematics and Computation 328 (2018) 26–44

and

σ = ϑy ϑz + (αx − iω0 )(ϑy (μz + iω0 ) + ϑz (μz + iω0 ) ) + τcr ϑy ϑz (αx − iω0 ).

(2.33)

Moreover

η 2 dx , αx μy μz + ηdx dy  −1 w3 = (αx + 2iω0 )(μy + 2iω0 )(μz + 2iω0 ) + ηdx dy e−2iω0 τcr , w2 =

(2.34)

and

ρc I , 1 = − ρc I , 2 =

η3 dx2  2

η 2 dx ϑx 2

3αx2 + ω02

ϑy ϑz (1 + αx τcr ) + 3αx (μz ϑy + μy ϑz )ϑx + ω02 ϑx (ϑz + ϑy ) + 2ω02 αx τcr ϑy ϑz ,

(ϑy ϑz (αx + τcr ϑx ) + ϑx (μy ϑz + μz ϑy )),

(2.35)

and

ρcII ,1 = αx ϑy ϑz + μy ϑx ϑz + μz ϑx ϑy + τcr η2 dx2 dy2 ,  ρcII ,3 = ηdx2 (1 + αx τcr )ϑy ϑz + αx (μz ϑy + μy ϑz ) + ω02 (ϑy + ϑz ) ,

(2.36)

and

κ11 = w3 (αx + iω0 )(μy + 2iω0 )(μz + 2iω0 ), κ12 = κ11 (αx + iω0 ), κ13 = w3 e

−2iω0 τcr

κ21 = κ13 (αx + iω0 ), κ22 = κ21 (αx + iω0 ),

(2.37)

κ23 = κ13 (αx + 2iω0 ).

,

Now we formulate the main result regarding the stability of the periodic solutions arising due to the Hopf bifurcation (i.e., when the assumptions of Theorem 2.5 (iii) are satisfied). Theorem 2.12. Let assumptions of Theorem 2.5 (iii) be fulfilled. If c˜ > 0, then the Hopf bifurcation is supercritical and if c˜ < 0, then it is subcritical, where



a c˜ = a0 + a1 fx (1 ) + a2 fy (1 ) + [ fx (1 ), fy (1 )] 11 a12 a0 = w2



a12 a22





fx (1 ) + a3 fx (1 ) + a4 fy (1 ), fy (1 )

ρcII ,1 2η2 αx μy μz dx3 − 2η2 αx dy dx2 ρcII ,3 + η4 dx4 (Re(κ11 σ (αx − iω0 )) − ηdy Re(κ21 σ ) ),

a1 = ρcI ,1 + w2 (ηdx (2ηαx dy dx − ϑx μy μz )ρcII ,1 + ηϑx dy ρcII ,3 ) + η4 dx3 dy Re(κ21 σ (αx − iω0 )) −



a2 = w2 1 2

+



2η2 αx μy μz dx3 − η3 dx4 dy

1 3 2 η dx (Re(κ12 σ (αx − iω0 )) − ηdy Re(κ22 σ ) ), 2

  2 2 ρcII ,1 − αx η dx ρcII ,3



2η4 dx4 Re(κ11 σ (αx + iω0 ) ) − η5 dx4 (dx Re(κ13 σ (αx − iω0 )) + Re(κ23 σ ) ) , dy

a 3 = ρc I , 2 , a4 =

η2 dx3 dy

ρc I , 2 ,

a11 = −w2 ηϑx dy ρcII ,1 − η3 dx dy Re(κ22 σ (αx − iω0 )), a12

w2 = 2 +

a22 = −

  dx ηϑx μy μz 2 2 αx η dx − ρcII ,1  2

η 3 dx 2

η3 dx4 dy

dy

ηdx Re(κ23 σ (αx − iω0 )) −

ρcII ,1 −

η5 dx5 dy



1 Re(κ12 σ (αx + iω0 )) , dy

Re(κ13 σ (αx + iω0 )).

M.J. Piotrowska et al. / Applied Mathematics and Computation 328 (2018) 26–44

35

Proof. We additionally substitute x˜=x − 1, y˜=y − μz , z˜=z − 1 obtaining the system which has a trivial steady state (0, 0, 0). This system rewritten in a functional form reads

x˜˙ = L1 (x˜t , y˜t , z˜t ) + G1 (x˜t , y˜t , z˜t ), y˜˙ = L2 (x˜t , y˜t , z˜t ) + G2 (x˜t , y˜t , z˜t ), z˜˙ = L3 (x˜t , y˜t , z˜t ) + G3 (x˜t , y˜t , z˜t ), where x˜t , y˜t , z˜t are x˜t (θ ) = x˜(t + θ ), y˜t (θ ) = y˜(t + θ ), z˜t (θ ) = z˜(t + θ ) for −τ ≤ θ ≤ 0 and Li and Gi are a linear and nonlinear parts of the right hand side of the shifted system (2.6), respectively, as follows

L1 (x˜t , y˜t , z˜t ) = −αx x˜t (0 ) − ηdx z˜t (0 ),

L2 (x˜t , y˜t , z˜t ) = dy x˜t (−τ ) − μy y˜t (0 ),

L3 (x˜t , y˜t , z˜t ) = y˜t (0 ) − μz z˜t (0 ),

   αx − μx (x˜t (0 ) + 1 ) − η fx (1 + z˜t (0 )) 1 + x˜t (0 ) − dx z˜t (0 ) ,  G2 (x˜t , y˜t , z˜t ) = α − μy μz + fy 1 + x˜t (−τ ) − dy x˜t (−τ ), G3 (x˜t , y˜t , z˜t ) = 0. G1 (x˜t , y˜t , z˜t ) =



To study the stability of new born periodic solutions i.e., the type of the Hopf bifurcation existing for τ = τcr one needs to investigate the sign of the expression

μ2 =

Re(c ) , Re(qD2 (iω0 , τcr )p )

(2.38)

where p, q are chosen such that  ,  = 1, and the characteristic matrix reads

 λ + αx

(λ, τ ) = −dy e−λτ 0

η dx

0

λ + μy −1

0

λ + μz



,

D2 denotes the derivative of a function with respect to the second variable, c = cI + cII + cIII ,

1 3 ¯ , qD G(0, τcr )(, , ) 2 1  cII = qD21 G(0, τcr ) ¯ (·, 0 ),  , cI =

cI I I =

 1 2 ¯ , qD1 G(0, τcr )  (·, 2iω0 ),  2

(2.39)

where Di1 , i = 2, 3 denote the derivatives of a considered function of the ith order with respect to the first variable and the function 1 : R × C → C3 is given by

1 (θ , a ) = eaθ (a, τcr )−1 D21 G(0, τcr )(, 1 ), ¯ or 1 =  for cII or cIII , respectively. with 1 =  The functions  ,  : R → C3 are defined as

(θ ) = qeiω0 θ

(θ ) = eiω0 θ p,

 ,  = qD1 (iω0 , τcr )p,

where p and q are eigenvectors of (iω0 , τ cr ), i.e., (iω0 , τcr )p = 0 and q(iω0 , τcr ) = 0. Clearly, because det (iω0 , τcr ) = 0 we have

−ηdx dy e−iω0 τcr = (αx + iω0 )(μy + iω0 )(μz + iω0 ).

(2.40)

The derivative of  with the respect to the first variable reads



D1 (λ, τ ) =

1

τ dy e−λτ 0

0 1 0



0 0 . 1

Thus, the left and right eigenvectors of (iω0 , τ cr ) reads



p=

−ηdx



(μz + iω0 )(αx + iω0 ) , αx + iω0



T (μy + iω0 )(μz + iω0 ) −ηdx q˜ = . −ηdx (μy + iω0 )

Now we calculate the normalisation constant ζ = q˜ D1 (iω0 , τcr )p such that q = q˜ ζ¯ /ζ 2 obtaining

  ζ¯ = −ηdx ν − τcr ηdx dy eiω0 τcr ,

where

ν = (μy − iω0 )(μz − iω0 ) + (αx − iω0 )(μz − iω0 ) + (αx − iω0 )(μy − iω0 ).

(2.41)

36

M.J. Piotrowska et al. / Applied Mathematics and Computation 328 (2018) 26–44





Now we calculate the denominator of (2.38) that is b = Re qD2 (iω0 , τcr )p , where



0 λdy e−λτ 0

D2 (λ, τ ) =



0 0 0



0 0 . 0



Therefore, b = dy Re iω0 e−iω0 τcr q2 p1 and because we are interested only in the sign of c we write

b=

η2 dx2 dy  η2 dx2 dy ω0  −iω0 τcr ¯ −iω0 τcr ¯ Re i ω e ζ = − Im e ζ . 0 ζ 2 ζ 2

Due to form of ζ¯ we have













Im e−iω0 τcr ζ¯ = −ηdx Im e−iω0 τcr ν − τcr ηdx dy = −ηdx Im e−iω0 τcr ν . Equality (2.40), the fact that ω0 , η, dx , dy > 0 and definition (2.41) imply the following





sign(b) = − sign Im (αx + iω0 )(μy + iω0 )(μz + iω0 )ν





= ω0 (μ2y + ω02 )(μ2z + ω02 ) + (αx2 + ω02 )(μ2z + ω02 ) + (αx2 + ω02 )(μ2y + ω02 ) = ω0





 ϑx ϑy + ϑx ϑz + ϑy ϑz > 0.

Thus, the sign of b does not depend on the values of the model parameters and for τ = τcr the inequality b < 0 holds. Hence it remains to calculate the numerator of (2.38), which strongly depends on the non-linear part of the system G = [G1 , G2 , G3 ]T . Let u, v, w : [−τ¯ , 0] → R3 be test functions, where τ¯ is sufficiently large i.e. τ¯ > τcr . Let φ = (x˜t , y˜t , z˜t ). We have

 αx − μx − η fx (1 + z˜t (0 )) u1 (0 )   − η (1 + x˜t (0 )) fx (1 + z˜t (0 )) + dx u3 (0 ),   D21 G1 (φ , τ )(u, v ) = −η fx (1 + z˜t (0 )) u1 (0 )v3 (0 ) + u3 (0 )v1 (0 ) D11 G1 (φ , τ )u =



− η (1 + x˜t (0 )) fx (1 + z˜t (0 ))u3 (0 )v3 (0 ),



D31 G1 (φ , τ )(u, v, w ) = −η fx (1 + z˜t (0 )) u1 (0 )v3 (0 )w3 (0 ) + u3 (0 )v1 (0 )w3 (0 )



+ u3 (0 )v3 (0 )w1 (0 ) − η (1 + x˜t (0 )) fx (1 + z˜t (0 ))u3 (0 )v3 (0 )w3 (0 ) and

D11 G2 (φ , τ )u = ( fy (1 + x˜t (−τ )) − dy )u1 (−τ ), D21 G2 (φ , τ )(u, v ) = fy (1 + x˜t (−τ ))u1 (−τ )v1 (−τ ), D31 G2 (φ , τ )(u, v, w ) = fy (1 + x˜t (−τ ))u1 (−τ )v1 (−τ )w1 (−τ ),

D1j G3 (φ , τ ) = 0, j = 1, 2, 3.

Because all terms cI , cII and cIII have a common positive factor 1/ζ 2 , we define

c˜I = ζ 2 cI ,

c˜II = ζ 2 cII ,

c˜III = ζ 2 cIII .

We calculate c˜I :





¯ = −η fx (1 ) p1 p3 p¯ 3 + p3 p1 p¯ 3 + p3 p3 p¯ 1 − η fx (1 )p3 p3 p¯ 3 , D31 G1 (0, τcr )(, , )



= η2 dx fx (1 ) 2ϑx +



αx2 − ω02 + 2iαx ω0



− ηϑx fx (1 )(αx + iω0 )

¯ = fy (1 )p1 p1 p¯ 1 e−iω0 τcr = −η3 dx3 fy (1 )e−iω0 τcr , D31 G2 (0, τcr )(, , )

1

Re(c˜I ) = Re

2

  η2 dx3   ¯ ζ¯ q˜ D31 G(0, τcr )(, , ) = ρcI ,1 fx (1 ) + ρcI ,2 fx (1 ) + f y (1 ) , dy

where ρcI ,1 and ρcI ,2 are defined by (2.35). Now, we proceed to calculate c˜II . First note that

(0, τcr )

−1

1 = αx μy μz + ηdx dy

 μy μz μz d y dy

−ηdx

αx μz αx



−ημy dx −ηdx dy .

αx μy

(2.42)

M.J. Piotrowska et al. / Applied Mathematics and Computation 328 (2018) 26–44

¯ obtaining Now we calculate D21 G(0, τcr )(, )

 

37





¯ = −η dx p1 p¯ 3 + p3 p¯ 1 + fx (1 )p3 p¯ 3 = 2η2 αx dx2 − η fx (1 )ϑx , D21 G1 (0, τcr )(, ) ¯ = fy (1 )p1 p¯ 1 = η2 dx2 fy (1 ). D21 G2 (0, τcr )(, ) ¯ = [κII,1 , κII,2 , κII,3 ]T , where Now, let define κII = (αx μy μz + ηdx dy )(0, τcr )−1 D21 G(0, τcr )(, )

  κII,1 = ημy μz 2ηαx dx2 − fx (1 )ϑx − η3 dx3 fy (1 ),   κII,3 = ηdy 2ηαx dx2 − fx (1 )ϑx + αx η2 dx2 fy (1 ).

(2.43)

Clearly, due to the form of the derivatives of the non-linear parts evaluated on test functions  we do not need κ II,2 . Moreover, κII,1 , κII,3 ∈ R and the function ¯ (θ , 0 ) = κII is constant. Direct calculations of D21 G(0, τcr ) ¯ (·, 0 ),  leads to













D21 G1 (0, τcr ) ¯ (·, 0 ),  = −η dx κII,1 + fx (1 )κII,3 (αx + iω0 ) + η2 dx2 κII,3 , D21 G2 (0, τcr ) ¯ (·, 0 ),  = −ηdx fy (1 )κII,1 e−iω0 τcr . Note that applying the norm to the both sides of (2.40) we obtain ϑx ϑy ϑz = η2 dx2 dy2 . Hence the following identities hold

  ζ¯ (μy + iω0 )(μz + iω0 ) = −ηdx ϑy ϑz 1 + τcr (αx − iω0 )   + (αx − iω0 ) ϑy (μz + iω0 ) + ϑz (μy + iω0 ) ,

ζ¯ e−iωτcr = −ηdx ν e−iω0 τcr + τcr η2 dx2 dy ν = (αx + iω0 )(μy + iω0 )(μz + iω0 ) + τcr η2 dx2 dy , dy

and



−ηζ¯ (μy + iω0 )(μz + iω )(αx + iω0 ) = η2 dx (αx + iω0 )ϑy ϑz + (μy + iω0 )ϑx ϑz



+ (μz + iω0 )ϑx ϑy + τcr η2 dx2 dy2 . Thus, we have

(c˜II ) = w2



   d κ ρcII ,1 κII,3 fx (1 ) + x II,1 fy (1 ) + dx κII,1 −κII,3 ρcII ,3 ,

(2.44)

dy

where ρcII , j , j = 1, 3 are given by (2.36) and w2 is given by (2.34). Now we proceed to the most complex part, namely c˜III . First, we calculate

D21 G1 (0, τcr )(, ) = −2ηdx p1 p3 − η fx (1 )p23 = 2η2 dx2 (αx + iω0 ) − η fx (1 )(αx + iω0 )2 , D21 G2 (0, τcr )(, ) = fy (1 )p21 e−2iω0 τcr = η2 dx2 fy (1 )e−2iω0 τcr . Now we calculate (2iω0 , τcr )−1 . We note that

det (2iω0 , τcr ) = (αx + 2iω0 )(μy + 2iω0 )(μz + 2iω0 ) + ηdx dy e−2iω0 τcr = w−1 3 yielding that

(2iω0 , τcr )−1 = w3



(μy + 2iω0 )(μz + 2iω0 ) dy (μz + 2iω0 )e−2iω0 τcr dy e−2iω0 τcr

−ηdx

(αx + 2iω0 )(μz + 2iω0 ) αx + 2iω0



−ηdx (μy + 2iω0 ) −ηdx dy e−2iω0 τcr . (αx + 2iω0 )(μy + 2iω0 )

(2.45)

Denote now κIII = [κIII,1 , κIII,2 , κIII,3 ]T = (2iω0 , τcr )−1 D21 G(0, τcr )(, ). Thus,  (θ , 2iω0 ) = e2iω0 θ ((2iω0 , τcr ))−1 D21 G(0, τcr )(, ) = e2iω0 θ κIII . Hence,

  κIII,1 = w3 2η2 dx2 (αx + iω0 ) − η fx (1 )(αx + iω0 )2 (μy + 2iω0 )(μz + 2iω0 ) − η3 dx3 fy (1 )e−2iω0 τcr ,    κIII,3 = w3 dy 2η2 dx2 (αx + iω0 ) − η fx (1 )(αx + iω0 )2 + η2 dx2 fy (1 )(αx + 2iω0 ) e−2iω0 τcr .

As before we do not need κ III,2 . However, this time κ III in not a real vector. We rewrite the expression for κ III, j as

κIII,1 = 2η2 dx2 κ11 − ηκ12 fx (1 ) − η3 dx3 κ13 fy (1 ), κIII,3 = 2η2 dx2 dy κ21 − ηdy κ22 fx (1 ) + η2 dx2 κ23 fy (1 ),

(2.46)

38

M.J. Piotrowska et al. / Applied Mathematics and Computation 328 (2018) 26–44

where κ ij are given by (2.37).  ¯ obtaining Now we calculate D21 G(0, τcr )  (·, 2iω0 ), 







¯ = −η dx (κIII,1 p¯ 3 + κIII,3 p¯ 1 ) + fx (1 )κIII,3 p¯ 3 D21 G1 (0, τcr )  (·, 2iω0 ), 

 







= −η dx (κIII,1 (αx − iω0 ) − ηdx κIII,3 ) + fx (1 )κIII,3 (αx − iω0 ) ,

¯ = fy (1 ) (−τ , 2iω0 )( ¯ −τ ) D21 G2 (0, τcr )  (·, 2iω0 ),  = fy (1 )κIII,1 p¯ 1 e−iω0 τcr = −ηdx fy (1 )κIII,1 e−iω0 τcr . Thus,

Re(c˜III ) = where

w3 2





fx (1 )ρcIII ,1 + fy (1 )ρcIII ,2 + ρcIII ,3 + ρcIII ,4 ,

(2.47)

  ρcIII ,1 = η2 dx Re κIII,3 σ (αx − iω0 ) ,  η2 dx2  ρc I I I , 2 = Re κIII,1 σ (αx + iω0 ) , dy   2 2 ρcIII ,3 = η dx Re κIII,1 σ (αx − iω0 ) ,   ρcIII ,4 = −η3 dx2 Re κIII,3 σ .

Now, using (2.42), (2.44), and (2.47) as well as definition (2.43) and (2.46) we obtain the formula for c˜ = Re(c ) · ζ 2 that have obviously the same sign as Re(c ), which completes the proof.  Clearly, due to the complexity of the obtained formula we are unable to find a reasonable analytic formula for the sign of Re(c ). However, we use direly the above calculations to investigate the type of the existing stability numerically in Section 3. 2.5. The type of the Hopf bifurcation in the case without delay Now we consider the p53-Mdm2 model without time delay. In order to analyse the direction of the Hopf bifurcation in this case we compute an analytical expression for the first Lyapunov coefficient l1 (see [34]). We observe that for τ = 0 system (2.6) reads

x˜˙ = L1 (x˜, y˜, z˜) + G1 (x˜, y˜, z˜), y˜˙ = L2 (x˜, y˜, z˜) + G2 (x˜, y˜, z˜),

(2.48)

z˜˙ = L3 (x˜, y˜, z˜) + G3 (x˜, y˜, z˜), where

L1 (x˜, y˜, z˜) = −αx x˜(t ) − ηdx z˜(t ), L2 (x˜, y˜, z˜) = dy x˜(t ) − μy y˜(t ), L3 (x˜, y˜, z˜) = y˜(t ) − μz z˜(t ),

   αx − μx (x˜(t ) + 1 ) − η fx (1 + z˜(t )) 1 + x˜(t ) − dx z˜(t ) ,  G2 (x˜, y˜, z˜) = α − μy μz + fy 1 + x˜(t ) − dy x˜(t ), G1 (x˜, y˜, z˜) =



G3 (x˜, y˜, z˜) = 0. Theorem 2.13. If

(αx + μy + μz )(μy μz + αx (μy + μz )) = αx μy μz + ηdx dy ,

(2.49)

holds, then the Hopf bifurcation occurs for system (2.48). Moreover, for l1 < 0 the Hopf bifurcation is supercritical, while for l1 > 0 it is subcritical, where l1 is given by

l1 =

1 dx ν2 ω0





a a0 + a1 fx (1 ) + a2 fy (1 ) + [ fx (1 ), fy (1 )] 11 a12

a12 a22



fx (1 ) fy (1 )





+ a3 fx (1 ) + a4 fy (1 ) ,

(2.50)

M.J. Piotrowska et al. / Applied Mathematics and Computation 328 (2018) 26–44

with

39

  αx η2 dx3 (μy + μz ) 2αx λ3 ω02 − ηdx dy (μy + μz ) + μy μz (3αx2 + ω02 ) λ3   2 2 3  η − dx (μy + μz )(αx2 + ω02 ) λ3 3αx − 2(μy + μz ) − λ3 ω02 2 2 2 ω0 λ3 + 4ω0   + αx2 (μy + μz ) − λ3 αx (μy + μz ) + 2αx ω02 (μy + μz ) + 2ω04 , 4

a0 = −

a1 = ω02 ηdx (μy + μz )(αx2 + ω02 )(μy + μz − 2αx )

  ηdx (μy + μz )(αx2 + ω02 ) ηdx dy (3αx − λ3 ) − μy μz (αx2 + ω02 ) λ3    2   2 η dx 2 2 +  ( μ 3λ3 3ω0 + αx (μy + μz ) + 2(2ω02 − λ23 )(μy + μz ) , y + μz ) αx + ω0 3 λ23 + 4ω02 2

−

      η3 dx4 6 λ23 + 4ω02 (μy + μz ) −2αx2 + αx λ3 − ω02 + 2αx μy μz 2 2 3λ3 (λ3 + 4ω0 )      2 2 2 2 2 2 + 3λ3 3αx (μy + μz ) + 3ω0 2αx + μy + μz + 2λ3 (2ω0 − λ3 ) 2αx (μy + μz ) − 3ω0 ,

a2 = −

a3 = (μy + μz )(αx2 + ω02 )2 ω02 , a4 = η3 dx4 ω02 and

2

η (μy + μz )(αx2 + ω02 )2 ,      2  1 η2 dx2 3 2 2 α a12 = − η2 dx2 ω02 − 2μy μz + 3 λ ( μ + μ ) + 3 ω + α ( μ + μ ) − ω α y z x y z 3 x x 0 0 λ3 6(λ23 + 4ω02 )   − 2ω02 (2ω02 − λ23 ) 3αx + μy + μz ,   η4 dx5 λ23 + 8ω02 . a22 = 2 2 λ3 (λ3 + 4ω0 ) a11 =

λ3

Proof. In the order to check the existence of the Hopf bifurcation for (0, 0, 0), we analyse the linear system

⎡ ⎤

 x˜˙ −αx ⎣y˜˙ ⎦ = dy 0 z˜˙

0 −μy 1

−ηdx 0 −μz

  x˜ y˜ z˜

(2.51)

for which the characteristic polynomial reads

W (λ ) = λ3 + λ2 (αx + μy + μz ) + λ(αx μy + αx μz + μy μz ) + αx μy μz + ηdx dy .

(2.52)

Clearly, W(λ) has a pair of purely imaginary roots ± iω0 (ω0 > 0) if condition (2.49) is satisfied and as there is a change of stability, the Hopf bifurcation occurs. To study the direction of the bifurcation for τ = 0 (see [34]) one needs to investigate the sign of the first Lyapunov expression

l1 =

1 2ω0



Re

    v, C (w, w, w¯ ) − 2 v, B(w, A−1 B(w, w¯ )) + v, B(w¯ , (2iω0 I3 − A )−1 B(w, w )) ,

where A is the Jacobian matrix of (2.48), B = [B1 , B2 , B3 ]T , C = [C1 , C2 , C3 ]T ,

Bi (x, y ) =

" " 3 ! ∂ 2 Gi ( ξ ) " ∂ 3 Gi ( ξ ) " x j yk , Ci (x, y, z ) = x y z, " ∂ ξ j ∂ ξk ξ =0 ∂ ξ j ∂ ξk ∂ ξl "ξ =0 j k l j,k=1 j,k,l=1 3 !

x = [x1 , x2 , x3 ]T , y = [y1 , y2 , y3 ]T , z = [z1 , z2 , z3 ]T and

    (μy − iω0 )(μz − iω0 ) −ηdx ν¯ −ηdx v=− , w = (μz + iω0 )(αx + iω0 ) ηdx ν2 −ηdx (μy − iω0 ) αx + iω0

with ν = −2ω02 + 2iω0 λ3 and λ3 = −(αx + μy + μz ) is the third eigenvalue of (2.52).

(2.53)

M.J. Piotrowska et al. / Applied Mathematics and Computation 328 (2018) 26–44

τ˜ = 15 1 mRNA

concentration (rescaled)

40

p53 Mdm2

0.5

1.5 1 0.4

0.5 0

0.8 0

200 400 time in minutes

600

1 p53

0.2

1.2

0

Mdm2

τ˜ = 32 2

p53 Mdm2

1.5

2 mRNA

concentration (rescaled)

Fig. 4. Trajectories of solutions of (2.6) as functions of time. Parameters given by (3.1), τ˜ = 15. Values of concentration rescaled as in (2.5) excluding time.

1 0.5 0

1 0.6

0

200

400 600 800 time in minutes

1,000

0.4 0.8

1

p53

1.2

1.4

0.2

Mdm2

Fig. 5. Trajectories of solutions of (2.6) as functions of time. Parameters given by (3.1), τ˜ = 32. Values of concentration rescaled as in (2.5) excluding time.

Using similar techniques as in the proof of Theorem 2.12 we conclude that the Lyapunov coefficient l1 has the form (2.50).  Remark 2.14. Note that (2.49) is equivalent to

2αx μy μz + ζ1 = αx μy μz ,

(2.54)

because using the equality (2.22) for τ = 0 we obtain condition (2.54). Remark 2.15. In the case of the p53-Mdm2 model without time delay it is more convenient to use formula (2.53) rather than (2.38). Formula (2.53) from Kuznetsov [34] allows as immediately determine the occurrence of stable oscillations while formula (2.38) from Diekmann et al. [33] requires additional calculations, i.e., the analysis of sign of the denominator in expression (2.38). It follows from the fact that the coefficient l1 of formula (2.53) is equal to c, which is given by (2.39), multiplied be a positive constant. 3. Numerical results and discussion To perform numerical simulations we take the parameters proposed, and partially estimated based on the biological data, in [14], that is:

μ p1 =

ln 2 −1 min , 100

μ p2 =

τ˜ = 15 min, α = 0.01,

93 ln 2 −1 min , 700 kp

αp

=

km

αm1 αm0

μm 1 = μm 2 = = 100,

n = 4.

ln 2 −1 min , 20

αm1 = 0.1, αm2 (3.1)

In [14] the particular value n = 4 was used to perform numerical simulations showing a good agreement with experimental data, see [14, Fig. 5]. Monk compared the model dynamics with the experimental data for irradiated mammalian cell lines, and claimed qualitative agreement between them although the experimental data were not shown, [14]. We decided to plot solutions of the rescaled system (2.6) and parameters corresponding to those considered by Monk, but without time rescaling in order not to loose the time scale. Thus, we give the delay values and the threshold for not rescale time that is τ˜ . In the left panels of Figs. 4 and 5 the trajectories in the space (x, y, z) are plotted. The value of the solution at time t = 0 is indicated by a dot.

M.J. Piotrowska et al. / Applied Mathematics and Computation 328 (2018) 26–44

41

Fig. 6. Trajectories of solutions of (2.6) as functions of time. Parameters given by (3.1), but with n = 10 and k p = 62.7059. In this case the steady state is unstable for τ = 0. Values of concentration rescaled as in (2.5) excluding time.

Fig. 7. Illustration of stability of the steady state and the stability of the Hopf bifurcation for model (1.1) with parameters given by (3.1) and varying Hill coefficient n.

For τ˜ = 15 minutes used by Monk [14] the amplitude of oscillations is decreasing, because time delay is sufficiently small and the steady state is asymptotically stable, see Fig. 4. For the parameters given by (3.1) the critical value of delay is τ˜cr ≈ 31.825 min and for that value the positive stable steady state loses its stability due to the Hopf bifurcation and the oscillations appear with a period around 187.65 min, see Fig. 5. In Fig. 5 (right panel) we show a limit cycle. As we have shown in Theorem 2.5, the steady state, for some parameters values, can be unstable for τ = 0. Due to the existence of an invariant set the solutions need to tend to some attractor. In our numerical experiments we get a limit cycle as it is shown in Fig. 6. Clearly, this does not exclude the possibility of the existence of strange attractors. This is a problem that we plan to address in a future. The plots presented in Fig. 6 were obtained for τ˜ = 0 and τ˜ = 100, n = 10 and k p = 62.7059. The rest of the parameters were as in (3.1). It is worth pointing out that in [14] no particular range for the Hill coefficient n was postulated. Our analysis summarized in Fig. 7 shows how essential the right choice of the degree of cooperativeness of nuclear import and binding of the protein to its promoter is. Nevertheless, derived conditions (2.22) for scaled model (2.6) in the case on non-scaled system (1.1) for fp and fr functions proposed by Monk become strongly implicit as it was shown in Proposition 2.8. However, for the particular parameters we are able to derive some results. For the parameters used by Monk (see (3.1)) and Corollary 2.11 we deduce that the steady state is unstable for large n. Now we estimate how large n should be. Using definition (2.29) of the function h we obtain

(μ p1 + μ p2 ) p¯ − α p 1 = = α p − μ p1 p¯ h2 ( p¯ )

μ p1 + μ p2 − αp − μ p1 p¯

Thus, we rewrite the function ( p¯ n ) as

nC 2 Ch2 αm 1

·

α

p

p¯

− μ p1



μ p1 + μ p2 −

αp   p¯

·

αp p¯

.

αm0 + αm1 − Ch h( p¯ )



1−

αm0  . Ch h( p¯ )

Clearly, the multiplication of the first two brackets decreases for p¯ > 2α p /(2μm1 + μm2 ) ≈ 19.24. On the other hand, we see that the multiplication of the last two brackets decreases for Ch2 h2 ( p¯ ) > αm0 (αm0 + αm1 ). Because the function h is decreasing, this means that this multiplication increases for p¯ > p˜ ≈ 22.3. Thus we may estimate the function ( p¯ n ) from the below, for n > 30 as

n·

α

p

p¯ ∞

− μ p1



μ p1 + μ p2 −

α p 

p¯ ∞

αm0 + αm1 − Ch h( p¯ 30 )



1−

αm0  > n · 1.69 · 10−4 . Ch h( p¯ 30 )

42

M.J. Piotrowska et al. / Applied Mathematics and Computation 328 (2018) 26–44

Fig. 8. Dependence of critical value of the Hill coefficient on the maximal Mdm2-induced p53 degradation rate (μ p2 ) and the p53 production rate (α p ). Value β represents a percentage change of α p according to its reference value given in (3.1).

On the other hand, the right-hand side of (2.24) is a decreasing function of p¯ n and thus we estimated it from above (for n > 30) by





μm 1 + μm 2



μm 1 μm 2 +

αp  αp p¯ 30

p¯ 30

+ μm 1 + μm 2



< 1.48 · 10−4 .

This shows, that for n > 30 condition (i) of Theorem 2.5 is fulfilled. Additionally, fixing the parameters as (3.1) (except n) and directly computing the stability conditions for each n from interval [1, 40] we conclude that for n = 1 the positive steady state of (1.1) is locally stable independently on the value of the bifurcation parameter (delay) and it is unstable for n ∈ {20, . . . , 40}. However, for arbitrary n from {2, . . . , 19} there always exists single Hopf bifurcation for positive the steady state, compare Fig. 7. Thus, for n ∈ {2, . . . , 19}, changing τ , one is always able to obtain damped oscillatory behaviour for sufficiently small delay and undamped oscillations otherwise, thus reflect the biologically observed phenomena. From the biological point of view the obtained results show how important the assumption regarding the shape of the function fr describing the Mdm2 mRNA production unregulated by p53 is. If we assume that fr looks like the Michaelis– Menten function, then we would not be able to obtain experimentally observed sustained oscillations of the protein level even if time delay in mRNA production is large. On the other hand, for fr being the Hill function with n > 1 the greater steepness (equivalent to higher n) of “switch” is the smaller delay in the Mdm2 mRNA production is required to obtain the oscillatory periodic behaviour. Clearly, the condition determining the type of the Hopf bifurcation presented in Theorem 2.12 is so complicated, that even for the rescaled system it is difficult to see any clear stability dependence on the model parameter. Thus, for the model parameter proposed by Monk (except n), we numerically check this condition obtaining supercritical Hopf bifurcation for all n ∈ {2, . . . , 19}. Additionally, again using the derived in Theorem 2.12 formula and parameters given by (3.1), we investigate dependence of the existence and the type of the Hopf bifurcation in (n, μ p2 ) and (n, α p ) parameter spaces, see Fig. 8, showing that this time for the considered ranges of parameters we always observe supercritical Hopf bifurcations. In general we distinguish two types of the Hopf bifurcation: sub- and supercritical (considered as a “safe” bifurcation). In the case of our model the difference between them would be reflected in the stability of the new born periodic solutions. For the subcritical bifurcation they would be unstable while for the supercritical new born periodic solutions would be stable. Clearly, whenever the stationary solution become unstable there is a risk that solutions of the system, i.e., the concentrations described by the model, would increase in unrestricted manner which is unrealistic from the biological point of view. However, we have shown that solutions to (2.6) would not leave the compact invariant set, thus such situation is impossible. Although such situation is not possible the instability of the stationary state (a point or a periodic solution) might lead to irregular, or even chaotic, behaviour of solutions. In Fig. 9 we present dependence of the critical value of the bifurcation (left) and periods of corresponding oscillatory solutions of system (1.1) varying parameters μ p2 and α p while rest of parameters were fixed as in (3.1). Value β represents a percentage change of parameter α p according to its reference value given in [14]. More precisely, β = +5% means the increase of the value by 5%, while β = −5% the decrease by 5%. Clearly, both τ cr and period are the decreasing functions of μ p2 parameter, see Fig. 9. Similarly, with the increase of α p the critical value of the bifurcation parameter and period decrease. That implies that the increase of α p leads to the reduction of the stability region of the positive steady state of system (1.1), while the increase of μ p2 enlarge that region. Thus, decrease of p53 production rate prevent oscillations of the concentrations of p53 and Mdm2 proteins, while the influence of the and Mdm2-induced p53 degradation rate is reverse i.e., its increase promotes oscillatory behaviour. Clearly, in our study we have used parameters proposed by Monk [14] and thus, our results are in agreement with experiments reported by Bar-Or et al. in [3]. On the other hand, results presented in Fig. 9 indicate that the considered model can also reproduce the solutions with period of 325 min reported by Geva-Zatorsky et al. [15], compare also Fig. 10

M.J. Piotrowska et al. / Applied Mathematics and Computation 328 (2018) 26–44

43

350

kp = 200

300 250

kp = 150

200

kp = 100 0.1

0.15 0.2 μp 2

0.25

80

kp = 200

τcr in minutes

period in minutes

Fig. 9. Dependence of the critical value of time delay τ cr and period of the oscillating solutions on the maximal Mdm2-induced p53 degradation rate (μ p2 ) and the p53 production rate (α p ). Value β represents a percentage change of α p according to its reference value given in (3.1). The red balls indicate values for the parameter set (3.1) originally considered by Monk in [14]. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

60 kp = 150

40

kp = 100 0.1

0.15 0.2 μp 2

0.25

Fig. 10. Dependence of the critical value of time delay τ cr and the period of the oscillating solutions on the maximal Mdm2-induced p53 degradation rate (μ p2 ).

where an increase of the parameter kp (chosen arbitrarily by Monk) up to 200 leads to increase of the oscillations period up to 5 h (compare [35]). In the considered model only the negative feedback loop mediated by Mdm2 protein is taken into account, while the positive feedback loop mediated by PTEN is neglected. For that reason the phenomenon (i.e., existence of to distinguishable thresholds of p53 level that determine the cell fate) described in the paper by Kracikova et al. [36] does not directly appear in our considerations. However, it should be noted that our findings are not in contradictory to the observations made by Kracikova et al. Clearly, changing model parameters, we change the values of the steady state. Thus, we may interpret that effect in the following way: an activation of the positive feedback loop mediated by PTEN leads to such change in the model parameters that solutions increase, and eventually an apoptotic threshold would be reached. Acknowledgements This work was partially supported by grant no. 2015/19/B/ST1/01163 of National Science Centre, Poland. References [1] C. Culmsee, X. Zhu, Q.-S. Yu, S.L. Chan, S. Camandola, Z. Guo, N.H. Greig, M.P. Mattson, A synthetic inhibitor of p53 protects neurons against death induced by ischemic and excitotoxic insults, and amyloid β -peptide, J. Neurochem. 77 (1) (2001) 220–228, doi:10.1046/j.1471-4159.2001.00220.x. [2] R.S. Morrison, Y. Kinoshita, M.D. Johnson, W. Guo, G.A. Garden, p53-dependent cell death signaling in neurons, Neurochem. Res. 28 (1) (2003) 15–27, doi:10.1023/a:1021687810103. [3] R. Bar-Or, M. Maya, L. Segel, U. Alon, A. Levine, M. Oren, Generation of oscillations by the p53-Mdm2 feedback loop: a theoretical and experimental study, PNAS 97 (20 0 0) 11250–11255. [4] D. Lane, p53, guardian of the genome, Nature 358 (1992) 15–16. [5] B. Vogelstein, D. Lane, A. Levine, Surfing the p53 network, Nature 408 (20 0 0) 307–310, doi:10.1038/35042675. [6] D. Michael, M. Oren, The p53-Mdm2 module and the ubiquitin system, Semin. Cancer Biol. 13 (1) (2003) 49–58, doi:10.1016/S1044-579X(02)0 0 099-8. [7] J. Piette, H. Neel, V. Marechal, Mdm2: keeping p53 under control, Oncogene 15 (1997) 1001–1010. [8] C. Prives, Signaling to p53: breaking the MDM2-p53 circuit, Cell 95 (1) (1998) 5–8, doi:10.1016/S0 092-8674(0 0)81774-2. [9] J. Momand, H.-H. Wu, G. Dasgupta, MDM2 - master regulator of the p53 tumor suppressor protein, Gene 242 (1–2) (20 0 0) 15–29, doi:10.1016/ S0378-1119(99)00487-4. [10] K. Ryan, A. Phillips, K.H. Vousden, Regulation and function of the p53 tumor suppressor protein, Curr. Opin. Cell Biol. 13 (3) (2001) 332–337, doi:10. 1016/S0955-0674(0 0)0 0216-7. ´ [11] K. Jonak, M. Kurpas, K. Szoltysek, P. Janus, A. Abramowicz, K. Puszynski, A novel mathematical model of ATM/p53/NF-κ B pathways points to the importance of the DDR switch-off mechanisms, BMC Syst. Biol. 10 (75) (2016) 1–12, doi:10.1186/s12918- 016- 0293- 0.

44

M.J. Piotrowska et al. / Applied Mathematics and Computation 328 (2018) 26–44

[12] X.-P. Zhang, F. Liu, W. Wang, Two-phase dynamics of p53 in the DNA damage response, PNAS 108 (22) (2011) 8990–8995, doi:10.1073/pnas. 110 060 0108. [13] S.L. Harris, A.J. Levine, The p53 pathway: positive and negative feedback loops, Oncogene 24 (17) (2005) 2899–2908. [14] N. Monk, Oscillatory expression of Hes1, p53, and NF-κ B driven by transcriptional time delays, Curr. Biol. 13 (2003) 1409–1413. [15] N. Geva-Zatorsky, N. Rosenfeld, S. Itzkovitz, R. Milo, A. Sigal, E. Dekel, T. Yarnitzky, Y. Liron, P. Polak, G. Lahav, U. Alon, Oscillations and variability in the p53 system, Mol. Syst. Biol. 2 (2006) 0033, doi:10.1038/msb4100068. [16] G. Lahav, N. Rosenfeld, A. Sigal, N. Geva-Zatorsky, A.J. Levine, M.B. Elowitz, U. Alon, Dynamics of the p53-Mdm2 feedback loop in individual cells, Nat. Genet. 36 (2004) 147–150, doi:10.1038/ng1293. [17] D.A. Ouattara, W. Abou-Jaoudé, M. Kaufman, From structure to dynamics: frequency tuning in the p53-Mdm2 network. II: Differential and stochastic approaches, J. Theor. Biol. 264 (4) (2010) 1177–1189, doi:10.1016/j.jtbi.2010.03.031. [18] K.E. Gordon, I.M. van Leeuwen, S. Lain, M.A. Chaplain, Spatio-temporal modelling of the p53-Mdm2 oscillatory system, Math. Model. Nat. Phenom. 4 (3) (2009) 97–116, doi:10.1051/mmnp/20094304. [19] M. Sturrock, A.J. Terry, D.P. Xirodimas, A.M. Thompson, M.A. Chaplain, Spatio-temporal modelling of the Hes1 and p53-Mdm2 intracellular signalling pathways, J. Theor. Biol. 273 (1) (2011) 15–31, doi:10.1016/j.jtbi.2010.12.016. [20] C.J. Proctor, D.A. Gray, Explaining oscillations and variability in the p53-Mdm2 system, BMC Syst. Biol. 2 (1) (2008) 75, doi:10.1186/1752- 0509- 2- 75. ´ [21] K. Puszynski, R. Bertolusso, T. Lipniacki, Crosstalk between p53 and nuclear factor-κ B systems: pro- and anti-apoptotic functions of NF-κ B, IET Syst. Biol. 3 (2009) 356–367, doi:10.1049/iet-syb.2008.0172. ´ [22] K. Puszynski, A. Gandolfi, A. d’Onofrio, The pharmacodynamics of the p53-Mdm2 targeting drug Nutlin: the role of gene-switching noise, PLoS Comput. Biol. 10 (12) (2014), doi:10.1371/journal.pcbi.1003991. ´ , B. Hat, T. Lipniacki, Oscillations and bistability in the stochastic model of p53 regulation, J. Theor. Biol. 254 (2) (2008) 452–465. [23] K. Puszynski [24] E. Batchelor, C.S. Mock, I. Bhan, A. Loewer, G. Lahav, Recurrent initiation: a mechanism for triggering p53 pulses in response to DNA damage, Mol. Cell 30 (3) (2008) 277–289, doi:10.1016/j.molcel.2008.03.016. [25] L. Ma, J. Wagner, J.J. Rice, W. Hu, A.J. Levine, G.A. Stolovitzky, A plausible model for the digital response of p53 to DNA damage, PNAS 102 (40) (2005) 14266–14271, doi:10.1073/pnas.0501352102. [26] G.I. Mihalas¸ , M. Neamt¸ u, D. Opris¸ , R. Horhat, A dynamic p53-Mdm2 model with time delay, Chaos Solitons Fractals 30 (4) (2006) 936–945. [27] A. Ciliberto, B. Novák, J.J. Tyson, Steady states and oscillations in the p53/Mdm2 network, Cell cycle 4 (3) (2005) 488–493. ´ [28] B. Hat, K. Puszynski, T. Lipniacki, Exploring mechanisms of oscillations in p53 and nuclear factor-κ B systems, IET Syst. Biol. 3 (5) (2009) 342–355, doi:10.1049/iet-syb.2008.0156. [29] T. Zhang, P. Brazhnik, J.J. Tyson, Exploring mechanisms of the DNA-damage response: p53 pulses and their possible relevance to apoptosis, Cell cycle 6 (1) (2007) 85–94. ´ ´ [30] K. Puszynski, R. Jaksik, A. Swierniak, Regulation of p53 by siRNA in radiation treated cells: simulation studies, Int. J. Appl. Math. Comput. Sci. 22 (4) (2012) 1011–1018, doi:10.2478/v10 0 06- 012- 0075- 9. [31] E. Kozłowska, K. Puszynski, Application of bifurcation theory and sirRNA-based control signal to restore the proper response of cancer cells to DNA damage, J. Theor. Biol. 408 (2016) 213–221, doi:10.1016/j.jtbi.2016.08.017. [32] K.L. Cooke, P. van den Driessche, On zeroes of some transcendental equations, Funkcj. Ekvacioj 29 (1986) 77–90. [33] O. Diekmann, S. van Giles, S. Lunel, Delay Equations, Springer-Verlag, New York, 1995. [34] Y.A. Kuznetsov, Elements of Applied Bifurcation Theory, Springer, New York, 2004, doi:10.1007/978- 1- 4757- 3978- 7. [35] E. Batchelor, A. Loewer, C. Mock, G. Lahav, Stimulus-dependent dynamics of p53 in single cells, Mol. Syst. Biol. 7 (488) (2011) 1–8, doi:10.1038/msb. 2011.20. [36] M. Kracikova, G. Akiri, A. George, R. Sachidanandam, S. Aaronson, A threshold mechanism mediates p53 cell fate decision between growth arrest and apoptosis, Cell Death Differ. 20 (2013) 576–588, doi:10.1038/cdd.2012.155.