Stability and Hopf bifurcation analysis for a Lac operon model with nonlinear degradation rate and time delay

Abstract: In this paper, we construct a discrete time delay Lac operon model with nonlinear degradation rate for mRNA, resulting from the interaction among several identical mRNA pieces. By taking a discrete time delay as bifurcation parameter, we investigate the nonlinear dynamical behaviour arising from the model, using mathematical tools such as stability and bifurcation theory. Firstly, we discuss the existence and uniqueness of the equilibrium for this system and investigate the effect of discrete delay on its dynamical behaviour. Absence or limited delay causes the system to have a stable equilibrium, which changes into a Hopf point producing oscillations if time delay is increased. These sustained oscillation are shown to be present only if the nonlinear degradation rate for mRNA satisfies specific conditions. The direction of the Hopf bifurcation giving rise to such oscillations is also determined, via the use of the so-called multiple time scales technique. Finally, numerical simulations are shown to validate and expand the theoretical analysis. Overall, our findings suggest that the degree of nonlinearity of the model can be used as a control parameter for the stabilisation of the system.


Introduction
Since the first mathematical model was developed by Brian in 1965 [10] and Joseph and co-workers [3,19,21] further developed the work, researchers have shown a growing interest in studying mathematical models of Genetic Regulatory Networks (GRNs) and/or gene expression.One of the reasons is that these models can be used to understand important aspects of living organisms and may help tackle illnesses influenced by genetic factors, such as genetic disorders, diabetes and cancer.Recent studies on gene expression and/or GRNs are rich and many significant results have been reported in the literature [21,23,28,37,39].These studies involve a number of different mathematical models and investigate multiple aspects of such networks, concentrating on biologically important themes such as network oscillations, effects of time delays, interactions between different constituents such as DNA, RNA and proteins, etc.
For example, the authors of Ref. [29] considered a typical GRNs model with a delay, and applied perturbation methods and centre manifold reduction to quantify how delay could drive oscillations in gene activity.Wang et al. [32] investigated the stability and Hopf bifurcation of a four-dimension GRNs model with double genes, whereas Zang and her collaborators [40] looked at the dynamics of a GRNs using normal form theory as proposed by Elowitz and Leibler [7].Further, Monk [21] developed a model for Hes1 with time delays, where transcription time delay describes the time needed to copy the information to a messenger RNA (mRNA) and translation time delay measures how long it takes to translate the protein from the ribosome.Other examples that constitute the background for the study we propose can be found, for instance, in Liu et al. [18], who introduced a model with nonlinear degradation using the interaction between mRNA and small non coding RNA (sRNA) and one protein, and showed how oscillations and bistability can be produced.This model was also investigated by Liu et al. in another paper [17], where they concluded that oscillations can be induced not just by interactions between RNAs and protein, but also by time delays during the diffusion or transportation of molecules in a cell.
Another interesting example is also given by the work of Yildirim and Mackey [39], who developed and analysed a nonlinear mathematical model of the Lac operon in bacterium Escherichia coli [27], which, in some cases, could be a strong contributor to colon carcinogenesis.Then, Nikolov et al. [23] proposed a model for the process of expression of B-galactosidase from a Lac Z gene, noting that a bifurcation associated with the presence of time delay causes the appearance of limit cycles in the model.Time delay is one the fundamental ingredients in these models, and will play a central role in the current work as well.This focus is justified because delays serve important functions in nature, and occur in a vast number of phenomena, some of which have been subjected to mathematical analysis.Examples are given by population dynamics [35], lake eutrophication ecosystems [44], insect outbreak dynamics [41], vegetation ecological system [12,43] and grazing ecosystem [42] to name a few.The key role of delay in all these cases is to be a driver of regime shifts between stable states or between oscillatory-non oscillatory regimes, with non trivial consequences for the behaviour of these systems.
A final, essential component of GRNs is mRNA, whose interactions with the other components of the system is essential to the dynamics of gene expression [1,4,8,9,20,36], and is included in a number of different types of genetic regulatory networks systems with time delay [24,30,31,34,45].The two major questions that we try and address in this work, and that are common to the formulations we have just referred to, are: 1. How do the degradation rates affect system dynamics?2. What will happen if mRNA breaks down into several parts?Because of these two questions, our modelling attempt includes an important new assumption describing the breakdown of mRNA into several identical parts.This, in turn, implies that a nonlinear degradation rate for mRNA is present, and has to be modelled accordingly.Our interest in a nonlinear degradation rate stems from the fact that the existing literature (for instance, see Refs.[5,23,38]) only partially describes the role of discrete time delay in destabilising the dynamics of GRNs.Most authors limit their analysis to the study of limit cycles (if present) when the delay is small.To expand the breadth of the analysis, we will instead (a) investigate the effect of the nonlinear degradation rate for mRNA for the stability of equilibrium solutions, and the amplitude and period of periodic solutions, and (b) provide insight into the behaviour of gene expression processes.
Our work starts with model development, which introduces a nonlinear degradation term for mRNA, and adds transcription and/or translation time delays to reflect the oscillatory phenomena observed experimentally.Mathematical tools, such as stability and bifurcation theory are then used to investigate the dynamics.Our results will show that different regimes exist and are influenced by time delay and nonlinear degradation.
This paper is organised as follows.The derivation and formulation of the Lac operon model is in Section 2. In Section 3, we study its steady state and the significance of the degree of nonlinearity for the model.In Section 4, we analyse the stability of the positive equilibrium with a nonlinear degradation rate under the effect of time delay, and give sufficient conditions to prove the existence and the direction of a Hopf bifurcation.In Section 5, some numerical simulations are discussed as a validation of the theoretical analysis and as a generalisation of the proven results.Discussion and a brief conclusion are given in Section 6, with Future directions being illustrated in Section 7.

