Linear Analysis of an Integro-Differential Delay Equation Model
International Journal of Differential Equations
Linear Analysis of an Integro-Differential Delay Equation Model
0 Department of Mathematics, California State University , Fullerton, CA , USA
1 Center for Computational and Applied Mathematics, California State University , Fullerton, CA , USA
This paper presents a computational study of the stability of the steady state solutions of a biological model with negative feedback and time delay. The motivation behind the construction of our system comes from biological gene networks and the model takes the form of an integro-delay differential equation (IDDE) coupled to a partial differential equation. Linear analysis shows the existence of a critical delay where the stable steady state becomes unstable. Closed form expressions for the critical delay and associated frequency are found and confirmed by approximating the IDDE model with a system of delay differential equations (DDEs) coupled to ordinary differential equations. An example is then given that shows how the critical delay for the DDE system approaches the results for the IDDE model as becomes large.
New genetic experiments [
] and mathematical approaches
] have been developed to help us better understand
how genes interact within a cell. Theoretically, the structure
of these interactions or networks are represented by the
various chemical reactions happening at a certain time. If
the reactions under consideration only involve a few genes,
then their dynamic behavior could be understood intuitively
and, most likely, confirmed with a biochemical experiment
]. However, if the system is formed of dozens of reactions,
then developing an intuitive understanding of the system?s
dynamics would be difficult. Nevertheless, current research in
the computational sciences [
] shows that the study of these
large gene networks is an important step which will help us
unravel some of the mysteries in the field of cell biology [
An important and popular modeling technique in the
applied sciences is based on differential equations in all its
various forms: linear [
], nonlinear [
], partial ,
], and delayed [
3, 6, 14
]. In this study we focus
our attention to a differential equation model with constant
delay, where the delay arises naturally as the time lag
associated with various intracellular processes, like movement
within the cell, synthesis of proteins, and transcription of
DNA, among many others. The model that motivated this
work was studied previously by [
] and is given by the
following set of delay differential equations (DDEs):
where the time dependent variables are the mRNA
concentration, (), and its associated protein concentration, (),
and the constants and are the decay rates of the mRNA
and protein molecules, respectively. The function (( ? ))
is generally a Hill equation representing the rate of
production of mRNA, where the delay, , is assumed to be a
positive constant. The associated biochemical representation
of the system is given in Figure 1(a) and the biological context
is the following: a gene is copied onto mRNA in the nucleus,
which is then translated into a protein in the cytoplasm of
the cell. T he protein then returns to the nucleus and acts
as a negative feedback regulator by repressing production of
mRNA (see [
] for more biological background).
In this paper, we analyze the steady state stability of a
model motivated by (1) and previously studied by the author
]. The model is given by an integro-delay differential
equation (IDDE) coupled to a partial differential equation
Production with delay, H((p ? T)) mRNA
Linear decay, ?
Proteins mRNA DNA
(PDE) and is characterized by an exponential ?weighting?
function that regulates the net repression effect on mRNA
based on protein synthesis location. T he model is given by
where = ( , ), = (, ), and ?|? | is a weighting
function that accounts for a ?stronger? mRNA repression for
proteins being synthesized closer to the nucleus than more
distant ones. The latter is due to the spatial distribution of
protein production within the cytoplasm, which occurs after
mRNA exits the nucleus and diffuses into the cytoplasm
where it is caught, read, and translated into a protein. The
exact location from the nucleus where this process occurs
is arbitrary and here we quantify it with a distance variable
0 ? ? 1 as explained in Figure 1(b). The latter yields
the integral term in (2) which represents the total sum of
the repression effect that newly synthesized proteins have on
The current work extends our previous study [
] in two
different ways. First, the biological setup explained above
sets our model (2)-(3) on firmer modeling grounds from
our previous study [
]. Here we assume the variable is
?distance? from the nucleus, as opposed to being a variable
that represents gene sites in the DNA as argued in [
This is a crucial difference that yields a better understanding
of our computational results. Second, the results from the
analysis of the steady state and its associated stability are
now confirmed via MATLAB?s dde23.m, which provides
more accurate approximations and numerical simulations for
the associated 2-dimensional system. The latter was not
accomplished in [
] and thus presented here for the first
The rest of the paper is organized as follows. In Section 2,
we present the associated linear stability analysis of (2)-(3)
and characterize the steady state solutions. Linear analysis
reveals the existence of a critical delay where the stable steady
state becomes unstable and thus closed form expressions for
the critical delay, cr, and associated frequency are found.
In Section 3, we construct a system of DDEs coupled
to ordinary differential equations (ODEs) and use these
to confirm the results obtained in Section 2. A numerical
example is then given in Section 4, which shows how the
critical delay for the DDE system approaches the results for
the IDDE model as becomes large. In Section 5, we discuss
our findings and conclusions.
2. Linear Stability Analysis
In this section, we consider the steady state behavior of (2)
and (3). Setting / = / = 0, we see from (3) that at
steady state ? = ?, where ( ?(), ?())represents the
steady state solution. Substituting the latter into (2) gives
The BVP (6)?(7) has a unique solution as long as the right
hand side (RHS) of (6) has bounded, positive, and continuous
partial derivatives with respect to ?. For the rest of this work,
which allows mathematical tractability for the stability
analysis presented below. Notice that since (8) satisfies all three
aforementioned BVP conditions, then we are guaranteed the
existence of a unique solution, which can be approximated
using a numerical technique for BVPs, such as finite
differences or a shooting method. See Section 4 for an example.
To study the stability of the steady state solution
( ?(), ?()), we set (, ) = ?() + (, ) and (, ) =
?() + (, ), substitute these into (2)-(3), and linearize the
resulting equations in (, ) and (, ) to obtain
= ? ? ?
= ? .
( ) ,
Setting (, ) = ()
and (, ) = ()
( + )
( + )
and differentiating twice, we obtain an equivalent
secondorder two-point boundary value problem (BVP) for the
equilibrium solution ? = ?()
then the eigenvalue problem (12) has real eigenvalues ? R.
To compute , we transform the integral equation (12) to the
following equivalent second-order BVP:
where the boundary conditions (BCs) are given by
where the RHS has a symmetric integral kernel, () is an
eigenfunction, and is the associated eigenvalue given by
Since (12) is a self-adjoint operator of the form
Dividing the expressions for sin( cr) and cos( cr) and
solving for cr, we obtain
which gives the value of the delay where the equilibrium
solution loses stability. The smallest value for cr is obtained
by setting = 0.73881 to obtain an expression in terms of the
decay rate . In Section 4, we present a numerical example to
confirm these results.
3. Approximating the IDDE-PDE Equations with a DDE-ODE System
To check our previous results, we ?discretize? the variables
(, ) and (, ) in (9) and (10) with a set of 2 variables
()and ()for = 1, 2, . . . , . This replaces the original
model (9) and (10) with a 2-dimensional system of DDEs
coupled to ODEs and replaces the integral in (9) with a sum
of terms as follows:
where = 1, 2, . . . , . By assuming solutions of the form =
and = and substituting them into (23), we obtain
? = ?
( + )
( + )
which yields the following eigenvalue problem:
where = and = ? ( + ) 2. For nontrivial solutions,
system (25) of equations must satisfy det( ? ) = 0,
where is the ? matrix = [ ?|?|/ ] and is its
associated eigenvalue. Since is a symmetric matrix, then all
of its eigenvalues are real. Furthermore, is positive definite
because det( ) > 0for = 1, 2, . . . , , where is the th
minor of along the main diagonal. Hence is a symmetric
positive definite matrix, which shows that is a positive real
number. The steady state stability results are thus summarized
(i) For = 0 , we have that = ? ? ??/, where
, , > 0. T his shows that Re () < 0 and so the
equilibrium solution with no delay is stable.
(ii) For = cr, we take the smallest value of for any
given and use (21) and (22) to obtain values for
and cr where we set = / . A numerical example
of this case is presented in the following section.
4. Numerical Example
In this section, we present a numerical example to compare
and confirm our previous results. From (6) and (8), we obtain
where = 1 + 2/ 2 > 0 gives the following solution:
( )= 1 sinh (?) +
2 cosh (?) +
which we have plotted in Figure 2 (solid). To confirm and
compare this result, we numerically integrate the system
for = 22, = 0.2, and = 0 using MATLAB?s
built-in function dde23.m. Figure 2 shows a summary of the
comparison between (29) and the 44-dimensional system
(30), where we can see that good agreement was found
between both systems as becomes large. In addition,
Figure 2 also presents a time course simulation for = 22 ,
where we exhibit the short transients to equilibrium for the
22 variables for = 1, 2, . . . , 22 .
Now we use (21) and (22) to compute the critical delay
and frequency where the steady state ?() loses its stability.
For the IDDE system, setting = 0.2 in (29) gives the values
= 0.83595 and cr = 0.56184, which we show as the
limiting value for the DDE system when becomes large.
Table 1 shows the results for = 0.2 for various values of
and Figure 3 presents the numerical simulations for = 0.2 ,
= 22 , and various delay values using MATLAB?s dde23.m.
For the case = 22 , Figure 3 shows that the equilibrium
solution is stable for = 0.55 < cr and unstable for
= 0.57 > cr. For cr = 0.56142, the system exhibits
oscillations with a frequency = 2/period = 2/(99.794 ?
92.281)= 0.83628 as predicted. Notice that in Figure 3 we
only plotted 1, which we use as one of the representatives
for the other 21 ?s, since they all exhibit the same time course
simulation. Table 1 also shows the approach to the limiting
value cr = 0.56184 which was approximately achieved in a
system of 2000 equations.
In this work, we investigated the equilibrium solutions and
their associated stability of a biological model with negative
feedback and time delay. T he model is formed of an IDDE
coupled to a PDE having time, , and distance, , as
independent variables. The study considers linear production and
degradation rates of mRNA and protein and an exponential
weighting function that models the net repression of all
proteins due to spatial distribution in the cytoplasm. Our
steady state analysis was accomplished by transforming the
steady state integral equation into a second-order two-point
boundary value problem, and it showed that the equilibrium
solution, ?, depends on the distance, . Stability analysis
then revealed that the nondelayed system is stable and that
there exists a critical value for the delay where the equilibrium
loses its stability.
We confirmed our results by ?discretizing? our original
model and approximating it with a system of DDEs
coupled to ODEs. This resulted in a 2-dimensional
system with delay where numerical evaluations for different
were performed and good agreement was found with
the ?continuous? IDDE counterpart as became large. In
particular, Table 1 shows that cdriscrete ? ccrontinuous = 0.56184
as ? ? , which was confirmed using MATLAB?s
builtin function dde23.m on the full DDE model in a system
of 2000 equations. Unfortunately, there are no numerical
routines available in MATLAB for the simulation of IDDEs,
but our results conf irm that it is possible to dissect and
understand the dynamics of such complicated equations via
a discretization approach, as the one presented in Section 3.
The current work corrects and extends our previous study [
via MATLAB?s dde23.m and thus providing more accurate
and reliable approximations (and numerical simulations) for
the associated 2-dimensional system used to confirm our
results. Table 1 summarizes and corrects our previous results
] by showcasing the = 22 numerical simulations and
their transition from stable to unstable behavior as seen in
Figure 3. It is thus hoped that our approach will be useful
to researchers in the field of computational mathematics and
gene networks trying to model physical or biological systems
characterized by IDDEs and PDEs. Future possible directions
for this work include choosing () nonlinear, multiple
delays, or a detailed bifurcation study proving that the system
undergoes a Hopf bifurcation when = cr.
Conflicts of Interest
The author declares that there are no conflicts of interest
regarding the publication of this paper.
Hindawi Publishing Corporation
Submit your manuscripts at
w w w.hindawi.com
International Journal of
 I. Palmeirim , D. Henrique , D. Ish-Horowicz , and O. Pourquie? , ? Avian hairy gene expression identifies a molecular clock linked to vertebrate segmentation and somitogenesis , ? Cell , vol. 91 , no. 5 , pp. 639 - 648 , 1997 .
 M. B. Elowitz and S. Leibier , ? A synthetic oscillatory network of transcriptional regulators , ? Nature , vol. 403 , no. 6767 , pp. 335 - 338 , 2000 .
 D. Bratsun , D. Volfson , L. S. Tsimring , and J. Hasty , ? Delayinduced stochastic oscillations in gene regulation , ? Proceedings of the National Acadamy of Sciences of the United States of America , vol. 102 , no. 41 , pp. 14593 - 14598 , 2005 .
 J. Lewis , ? Autoinhibition with transcriptional delay: a simple mechanism for the zebrafish somitogenesis oscillator , ? Current Biology , vol. 13 , no. 16 , pp. 1398 - 1408 , 2003 .
 N. A. M. Monk , ? Oscillatory expression of Hes1, p53, and NFB driven by transcriptional time delays , ? Current Biology , vol. 13 , no. 16 , pp. 1409 - 1413 , 2003 .
 A. Verdugo and R. Rand , ? Hopf bifurcation in a DDE model of gene expression , ? Communications in Nonlinear Science and Numerical Simulation , vol. 13 , no. 2 , pp. 235 - 242 , 2008 .
 H. de Jong, ? Modeling and simulation of genetic regulatory systems: a literature review , ? Journal of Computational Biology , vol. 9 , no. 1 , pp. 67 - 103 , 2002 .
 T. Schlitt and A. Brazma , ? Current approaches to gene regulatory network modelling,? BMC Bioinformatics , vol. 8 , no. 6 , article no. S9 , 2007 .
 H. De Jong , J.-L. Gouze?, C. Hernandez , M. Page , T. Sari , and J. Geiselmann , ? Qualitative simulation of genetic regulatory networks using piecewise-linear models,? The Bulletin of Mathematical Biology , vol. 66 , no. 2 , pp. 301 - 340 , 2004 .
 A. Mochizuki , ? Structure of regulatory networks and diversity of gene expression patterns , ? Journal of Theoretical Biology , vol. 250 , no. 2 , pp. 307 - 321 , 2008 .
 S. Turner , J. A. Sherratt , K. J. Painter , and N. J. Savill , ? From a discrete to a continuous model of biological cell movement,? Physical Review E , vol. 69 , 2004 .
 J. Goutsias and S. Kim , ? Stochastic transcriptional regulatory systems with time delays: A mean-field approximation , ? Journal of Computational Biology , vol. 13 , no. 5 , pp. 1049 - 1076 , 2006 .
 A. Ribeiro , R. Zhu , and S. A. Kauffman , ? A general modeling strategy for gene regulatory networks with stochastic dynamics , ? Journal of Computational Biology , vol. 13 , no. 9 , pp. 1630 - 1639 , 2006 .
 R. Edwards , P. Van Den Driessche , and L. Wang , ? Periodicity in piecewise-linear switching networks with delay , ? Journal of Mathematical Biology , vol. 55 , no. 2 , pp. 271 - 298 , 2007 .
 A. Verdugo and R. H. Rand , ? DDE Model of Gene Expression: A Continuum Approach ,? in ASME Proceedings of IMECE , pp. 1112 - 1120 , 2008 .