Bifurcation Analysis of a Nonlinear Genetic Network Model with Time Delay ()
1. Introduction
The basic functions of life are carried out by millions of cells, all of which are formed of thousands of proteins and genes interacting [1] . An important goal in the areas of applied mathematics and biology is trying to understand these network of interactions. Theoretically, the network represents how certain genes affect other genes [2] and, in the case when the network involves a small number of genes, its dynamic properties can be studied directly in a laboratory [3] [4] . However, in the case when the network is formed of a large number of genes, then understanding its dynamic behavior might be difficult and nonintuitive [5] [6] . In recent years, various computational techniques have been developed that help us understand some of the intrinsic properties of these networks. For an extensive review see [7] - [14] .
The biology behind the model in this work can be roughly summarized in a three step process (see Chapter 2 for more details): 1) in the nucleus a gene is copied onto messenger RNA (mRNA), 2) the mRNA then gets translated into a protein, and 3) the protein then returns to the nucleus where it represses production of mRNA (see Figure 1). This mRNA-protein feedback mechanism has been studied previously in several works [14] [15] . One of the more promising and relevant mRNA-protein models with feedback was first studied by Monk [15] and subsequently analyzed in more detail by a few others [14] [16] .
In this work we extend the author’s previous study [14] of a single mRNA-protein model to a network formed of N coupled mRNA-protein systems. The network structure in this work is described in Section 4 and the model is given by a system of N nonlinear DDEs coupled to N linear ODEs. Section 2 gives a summary of the biological background of our model. In Section 3 we present the linear analysis of a network of two coupled mRNA-protein systems. The network of N coupled mRNA-protein systems and its associated linear analysis is presented in Section 4 followed by its nonlinear analysis in Section 5. Section 6 presents a computational analysis and comparison of MATLAB’s numerical results with our perturbation results. In Section 7 we discuss our conclusions and future directions.
2. Biological Background
Transcription and translation are the main processes by which a cell expresses the instructions encoded in its genes. Transcription is the first step in gene expression and it includes the replication of a gene into messenger RNA (mRNA). The second step is the translation process, where the mRNA diffuses out of the nucleus and is captured, read, and translated by a ribosome yielding as a final product an associated protein. The final step in this feedback mechanism happens when the protein then diffuses back into the nucleus and represses production of its own mRNA. From these mechanisms, mRNA and protein concentrations arise naturally as the main intracellular regulatory agents in our model [1] [2] .
Figure 1. Diagram of the mRNA-protein feedback inhibition mechanism. The gene is copied onto mRNA, which then diffuses out of the nucleus and gets translated into a protein. The protein then diffuses back into the nucleus where it represses the transcription of its own mRNA.
There are several mechanisms that the cell uses to regulate levels of mRNA and protein concentrations. One important mechanism is the cell’s ability to adjust and control the breakdown rate (degradation rate) of old mRNA or proteins by increasing or decreasing the concentration of enzymes that degrade these. Another important example is the cell’s ability to regulate the production rate (transcription rate) of a specific mRNA. This is usually accomplished by a special type of proteins called transcription factors, which either promote or repress mRNA production by binding directly to the gene. Although other cellular mechanisms happen parallel during this transcription-translation process, in Figure 2 we only present the main mechanisms relevant to the description of our model. Here feedback arises when the protein product returns to the nucleus as a transcription factor and slows down the transcription of its own mRNA by binding to the gene’s promoter site [3] [15] .
Previous findings [15] show that there are time delays associated with the feedback mechanism presented in Figure 2. These delays arise naturally as transcriptional delays (time it takes the gene to get copied into mRNA) and translational delays (time it takes the ribosome to translate mRNA into protein). Furthermore, recent studies [14] [15] have shown that, due to the different time scales in which these two processes occur, we can have an accurate dynamic model by considering only the transcriptional delays. As explained in [14] and [15] , the ordinary differential equations (ODEs) associated to this system are given as follows
(1)
(2)
where the time dependent variables are mRNA concentration,
, and its associated protein concentration,
. The constants
and
are the degradation rates of mRNA and protein,
and
the synthesis rates,
is the repression threshold of protein concentration, n is the hill coefficient, and
is the Hill term representing the rate of delayed mRNA synthesis, where the delay is assumed to be positive and constant,
.
3. Linear Analysis of Two Coupled mRNA-Protein Systems
The single mRNA-protein system was studied in detail by the author in [14] . Here we present a network of two mRNA-protein systems coupled via protein repression. The protein of the first system,
, represses production of its own mRNA,
, and the second system’s mRNA,
. The geometric representation
Figure 2. (a) Diagram representing the synthesis and degradation mechanisms of mRNA and protein. Solid and dashed lines represent chemical reactions and indirect regulatory signals, respectively. The synthesis and degradation rates are given by
and
, respectively, and five small circles represent degradation byproducts. (b) Compact notation for the same mRNA-protein system presented in (a). Here the arrow (
) represents production (upregulation) and the perpendicular symbol (
) represents repression (downregulation).
of the network is depicted in Figure 3 and the associated 4-dimensional mathematical model is given by
(3)
(4)
(5)
(6)
where we assume
,
,
, and
for both systems. Setting
and
sacrifices biological flexibility in regards to synthesis and degradation, but it allows us to gain mathematical and computational tractability for the analysis of the general network presented in Section 4.
Figure 3. Network diagram of two mRNA-protein systems coupled through protein repression. Protein from the first system,
, represses its own mRNA production,
, and the second’s system’s mRNA,
. Each system produces it’s own protein and they all share the same degradation rates
with production rates assumed as
.
The first step in the linear analysis is to set
in (3)-(6) and find expressions for the equilibrium point
. We substitute the latter and eliminate
and
to obtain
(7)
(8)
From (7)-(8) we have that
and thus
. This gives the following expression for the equilibrium solution
(9)
which can be numerically approximated using a root finding method.
Next we define
,
,
, and
to be deviations from equilibrium:
,
,
, and
. Substituting these into (3)-(6) results
(10)
(11)
(12)
(13)
Expanding for small
Equations (10)-(11) become
(14)
(15)
where
is given by
(16)
(17)
and thus the associated linearized system coming from Equations (12)-(15) is given by
(18)
(19)
(20)
(21)
For
(no delay) the system (18)-(21) has the following Jacobian at the origin
(22)
with eigenvalues
(23)
and since
and
then all real parts are negative and the origin is stable.
To find the critical delay,
, we assume solutions of the form
and
for
and substitute into (18)-(21)
(24)
(25)
(26)
(27)
Substituting (26)-(27) into (24)-(25), setting
, and
we obtain
(28)
(29)
The algebraic Equations (28)-(29) will have solutions if
such that
(30)
where
(31)
Equation (30) yields
which implies Im(c) = 0, that is
(32)
(33)
Also,
yields
(34)
(35)
where we have used the identities
(36)
(37)
The stability analysis presented above shows that the stable equilibrium point
of the system (10)-(13) loses its stability when
giving rise to a pair of pure imaginary eigenvalues
corresponding to the solutions
(38)
(39)
where
and
are the amplitudes of the
and
oscillations, and where
is a phase angle. Note that without loss of generality we may set the phase
and thus for values of delay T close to
we will have
where the nonlinear system (3)-(6) is expected to exhibit a periodic solution for small
. The nonlinear analysis for the two coupled mRNA-protein system follows the same computational steps as the N coupled mRNA-protein system. Thus to avoid redundancy we will postpone its presentation until Section 5.
4. Linear Analysis of N Coupled mRNA-Protein Systems
We now extend our previous linear results to a system of N coupled mRNA-protein systems with feedback and delay. The present analysis follows the steps as the one presented in Section 3, hence, and to avoid repetition, we will present only the main results. The system is depicted in Figure 4 and modeled by the following equations
(40)
(41)
where
and
. Setting
on (40) and (41) we can show that
is the equilibrium point, where
(42)
Figure 4. Network diagram representing N coupled mRNA-protein systems. Here each mRNA,
, is repressed through the interaction with every protein,
, via a nonlinear Hill function as given by Equatuon (40). Protein production is affected linearly only by its own mRNA,
, and by the associated intrinsic degradation rate
as modeled by Equatuon (41). The network exhibits negative feedback, nonlinear interactions, and constant delay, all of which are important components for oscillations.
can be multiplied by the product
, and thus be made into a polynomial of degree (
). The roots of the resulting polynomial can then be approximated using a numerical root finding technique.
Next we define
and
to be deviations from equilibrium, substitute them into Equations (40)-(41), and expand for small
to obtain
(43)
(44)
where the Taylor coefficients,
, are given as follows
(45)
(46)
(47)
To show the steady state is stable when
(no delay), we consider the 2N-dimensional linearized system coming from Equations (43)-(44)
(48)
(49)
with its associated Jacobian matrix
(50)
where I is the
identity matrix,
, and
(51)
To find the eigenvalues,
, of the block Jacobian matrix, J, we compute
(52)
which gives the following eigenvalues
(53)
where
(multiplicity
) and since
then the square root term in Equation (53) will be imaginary and thus
for all
. This shows that the equilibrium solution is stable for
.
To find the critical delay,
, we assume solutions of the form
and
, substitute them into (48)-(49), and set
and
to obtain the linear system
(54)
where
. System (54) will have nontrivial solutions when
which gives
and thus
(55)
(56)
where we have used the identities (36)-(37) and where the nonlinear system (40)-(41) is expected to exhibit periodic solutions for
when
.
5. Nonlinear Analysis of the N Coupled System
We use the Poincare-Lindstedt’s perturbation method to find closed form approximate expressions for the amplitude of the limit cycle born at the Hopf bifurcation. This will be accomplished by considering the full nonlinear system (43)-(44) and perturbing off of the critical delay,
. We start by combining Equations (43)-(44) into a system of second order DDEs
(57)
where
, and
,
, and
are the Taylor coefficients given by Equations (45)-(47), respectively. Next we introduce a small parameter
via the following scaling
(58)
stretch time by defining a new independent variable
(59)
and expand
in a power series in
as follows
(60)
where we omit the
term since it turns out to be zero. We substitute Equations (58) and (59) into (57) to obtain
(61)
which will be used after perturbing the delay about
and Taylor expanding
and
. We accomplish the latter by scaling the detuning
as follows
(62)
and expanding
in a power series in
(63)
Substituting Equations (60) and (62) into
and Taylor expanding about
we obtain
(64)
Substituting Equations (60), (63) and (64) into (61), and collecting like powers of
we find
(65)
(66)
(67)
where
stands for terms in
and
, omitted here for brevity.
The next step is to solve Equations (65)-(67). We start with Equation (65) by setting the solutions,
, as
(68)
where by Equation (58) we know
. Next we substitute (68) into (66) and obtain the following expression for
(69)
where
is given by the equation:
(70)
where
(71)
and where we omit the expressions for
and
for brevity.
Next we substitute Equations (68) and (69) into (67), and after trigonometric simplifications have been performed, we equate to zero the coefficients of the resonant terms
and
. This yields the amplitude, A, of the limit cycle that was born in the Hopf bifurcation
(72)
where
(73)
(74)
and
(75)
(76)
and where the
’s are given by Equation (71) and Equations (45)-(47). Equation (72) gives the amplitude of the limit cycle born at the Hopf bifurcation in terms of the original system parameters. The numerical values are calculated via a computer routine and the outcome for different N will be presented in Section 6.
6. Numerical Results
In this section we use our closed form expressions to investigate the nonlinear behavior of Equations (40)-(41) close to the Hopf bifurcation. Setting
,
,
, and
, we compute the amplitude of the limit cycle born at the Hopf by using Equation (72). Our results for
, summarized in Figure 5, show that the Hopf occurs at
and the amplitude of the limit cycle grows as
becomes larger. Several points along the amplitude curve were tested using MATLAB’s built-in function dde23.m on the original DDE system (40)-(41). Figure 5 shows a comparison between our perturbation curve and eight points computed with MATLAB, all of which were in agreement. In addition, we also present one dde23.m numerical simulation, point a5, which shows that the amplitude is approximately 4.350 when
.
Figure 5. Amplitude VS Detuning graph.
was calculated using Equatuon (55) and the amplitudes calculated using various detuning (Delta) values as described in the table. Here we present one explicit numerical simulation result for
, labeled as point a5, which was obtained using MATLAB’s built-in function dde23.m. The rest of the amplitude values were confirmed using MATLAB on the full nonlinear system (40)-(41).
Figure 6. Amplitude VS Detuning plots for N = 2, 7, 15, and 30. Other parameter values were set as
and
for
. Table shows numerical values for the different amplitudes obtained for various N.
Figure 6 shows the same analysis described above for
and 30. Different values of
were found for each N and the detuning was set as
for comparison purposes. These results were all confirmed with MATLAB and they show the accuracy of our closed form expressions in predicting the critical delay, stability, and amplitude of the oscillation for large values of N. In addition, they also show the influence and significance of the original model parameters on the system’s dynamics, where (regardless of N) our model will always exhibit a transition between equilibrium and oscillation due to a Hopf bifurcation. We point out, however, that our perturbation results did not match MATLAB’s numerical findings for
. This is attributed to the nature of our perturbation technique, where the results are guaranteed for small
as given by Equation (62). Thus, the limits of this approach rely on
, that is, they are limited to a small neighborhood about
, which is commonly the case with any kind of perturbation technique.
7. Conclusions
In this paper we investigated the periodic solutions of a mRNA-protein nonlinear model with negative feedback and delay. Nonlinear models exhibiting these features are usually good candidates for periodic behavior. In this work we show that, regardless of the size of the network, having these features yields a dynamical system that always transitions from a stable steady state to a limit cycle oscillation via a Hopf bifurcation. The standard approach to study Hopf bifurcations for DDEs is to use the more conventional center manifold approximation [17] [18] . However, for a network of this size, this would be a computationally intractable endeavor and thus our approach presents a new way on how these perturbation techniques are a useful tool to better understand the dynamics of large networks with nonlinear interactions and delays.
Our closed form approximate expressions for the amplitude of the limit cycle provide a useful way to understand the different roles that the system’s parameters play in the dynamics of the periodic motion. We point out, however, that our linear and nonlinear computational analyses are possible due to the symmetric nature of our network: where we have assumed that all the systems share the same degradation rate
, every mRNA is equally repressed by every protein, and biological location within the cell is not taken into account. Regardless of these simplifications, this study shows the importance of the nonlinearities arising from the Hill term in Equation (40) and the delay, both of which give rise to the periodic dynamics of the system. Future research efforts should be directed towards relaxing these symmetric conditions. For example, including a “weighting” function on the Hill term of Equation (40) would take into account the net effect that a protein has on each mRNA. The latter would be more significant, from a biological viewpoint, because not every protein represses equally each mRNA. Another direction would be to assume different degradation rates
’s for each protein and mRNA. Finally, considering different synthesis rates,
and
, would be another interesting direction that could allow us to classify the significance of each of these biological parameters with regards to the system’s periodic behavior. The analysis of the system in these new directions proves to be computationally cumbersome from a perturbation analysis point of view. However, new findings in the fields of DDEs and perturbation theory will prove that these research directions will provide more insight into the rich dynamics exhibited by these networks.