Model derivation and formulation
A number of diseases generally indicated as dynamical diseases exist.These diseases are characterised by a dramatic change in the dynamics of a physiological variable [2].It is thus essential to understand the underlying mechanisms of genetic regulatory dynamics and account for different regulation scenarios mediated by several mRNA, which may drive such change.In particular and as mentioned above, oscillatory phenomena revealed by several mRNA regulation mechanisms can help understand the crucial role of mRNA in gene regulation.
Given all this, we introduce a new assumption to the Lac operon model first described in Ref. [23] and require that mRNA breakdown results in several identical parts and that these parts translate to enzyme B-galactosidase, as sketched in Figure 1.The process of transcription and translation of information is assumed to be regulated and represented by the following chemical reactions where H(p) is the Hill function, and the other variables will be introduced shortly.If the transcription factor is an activator, this function can be represented as where n is coefficient of the Hill function or the cooperativity, n 1 represents transcription factors and k = k 1 k 2 .On the other hand, if the transcription factor is a repressor, we have that  where n measures the cooperativity of the end product repression.Mathematically, the regulation can be modelled with three equations, as follow: where the temporal evolution of mRNA (m), enzyme B-galactosidase (y) and repressor protein (p) are described, and G(m) = k 3 m N , is the nonlinear degradation rate.The important constant N describes the degree of nonlinearity and we always have N > 0. As it is customary, the following dimensionless variables are introduced m = γ 1 m , y = γ 2 y , p = γ 3 p and t = γ 4 t, , Then, after dropping the prime and the hat, the system (2.2) becomes adimensional and the new equations read where σ 1 = k 3 γ N 4 , σ 2 = k 4 γ 4 , and σ 3 = k 5 γ 4 .As mentioned, time delay is one of the key features of the present paper.If we introduce transcription and translation delays in the process, we arrive at the final system where τ j , j = 1, 2, 3 are time delays.It is important to note that if N = 1, namely a linear degradation rate is considered, the model is the same as the one studied in Ref. [38], foe which the authors have shown that a Hopf bifurcation exists.

Solving for the steady state
The steady state of the model with no delay, i.e. equations (2.3), is determined by equating the right-hand sides of the equations to zero: Then, we obtain solutions y * = σ 3 p * and m * = σ 2 σ 3 p * , where p * is the root of equation < 0 and lim p→+∞ H(p) = +∞, implying that at least one positive root exists.A straightforward calculation shows that meaning that H(p) is monotonically increasing in the half-plane (0, +∞).As a result, the positive root of equation (3.1) is unique.This means that the uniqueness of the positive equilibrium of (2.3) has been proven.

Steady state of the model and significance of the parameter N
The differential equations of the model (2.3) could be numerically solved to explore the effect of the nonlinear degradation rate for mRNA.Looking at how steady solutions change as N is varied, we note that the positive equilibrium for variables m, y, and p, increases or decreases in a common fashion for all variables.So, the steady state of p, as a function of the parameter N, for example, for various values of n, is shown in Figure 2. From this diagram we can see that there is a range of N values i.e.N ∈ (0, 2], at which a decrease in protein level occurs until a larger value for the steady state is attained at larger N. Thinking about N as a bifurcation parameter, the plot shows that a regime switch, with the higher value of the steady state as a "on" position and the lower as a "off" position, can be driven by N. The diagram also suggests that the decrease occurs for both non-cooperative repression (n = 1), and cooperative repression (n > 1), in the end product repression.An essential prediction of the model is that this switching mechanism is only active when the degree nonlinearity is in the interval N ∈ (0, 2].When N > 2, then there is no role for the cooperativity of the end product repression and concentration of protein is stable, with no regime changes.The degree of nonlinearty N

