Skip to main content


Periodicity and dosage optimization of an RNAi model in eukaryotes cells

Article metrics



As a highly efficient and specific gene regulation technology, RNAi has broad application fields and good prospects. The effect of RNAi enhances as the dosage of siRNA increases, while an exorbitant siRNA dosage will inhibit the RNAi effect. So it is crucial to formulate a dose-effect model to describe the degradation effects of the target mRNA at different siRNA dosages.


In this work, a simple RNA interference model with hill kinetic function (Giulia Cuccato et al. (2011)) is extended. Firstly, by introducing both the degradation time delay τ1 of mRNA caused by siRNA and the transportation time delay τ2 of mRNA from the nucleus to the cytoplasm during protein translation, one acquires a novel delay differential equations (DDEs) model with physiology lags. Secondly, qualitative analyses are executed to identify regions of stability of the positive equilibrium and to determine the corresponding parameter scales. Next, the approximate period of the limit cycle at Hopf bifurcation points is computed. Furthermore we analyze the parameter sensitivity of the limit cycle. Finally, we propose an optimal strategy to select siRNA dosage which arouses significant silencing efficiency.


Our researches indicate that when the dosage of siRNA is large, oscillating periods are identical for disparate number of siRNA target sites even if it greatly impacts the critical siRNA dosage which is the switch of oscillating behavior. Furthermore, parametric sensitivity analyses of limit cycle disclose that both of degradation lag and maximum degradation rate of mRNA due to RNAi are principal elements on determining periodic oscillation. Our explorations will provide evidence for gene regulation and RNAi.


The mechanism for sequence-specific post-transcriptional gene silencing that is induced by double-stranded RNA (dsRNA), leading to the regression of the target messenger RNA (mRNA) [1]. This common phenomenon in many eukaryotes, including insects, is named RNA interference (RNAi) by Fire et al. RNAi in animals [2] and in plants [3], is an evolutionarily conservative defense against transgenic or exotic virus infringement mechanism [4]. The process of RNAi can be divided into four stages:

  • Step 1. Double stranded RNA (dsRNA) expressed in or introduced into the cell is cleaved into fragments of 21-23 base pairs (called small interfering RNA, abbreviated as siRNA) by the Dicer enzyme.

  • Step 2. siRNAs are firstly adhered to RNA Induced Silencing Complex (RISC), and whereafter split into the sense strands which are deserted [5], and the antisense strands which are still roped to RISC.

  • Step 3. An available siRNA-RISC complex, includes the siRNA loaded to the Ago protein, is packaged by antisense strand. Then it identifies and unites target mRNAs via the principle of complementary base pairing.

  • Step 4. The antisense strand commands a endonuclease bound to RISC (an Argonaute protein called ‘slicer’) to operate the degradation of the target mRNA. And next the complex is liberated to dispose further mRNA targets.

In recent years, some studies have shown that synthetic siRNAs can effectively trigger RNAi in eukaryotes [6] and the siRNAs seem to avoid off-target effects prompted by longer double-stranded RNAs in mammalian cells [7]. This discovery makes the application of RNAi technology more convenient. The high efficiency and specificity of RNAi make it become a powerful tool for researching gene function. RNAi also provides a novel idea to schedule synthetic biological circuits for synthetic biology [8]. In the treatment of certain genetic diseases, for example viral infections [9], cancer [10] and inherited genetic disorders [11], RNAi has the potential to become a new type of therapeutic tool. In the field of pest management, RNAi also shows its talents [12]. And RNAi technology first approved by the US Environmental Protection Agency as a pesticide in 2017.

Because excessive siRNAs not only affect its efficiency [13], but attract off-target effect [7]. For RNAi application, it is necessary to find a quantitative mathematical model that can describe the relationship between the dosage of siRNA and the RNAi effect. Giulia Cuccato et al. (2011), according to vitro experimental data and squared error measure, capture the most efficient mathematical model of RNA interference in [13].

However, for the model, we consider that there are two important time delays that cannot be ignored during the entire RNAi process. First, degradation of mRNA due to RNAi. Here, we use τ1 to describe this time delay. Next, carriage of mRNA from nucleus to cytoplasm. Thus, we introduce τ2 to represent this time delay. In our work, we start from the model proposed in [13] and then modify it. First, we conduct a qualitative analysis of the model with delay. Our result show that the stability of the only positive equilibrium has changed: it is stable while the original model without time delays, as the time delay increases, it will turn into damped oscillation and lose its stability via a Hopf Bifurcation. Therefore, time delay plays an important role in dynamics of RNAi model and should not be ignored in the modeling of genetic regulation. Next, we introduce the solution to the periodic value of the periodic solution of the system with the limit cycle. And we analyze the parameter sensitivity of the amplitude and period of a periodic solution for a system with a limit cycle. Finally, we give optimal control for quantitative RNAi model by optimization theory.


Qualitative analysis

When the delays are finite, the characteristic equations are functions of delays. As values of the delays change, the stability of the trivial solution may also changes. Such phenomena is often refereed to as stability switches. Next the qualitative analysis of model (20) will be conducted.

Stability and Hopf Bifurcation

In this section, we discuss the local asymptotic stability of the unique positive equilibrium \(Q^{*}(\tilde M, \tilde P)\) and the existence of Hopf bifurcation. Setting β=rSn/(θn+Sn), \(\tilde M \) and \(\tilde P\) are denoted by

