# Linear Analysis of an Integro-Differential Delay Equation Model

Authors: Anael Verdugo
Journal: International Journal of Differential Equations
 (2018)

---

## Abstract

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 ofN delay differential equations (DDEs) coupled to N 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 N becomes large.


---

## Body

## 1. Introduction

New genetic experiments [1–3] and mathematical approaches [4–6] 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 [2, 3]. 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 [7, 8] 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 [5, 6].An important and popular modeling technique in the applied sciences is based on differential equations in all its various forms: linear [9], nonlinear [6, 10], partial [11], stochastic [12, 13], 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 [4–6] and is given by the following set of delay differential equations (DDEs): (1)dmdt=-μmmt+Hpt-T,dpdt=mt-μppt,where the time dependent variables are the mRNA concentration, mt, and its associated protein concentration, pt, and the constants μm and μp are the decay rates of the mRNA and protein molecules, respectively. The function Hpt-T is generally a Hill equation representing the rate of delayed production of mRNA, where the delay, T, 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. The protein then returns to the nucleus and acts as a negative feedback regulator by repressing production of mRNA (see [4–6] for more biological background).Figure 1
(a) Biological circuit diagram of (1). Protein production, protein decay, and mRNA decay are assumed to be linear processes. Production of mRNA is considered as a process affected by a delayed response of protein repression. Here the arrow (↑) represents activation and the perpendicular symbol (⊥) represents repression. Solid and dashed lines represent direct chemical reactions and indirect regulatory signals, respectively. Five small circles represent degradation byproducts. (b) Spatial distribution of protein production in the cytoplasm. Protein synthesis happens at various locations from the nucleus. The distance from the nucleus is represented here by the variable x, where 0≤x≤1.
(a) (b)In this paper, we analyze the steady state stability of a model motivated by (1) and previously studied by the author in [15]. The model is given by an integro-delay differential equation (IDDE) coupled to a partial differential equation (PDE) and is characterized by an exponential “weighting” function that regulates the netrepression effect on mRNA based on protein synthesis location. The model is given by(2)dmdt=-μm+∫01e-x-x-Hpx-,t-Tdx-,(3)dpdt=m-μp,where m=mx,t, p=p(x,t), and e-x-x- 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 adistance variable 0≤x≤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 mRNA.The current work extends our previous study [15] in two different ways. First, the biological setup explained above sets our model (2)-(3) on firmer modeling grounds from our previous study [15]. Here we assume the variable x is “distance” from the nucleus, as opposed to x being a variable that represents gene sites in the DNA as argued in [15]. 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’sdde23.m, which provides more accurate approximations and numerical simulations for the associated 2N-dimensional system. The latter was not accomplished in [15] and thus presented here for the first time.The rest of the paper is organized as follows. In Section2, 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, Tcr, and associated frequency ω are found. In Section 3, we construct a system of N DDEs coupled to N 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 N 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 dm/dt=dp/dt=0, we see from (3) that at steady state m∗=μp∗, where (m∗(x),p∗(x)) represents the steady state solution. Substituting the latter into (2) gives (4)μ2p∗x=∫01e-x-x-Hp∗x-dx-.Splitting the integration limits (5)μ2p∗x=e-x∫0xex-Hp∗x-dx-+ex∫x1e-x-Hp∗x-dx-,and differentiating twice, we obtain an equivalent second-order two-point boundary value problem (BVP) for the equilibrium solution p∗=p∗(x)(6)d2p∗dx2=p∗-2μ2Hp∗,where the boundary conditions (BCs) are given by (7)p∗0=1μ2∫01e-x-Hp∗x-dx-,p∗1=1eμ2∫01ex-Hp∗x-dx-.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 p∗. For the rest of this work, we let (8)Hpx,t=1-px,t,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 solutionm∗x,p∗x, we set px,t=p∗x+ηx,t and mx,t=m∗x+ξx,t, substitute these into (2)-(3), and linearize the resulting equations in ηx,t and ξx,t to obtain (9)dξdt=-μξ-∫01e-x-x-ηdx-dx-,(10)dηdt=ξ-μη.Setting ξx,t=ϕxeλt and ηx,t=ψxeλt gives (11)-eλTλ+μϕx=∫01e-x-x-ψx-dx-,λ+μψx=ϕx,which yields (12)rψx=∫01e-x-x-ψx-dx-,where the RHS has asymmetric integral kernel, ψx is an eigenfunction, and r is the associated eigenvalue given by (13)r=-eλTλ+μ2.Since (12) is a self-adjoint operator of the form (14)L·=∫01Kx,x-·dx-,then the eigenvalue problem (12) has real eigenvalues r∈R. To compute r, we transform the integral equation (12) to the following equivalent second-order BVP: (15)d2ψdx2+ρψ=0,with solutions (16)ψx=c1sin⁡ρx+c2cos⁡ρx,where c1 and c2 are constants and ρ=2/r-1. The endpoint BCs are obtained from (12) as follows: (17)ψ0=ρ2+12∫01e-x-ψx-dx-,ψ1=ρ2+12e∫01ex-ψx-dx-.Substituting the solution (16) into the BCs (17) gives the system of equations (18)ρsin⁡ρ-cos⁡ρ-e-sin⁡ρ-ρcos⁡ρ+eρeρsin⁡ρ-ecos⁡ρ-1-esin⁡ρ-eρcos⁡ρ+ρc1c2=0,which yields the condition on ρ for nontrivial solutions (19)ρ2-1sin⁡ρ-2ρcos⁡ρ=0.Using a numerical root finding technique on (19), we obtain ρ=1.30654,3.67319,6.58462,… which gives the corresponding values for r=2/1+ρ2=0.73881,0.13800,0.04509,…. Thus to determine λ from r we have two cases:(i)
ForT=0, (13) gives λ=-μ±-r and since r=2/(1+ρ2)>0 then Reλ=-μ<0 for μ>0. The latter shows that the equilibrium solution is stable when there is no delay.(ii)
ForT=Tcr and λ=iω, (13) becomes r=-eiωTcriω+μ2 which gives the two real equations (20)r=2μωsin⁡ωTcr+ω2-μ2cos⁡ωTcr,0=ω2-μ2sin⁡ωTcr-2μωcos⁡ωTcr.Solving (20) for sinωTcr and cosωTcr, and using the identity sin2ωTcr+cos2ωTcr=1, we obtain (21)ω=r-μ2.Dividing the expressions for sinωTcr and cosωTcr and solving for Tcr, we obtain (22)Tcr=1r-μ2arctan⁡2μr-μ2r-2μ2,which gives the value of the delay where the equilibrium solution loses stability. The smallest value for Tcr is obtained by setting r=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ξx,t and ηx,t in (9) and (10) with a set of 2N variables ξit and ηit for i=1,2,…,N. This replaces the original model (9) and (10) with a 2N-dimensional system of N DDEs coupled to N ODEs and replaces the integral in (9) with a sum of N terms as follows: (23)ξ˙i=-μξi-1N∑j=1Ne-i-j/Nηjt-T,η˙i=ξi-μηi,where i=1,2,…,N. By assuming solutions of the form ξi=ϕieλt and ηi=ψieλt and substituting them into (23), we obtain (24)-eλTλ+μϕi=1N∑j=1Ne-i-j/Nψj,λ+μψi=ϕi,which yields the following eigenvalue problem: (25)cψi=∑j=1Ne-i-j/Nψj,where c=Nr and r=-eλTλ+μ2. For nontrivial solutions, system (25) of N equations must satisfy detK-cI=0, where K is the N×N matrix K=e-i-j/N and c is its associated eigenvalue. Since K is a symmetric matrix, then all of its eigenvalues are real. Furthermore, K is positive definite because detMii>0 for i=1,2,…,N, where Mii is the ith minor of K along the main diagonal. Hence K is a symmetric positive definite matrix, which shows that c is a positive real number. The steady state stability results are thus summarized as follows:(i)
ForT=0, we have that λ=-μ±-c/N, where μ,N,c>0. This shows that Re(λ)<0 and so the equilibrium solution with no delay is stable.(ii)
ForT=Tcr, we take the smallest value of c for any given N and use (21) and (22) to obtain values for ω and Tcr where we set r=c/N. 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 (26)d2p∗dx2-γp∗=1-γ,where γ=1+2/μ2>0 gives the following solution: (27)p∗x=c1sinh⁡γx+c2cosh⁡γx+2μ2γ.Substituting (27) into the BCs (7), we obtain (28)c1=1-eγ1-γ-1+γeγγμ2γ+μ2+1e2γ+μ2γ-μ2-1,c2=1+eγ1-γ-1+γeγγμ2γ+μ2+1e2γ+μ2γ-μ2-1.Letting μ=0.2, we obtain (29)p∗x=0.12sinh⁡7.14x-0.12cosh⁡7.14x+0.98,which we have plotted in Figure 2 (solid). To confirm and compare this result, we numerically integrate the system (30)m˙i=-μmi+1N∑j=1Ne-i-j/N1-pjt-T,p˙i=mi-μpi,for N=22, μ=0.2, and T=0 using MATLAB’s built-in functiondde23.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 N becomes large. In addition, Figure 2 also presents a time course simulation for N=22, where we exhibit the short transients to equilibrium for the 22 variables pi for i=1,2,…,22.Figure 2
Numerical comparison between the steady state solutions for (29) and the 44-dimensional system given by (30). Here we also present a time course simulation for N=22, which shows the short transient response to equilibrium. The asterisk marks (∗) are the numerical values extracted from t=40 in the time course simulation for the 22 variables pi for i=1,2,…,22.Now we use (21) and (22) to compute the critical delay and frequency where the steady state p∗x loses its stability. For the IDDE system, setting μ=0.2 in (29) gives the values ω=0.83595 and Tcr=0.56184, which we show as the limiting value for the DDE system when N becomes large. Table 1 shows the results for μ=0.2 for various values of N and Figure 3 presents the numerical simulations for μ=0.2, N=22, and various delay values using MATLAB’sdde23.m. For the case N=22, Figure 3 shows that the equilibrium solution is stable for T=0.55<Tcr and unstable for T=0.57>Tcr. For Tcr=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 p1, which we use as one of the representatives for the other 21 pi’s, since they all exhibit the same time course simulation. Table 1 also shows the approach to the limiting value Tcr=0.56184 which was approximately achieved in a system of 2000 equations.Table 1
Numerical results forμ=0.2 and various N.
N c ω Tcr 2 1.6065 0.87365 0.51518 5 3.7453 0.84206 0.55386 10 7.4137 0.83748 0.55982 15 11.0992 0.83663 0.56094 22 16.2655 0.83627 0.56142 30 22.1729 0.83612 0.56161 50 36.9457 0.83601 0.56175 100 73.8836 0.83596 0.56181 200 147.7634 0.83595 0.56183 1000 738.8111 0.83595 0.56184Figure 3
Numerical simulations forμ=0.2, N=22, and different delay values. Here we plotted p1, which serves as a representative for the other 21 pi’s, since they all exhibit the same time course simulations. We can see that the equilibrium solution is stable for T=0.55<Tcr and unstable for T=0.57>Tcr. For Tcr=0.56142 the system exhibits oscillations with a frequency ω=2π/period = 2π/(99.794-92.281) = 0.83628. These are the simulations associated with the N=22 case in Table 1.


## 5. Conclusions

In this work, we investigated the equilibrium solutions and their associated stability of a biological model with negative feedback and time delay. The model is formed of an IDDE coupled to a PDE having time,t, and distance, x, 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, p∗, depends on the distance, x. 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 ofN DDEs coupled to N ODEs. This resulted in a 2N-dimensional system with delay where numerical evaluations for different N were performed and good agreement was found with the “continuous” IDDE counterpart as N became large. In particular, Table 1 shows that Tcrdiscrete→Tcrcontinuous=0.56184 as N→∞, which was confirmed using MATLAB’s built-in functiondde23.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 confirm 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 [15] via MATLAB’sdde23.m and thus providing more accurate and reliable approximations (and numerical simulations) for the associated 2N-dimensional system used to confirm our results. Table 1 summarizes and corrects our previous results [15] by showcasing the N=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 H(p) nonlinear, multiple delays, or a detailed bifurcation study proving that the system undergoes a Hopf bifurcation when T=Tcr.



---

*Source: 5035402-2018-05-02.xml*