Dynamical behaviour of model's equilibrium and delay
The dynamical behaviour of the equilibrium is now studied, in the absence and presence of delay, via an analysis of the stability and the Hopf bifurcation in the system.As customary practice, we linearize the model around its equilibrium and check the eigenvalues of the corresponding characteristic equation.Furthermore, a Hopf bifurcation can be observed from genetic expression as for previous Refs.[7,10,21], yielding helpful information about the stability of the periodic solution and its character.

Stability in the absence of delay
Without delay, i.e. with τ = 0, the characteristic equation ( 4.3) becomes To determine the sign of the real parts of roots of equation (4.4) and find their stability, we make the following, important hypothesis: (H 1 ) : The above assumption leads to the following lemma.Proof.The principal diagonal minors, ∆ i , of the Hurwitz matrix for equation (4.4) are: Notice that q i > 0 (i = 1, 2, and 3) and F 1 > 0. If hypothesis (H 1 ) holds, then all leading principal minors of the Hurwitz matrix are positive.Therefore, we conclude that the roots of (4.4) have negative real part, meaning that the unique positive equilibrium is stable.
Next, we will study the effect of discrete delay in the dynamical behaviour of equation (2.4), always assuming that condition (H 1 ) holds.

Time delay and oscillation
Assume that λ = iω(ω > 0) is a pure imaginary root of equation (4.3).If we substitute it into equation (4.3), we obtain Then we have from which we obtain If we now let u = ω 2 , we obtain where Noticing that p 1 > 0 and q > 0, it can be verified that, when r ≥ 0, equation (4.8) has no positive solution.Contrary, if r < 0, the equation has precisely one positive solution u 0 as G(0) = r < 0 and G (u) > 0 for u ∈ [0, +∞).Note also that r < 0 if and only if q 3 < F 1 .
For the root u 0 of equation (4.8), we now let ω * = √ u 0 , which is also known as the critical frequency, and obtain the corresponding critical time delays τ = τ * j , i.e.
At critical time delays τ * j , we have that the following lemma holds, which will help determine the conditions for oscillatory solutions.Then substituting λ = iω * , one has Re dλ dτ 1 ω * > 0, since G (u) > 0. This is completes the proof.
From the above analysis, Lemma 4.1 and Lemma 4.2, the following important result on the existence of oscillations holds.(i) all roots have negative real part if one of the following holds 2 and τ = τ * j , with j = 0, 1, . .., only purely imaginary roots λ = ±iω * exist; (iii) system (4.3) has at least one root with positive real part if τ ∈ (τ * 0 , +∞) and σ 1 < From the above Theorem, we can draw the the following conclusions.
Theorem 4.4.For system (2.4), (i) the equilibrium E * = (m * , y * , p * ) is stable if one of the following conditions holds Given these results, we will now determine the stability of the bifurcating periodic solution by using the so-called multiple time scale technique.

The direction of the Hopf bifurcation
The direction of the Hopf bifurcation is studied by applying the so-called multiple time scales technique.The direction (or type) can be either supercritical and subcritical [14]: a supercritical Hopf bifurcation leads to a stable branch of limit cycle with sustained oscillations, whereas the opposite occurs for the subcritical case.For simplicity of discussion, we assume τ 1 = τ 2 = τ 3 = τ and let t = t * τ, so that only a single, normalised time delay is present.After dropping the asterisks for notational simplicity, system (4.1)becomes where 3  , To determine the normal form of the centre manifold near the equilibrium of system (4.10),we introduce time scales T 0 = t and T 2 = 2 t, where | | 1 is a non-dimensional bookkeeping parameter [22,40].The derivative with respect to t can be written as and the solution of the system can be expanded in terms of as follows So, At τ * = τ * j , let τ = τ * + 2 δ, where δ > 0 is known as the detuning parameter [22,40], and substitute equations (4.11)-(4.13)into (4.10).After equating coefficients with same powers of , we obtain for where T and the overhead bar indicates complex conjugation.Then, using equation (4.17), we obtain A straightforward calculation shows that equation (4.15) has a particular solution of the following form: Substituting equation (4.18) into (4.15) and equating coefficients for the term e 2iω * τ * T 0 , we find where Also, equating coefficients for term BB, allows us to find that Substituting equations (4.17) and (4.18) into (4.16),we have where CC refers to the complex conjugate of the secular terms and NST represents the non-secular terms.As shown in Refs.[22,45], if a so-called "solvability condition" is satisfied, then equation (4.20) has a solution.Let us now assume that the particular solution of the above equation is of the form X 3 (T 0 , T 2 ) = φ(T 2 )e iω * τ * T 0 + CC.