$$\tilde M=\frac{k_{m}}{d_{m}+\beta},~~~~~~\tilde P=\frac{k_{p}}{d_{p}}\tilde{M}. $$

For τ1,2>0, characteristic equation of model (20) is given by

$$ \begin{array}{l} (d_{m}+\beta e^{-\lambda\tau_{1}}+\lambda)(d_{p}+\lambda)=0. \end{array} $$

Obviously, λ1=−dp is a negative root of the Eq. (1). Next let the first item of the left side of the Eq. (1) be

$$ \begin{array}{l} f(\lambda)=d_{m}+\beta e^{-\lambda\tau_{1}}+\lambda. \end{array} $$

Lemma 1

For ω[π/(2τ1),π/τ1], let \(\phantom {\dot {i}\!}\beta _{0}=e^{-d_{m}\tau _{1}-1}/\tau _{1},\) β1=−dm/ cos(ωτ1)>0. Then the following results hold.

(a) If β<β0, f(λ) has two real negative roots.

(b) If β=β0, f(λ) has one real negative root.

(c) If β0<β<β1, f(λ) has two complex conjugate roots with Re(λ)<0.

(d) If β=β1, f(λ) has two complex conjugate roots with Re(λ)=0.

(e) If β>β1, f(λ) has two complex conjugate roots with Re(λ)>0.


Function (2) implies f(−)=+, f(+)=+ and f(0)=dm+β>0. Then, letting \(\phantom {\dot {i}\!}f^{\prime }(\lambda)=1-\beta \tau _{1}e^{-\lambda \tau _{1}}=0\) yields

$$\lambda^{*}=\frac{1}{\tau_{1}}\ln(\beta\tau_{1}). $$

So, f(λ) maybe has negative root only if βτ1<1. In addition, because f(λ) is minimum of f(λ) for every \(\lambda \in \mathbb {R}\), thus the function (2) has one real negative root λ if f(λ)=0, namely, β=β0. Hence (b) is proved. If β<β0, we obtain f(λ)<0 and the function (2) has two negative real roots, then (a) is proved too.

For β>β0, the Eq. (2) may has two complex roots. For that, we assume that there exists a solution of the characteristic equation of the form λ=iω(ω>0). Putting it into f(λ), it follows

$$\begin{array}{l} d_{m}+\beta\cos(\omega\tau_{1})-\beta i\sin(\omega\tau_{1})+i\omega=0. \end{array} $$

Comparing real and imaginary parts we get,

$$ \begin{array}{ll} \cos(\omega\tau_{1})=-\frac{d_{m}}{\beta},~~~ \sin(\omega\tau_{1})=\frac{\omega}{\beta}.\vspace{0.5em} \end{array} $$

Squaring and adding the first and the second of (3), we get \((d_{m}^{2}+\omega ^{2})/\beta ^{2}=1\), that is, \(\omega ^{2}=\beta ^{2}-d_{m}^{2}\). Hence positive solution \(\omega _{0}=\sqrt {\beta ^{2}-d_{m}^{2}}\) exists if β>dm. And, corresponding to λ=iω0 and the first equation of (3), there exists \(\tau _{1}^{*}>0\) such that,

$$\begin{array}{l} \tau_{1}^{*}=\frac{1}{\omega_{0}}\arccos\left(-\frac{d_{m}}{\beta}\right),\\ \beta_{1}=-\frac{d_{m}}{\cos(\omega_{0}\tau_{1})}, \qquad \omega_{0}\in[\pi/(2\tau_{1}),\pi/\tau_{1}], \end{array} $$

and (c) and (d) are proved. When β>β1, Re(λ)>0, so (e) is proved. □

Theorem 1

For model (20), the following results hold.

(a) If ββ0, then the equilibrium Q is asymptotically stable.

(b) If β0<β<β1, then the equilibrium Q is oscillatory stable.

(c) If β>β1, then the equilibrium Q is unstable. Furthermore if β=β1, Hopf bifurcation occurs.


(a) and (b) are apparently valid by Lemma 1. Now differentiating (1) with respect to τ1 gives

$$ \begin{array}{lll} \left(\frac{d\lambda}{d\tau_{1}}\right)^{-1}&=\frac{2\lambda+\beta d_{p}e^{-\lambda\tau_{1}}(-\tau_{1})+\beta e^{-\lambda\tau_{1}}+\beta\lambda e^{-\lambda\tau_{1}}(-\tau_{1})+d_{m}+d_{p}}{\beta\lambda e^{-\lambda\tau_{1}}(d_{p}+\lambda)} \vspace{0.5em} \\ &=\frac{2\lambda+d_{m}+d_{p}}{\beta\lambda e^{-\lambda\tau_{1}}(d_{p}+\lambda)}-\frac{\tau_{1}}{\lambda}+\frac{1}{\lambda(d_{p}+\lambda)} \vspace{0.5em} \\ &= \frac{2\lambda+d_{m}+d_{p}}{\lambda(-\lambda^{2}-(d_{m}+d_{p})\lambda-d_{m}d_{p})}-\frac{\tau_{1}}{\lambda}+\frac{1}{\lambda(d_{p}+\lambda)}. \end{array} $$

Then at λ=iω0, one gets

