Adv Syst Sci Appl 2020; 03:73–90 Published online at https://ijassa.ipu.ru. Delay Differential Equation Model of Gene Expression Amit Sharma1*, Neeru Adlakha1 1Applied Mathematics & Humanities Department, S.V. National Institute of Technology, Surat-395007, Gujarat, India. Abstract: In this study, a delay differential equation model of gene expression for both retroviruses and normal cell is proposed to study the dynamics of functional gene products. The model is categorised into two sub-models to understand the characteristics of a cell by incorporating time delays in the processes of gene expression. The first model which is for retroviruses, involves time delay in replication, transcription, reverse transcription and translation processes taking place in the cell, while in the second model which is for normal cell, the time delay in transcription and translation processes are incorporated. A numerical solution is obtained using semi-temporal data set. The impact of time delays on temporal concentration profile of DNA, mRNA and proteins have been analysed which gives better insight into the normal cell as well as retroviruses. Further, sensitivity analysis has been performed for both models to study the behaviour of gene expression in the cell. The results obtained from such models can be useful for biomedical applications. Keywords: DNA, RNA, proteins, delay differential equations 1. INTRODUCTION The cell has a very sophisticated control system which dynamically initiates, sustains and terminates the processes of replication, transcription, reverse transcription and translation in the cell. Another vital characteristic of a control system of the cell is its ability to dynamically delay the processes of replication, transcription, reverse transcription and translation in order to regulate the concentrations of DNA, mRNA and proteins. Apart from this, the delay may be caused in these processes of the cell due to some external influence, noise in the cell signal, any abnormality in the environment or due to any disease. The delay in these processes may either be useful to the cell or may have some harmful effects on the cell and the organism. Thus, these delays in the processes may contribute positively to the control system of the cell and its gene expression or it may have a negative impact on the control system of the cell. Experimental investigations have been performed by a number of research workers to study the expression [2, 3, 6, 8]. The experimental investigations are quite expensive and time-consuming, and therefore the scientists have also explored the theoretical approaches to study the gene expression [5, 15, 18, 20, 21, 22]. Further, mathematical modeling plays a vital role for better understanding of the complex real-world problems in every area, for example, supply chain management, inventory theory etc [23]. In terms of differential equations, delay differential equations are more appropriate for modeling of the complex real-world problem, for example, prey-predator model [12, 17], chemostat models [31], circadian rhythms [24], epidemiology [7], the respiratory system [27], tumor growth [28] and neural networks [4]. Thus, delay differential equation is one of the possible approaches to deal and unravel the complexities of real-world problems, and gene ∗Corresponding author: amitsharmajrf@gmail.com 74 A. SHARMA, N. ADLAKHA expression entities are also one of them. According to experimental data, there are some delays in the processes of replication, transcription, reverse transcription and translation to be completed, for example, the transcription process takes about 20 minutes to complete and a model is reported having both transcriptional and translational time delays [13]. Some delay differential equation models are reported in the literature for the study of gene expression [1, 10, 11]. A model with discrete time delay, by using Lindstedt’s method and the Hopf bifurcation is reported to study the gene expression [26]. Also, a temporal model, which demonstrates the intracellular signalling using delay differential equations is reported in the literature [25]. Also, some research workers have developed a model with distributed time delay and with translational time delay using the fourth order Runge-Kutta method to study gene expression [9, 16]. Further, in spite of experimental and theoretical investigations, the gene expression has still not been well understood. Thus, in order to have a better understanding of gene expression and control system of the cell, it is necessary to develop models of gene expression involving delays in the processes of a cell. Also, from the literature survey, it is evident that very few attempts are reported in the literature for modeling of gene expression using delay differential equation, and no attempts are reported in the literature for the development of Michaelis-Menten’s mechanism-based delay differential equations model to study gene expression. In this paper, two models based on delay differential equation are proposed to study the gene expression. The Michaelis-Menten’s mechanism is incorporated in these models. In the first model, the time delay in all four processes replication, transcription, reverse transcription and translation has been considered. In the second model, transcription and translation processes are considered with different time delay. The impact of time delays on the temporal concentration profile of DNA, mRNA and proteins have been analysed with the help of numerical results. The mathematical model is presented in the next section. 2. MATHEMATICAL MODEL In this section, two different delay differential equation models are proposed for the gene expression in retroviruses and a normal cell, respectively. In order to develop the model, the following assumptions are made: 1. There can be time delays in the processes of replication, transcription, reverse transcription and translation. 2. The variation in the temperature is constant throughout the model, so that, there is no effect of temperature on replication, transcription, reverse transcription and translation processes in the gene. 3. The rates of replication, transcription, reverse transcription and translation processes in the cell lie in the interval [0, 1] . 4. A single cell is considered throughout the model, and as the replication of DNA occurs at the time of cell division in the baby cell from the mother cell, it is assumed that the replication occurs in the cell. 5. A single protein is synthesized in a single mRNA transcript. In the present study, two mathematical models have been proposed using Michaelis- Menten’s mechanism. In the first model, the processes of replication, transcription, reverse transcription and translation are incorporated as these processes take place in retroviruses. The first model called “Model-I” is described by the system of delay differential equations (2.1)− (2.3): dw dt = k0w(t− τ1)− k1w(t− τ2) + k2x(t− τ3), (2.1) dx dt = k1w(t− τ2)− k2x(t− τ3)− k3x(t− τ4), (2.2) Copyright © 2020 ASSA. Adv Syst Sci Appl (2020) DELAY DIFFERENTIAL EQUATION MODEL OF GENE EXPRESSION 75 dp dt = k3x(t− τ4), t ≥ 0, τi > 0, i = 1, 2, 3, 4. (2.3) Here k0, k1, k2 and k3 are constants. In the second model, it is assumed that the replication and reverse transcription processes are absent, and transcription and translation processes are taking place in the cell. In general, the replication process stops taking place in a mature cell and the reverse transcription is absent in a normal cell. In view of the above, the second model is proposed incorporating transcription and translation processes only. The second model called “Model-II” is described by the system of delay differential equations (2.4)− (2.6): dw dt = −k1w(t− τ2), (2.4) dx dt = k1w(t− τ2)− k3x(t− τ4), (2.5) dp dt = k3x(t− τ4), t ≥ 0, τi > 0, i = 2, 4. (2.6) We consider systems of differential equations (2.1)− (2.3) and (2.4)− (2.6) with the initial condition function φ : [−τ, 0]→ R3, where τ represents time delay. The time delays are taken as positive constants [30]. Further, if all τi = 0, then there is no delay in the system and for without delay system, the initial condition is consider as w(t) = w0, x(t) = x0, p(t) = p0. Here —————————————————————————————- w(t) Concentration of DNA in the cell at time t (in second). x(t) Concentration of mRNA in the cell at time t (in second). p(t) Concentration of protein in the cell at time t (in second). k0 Rate of replication (microgram/second). k1 Rate of transcription (microgram/second). k2 Rate of reverse transcription (microgram/second). k3 Rate of translation (microgram/second). φ Initial history function. τ Time delay. τ1 Delay in replication process. τ2 Delay in transcription process. τ3 Delay in reverse transcription process. τ4 Delay in translation process. —————————————————————————————- For numerical solution, we use the built-in MATLAB programme in the optimization toolbox [19] for both Model-I and Model-II, which is presented in the next section. 3. RESULTS The rates of replication, transcription, reverse transcription and translation processes will depend on the circumstances and capacity of the cell, and therefore k0, k1, k2 and k3 can be assigned different values. The data of the strain TJK16 is used to compute the results [6]. Generally, the values of k0, k1, k2 and k3 lie between [0, 1] [14, 29]. The results obtained for the above system of delay differential equations of both Model-I and Model-II, Table 3.1, along with the comparison study between with time delay and without time delay of both models, are shown in Fig. 3.1 and Fig. 3.2. Here, four cases are discussed regarding the time lag of replication, transcription, reverse transcription and translation processes, which are shown in Table 3.2. In the first case, Copyright © 2020 ASSA. Adv Syst Sci Appl (2020) 76 A. SHARMA, N. ADLAKHA 0 20 40 60 80 100 120 140 160 180 200 −50 0 50 100 150 Time C on ce nt ra tio n Model−I, with delay (initial condition are 25, 0, 0) 0 20 40 60 80 100 120 140 160 180 200 0 10 20 30 40 50 Time C on ce nt ra tio n Model−I, without delay (initial condition are 25, 0, 0) 0 20 40 60 80 100 120 140 160 180 200 −50 0 50 100 150 Time C on ce nt ra tio n Model−I, with delay (initial condition are 25, 25, 25) 0 20 40 60 80 100 120 140 160 180 200 0 20 40 60 80 100 Time C on ce nt ra tio n Model−I without delay (initial condition are 25, 25, 25) DNA mRNA Protein Fig. 3.1. Graphical representation of the solution of the concentrations of DNA, mRNA and protein with time delay in the replication, transcription, reverse transcription and translation processes (both left figures where the initial conditions of DNA, mRNA and protein are (25, 0, 0) and (25, 25, 25) respectively), and without time delay in the replication, transcription, reverse transcription and translation processes are shown (both right figures where the initial conditions of DNA, mRNA and protein are (25, 0, 0) and (25, 25, 25) respectively, of Model-I). The values of parameters are spatio-temporal which exhibit the stability of the system. 0 20 40 60 80 100 120 140 160 180 200 −20 −10 0 10 20 30 40 50 Time C on ce nt ra tio n Model−II, with delay (initial condititons are 25, 0, 0) 0 20 40 60 80 100 120 140 160 180 200 0 5 10 15 20 25 Time C on ce nt ra tio n Model−II, without delay (initial condititons are 25, 0, 0) DNA mRNA Protein 0 20 40 60 80 100 120 140 160 180 200 −20 0 20 40 60 80 100 120 Time C on ce nt ra tio n Model−II, with delay (initial condititons are 25, 25, 25) 0 20 40 60 80 100 120 140 160 180 200 0 10 20 30 40 50 60 70 Time C on ce nt ra tio n Model−II, without delay (initial condititons are 25, 25, 25) Fig. 3.2. Graphical representation of the solution of the concentrations of DNA, mRNA and protein with time delay in the replication, transcription, reverse transcription and translation processes (both left figures where the initial conditions of DNA, mRNA and protein are (25, 0, 0) and (25, 25, 25) respectively), and without time delay in the replication, transcription, reverse transcription and translation processes are shown (both right figures where the initial conditions of DNA, mRNA and protein are (25, 0, 0) and (25, 25, 25) respectively, of Model-II). The parameter values of transcription and translation rates K1 = 0.42, k3 = 2.59 respectively, are taken from the literature [29]. Copyright © 2020 ASSA. Adv Syst Sci Appl (2020) DELAY DIFFERENTIAL EQUATION MODEL OF GENE EXPRESSION 77 Table 3.1. Parameters related to functional gene products DNA, mRNA and protein for strain TJK16 of both models I and II, and the time is t=200 sec. Model-I Model-II Replication lag τ1 (in second) 1.7 - Transcription lag τ2 (in second) 2.3 3.7 Reverse transcription lag τ3 (in second) 5.2 - Translation lag τ4 (in second) 8.6 0.6 Replication rate k0 (microgram/second) 0.1425 - Transcription rate k1 (microgram/second) 0.4785 0.42 Reverse transcription rate k2 (microgram/second) 0.2568 - Translation rate k3 (microgram/second) 0.3691 2.59 Time t (in second) 200 200 Table 3.2. Parameters related to the time lag in the replication, transcription, reverse transcription and translation processes of functional gene products DNA, mRNA and protein for strain TJK16 at time t=200 sec. Model-I τ1 < τ2 < τ3 < τ4 τ1 > τ2 > τ3 > τ4 τ1 = τ2 > τ3 = τ4 τ1 > τ2 < τ3 < τ4 Replication lag τ1 (in second) 1.7 7.0 4.2 1.7 Transcription lag τ2 (in second) 2.3 6.2 4.2 1.2 R. transcri. lag τ3 (in second) 5.2 3.1 1.9 5.2 Translation lag τ4 (in second) 8.6 1.9 1.9 8.6 Time t (in sec- ond) 200 200 200 200 the time lag of replication, transcription, reverse transcription and translation processes are taken in increasing order as τ1 < τ2 < τ3 < τ4 respectively. In the second case, the time lag of replication, transcription, reverse transcription and translation processes are taken in decreasing order as τ1 > τ2 > τ3 > τ4 respectively. In the third case, the time lag of replication process is equal to the transcription process, and the time lag of reverse transcription process is equal to the translation process, which is taken as τ1 = τ2 > τ3 = τ4. Finally, in the fourth case, the time lag of the replication process is greater than that of the transcription process and in turns lesser than both the reverse transcription and translation processes, which are taken as τ1 > τ2 < τ3 < τ4. The above mentioned conditions of time lag for Model-I are taken in the replication, transcription, reverse transcription and translation processes, respectively and the results are shown in Fig. 3.3. Copyright © 2020 ASSA. Adv Syst Sci Appl (2020) 78 A. SHARMA, N. ADLAKHA 0 50 100 150 200 −200 0 200 Time C on ce nt ra tio n tau 1 < tau 2 < tau 3 < tau 4 0 50 100 150 200 −2000 0 2000 Time C on ce nt ra tio n tau 1 > tau 2 > tau 3 > tau 4 0 50 100 150 200 −100 0 100 Time C on ce nt ra tio n tau 1 > tau 2 > tau 3 > tau 4 0 50 100 150 200 −1 0 1 x 10 11 Time C on ce nt ra tio n tau 1 < tau 2 < tau 3 < tau 4 0 50 100 150 200 −100 0 100 Time C on ce nt ra tio n tau 1 = tau 2 > tau 3 = tau 4 0 50 100 150 200 −1 0 1 x 10 9 Time C on ce nt ra tio n tau 1 = tau 2 < tau 3 = tau 4 0 50 100 150 200 −200 0 200 Time C on ce nt ra tio n tau 1 > tau 2 < tau 3 < tau 4 DNA mRNA Protein 0 50 100 150 200 −100 0 100 Time C on ce nt ra tio n tau 1 > tau 2 > tau 3 < tau 4 Fig. 3.3. Change in the concentrations profiles of functional gene products with the change in time lags of Model-I. The above graphical representation shows the change in the concentrations of DNA, mRNA and protein with different time lags taken in the replication, transcription, reverse transcription and translation processes, respectively. The left figures show the stability of the model while the right figures except for the fourth one (right lower figure) indicate that the system is not stable and the system bursts. The right figures are generated after reversing the order of time lags corresponding to each left figures, respectively. 4. SENSITIVITY ANALYSIS The model can be demonstrated precisely using sensitivity analysis. Based on sensitivity, it can be revealed that the particular model is stable, semi-stable or unstable. Regarding this, we consider both models I and II for sensitivity analysis. First, we analyze the sensitivity of both the models based on the rates of replication, transcription, reverse transcription and translation processes. For Model-I, the results are given in Table 4.3 and the graphical solutions are shown in Fig. 4.4, Fig. 4.5, Fig. 4.6 and Fig. 4.7, while for Model-II, the results are shown in Table 4.4, and the graphical solutions are shown in Fig. 4.8 and Fig. 4.9. Secondly, we analyze the sensitivity of both models based on lags or delays in replication, transcription, reverse transcription and translation processes. For Model-I, the results are given in Table 4.5, and shown in Fig. 4.10, 4.11, 4.12 and 4.13, while for Model-II, the results are given in Table 4.6, and shown in Fig. 4.14 and 4.15. 5. DISCUSSION In the left parts of Fig. 3.1 we observe the variation in the concentration of DNA, mRNA and protein with delay in the processes and without delay in the processes for the initial values of DNA, mRNA and protein as (25, 0, 0) for a newly born cell and (25, 25, 25) for mature cell, respectively. The oscillations are observed in DNA, mRNA and protein concentrations due to the delays in the processes in replication, transcription, reverse transcription and translation (both left parts of Fig. 3.1), and in case of without delays, the smooth curves are seen in DNA, mRNA and proteins concentrations when the delay is absent in replication, transcription, reverse transcription and translation processes (both right parts of Fig. 3.1). Due to change in the initial values of DNA, mRNA and proteins concentrations from (25, 0, 0) to (25, 25, 25), the small oscillation is observed (both left parts of Fig. 3.1). In case of processes Copyright © 2020 ASSA. Adv Syst Sci Appl (2020) DELAY DIFFERENTIAL EQUATION MODEL OF GENE EXPRESSION 79 Table 4.3. Sensitivity analysis of Model-I with the parameter k0, k1, k2 and k3, respectively, and the initial condition of each DNA, mRNA and protein concentration as 25 micrograms. The value of delays τ1, τ2, τ3 and τ4 is 1.7, 2.3, 5.2 and 8.6, respectively. Model-I Replication rate k0 (microgram/second) 0.1425 +20% +50% -20% -50% Replication rate k0 (microgram/second) 0.1710 0.2138 0.1140 0.0713 DNA concentration (in micro- gram) -52.3469 2.4164e+03 -0.3831 -5.3790 mRNA concentration (in micro- gram) 35.1480 -2.6285e+03 0.2160 -7.9647 Protein concentration (in micro- gram) 144.1702 34.4349 87.6179 96.5009 Transcription rate k1 (microgram/second) 0.4785 Transcription rate k1 (microgram/second) 0.5742 0.7178 0.3828 0.2393 DNA concentration (in micro- gram) 126.7853 5.6552e+05 1.8871e+03 -1.7556e+06 mRNA concentration (in micro- gram) -93.0348 -4.6064e+05 -1.1285e+03 9.6616e+06 Protein concentration (in micro- gram) 8.8427 -1.9215e+05 -1.4282e+03 -1.2712e+07 Reverse transcription rate k2 (microgram/second) 0.2568 Reverse transcription rate k2 (microgram/second) 0.3082 0.3852 0.2054 0.1284 DNA concentration (in micro- gram) 1.4106e+03 -1.9018e+03 0.2606 488.0431 mRNA concentration (in micro- gram) -1.3330e+03 1.6340e+06 -0.0292 -8.5203 Protein concentration (in micro- gram) -32.0213 7.4441e+05 90.3924 -416.6525 Translation rate k3 (microgram/second) 0.3691 Translation rate k3 (microgram/second) 0.4430 0.5537 0.2953 0.1846 DNA concentration (in micro- gram) -0.4949 -6.2232e+03 30.9181 -750.3591 mRNA concentration (in micro- gram) 0.4716 -5.5662e+03 -31.0534 742.9808 Protein concentration (in micro- gram) 88.9815 1.5255e+03 103.0534 295.4582 without delay, when the initial values of DNA, mRNA and proteins concentrations changes from (25, 0, 0) to (25, 25, 25), the oscillation in the concentration of mRNA has vanished (both right parts of Fig. 3.1). For the proposed Model-II, there are only transcription and translation processes with delays. In Fig. 3.2, the same behaviour of the DNA, mRNA and Copyright © 2020 ASSA. Adv Syst Sci Appl (2020) 80 A. SHARMA, N. ADLAKHA 0 50 100 150 200 −100 −50 0 50 100 150 Time C on ce nt ra tio n +20% 0 50 100 150 200 −50 0 50 100 150 Time C on ce nt ra tio n −20% 0 50 100 150 200 −3000 −2000 −1000 0 1000 2000 3000 Time C on ce nt ra tio n +50% 0 50 100 150 200 −50 0 50 100 150 Time C on ce nt ra tio n −50% DNA mRNA Protein Fig. 4.4. Sensitivity analysis of the parameter replication rate of Model-I with an increase of 20% and 50%, and decrease of 20% and 50% in the base value. 0 50 100 150 200 −200 −100 0 100 200 Time C on ce nt ra tio n +20% 0 50 100 150 200 −2000 −1000 0 1000 2000 3000 Time C on ce nt ra tio n −20% 0 50 100 150 200 −6 −4 −2 0 2 4 6 x 10 5 Time C on ce nt ra tio n +50% 0 50 100 150 200 −2 −1 0 1 2 x 10 7 Time C on ce nt ra tio n −50% DNA mRNA Protein Fig. 4.5. Sensitivity analysis of the parameter transcription rate of Model-I with an increase of 20% and 50%, and decrease of 20% and 50% in the base value. protein concentration profiles are observed as shown in Fig. 3.1 where oscillations occur due to delays in transcription and translation processes and smooth curves are observed due to no delay in transcription and translation processes. The values of concentrations profiles of DNA, mRNA and protein are different in both Fig. 3.1 and Fig. 3.2. Overall, both figures depicted that Model-I and Model-II are stable. Further, four cases related to delays in the replication, transcription, reverse transcription and translation processes of Model-I are shown in Fig. 3.3. The concentration profiles of DNA, mRNA and protein vary with the change in delays of replication, transcription, reverse Copyright © 2020 ASSA. Adv Syst Sci Appl (2020) DELAY DIFFERENTIAL EQUATION MODEL OF GENE EXPRESSION 81 0 50 100 150 200 −2000 −1000 0 1000 2000 Time C on ce nt ra tio n +20% 0 50 100 150 200 −50 0 50 100 150 Time C on ce nt ra tio n −20% 0 50 100 150 200 −2 −1 0 1 2 x 10 6 Time C on ce nt ra tio n +50% 0 50 100 150 200 −1500 −1000 −500 0 500 1000 1500 Time C on ce nt ra tio n −50% DNA mRNA Protein Fig. 4.6. Sensitivity analysis of the parameter reverse transcription rate of Model-I with an increase of 20% and 50%, and decrease of 20% and 50% in the base value. 0 50 100 150 200 −50 0 50 100 150 Time C on ce nt ra tio n +20% DNA mRNA Protein 0 50 100 150 200 −50 0 50 100 150 Time C on ce nt ra tio n −20% 0 50 100 150 200 −2 −1 0 1 2 x 10 4 Time C on ce nt ra tio n +50% 0 50 100 150 200 −1000 −500 0 500 1000 Time C on ce nt ra tio n −50% Fig. 4.7. Sensitivity analysis of the parameter translation rate of Model-I with an increase of 20% and 50%, and decrease of 20% and 50% in the base value. transcription and translation processes, respectively. The left four figures along with right lower figure indicate the stability of the model while the figures on right hand side except for the fourth one(right lower figure) indicate that the model is unstable and the system bursts. The sensitivity analysis of Model-I, Fig. 4.4 shows the concentration profiles of DNA, mRNA and protein when the replication rate increases by 20% and 50%, and decreases by 20% and 50% in the base value. The oscillation decreases with the decrease in the replication rate (both right parts of Fig. 4.4). It is observed from the figure that when the replication rate is lower, the model is stable (both right parts of the Fig. 4.4) while the increase in replication Copyright © 2020 ASSA. Adv Syst Sci Appl (2020) 82 A. SHARMA, N. ADLAKHA Table 4.4. Sensitivity analysis of Model-II with the parameter k1 and k3, respectively, and the initial condition of each DNA, mRNA and protein concentration as 25 micrograms. The value of delays τ2 and τ4 is 3.7 and 0.6, respectively. Model-II Transcription rate k1 (microgram/second) 0.42 +20% +50% -20% -50% Transcription rate k1 (microgram/second) 0.50 0.63 0.34 0.21 DNA concentration (in micro- gram) 1.8443e+03 1.3452e+08 -0.0028 5.3398e-11 mRNA concentration (in micro- gram) 2.7789e+03 3.9139e+06 0.8346 0.8983 Protein concentration (in micro- gram) -4.5482e+03 -1.3844e+08 74.1682 74.1017 Translation rate k3 (microgram/second) 2.59 Translation rate k3 (microgram/second) 3.11 3.89 2.07 1.30 DNA concentration (in micro- gram) -16.4334 -16.4334 -16.4106 -16.4070 mRNA concentration (in micro- gram) -1.3011e+19 3.2367e+42 -1.5741 -1.8137 Protein concentration (in micro- gram) 1.3011e+19 -3.2367e+42 92.9847 93.2207 0 50 100 150 200 −1 −0.5 0 0.5 1 x 10 4 Time C on ce nt ra tio n +20% 0 50 100 150 200 −20 0 20 40 60 80 100 Time C on ce nt ra tio n −20% DNA mRNA Protein 0 50 100 150 200 −1 −0.5 0 0.5 1 x 10 8 Time C on ce nt ra tio n +50% 0 50 100 150 200 −20 0 20 40 60 80 100 Time C on ce nt ra tio n −50% Fig. 4.8. Sensitivity analysis of the parameter transcription rate of Model-II with an increase of 20% and 50%, and decrease of 20% and 50% in the base value. rate causes bursts in the system and the model will be unstable (both left parts of the Fig. Copyright © 2020 ASSA. Adv Syst Sci Appl (2020) DELAY DIFFERENTIAL EQUATION MODEL OF GENE EXPRESSION 83 0 50 100 150 200 −1 −0.5 0 0.5 1 x 10 19 Time C on ce nt ra tio n +20% 0 50 100 150 200 −20 0 20 40 60 80 100 Time C on ce nt ra tio n −20% 0 50 100 150 200 −3 −2 −1 0 1 2 3 x 10 42 Time C on ce nt ra tio n +50% 0 50 100 150 200 −20 0 20 40 60 80 100 Time C on ce nt ra tio n −50% DNA mRNA Protein Fig. 4.9. Sensitivity analysis of the parameter translation rate of Model-II with an increase of 20% and 50%, and decrease of 20% and 50% in the base value. 0 50 100 150 200 −50 0 50 100 150 Time C on ce nt ra tio n +20% 0 50 100 150 200 −50 0 50 100 150 Time C on ce nt ra tio n −20% 0 50 100 150 200 −50 0 50 100 150 Time C on ce nt ra tio n +50% 0 50 100 150 200 −1500 −1000 −500 0 500 1000 Time C on ce nt ra tio n −50% DNA mRNA Protein Fig. 4.10. Sensitivity analysis of the parameter replication with time delay for Model-I with an increase of 20% and 50%, and decrease of 20% and 50% in the base value. 4.4). In Fig. 4.5, it can be seen that the increase or decrease in transcription rate leads to system burst and makes the model unstable. Thus in Fig. 4.5, the model is highly sensitive as the variation in transcription rate leads to instability. In Fig. 4.6, it is seen that for the small decrease in reverse transcription rate, the model remains stable (right upper part of the Fig. 4.6) but for a substantial and larger amount of decrease in reverse transcription rate, the model becomes unstable (right lower part of Fig. 4.6). Also, the model becomes unstable on increasing the reverse transcription rate (both left parts of Fig. 4.6). In Fig. 4.7, it is observed that the model is stable for a small increase in translation rate (left upper part of Fig. 4.7) Copyright © 2020 ASSA. Adv Syst Sci Appl (2020) 84 A. SHARMA, N. ADLAKHA Table 4.5. Sensitivity analysis of Model-I with the parameter τ1, τ2, τ3 and τ4, respectively, and the initial condition of each DNA, mRNA and protein concentration is 25 micrograms. The value of rates k0, k1, k2 and k3 are 0.1425, 0.4785, 0.2568 and 0.3691, respectively. Model-I Replication lag τ1 (second) 1.7 +20% +50% -20% -50% Replication lag τ1 (second) 2.04 2.55 1.36 0.85 DNA concentration (in micro- gram) -0.3808 0.0017 35.9046 -1.0507e+03 mRNA concentration (in micro- gram) 0.2621 2.1784e-04 -31.9812 815.3509 Protein concentration (in micro- gram) 96.1922 99.6338 89.4417 348.8821 Transcription lag τ2 (second) 2.3 Transcription lag τ2 (second) 2.76 3.45 1.84 1.15 DNA concentration (in micro- gram) -4.3923e+05 -6.3182e+07 0.00029 7.2881 mRNA concentration (in micro- gram) 4.4777e+05 5.6421e+08 7.0222e-04 -24.9418 Protein concentration (in micro- gram) 8.0491e+03 -1.5083e+09 96.8251 125.2183 Reverse transcription lag τ3 (second) 5.2 Reverse transcription lag τ3 (in second) 6.24 7.8 4.16 2.6 DNA concentration (in micro- gram) -211.9863 6.6169e+05 0.3168 -2.2200 mRNA concentration (in micro- gram) 209.0774 -7.9784e+05 -0.1869 -1.6293 Protein concentration (in micro- gram) 127.0256 -4.6294e+04 89.3957 87.2496 Translation lag τ4 (second) 8.6 Translation lag τ4 (second) 10.32 12.9 6.88 4.3 DNA concentration (in micro- gram) -10.3475 -1.1970e+04 -3.2942e+05 1.5779e+10 mRNA concentration (in micro- gram) -8.1428 5.9625e+03 1.1022e+05 -3.3152e+10 Protein concentration (in micro- gram) 109.7224 8.6237e+03 4.2001e+05 1.7112e+10 but for the larger increase in translation rate, the model becomes unstable (left lower part of Fig. 4.7). Also, with a decrease in translation rate, the increase in the size of oscillations is observed and it leads to instability (both right parts of Fig. 4.7). For Model-II, in Fig. 4.8, it can be seen that the size of oscillations increases with the increase in the transcription rate causing the system burst, and thus, the model will be unstable (both left parts of Fig. 4.8). Also, the size of oscillations decreases with a decrease Copyright © 2020 ASSA. Adv Syst Sci Appl (2020) DELAY DIFFERENTIAL EQUATION MODEL OF GENE EXPRESSION 85 0 50 100 150 200 −5 0 5 x 10 5 Time C on ce nt ra tio n +20% 0 50 100 150 200 −50 0 50 100 150 Time C on ce nt ra tio n −20% DNA mRNA Protein 0 50 100 150 200 −3 −2 −1 0 1 2 3 x 10 9 Time C on ce nt ra tio n +50% 0 50 100 150 200 −50 0 50 100 150 Time C on ce nt ra tio n −50% Fig. 4.11. Sensitivity analysis of the parameter transcription with time delay for Model-I with an increase of 20% and 50%, and decrease of 20% and 50% in the base value. 0 50 100 150 200 −300 −200 −100 0 100 200 300 Time C on ce nt ra tio n +20% 0 50 100 150 200 −50 0 50 100 150 Time C on ce nt ra tio n −20% DNA mRNA Protein 0 50 100 150 200 −50 0 50 100 150 Time C on ce nt ra tio n −50% 0 50 100 150 200 −1 −0.5 0 0.5 1 x 10 6 Time C on ce nt ra tio n +50% Fig. 4.12. Sensitivity analysis of the parameter reverse transcription with time delay for Model-I with an increase of 20% and 50%, and decrease of 20% and 50% in the base value. in transcription rate and the oscillations converge to make the system stable (both right parts of the Fig. 4.8). In Fig. 4.9, it can be seen that with the decrease in translation rate, the oscillations converge thereby making the model stable (both right parts of Fig. 4.9). But the increase in translation rate leads to an increase in the size of oscillations, thereby making the model unstable (both right parts of the Fig. 4.9). In Fig. 4.10, it can be seen that the decrease in the time delay in replication process makes the model unstable (both right parts of Fig. 4.10), while model attains stability with the increase in the time delay in replication process (both left parts of Fig. 4.10). Copyright © 2020 ASSA. Adv Syst Sci Appl (2020) 86 A. SHARMA, N. ADLAKHA 0 50 100 150 200 −50 0 50 100 150 Time C on ce nt ra tio n +20% 0 50 100 150 200 −10 −5 0 5 x 10 5 Time C on ce nt ra tio n −20% 0 50 100 150 200 −1.5 −1 −0.5 0 0.5 1 1.5 x 10 4 Time C on ce nt ra tio n +50% 0 50 100 150 200 −4 −3 −2 −1 0 1 2 x 10 10 Time C on ce nt ra tio n −50% DNA mRNA Protein Fig. 4.13. Sensitivity analysis of the parameter translation with time delay for Model-I with an increase of 20% and 50%, and decrease of 20% and 50% in the base value. Table 4.6. Sensitivity analysis of Model-II with the parameter τ2 and τ4, respectively. Here the initial values of each DNA, mRNA and protein concentration are 25 micrograms. The value of rates k1 and k3 is 0.42 and 2.59, respectively. Model-II Transcription lag τ2 (second) 3.7 +20% +50% -20% -50% Transcription lag τ2 (second) 4.44 5.55 2.96 1.85 DNA concentration (in micro- gram) 5.7725e+03 7.6598e+05 -2.1140e-05 7.1306e-23 mRNA concentration (in micro- gram) -604.6839 3.2312e+04 0.8043 0.8564 Protein concentration (in micro- gram) -5.0929e+03 -7.9822e+05 74.1957 74.1436 Translation lag τ4 (second) 0.6 Translation lag τ4 (second) 0.72 0.9 0.48 0.3 DNA concentration (in micro- gram) -16.4332 -16.4329 -16.4184 -16.4080 mRNA concentration (in micro- gram) 7.9425e+15 -4.8743e+28 -1.2182 -0.9982 Protein concentration (in micro- gram) -7.9425e+15 4.8743e+28 92.6366 92.4061 In Fig. 4.11, it can be seen that the increase in the time delay in transcription process makes the model unstable (both left parts of Fig. 4.11), while a small decrease in the time delay in transcription process, oscillations converge and lead to the stability of the system (right upper part of Fig. 4.11), and further, for a larger decrease in the time delay Copyright © 2020 ASSA. Adv Syst Sci Appl (2020) DELAY DIFFERENTIAL EQUATION MODEL OF GENE EXPRESSION 87 0 50 100 150 200 −4000 −2000 0 2000 4000 6000 Time C on ce nt ra tio n +20% 0 50 100 150 200 −20 0 20 40 60 80 100 Time C on ce nt ra tio n −20% DNA mRNA Protein 0 50 100 150 200 −5 0 5 x 10 5 Time C on ce nt ra tio n +50% 0 50 100 150 200 −20 0 20 40 60 80 100 Time C on ce nt ra tio n −50% Fig. 4.14. Sensitivity analysis of the parameter transcription with time delay for Model-II with an increase of 20% and 50%, and decrease of 20% and 50% in the base value. 0 50 100 150 200 −1 −0.5 0 0.5 1 x 10 16 Time C on ce nt ra tio n +20% 0 50 100 150 200 −20 0 20 40 60 80 100 Time C on ce nt ra tio n −20% DNA mRNA Protein 0 50 100 150 200 −4 −2 0 2 4 x 10 28 Time C on ce nt ra tio n +50% 0 50 100 150 200 −20 0 20 40 60 80 100 Time C on ce nt ra tio n −50% Fig. 4.15. Sensitivity analysis of the parameter translation with time delay for Model-II with an increase of 20% and 50%, and decrease of 20% and 50% in the base value. in transcription process, model is semi-stable (right lower part of Fig. 4.11). In Fig. 4.12, for the decrease in the time delay in reverse transcription process, the oscillations converge and lead to stability of the system (both right parts of Fig. 4.12), while for the increase in the time delay in reverse transcription process, the model becomes unstable (both left parts of Fig. 4.12). In Fig. 4.13, for the small increase in the time delay in translation process, the oscillations converge and the model becomes stable (left upper part of Fig. 4.13), while a larger increase in the time delay in translation process, leads to instability as oscillations do not converge (left lower part of Fig. 4.13). The decrease in the time delay in translation Copyright © 2020 ASSA. Adv Syst Sci Appl (2020) 88 A. SHARMA, N. ADLAKHA process makes the model unstable as oscillations do not converge (both right parts of Fig. 4.13). For Model-II, in Fig. 4.14, it can be seen that the model is stable as oscillations converge when time delays in transcription process decrease (both right parts of Fig. 4.14), while the model is unstable when the time delays in transcription process increase (both left parts of Fig. 4.14). In Fig. 4.15 it can be seen that when the time delay in translation process increases then the model becomes unstable and due to oscillations in an increasing manner the system bursts (both left parts of Fig. 4.15), while the decrease in the time delay in the translation process causes a decrease in size of oscillations and it makes the system stable (both right parts of Fig. 4.15). 6. CONCLUSION The proposed model, categorised into two models for retroviruses and normal cell, is employed to study the effect of time delays in the processes of replication, transcription, reverse transcription and translation of the gene expression. It is concluded from the results that the delay in these processes like replication, transcription, reverse transcription and translation causes more dynamic changes in concentrations of DNA, mRNA and proteins in comparison to the system without delay in these processes. The concentration profiles of DNA, mRNA and proteins in the cell for the case without delay are smooth. But the delay in the processes causes the disturbance in the system and therefore the control system of the cell exerts control on these processes to coordinate with each other to regulate the concentration profiles of DNA, mRNA and proteins and thus leads to oscillations. If the delay in these processes is larger, it causes larger disturbances which go beyond the control of the cell thereby making the system unstable. It is concluded that the system is highly sensitive to the rates of these processes replication, transcription, reverse transcription and translation as well as time delay in these processes. The impact of time delays on temporal concentration profile of DNA, mRNA and proteins have been analysed which gives better insight into the normal cell as well as retroviruses. Thus, these two models give us interesting and useful information about the impact of rates of replication, transcription, reverse transcription and translation in presence or absence of time delay in these processes of gene expression which can be useful for various biomedical applications. ACKNOWLEDGEMENTS The first author is grateful to Council of Scientific & Industrial Research (CSIR), New Delhi, India, award no.- 09/1007(0002)/2009 for giving financial assistance as JRF/SRF. Also, the authors are thankful to Department of Biotechnology (DBT), New Delhi, India for providing Bioinformatics Infrastructure Facility at SVNIT, Surat to carry out this work. REFERENCES [1] Bernard, S., Čajavec, B., Pujo-Menjouet, L., Mackey, M. C., and Herzel, H. (2006). Modelling transcriptional feedback loops: the role of gro/tle1 in hes1 oscillations. Philosophical Transactions of the Royal Society of London A: Mathematical, Physical and Engineering Sciences, 364(1842):1155–1170. [2] Blake, W., Kaern, M., Cantor, C., and Collins, J. (2003). Noise in eukaryotic gene expression. Nature, 422(6932):633–637. Copyright © 2020 ASSA. Adv Syst Sci Appl (2020) DELAY DIFFERENTIAL EQUATION MODEL OF GENE EXPRESSION 89 [3] Bremer, H. and Dennis, P. (1996). Escherichia coli and Salmonella: cellular and molecular biology, chapter Modulation of chemical composition and other parameters of the cell by growth rate, pages 1553–1569. American Society for Microbiology, Washington (DC). [4] Campbell, S. A., Edwards, R., and van den Driessche, P. (2004). Delayed coupling between two neural network loops. SIAM Journal on Applied Mathematics, 65(1):316– 335. [5] Chen, T., He, H. L., Church, G. M., et al. (1999). Modeling gene expression with differential equations. In Pacific symposium on biocomputing, volume 4, pages 29–40. World Scientific. [6] Churchward, G., Bremer, H., and Young, R. (1982). Transcription in bacteria at different dna concentrations. Journal of bacteriology, 150(2):572–581. [7] Cooke, K., Van den Driessche, P., and Zou, X. (1999). Interaction of maturation delay and nonlinear birth in population and epidemic models. Journal of Mathematical biology, 39(4):332–352. [8] Crick1958 (1958). On protein synthesis. Symposia of the Society of Experimental Biology, 12:138–163. [9] Fu, G., Wang, Z., Li, J., and Wu, R. (2011). A mathematical framework for functional mapping of complex phenotypes using delay differential equations. Journal of theoretical biology, 289:206–216. [10] Jensen, M., Sneppen, K., and Tiana, G. (2003). Sustained oscillations and time delays in gene expression of protein hes1. Febs Letters, 541(1):176–177. [11] Lewis, J. (2003). Autoinhibition with transcriptional delay: a simple mechanism for the zebrafish somitogenesis oscillator. Current Biology, 13(16):1398–1408. [12] Mohr, M., Barbarossa, M. V., and Kuttler, C. (2013). Predator-prey interactions, age structures and delay equations. arXiv preprint arXiv:1308.2532. [13] Monk, N. A. M. (2003). Oscillatory expression of hes1, p53, and nf-κb driven by transcriptional time delays. Current Biology, 13(16):1409–1413. [14] Ozbudak, E. M., Thattai, M., Kurtser, I., Grossman, A. D., and van Oudenaarden, A. (2002). Regulation of noise in the expression of a single gene. Nature Genetics, 31(1):69– 73. [15] Paulsson, J. (2005). Models of stochastic gene expression. Physics of life reviews, 2(2):157–175. [16] Rateitschak, K. and Wolkenhauer, O. (2007). Intracellular delay limits cyclic changes in gene expression. Mathematical biosciences, 205(2):163–179. [17] Ruan, S. (2009). On nonlinear dynamics of predator-prey models with discrete delay. Mathematical Modelling of Natural Phenomena, 4(02):140–188. [18] Rué, P. and Garcia-Ojalvo, J. (2013). Modeling gene expression in time and space. Annual review of biophysics, 42:605–627. [19] Shampine, L. F., Thompson, S., and Kierzenka, J. (2002). Solving delay differential equations with dde23. Manuscript, available at http://www. mathworks. com/dde tutorial. Copyright © 2020 ASSA. Adv Syst Sci Appl (2020) 90 A. SHARMA, N. ADLAKHA [20] Sharma, A. and Adlakha, N. (2014). Markov chain model to study the gene expression. Advances in Applied Science Research, 5(2):387–393. [21] Sharma, A. and Adlakha, N. (2015). A computational model to study the concentrations of dna, mrna and proteins in a growing cell. Journal of Medical Imaging and Health Informatics, 5(5):945–950. [22] Sharma, A. and Adlakha, N. (2018). Fuzzy system model for gene expression. Egyptian Journal of Medical Human Genetics, 19(4):301–306. [23] Sharma, A., Goel, R., and Dua, N. K. (2012). Optimal policy for eoq model with two level of trade credits in one replenishment cycle. American Journal of Operations Research, 2(1):51–58. [24] Smolen, P., Baxter, D. A., and Byrne, J. H. (2002). A reduced model clarifies the role of feedback loops and time delays in the drosophila circadian oscillator. Biophysical Journal, 83(5):2349–2359. [25] Sturrock, M., Terry, A. J., Xirodimas, D. P., Thompson, A. M., and Chaplain, M. A. (2011). Spatio-temporal modelling of the hes1 and p53-mdm2 intracellular signalling pathways. Journal of Theoretical Biology, 273(1):15–31. [26] Verdugo, A. and Rand, R. (2008). Hopf bifurcation in a dde model of gene expression. Communications in Nonlinear Science and Numerical Simulation, 13(2):235–242. [27] Vielle, B. and Chauvet, G. (1998). Delay equation analysis of human respiratory stability. Mathematical biosciences, 152(2):105–122. [28] Villasana, M. and Radunskaya, A. (2003). A delay differential equation model for tumor growth. Journal of Mathematical Biology, 47(3):270–294. [29] Xie, P. (2014). An explanation of biphasic characters of mrna translocation in the ribosome. Biosystems, 118:1–7. [30] Zhang, T., Song, Y., and Zang, H. (2012). The stability and hopf bifurcation analysis of a gene expression model. Journal of Mathematical Analysis and Applications, 395(1):103– 113. [31] Zhao, T. (1995). Global periodic-solutions for a differential delay system modeling a microbial population in the chemostat. Journal of mathematical analysis and applications, 193(1):329–352. Copyright © 2020 ASSA. Adv Syst Sci Appl (2020) Introduction Mathematical Model Results Sensitivity Analysis Discussion Conclusion