Periodicity and dosage optimization of an RNAi model in eukaryotes cells

Background 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. Results 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. Conclusions 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.


Background
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: 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 * (M,P) and the existence of Hopf bifurcation. Setting β = rS n /(θ n + S n ),M andP are denoted bỹ For τ 1,2 > 0, characteristic equation of model (20) is given by Obviously, λ 1 = −d p is a negative root of the Eq. (1). Next let the first item of the left side of the Eq. (1) be So, f (λ) maybe has negative root only if βτ 1 < 1. In addition, because f (λ * ) is minimum of f (λ) for every λ ∈ 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 d m + β cos(ωτ 1 ) − βi sin(ωτ 1 ) + iω = 0.
Proof (a) and (b) are apparently valid by Lemma 1. Now differentiating (1) with respect to τ 1 gives Then at λ = iω 0 , one gets So, Re(λ) > 0 providing τ 1 > τ * 1 . By the Hopf bifurcation theorem, the condition with τ 1 > τ * 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 k m = 10, d m = 0.05, k p = 1, d p = 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.

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, Linearising (4) at the first component of the steady state, By looking for solutions m(t) in the form m(t) = ce λt , we get 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, We knew that the steady state Q * is stable if 0 < τ 1 < τ * 1 and the delay Eq. (20) has an stable periodic solution for τ 1 = τ * 1 . In the latter case we expect the solution to exhibit stable limit cycle behaviour. The critical value τ 1 = τ * 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 τ * 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 The solution λ = μ+iω, of (7), with the Re(λ) = 0 when For ε small we expect μ and ω to differ from μ = 0 and ω = ω 0 also by small quantities so let where δ and σ are to be determined. Substituting these into the second of (7) and expanding for small δ, σ and ε gives while the first of (7) gives Considering that the imaginary part of λ is the root cause of m(t) periodicity, near the bifurcation, the first of (5) with (6) gives This shows that the delay Eq. (20) has an stable periodic solution due to the occurrence of the Hopf bifurcation with period 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 S n which is the switch of oscillating behavior.
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 d m : In similar way, we get the sensitivity equation with respect to τ 1 : 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. 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 where P s is the cost of per unit siRNA, T is the terminal time. The first part of (10) represents the cost of siRNA  [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 Introduce a Hamiltonian function H as follows: Then corresponding costate equations is determined by ⎧ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎩λ and the corresponding gradients of the cost function (11) with respect to S is
Rewriting the cost function (10) in the same way, yields Define a Hamiltonian function H by The corresponding costate equations is dominated by ⎧ ⎨ ⎩λ with jump conditions and boundary conditions The corresponding gradients of the cost function (14) with respect to S is governed by According to [14], we take k m =10, d m =0.05, k p =1, d p =0.01, and set other parameter values r=0.02, θ=10, n=4, P s = 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).

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.

Discussion
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.

Conclusions
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.

Methods
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: where M is transcribed at a rate k m from the promoter; d m and d p are the degradation rates of M and P, respectively. P is translated at a rate k p 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): 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.