Mathematical Biosciences and Engineering Volume 16, Issue 3, xxx-xxx
We can then obtain the following: (4.21) The solvability condition requires to find a vector d such that the next adjoint homogeneous equation is satisfied where Here, we now apply the following condition which makes d 1 unique and yields where the following coefficients have been found Finally, if we now multiply d T from the left to each side of equation (4.21), we obtain where and This allows us to arrive at the final result if we let B = ze −iω * τ * t with z = x + iy.In fact, we can now convert equation (4.24) into a more familiar normal form term, i.e.
where we used real coordinates.Using x = ρ cos θ and y = ρ sin θ, provides the polar form of (4.25) This last equation allows us to arrive at the stability of the bifurcating periodic solution and the direction of the Hopf bifurcation of system (4.1),studying the sign of Re(m 1 )Re(m 2 ) [14].We can summarise the final result in the following theorem.

Numerical simulations
In this section, numerical simulations are carried out to expand on our theoretical analysis, using Matlab c code (dde23) and the DDE-Biftool software.Based on published experimental and theoretical results such as those of Ref. [38], parameters are chosen as follows: σ 1 = 0.25, σ 2 = 0.1, σ 3 = 0.36 and a Hill coefficient of n = 2.Then, straightforward calculation shows that Given the results discussed in the previous section, the Hopf bifurcation induced by time delay occurs when N ∈ (0, 2) and, as a result, in our simulation we probe different values of N within this interval and outside, i.e.N = 0.5, 1, 1.5, 1.75, 2 and 10.The effect of the degree of nonlinearity N on the values of the critical delay is shown in Table 2.It is found that the critical value τ * increases with N, showing that τ * is highly sensitive to the change of the degree of nonlinearity.However, when N ≥ 2, the equilibrium E * remains stable for different values of delay.This is also understood by virtue of Theorem 4.4, since the positive equilibrium is always stable when τ = 0.In Figure 3, time series of system variables (and no time delay) for cases N = 0.5, 1.5, 2 and 10 are shown.The system reaches a stable equilibrium also for a time delay smaller than the critical value, namely τ ∈ (0, τ * ) with depicting a phase portrait for the case N = 1.75.In both cases, it is clear that after a short oscillating transient, the system is damped down to a stable plateau.When the delay is particularly large, on the other hand, the equilibrium is still stable.In particular, if , then even a very large delay gives rise to a stable equilibrium.Figure 5 shows a typical behaviour of this kind, for the case N = 2, τ 1 + τ 2 + τ 3 = 100.Note the different time scale from Figure 4 and the longer transient.
When the previously mentioned condition is not met, if time delay exceeds its critical value, the system shows a periodic solution, as depicted in Figure 6.Notably, the amplitude and the period of the periodic solution get larger as N is increased, as reported in Table 2, and shown in Figure 6.As expected, the Hopf bifurcation is supercritical and, as shown in Figure 7 where periodic solutions for N = 1.5, and N = 1.75 are plotted, the basin of attraction of the periodic branch includes parameter that can be less than those at which the Hopf point occurs.Studying the distribution of the roots of the characteristic equation using DDE-Biftool, we can find when the complex conjugate characteristic root pair crosses the imaginary axis and obtain numerical values for the parameters at which the Hopf bifurcation occurs (see Figure 8).From all the above results, we conclude that sustained oscillations are induced by time delay and its destabilising role in the dynamic of GRNs non trivially depends on the the degree of nonlinearity N. In other words, N can be used as a control parameter for the given GRNs.