$$ \begin{array}{llll} \left(\frac{d\lambda}{d\tau_{1}}\right)^{-1}\bigg|_{\lambda=i\omega_{0}}&=\frac{2i\omega_{0}+d_{m}+d_{p}}{i\omega_{0}(\omega_{0}^{2}-(d_{m}+d_{p})i\omega_{0}-d_{m}d_{p})}\,-\,\frac{\tau_{1}}{i\omega_{0}}\,+\,\frac{1}{d_{p}i\omega_{0}-\omega_{0}^{2}} \\ &=\frac{(d_{m}+d_{p})^{2}\omega_{0}^{2}+2\omega_{0}(\omega_{0}^{3}-d_{m}d_{p}\omega_{0})}{((d_{m}+d_{p})\omega_{0}^{2})^{2}+(\omega_{0}^{3}-d_{m}d_{p}\omega_{0})^{2}}-\frac{1}{d_{p}^{2}+\omega_{0}^{2}} \\ &\quad+\frac{2(d_{m}+d_{p})\omega_{0}^{3}+(d_{m}+d_{p})(\omega_{0}^{3}-d_{m}d_{p}\omega_{0})}{((d_{m}+d_{p})\omega_{0}^{2})^{2}+(\omega_{0}^{3}-d_{m}d_{p}\omega_{0})^{2}}i \\ &\quad+\frac{\tau_{1}}{\omega_{0}}i-\frac{d_{p}}{\omega_{0}(d_{p}^{2}+\omega_{0}^{2})}i. \end{array} $$


$${\begin{aligned} \begin{array}{lll} \left(\frac{dRe\lambda(\tau_{1}^{*})}{d\tau_{1}}\right)^{-1}\bigg|_{\lambda=i\omega_{0}}&=\frac{(d_{m}+d_{p})^{2}\omega_{0}^{2}+2\omega_{0}(\omega_{0}^{3}-d_{m}d_{p}\omega_{0})}{((d_{m}+d_{p})\omega_{0}^{2})^{2}+(\omega_{0}^{3}-d_{m}d_{p}\omega_{0})^{2}}-\frac{1}{d_{p}^{2}+\omega_{0}^{2}} \\ &=\frac{\omega_{0}^{2}((d_{m}+d_{p})^{2}d_{p}^{2}+(\omega_{0}^{2}-d_{m}d_{p})(2d_{p}^{2}+\omega_{0}^{2}+d_{m}d_{p}))}{(((d_{m}+d_{p})\omega_{0}^{2})^{2}+(\omega_{0}^{3}-d_{m}d_{p}\omega_{0})^{2})(d_{p}^{2}+\omega_{0}^{2})} \vspace{0.5em} \\ &=\frac{\omega_{0}^{2}(\omega_{0}^{4}+2\omega_{0}^{2}d_{p}^{2}+d_{p}^{4})}{(((d_{m}+d_{p})\omega_{0}^{2})^{2}+(\omega_{0}^{3}-d_{m}d_{p}\omega_{0})^{2})(d_{p}^{2}+\omega_{0}^{2})}>0\vspace{0.5em}, \end{array} \end{aligned}} $$
$$\begin{array}{l} sign\left(\frac{dRe\lambda(\tau_{1}^{*})}{d\tau_{1}}\right)=sign\left[\left(\frac{dRe\lambda(\tau_{1}^{*})}{d\tau_{1}}\right)^{-1}\right]=1. \end{array} $$

So, Re(λ)>0 providing \(\tau _{1}>\tau _{1}^{*}\). By the Hopf bifurcation theorem, the condition with \(\tau _{1}>\tau _{1}^{*}\) and Re(λ)>0 guarantees that the Hopf bifurcation at β=β1 is supercritical. The result (c) is proved.

The bifurcation diagram of Eq. (20) as a function of the delay τ1 and of the parameter β is shown in Fig. 1a. Using some parameter values suggested in [14], we take km=10, dm=0.05, kp=1, dp=0.01, other parameter values are set by θ=10, n=4, S=30. It is always possible to choose values of r, θ, n and S such that β>β1(τ1), where β1(τ1) is the parameter that determines a supercritical Hopf Bifurcation. In this case, the delay model (20) has asymptotically stable oscillatory solutions (limit cycle solutions) in Fig. 1a, and the time evolution of protein is shown in Fig. 1c. The phase diagrams for system (20) in damped oscillation region and at the limit cycle are shown in Fig. 1b and d, respectively.

Fig. 1

(a) Bifurcation diagram for delay model. The functions β0(τ1) and β1(τ1) are defined in “Qualitative analysis” section. In the region marked ‘stable’, the positive equilibrium Q is stable. In the region marked ‘damped oscillations’, the solutions of the delay model converge steadily to the positive equilibrium as t is large enough. At Hopf Bifurcation value, the solutions of the delay model converge consistently to limit cycle. (b) Phase diagram of delay model with parameters r=0.5062, τ1=2.5 and τ2=1 corresponding to ‘damped oscillations’ region. (c) Sustained oscillations of protein at Hopf Bifurcation value τ1=3.3594. (d) Phase diagram of delay model at Hopf Bifurcation value τ1=3.3594

Period of the bifurcating oscillatory solution

We knew in the above subsection how the delay differential equation model was capable of generating limit cycle periodic solutions. One indication of their existence is if the steady state is unstable by growing oscillations, although this is certainly not conclusive. From the analysis of the previous section, we found that τ2 has no effect on the stability of the system, and the first equation of (20) is independent. So, resorting to the periodicity of the first equation, we intent to analyze the period of the whole system as the delay occurs. Hence pull out the first equation of (20) separately,

$$ \begin{array}{l} \frac{dM(t)}{dt}=k_{m}-d_{m}M(t)-\frac{rS^{n}}{\theta^{n}+S^{n}}M(t-\tau_{1}). \end{array} $$

Linearising (4) at the first component of the steady state, \(\tilde {M}=k_{m}/(d_{m}+\beta)\), that is, writing \(M(t)-\tilde {M}=m(t)\), yields

$$ \begin{array}{l} \frac{dm(t)}{dt}=-d_{m}m(t)-\beta m(t-\tau_{1}). \end{array} $$

By looking for solutions m(t) in the form m(t)=ceλt, we get

$$ \begin{array}{l} \lambda=-d_{m}-\beta e^{-\lambda\tau_{1}}, \end{array} $$

where c is a constant and the eigenvalues λ are solutions of (6), a transcendental equation in which τ1>0. It is not easy to find the analytical solutions of (6). However, all we really want to know from a stability point of view is whether there are any solutions with Re(λ)>0 which from the form of m(t) implies instability since in this case m(t) grows exponentially with time.

Putting λ=μ+iω, in (6), and now take the real and imaginary parts of the transcendental equation in (6), namely,

$$ \begin{array}{l} \mu=-d_{m}-\beta e^{-\mu\tau_{1}}\cos\omega\tau_{1},\qquad\omega=\beta e^{-\mu\tau_{1}}\sin\omega\tau_{1}. \end{array} $$

We knew that the steady state Q is stable if \(0<\tau _{1}<\tau _{1}^{*}\) and the delay Eq. (20) has an stable periodic solution for \(\tau _{1}=\tau _{1}^{*}\). In the latter case we expect the solution to exhibit stable limit cycle behaviour. The critical value \(\tau _{1}=\tau _{1}^{*}\) is the bifurcation value. The effect of delay in models is usually to increase the potential for instability. Here as τ1 is increased beyond the bifurcation value \(\tau _{1}^{*}\), the steady state becomes unstable.

Near the bifurcation value we can get a estimate of the period of the bifurcating oscillatory solution as follows. Consider the dimensionless form and let

$$ \begin{array}{l} \tau_{1}=\tau_{1}^{*}+\varepsilon,\qquad 0<\varepsilon\ll1. \end{array} $$

The solution λ=μ+iω, of (7), with the Re(λ)=0 when \(\tau _{1}=\tau _{1}^{*}\) is μ=0, \(\omega =\omega _{0}=\sqrt {\beta ^{2}-d_{m}^{2}}\). For ε small we expect μ and ω to differ from μ=0 and ω=ω0 also by small quantities so let

$$ \begin{array}{l} \mu=\delta,\quad\omega=\omega_{0}+\sigma,\quad0<\delta\ll1,\qquad\left|\sigma\right|\ll1, \end{array} $$

where δ and σ are to be determined. Substituting these into the second of (7) and expanding for small δ, σ and ε gives

$$\begin{array}{l} \omega_{0}+\sigma=\beta e^{-\delta(\tau_{1}^{*}+\varepsilon)}\sin[(\omega_{0}+\sigma)(\tau_{1}^{*}+\varepsilon)]\Rightarrow\ \\ \qquad\,\,\sigma\approx -d_{m}\sigma\tau_{1}^{*}-d_{m}\varepsilon\omega_{0}-\omega_{0}\delta\tau_{1}^{*}, \end{array} $$

while the first of (7) gives

$$ \begin{array}{l} \delta=-d_{m}-\beta e^{-\delta(\tau_{1}^{*}+\varepsilon)}\cos[(\omega_{0}+\sigma)(\tau_{1}^{*}+\varepsilon)]\Rightarrow\ \\ \delta\approx \omega_{0}\sigma\tau_{1}^{*} +\omega_{0}^{2}\varepsilon-d_{m}\delta\tau_{1}^{*}. \end{array} $$

Thus on solving these simultaneously

$$\begin{array}{l} \sigma\approx \frac{-(1+d_{m}\tau_{1}^{*})d_{m}\varepsilon\omega_{0}-\omega_{0}^{3}\tau_{1}^{*}\varepsilon}{(1+d_{m}\tau_{1}^{*})^{2}+\omega_{0}^{2}{\tau_{1}^{*}}^{2}}. \end{array} $$

Considering that the imaginary part of λ is the root cause of m(t) periodicity, near the bifurcation, the first of (5) with (6) gives

$${\begin{aligned} \begin{array}{ll} M(t)&=\tilde{M}+Re\left\{c\exp[\delta t+i(\omega_{0}+\sigma)t]\right\} \vspace{0.5em} \\ &\approx \tilde{M}\!+Re\left\{c\exp(\delta t)\exp[it(\omega_{0}\,-\,\frac{(1+d_{m}\tau_{1}^{*})d_{m} \varepsilon\omega_{0}+\omega_{0}^{3}\tau_{1}^{*}\varepsilon}{(1+d_{m}\tau_{1}^{*})^{2}+\omega_{0}^{2}{\tau_{1}^{*}}^{2}})]\!\right\}. \end{array} \end{aligned}} $$