Discussion and conclusion
The goal of this study is to explore the equilibrium and oscillatory dynamics of a chosen GNRs architecture, mediated by the interaction among several identical parts of mRNA.The main biological hypothesis in our approach is that the breakdown of the original mRNA in the Lac operon model is coherent with the framework introduced by Ref. [23].It should be noted that, in that work, the authors focussed on the role of cooperativity of the end product repression, and concluded that cooperativity has an essential role in the stability of equilibrium of the system.For example, they found that the magnitude of oscillations enlarges as the value of cooperativity is increased.In addition, authors in Ref. [5] concluded that discrete delay has a destabilising role in the dynamics of the Lac operon model.
Differently than the existing literature, in this model we have shown that the degree of nonlinearity plays also an essential role in the stability of the equilibrium.For example, we have shown that a degree of nonlinearity N ≥ 2, suppresses any effect of cooperativity and the equilibrium of the model remains stable.However, when N ∈ (0, 2), there is a notable decrease in the concentration of the protein, which can be further depressed as cooperativity is enhanced.An important key observation can be made from Figure 2 with regards to the nonlinear degradation rate N, for mRNA: an increase in the end product repression measure when N ∈ (0, 2), reduces the concentration of protein.This observation may suggest a possible way to monitor a change in the status of GRNs networks that can be indicative of an occurrence of genetic illnesses.
On the other hand, the use of hypothesis (H 1 ) warrants a stable, positive equilibrium for the system in the absence of time delays.If, on the contrary, a time delay is present, conditions have been found for the existence of a Hopf bifurcation, which produces oscillatory solutions.In the case of a homogeneous delay, the Hopf bifurcation occurs when τ is larger than a critical value and its direction has been found to be supercritical, implying a sustained oscillatory behaviour for the model.
Another interesting characteristic of our model is that limit cycles increase their amplitude as the degree of nonlinearity is increased, as Figure 6 shows, up until the limit value N = 2.When N ≥ 2, the oscillations are suppressed and the equilibrium of system (2.4) regains its stability.According to Ref. [1], wide-amplitude protein oscillations are estimated to coexist with the steady state solutions, indicating that GNRs tend to respond to significantly large perturbation with apoptosis.Our model may show that the occurrence of some diseases such as cancer could be related to proteins being in a so-called "off" positions, which dynamically correspond to the lowest protein steady states.

Further research
The model proposed in this work can be expanded to tackle more questions and increase the level of realism of the system.One of the interesting topics that has not been discussed in this work is the role of noise in the processes and results we have illustrated.This is particularly relevant as noise in biological systems is ubiquitous, and, as far as gene expression dynamics is concerned, most processes are stochastic.Both intrinsic and extrinsic noise [46] can play a part in our GRNs: intrinsic noise results from the probabilistic nature of chemical dynamics, whereas extrinsic noise occurs due to random fluctuations in environmental parameters, such as temperature and pressure.
In particular, there have been a number of interesting recent papers on gene dynamics [6,11,25] that describe the effect of noise on oscillations.In Ref. [11] the authors found that the cooperative binding of repressor molecules is not essential for the oscillatory behaviour and provided new insights into the nature of network oscillations.Cruz et al. [6] investigated different escape dynamics from a stable limit cycle in an oscillatory biochemical state, whereas Potoyan and Wolynes [25] contributed two new mathematical approaches to gain understanding in dephasing phenomena in gene oscillators, showing that essential behaviour arising from model bifurcations results from dichotomic noise from gene-state switching.In conclusion and in line with these studies, further research will consider the effect of noise onto our model and investigate the dynamical behaviour of the network and its oscillations.

Figure 1 .
Figure 1.Major biological assumption of the present model: LacZ gene transcribes to mRNA j , j = 1, 2, . . ., N, when the promoter is active.Then, mRNA j is translated to protein B-galactosidase, which has an effect on the synthesis of repressor protein.

Figure 2 .
Figure 2. Values for the steady state solution of variable p as a function of parameter N, for various values of the parameter n.

Lemma 4 . 1 .
If (H 1 ) holds, then the roots of the characteristic equation (4.4) have negative real part.

Theorem 4 . 3 .
Assume hypothesis (H 1 ) holds.Then, for the characteristic equation of system (4.3)we have the following:

Figure 7 .
Figure 7. Direction of the Hopf bifurcation: phase trajectories for N = 1.5 and N = 1.75, showing that the bifurcating periodic solution is stable and the Hopf point is supercritical.

Figure 8 .
Figure 8. Numerical proof of the existence of the positive roots for the characteristic equation of the system for different values of N. Proceeding from the top left panel, clockwise, plots are for N = 0.5, N = 1.0,N = 2 and N = 1.5.For every value N ≥ 2 note that no intersection is present, indicating that no Hopf bifurcation exists.

)
Mathematical Biosciences and EngineeringVolume 16, Issue 3, xxx-xxx where e 1 is the canonical basis.Since equation (4.14) is homogeneous, one can easily obtain the general solution as

Table 1 .
Bifurcation point τ * versus the degree of nonlinearity and amplitude and period of bifurcating periodic solution.