This shows that the delay Eq. (20) has an stable periodic solution due to the occurrence of the Hopf bifurcation with period

$$ \begin{array}{l} T= \frac{2\pi}{\omega_{0}-\frac{(1+d_{m}\tau_{1}^{*})d_{m}\varepsilon\omega_{0}+\omega_{0}^{3}\tau_{1}^{*}\varepsilon}{(1+d_{m}\tau_{1}^{*})^{2}+ \omega_{0}^{2}{\tau_{1}^{*}}^{2}}} \approx \frac{2\pi}{\omega_{0}}=\frac{2\pi}{\arccos(-\frac{d_{m}(\theta^{n}+S^{n})}{rS^{n}})}\tau_{1}^{*}. \end{array} $$

T describes the relations of the Hopf bifurcation period with the maximal regression rate r, siRNA dosage S, the half saturation coefficient θ and the number of siRNA target sites n. Here, Fig. 2 reveals that the levels of mRNA and the protein are persistence when the dosage of siRNA are small, otherwise the periodic oscillating happens. Meanwhile, it indicates that when the dosage of siRNA is large, oscillating periods are identical for disparate number of siRNA target sites even though it greatly impacts the critical siRNA dosage Sn which is the switch of oscillating behavior.

Fig. 2

Relationship among period T, siRNA dosage S and the number of siRNA target sites n

Our delayed differential equations are applied to model gene regulatory network due to RNAi. The periodic solutions of delayed differential equations are subjected to parameters. So, it is necessary that a parametric sensitivity analysis for amplitude and period of periodic solutions.

Parametric Sensitivity

In this section, we use a sensitivity analysis method proposed in [15], and focus on sensitivities of amplitude and of the period when our delay model (20) possesses a periodic solution. Then define the sensitivity equation for parameter dm:

$$\left\{ \begin{array}{ll} \frac{dR_{d_{m}}^{M}(t)}{dt}=-d_{m}R_{d_{m}}^{M}(t)-\frac{rS^{n}}{S^{n}+\theta^{n}}R_{d_{m}}^{M}(t-\tau_{1})-M(t), \\ \frac{{dR}_{d_{m}}^{P}(t)}{dt}=-d_{p}R_{d_{m}}^{P}(t)+k_{p}R_{d_{m}}^{M}(t-\tau_{2}) \vspace{0.5em}. \end{array}\right. $$

In similar way, we get the sensitivity equation with respect to τ1:

$${\begin{aligned} \left\{\begin{array}{lll} \frac{{dR}_{\tau_{1}}^{M}(t)}{dt}&= -d_{m}R_{\tau_{1}}^{M}(t)-\frac{rS^{n}}{S^{n}+\theta^{n}}[R_{\tau_{1}}^{M}(t-\tau_{1})-(k_{m}-d_{m}M(t-\tau_{1}) \vspace{1em}\\ &\quad-\frac{rS^{n}}{S^{n}+\theta^{n}}M(t-2\tau_{1}))], \vspace{0.5em} \\ \frac{{dR}_{\tau_{1}}^{P}(t)}{dt}&= -d_{p}R_{\tau_{1}}^{P}(t)+k_{p}R_{\tau_{1}}^{M}(t-\tau_{2}).\vspace{0.5em} \end{array}\right. \end{aligned}} $$

Analogously, the other sensitivity equations on the rest parameters can be captured, I won’t list them here. Solving the there equations, and according to the circumscription of sensitivities of the limit cycle in [15], we obtain the relative sensitivities of the amplitude and of the period shown in Fig. 3.

Fig. 3

Relative sensitivities of the amplitude (a) and of period (b) in the dosage of protein. Nominal parameter values: km=10, dm=0.05, r=0.5062, S=30, θ=10, n=4, kp=1, dp=0.01, τ1=3.3594, τ2=1

We observe that τ1, RNAi process delay, has a effective impact on both amplitude and period, while τ2, the mRNA translational delay, has inappreciable influence. Because, the occurrence of the limit cycle is only related to the value of τ1, and τ2 does not affect the stability of the equilibrium point of model (20). Moreover, parameter r, the maximal degradation rate of the mRNA due to RNAi, has a important affection on period too. This is because the value of r determines the satisfaction of β=β1. When β=β1, the system (20) will have a limit cycle, where β is the degradation rate of mRNA due to RNAi. In other words, in eukaryotic cells, if the rate of degradation of mRNA due to RNAi is greater than the rate of degradation of the mRNA itself, τ1 and r will be important parameters in the quantitative delay system (20).

Optimizing the dosage of siRNA in RNAi

During the RNAi, excessive siRNAs not only affect the efficiency, but also attracts off-target effect. So the rational dosage of siRNA is crucial for both enhancing RNAi efficiency and reducing cost.

Optimal control for model without delay

Define a cost function as

$$ \begin{array}{l} J=P_{s}ST+\int_{0}^{T} P(t) \,dt, \end{array} $$

where Ps is the cost of per unit siRNA, T is the terminal time. The first part of (10) represents the cost of siRNA consumed in [0,T], and the second part shows the accumulation of protein (denoted by PA) in [0,T]. The aim of this work is to minimize the cost of siRNA and the accumulation of protein.

Problem(Q1). For model (19), choose S[0,200] (according to the experiments in [13]) to minimize the cost function (10).

Since the constraint of the cost function (10) is only the state equation, it must be observed during [0,T]. Then the Lagrange multiplier vector can be used to introduce the equality constraint into the integrand part of the definite integral, thus transforming the constrained optimization problem into an unconstrained optimization problem. Then there is

$$ {\begin{aligned} \begin{array}{ll} \tilde{J}&= P_{s}ST+\int_{0}^{T}P(t)+\lambda_{1}\left(k_{m}-d_{m}M(t)-r\frac{S^{n}}{\theta^{n}+S^{n}}M(t)-\dot{M}(t)\right) \\ &\quad+\lambda_{2}\left(k_{p}M(t)-d_{p}P(t)-\dot{P}(t)\right)dt \vspace{0.5em}. \end{array} \end{aligned}} $$

Introduce a Hamiltonian function H as follows:

$$\begin{array}{l} H=P(t)+\lambda_{1}\left(k_{m}-d_{m}M(t)-r\frac{S^{n}}{\theta^{n}+S^{n}}M(t)\right)\\ \quad\quad+\lambda_{2}\left(k_{p}M(t)-d_{p}P(t)\right). \end{array} $$

Then corresponding costate equations is determined by

$$ {\left\{ \begin{array}{lll} \dot{\lambda}_{1}(t)=-\frac{\partial H}{\partial M(t)}=\lambda_{1}(t)(d_{m}+r\frac{S^{n}}{\theta^{n}+S^{n}})-\lambda_{2}k_{p}, \vspace{0.5em} \\ \dot{\lambda}_{2}(t)=-\frac{\partial H}{\partial P(t)}=\lambda_{2}(t)d_{p}-1, \vspace{0.5em} \\ \lambda_{1}(T)=0,\quad \lambda_{2}(T)=0, \vspace{0.5em} \end{array}\right.} $$

and the corresponding gradients of the cost function (11) with respect to S is

$$ \begin{array}{l} \frac{\partial\tilde{J}(S)}{\partial S}=P_{s}T-\int_{0}^{T} \lambda_{1}M(t)r\theta^{n}\frac{nS^{n-1}}{(\theta^{n}+S^{n})^{2}}dt. \end{array} $$

Optimal control for model with delay

Problem (Q2) For model (20), opt for S[0,200] (according to the experiments in [13]) to minimize the cost function (10).

Rewriting the cost function (10) in the same way, yields

$$ \begin{aligned} \begin{array}{ll} \tilde{J}=& P_{s}ST+\int_{0}^{T}P(t)+\lambda_{1}\\&\left(k_{m}-d_{m}M(t)-r\frac{S^{n}}{\theta^{n}+S^{n}}M(t-\tau_{1})-\dot{M}(t)\right) \vspace{0.5em} \\ &+\lambda_{2}\left(k_{p}M(t-\tau_{2})-d_{p}P(t)-\dot{P}(t)\right)dt \vspace{0.5em}. \end{array}\end{aligned} $$

Define a Hamiltonian function H by

$$\begin{array}{ll} H&=P(t)+\lambda_{1}\left(k_{m}-d_{m}M(t)-r\frac{S^{n}}{\theta^{n}+S^{n}}M(t-\tau_{1})\right)\vspace{0.5em}\\ &\quad+\lambda_{2}\left(k_{p}M(t-\tau_{2})-d_{p}P(t)\right). \end{array} $$

The corresponding costate equations is dominated by

$$ \left\{ \begin{array}{ll} \dot{\lambda_{1}}(t)= d_{m}\lambda_{1}(t)+r\frac{S^{n}}{\theta^{n}+S^{n}}\lambda_{1}(t+\tau_{1})-k_{p}\lambda_{2}(t+\tau_{2}), \vspace{0.5em} \\ \dot{\lambda_{2}}(t)=d_{p}\lambda_{2}(t)-1, \vspace{0.5em} \end{array}\right. $$

with jump conditions

$$ \begin{array}{l} \lambda_{1}(T^{-})=\lambda_{1}(T^{+}),~~~\lambda_{2}(T^{-})=\lambda_{2}(T^{+}), \end{array} $$

and boundary conditions

$$ \begin{array}{l} \lambda_{1}(t)=0,~~~\lambda_{2}(t)=0,~~~t\geq T. \end{array} $$

The corresponding gradients of the cost function (14) with respect to S is governed by

$$ \begin{array}{l} \frac{\partial\tilde{J}(S)}{\partial S}=P_{s}T-\int_{0}^{T} \lambda_{1}M(t-\tau_{1})r\theta^{n}\frac{nS^{n-1}}{(\theta^{n}+S^{n})^{2}} \,dt. \end{array} $$

According to [14], we take km=10, dm=0.05, kp=1, dp=0.01, and set other parameter values r=0.02, θ=10, n=4, Ps=1, T=60, M(0)=160, P(0)=10000. Then, we solve two optimal problems with these parameter values by using Matlab programs.

Simulation 1. Comparison of the optimal value and not optimal value about model without delay

The solution obtained by the optimizer is S=43.01. Substituting S into β, one gets β=0.0199≈r=0.02. This shows that the degradation rate of mRNA due to RNAi has reached the maximum. When S>43.01, correspondingly, β is almost unanimously close to r. It implies that the amounts of protein accumulated are identical and the degradation of mRNA is subject to saturation effects when siRNA dosage is larger. Meanwhile, by employing the optimal result S=43.01 together with a no-optimal value S=10 of siRNA dosage and remaining parameters are as above, we make a comparison about the time evolution of mRNA and protein dosages of model without delay under different siRNA dosages controls (see Fig. 4).

Fig. 4

Comparison chart of the time evolution of mRNA and protein dosages of model without delay under different siRNA dosages controls. The blue dotted line corresponds to the parameter S=10, the red solid line corresponds to the parameter S=43.01

Simulation 2. Comparison of the optimal value and not optimal value about model with delay

For this case, take τ1=2.5, τ2=1 and the other parameters are taken as those in simulation 1. The solution obtained by the optimizer is S=30.37. After substituting, one gets β=0.0198≈r=0.02. Although the two results differ by 12.64, the corresponding β is almost the same. This shows that the degradation rate of mRNA due to RNAi is almost the highest under the optimal conditions. Similarly, when S>30.37, the RNAi-mediated degradation of mRNA is subject to saturation effects, we also make a comparison like simulation 1 at S=30.37 and S=10 (see Fig. 5). In addition, Table 1 gives the value of the optimal siRNA dosage, protein accumulation (PA) and the cost function value J. Obviously, with the participation of time delays, the accumulation of protein is much lower than when there is no time delays, although both of S are taken at the best value.

Fig. 5

Comparison chart of the time evolution of mRNA and protein dosages of model with delay under different siRNA dosages controls. The blue dotted line corresponds to the parameter S=10, the red solid line corresponds to the parameter S=30.37

Table 1 Comparison of effects with and without delays


What we interest in is a mathematical model that reflects the relationship between the RNAi effect and the siRNA dose, which is called the dose-effect model. The study had three primary goals. The first was to depict and forecast the evolution rules of mRNA and protein by the dynamic analysis. The second was to study the effect of parameters on periodic oscillation. The third was to explore the optimal dosage for the significant silencing efficiency. Our work provides a theoretical basis for more precise and economical RNAi experiments and applications. Even so, there are some questions worth exploring further. One is that the degradation and amplification process of siRNA should be considered in RNAi model. The second is that the stochastic effects and variable siRNA dosage should be involved in our model. These factors will result in more complicated dynamic behaviors and reveal more mechanisms of RNAi.


In this paper, we reference a simple Hill kinetic model proposed by [13] and consider the potential effect of two time delays. One is degradation of mRNA due to RNAi, other one is carriage of mRNA from nucleus to cytoplasm. For the improved time-delay system, the role of time delays and the dynamic behavior of system are discussed. Qualitative analyses indicate that the introduction of time delays changes the dynamic behaviors of the system. In detail, as delays increase, the unique positive equilibrium firstly is oscillatory stable and then loses its stability via a Hopf Bifurcation. Furthermore, we give the corresponding parameter scales for these results. Meanwhile, the period of the oscillation solution shows that when the dosage of siRNA is large, oscillating periods are identical for disparate number of siRNA target sites in spite of it greatly impacts the critical siRNA dosage which is the switch of oscillating behavior. And then, parametric sensitivities of the limit cycle is determined. The results indicate that both of degradation lag and maximum degradation rate of mRNA due to RNAi are principal elements on determining periodic oscillation. After that, we propose and solve a simple optimization problem for ODEs model (19) and DDEs model (20) based on the optimization theory. The rational dosage of siRNA is given for both enhancing RNAi efficiency and reducing cost by a Matlab program. The results imply that the optimal dosage of siRNA with delay effects is less than one without time delay.


In this section, we apply and expand the model recommended in [13]. This model well describes the mRNA and protein level in RNA interference process for different dosages of siRNA in mammalian cells in vitro, and great predicts the saturation effect observed experimentally of the RNAi process [13]. The RNAi process caused by siRNA (S) is encapsulated into a whole, and the degradation of the target mRNA (M) due to RNAi is expressed in the form of a functional reaction. In addition, the protein corresponding to the target mRNA is denoted as P. The time evolution of the dosages of mRNA and protein can be described by the ordinary differential equations (ODEs) as follows:

$$ \left\{ \begin{array}{ll} \frac{dM(t)}{dt}=k_{m}-d_{m}M(t)-\frac{rS^{n}}{\theta^{n}+S^{n}}M(t),\\ \frac{dP(t)}{dt}=k_{p}M(t)-d_{p}P(t), \end{array}\right. $$

where M is transcribed at a rate km from the promoter; dm and dp are the degradation rates of M and P, respectively. P is translated at a rate kp form M. The extra degradation rate of M as a result of RNAi is the third segment of the first equation of (19), which is a Hill-kinetic model. Positive integer n is a Hill coefficient, representing the number of siRNA bounded on the target mRNA (or the number of siRNA target sites). r and θ tie to the potency of RNAi induced by siRNA [16]: r denotes the maximal regression rate of M because of RNAi, θ is the dosage of S required to reach half of the maximal degeneration rate r.

Time delay plays an important role in many biological dynamical systems. There are two important biological delays that must be considered when modeling RNAi. One is the RNAi process caused by siRNA, using τ1 to describe it. The other one is the transportation process of mRNA from nucleus to cytoplasm, introducing τ2 to represent it. Then, the time evolution of the dosages of mRNA and protein can be described by the following delay differential equations (DDEs):

$$ \left\{ \begin{array}{ll} \frac{dM(t)}{dt}=k_{m}-d_{m}M(t)-\frac{rS^{n}}{\theta^{n}+S^{n}}M(t-\tau_{1}),\vspace{0.5em} \\ \frac{dP(t)}{dt}=k_{p}M(t-\tau_{2})-d_{p}P(t),\vspace{0.5em}\\ \end{array} \right. $$

with the initial condition: M(t)=M(0) and P(t)=P(0) for −max{τ1,τ2}≤t≤0. It is assumed that all the parameters of model (20) are positive.

In real RNAi experiments and applications, the biological time delays are ubiquitous, such as inhibiting the expression of chitinase of migratory locust, gene knockout in animal and inhibiting cancer proliferation. Therefore, our improved time delay model is more convincing in describing the relationship between siRNA measurement and RNAi efficiency in eukaryotic cells.

Availability of data and materials

Data sharing is not applicable to this article as no datasets were generated or analysed during the current study.



Delay differential equations


Double-stranded RNA


Messenger RNA


Ordinary differential equations


RNA induced silencing complex


RNA interference


Small interfering RNA


  1. 1

    Fire A, Xu S, Montgomery MK, Kostas SA, Driver SE, Mello CC. Potent and specific genetic interference by double-stranded rna in caenorhabditis elegans. Nature. 1998; 391(6669):806–811.

  2. 2

    Hannon GJ. Rna interference. Nature. 2002; 418(6894):244–251.

  3. 3

    Baulcombe DC. Rna silencing in plants. Nature. 2004; 431(7006):356–363.

  4. 4

    Chen YMZ. Rna interference. Chinese Journal of Biological Engineering. 2003; 23(3):39–43. in Chinese.

  5. 5

    Filipowicz W. Rnai: the nuts and bolts of the risc machine. Cell. 2005; 122(1):17–20.

  6. 6

    Brummelkamp TR, Bernards R, Agami R. A system for stable expression of short interfering rnas in mammalian cells. Science. 2002; 296(5567):550–553.

  7. 7

    Caplen NJ, Parrish S, Imani F, Fire A, Morgan RA. Specific inhibition of gene expression by small double-stranded rnas in invertebrate and vertebrate systems. Proceedings of the National Academy of Sciences of the United States of America. 2001; 98(17):9742–9747.

  8. 8

    Deans TL, Cantor CR, Collins JJ. A tunable genetic switch based on rnai and repressor proteins for regulating gene expression in mammalian cells. Cell. 2007; 130(2):363–372.

  9. 9

    Barik S, Bitko V. Prospects of rna interference therapy in respiratory viral diseases: update 2006. Expert Opinion on Biological Therapy. 2006; 6(11):1151–1160.

  10. 10

    Takeshita F, Ochiya T. Therapeutic potential of rna interference against cancer. Cancer Science. 2006; 97(8):689–696.

  11. 11

    Aagaard L, Rossi JJ. Rnai therapeutics: Principles, prospects and challenges. Advanced Drug Delivery Reviews. 2007; 59(2):75–86.

  12. 12

    Price DRG, Gatehouse JA. Rnai-mediated crop protection against insects. Trends in Biotechnology. 2008; 26(7):393–400.

  13. 13

    Cuccato G, Polynikis A, Siciliano V, Graziano M, Bernardo MD, Bernardo DD. Modeling rna interference in mammalian cells. BMC Systems Biology. 2011; 5(1):19–19.

  14. 14

    Zhou P, Cai S, Liu Z, Wang R. Mechanisms generating bistability and oscillations in microrna-mediated motifs. Physical Review E. 2012; 85(4):041916.

  15. 15

    Ingalls B, Mincheva M, Roussel MR. Parametric sensitivity analysis of oscillatory delay systems with an application to gene regulation. Bulletin of Mathematical Biology. 2017; 79(7):1539–1563.

  16. 16

    Khanin R, Vinciotti V. Computational modeling of post-transcriptional gene regulation by micrornas. Journal of Computational Biology. 2008; 15(3):305–316.

Download references


The authors thank the referees for their careful reading of the original manuscript and many valuable comments and suggestions, which greatly improved the presentation of this paper.


Funding bodies did not play any role in the design of the study and in writing this manuscript.

Author information

YP presented the ideas and designed the frame of this paper; TM and MZ finished the proofs, computes and writing of the first draft; CL polished, revised the last draft. All authors read and approved the final manuscript.

Correspondence to Yongzhen Pei.

Ethics declarations

Ethics approval and consent to participate

Not applicable.

Consent for publication

Not applicable.

Competing interests

The authors declare that they have no competing interests.

Additional information

Publisher’s Note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Rights and permissions

Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (, which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The Creative Commons Public Domain Dedication waiver( applies to the data made available in this article, unless otherwise stated.

Reprints and Permissions

About this article

Verify currency and authenticity via CrossMark


  • RNA interference
  • Delay
  • Oscillation period
  • Sensitivity analyse
  • Optimal control