In this study, we formulate a computational reaction model following a chemical kinetic theory approach to predict the binding rate constant for the siRNA-RISC complex formation reaction. The model allowed us to study the potency difference between 2-nt 3' overhangs against blunt-ended siRNA molecules in an RNA interference (RNAi) system. The rate constant predicted by this model was fed into a stochastic simulation of the RNAi system (using the Gillespie stochastic simulator) to study the overall potency effect. We observed that the stochasticity in the transcription/translation machinery has no observable effects in the RNAi pathway. Sustained gene silencing using siRNAs can be achieved only if there is a way to replenish the dsRNA molecules in the cell. Initial findings show about 1.5 times more blunt-ended molecules will be required to keep the mRNA at the same reduced level compared to the 2-nt overhang siRNAs. However, the mRNA levels jump back to saturation after a longer time when blunt-ended siRNAs are used. We found that the siRNA-RISC complex formation reaction rate was 2 times slower when blunt-ended molecules were used pointing to the fact that the presence of the 2-nt overhangs has a greater effect on the reaction in which the bound RISC complex cleaves the mRNA.
RNA interference (RNAi) refers to a post-transcriptional gene silencing mechanism with potential therapeutic application for the treatment of various diseases including cancer, viral infections, and neurodegenerative disorders . The RNAi pathway involves the introduction of a small interfering double stranded (ds) RNA, siRNA, promoting degradation of a target mRNA . RNAi molecules have sequence complementation to that of the siRNA antisense (or guide) strand ultimately inhibiting the translation of the encoded protein. Currently, RNAi is being extensively used to study the functions of individual genes based on its intrinsic property of regulating the expression of a distinct mRNA species in mammalian cells [3,4] as well as offering additional improvements over alternative technologies including knock-out by homologous recombination and antisense. In Ref. , the authors performed a comparative analysis of the suppressive effects of three knockdown methods, namely, methods based on RNA interference (RNAi), antisense ODNs, and ribozymes, using a luciferase reporter system. Their dose-response experiments revealed that the IC50 value for the siRNA was about 100-fold lower than that of the antisense ODN besides providing useful information about the positional effects in RNAi.
An important consideration towards selecting an effective siRNA-based gene silencing tool is the duration of effect and efficacy of a candidate molecule. Extensive studies determining the intensity of gene silencing mediated by siRNA were reported in Ref.  & by finding the optimal sequence of siRNA. Upon introduction of siRNAs with an optimized sequence into cells, the target mRNA is degraded resulting in a lowering of the corresponding protein. It has been reported that 21-nt siRNAs with 2-nt 3' overhangs are more potent than other types of siRNA molecules in terms of reduction in protein levels . Recent studies have reported a higher potency effect due to 27-mer RNA duplexes which we will consider in the Discussions section.
Upon introduction into a cell, a siRNA molecule will be diluted over time due to its degradation and cellular proliferation resulting in a decrease in its effective concentration. Consequently, the expression level of the target gene will return to a normal level after the gene silencing period, dependent upon the number of siRNA molecules actually entering the cell. To use siRNA for silencing target gene expression, it is important to understand how long the target mRNA or protein is suppressed by the siRNA. Maximal inhibitory efficiency of siRNA, a parameter that has frequently been used to express the potency of each siRNA, should be discussed in context with the duration or persistence of its effect.
Figure 1 illustrates the RNAi system which consists of the following major steps . Long dsRNAs (e.g. >200 base pairs) get processed into 19–25 nucleotide (nt) small interfering RNAs (siRNAs) by an RNase III-like enzyme called Dicer (initiation step) or in the case of mammalian RNAi are introduced to cells directly as "active"(19–25 nt) siRNA molecules. The siRNAs assemble into endoribonuclease containing complexes known as RNA-induced silencing complexes (RISCs), unwinding in the process. Activated RISC then binds to complementary transcript by base pairing interactions between the siRNA antisense strand and the mRNA. The bound mRNA is cleaved and sequence specific degradation of mRNA results in gene silencing. We refer to the aberrant pieces of RNA after cleavage as "garbage RNA" (referred to as gRNA in this paper). It has been reported that a substantial fraction of the siRNAs in Caenorhabditis elegans is not derived directly from the introduced dsRNA . To explain this, two amplification routes have been proposed: primed and unprimed amplification [11-13]. In both cases, RNA-dependent RNA polymerase (RDR) synthesizes dsRNA: in the case of primed amplification siRNA binds to mRNA to initiate dsRNA synthesis, whereas in the case of unprimed amplification the mere presence of aberrant garbage RNA is sufficient to trigger RDR. In short, the generally accepted pathway of RNA silencing consists of the degradation of mRNA via RISC and an amplification pathway through RDR.
Figure 1. RNAi system overview.
In this study, we use computational modeling to investigate the potency effects of the widely used 21-nt siRNAs with 2-nt 3' overhangs as compared to blunt-ended siRNA molecules to assess the effectiveness of the latter as an alternative structural entity. We have developed a simple systems biology simulation to quantitatively assess both the intensity and duration of gene silencing by siRNA. The stochastic simulation framework presented here allows us to predict some quantitative and qualitative aspects of the RNAi system for the two different types of siRNAs.
Materials and methods
The reaction model was built to investigate the potency difference between 2-nt 3' overhangs and blunt-ended siRNA molecules. We identified that the siRNA-RISC complex formation reaction rate was altered due to the different siRNA molecular structures. We assume that once the bound RISC complex is formed, it will cleave the mRNA with the same rate irrespective of the presence of the overhangs on the siRNA. Findings in Ref  suggest that the direction of Dicer processing confers some kind of functional polarity, possibly at the level of RISC loading. However, the mRNA cleavage reaction rate has not been explicitly studied for different types of siRNAs as yet.
Hence, our first contribution is a chemical kinetic theory based analytical model for measuring the rate constant of the siRNA-RISC complex formation reaction. The two different rate constants predicted by this model for the different siRNA structures were then fed into a stochastic simulation (using the widely used Gillespie stochastic simulator ) to study the potency effects. This gave us a relative quantification between the gene silencing period and the concentration of the two types of siRNA molecules required to keep the mRNA count at the same level. We also incorporated the gene transcription/translation reactions to study whether a burst in mRNA production  has a role to play in the RNAi system.
Consider the elementary reaction with three types of molecules siRNA, RISC and the siRNA-RISC complex:
We divide the reaction event into two independent micro-events as follows; 1) Random collisions between the reactants; this allows us to compute the probability of collision (pc) between the reactant molecules. 2) A reaction will occur only when the kinetic energy of the colliding reactant exceeds the activation energy requirement for the reaction. Using these two events, allows us to compute the probability of reaction (pr).
The total probability for reaction after a collision is hence the joint probability of these two events. To model this reaction analytically in the time domain, we first assume that the siRNA molecules enter the cell one at a time. Note that siRNAs are typically delivered via transfection thereby introducing a bolus of molecules into the cell. Thus, we need to consider the effective number of siRNAs in the cell while computing the binding rate. We will show how the computations change while we consider a certain concentration of siRNA molecules for deriving the overall binding rate subsequently in this section. We also assume that the cell contains a fixed number, n2, of RISC molecules. Note that while modeling the reaction, it is not necessary to consider the fact that the siRNAs or RISC complexes can also take part in other reactions with other reactants or can degrade independently. The time domain model is based solely on the current instance of these two reactants in the cell as the time taken to complete this reaction will generally be less in comparison to the time taken to degrade a siRNA or a RISC molecule. The idea here is to discretize this reaction from other competing reactions that can change the concentrations of the siRNA/RISC molecules. Though this approximation might lead to less accurate predictions of the binding rate for this reaction, we can still consider the effects of such competing reactions by the system simulation of the RNAi pathway (as shown later).
We follow the principles of collision theory for hard spheres  to model the chemical kinetics of the reaction. As shown in Figure 2, we model each reactant molecule as a rigid sphere with radii r1 and r2 for the siRNA and RISC molecules respectively. We define our coordinate system such that the RISC molecule is stationary with respect to the siRNA molecule for the reaction, so that the siRNA moves towards the RISC molecule with a relative velocity U12. The siRNA molecule moves through the cell cytoplasm to sweep out a collision cross section A = π r122 (as illustrated in Figure 3), where r12 is the collision radius given by: r12 = r1 + r2. If the center of a RISC molecule comes within a distance of r12 of the center of a siRNA molecule, they will collide which might result in a successful reaction.
Figure 2. Schematic diagram of siRNA and RISC molecules.
Figure 3. Volume swept out by the siRNA molecule in time Δt.
To discretize the system, we consider the dynamics of this process within a small time Δt. We assume that the temporal reaction process is an independent sequence of events separated by Δt. In time Δt, the siRNA molecule sweeps out a volume ΔV given by: ΔV = π r122 U12 Δt. Note that, in Figure 3 the siRNA molecule actually sweeps out a cylindrical volume that allows us to estimate ΔV as the length of the cylinder in time Δt, thus is given by U12 Δt.
Now, the probability of a siRNA molecule being present in the collision volume ΔV is psiRNA = 1. This is because we have already assumed that one siRNA molecule entered the cell creating a collision volume of ΔV.
Probability of at least one molecule of RISC being present in an arbitrary uniformly distributed ΔV in V is pRISC = ΔV.n2/V, where V denotes the cell volume. Ideally V should be the volume of the cytoplasm, which can be approximated by the entire cell volume. The probability that a siRNA molecule collides with a RISC molecule in Δt is given by:
Thus we have a stochastic sequence of events characterized by the probability of collision, and it is important to determine whether the collision will create the reaction. To complete the reaction, the molecules have to bind to each other. Different types of bonds (ionic, covalent, hydrogen etc.) require different activation energies for binding. We next assume that the colliding molecules must cross an energy threshold, defined by the free energy EAct, to provide the energy to react. Also, we assume that only the kinetic energy directed along the line of centers of the two reacting molecules contribute to the reaction as the effects of other forces (e.g. coulomb force) can been captured by the velocity distribution of the siRNA molecules in the cell.
These two assumptions define the probability of another independent event: successful reaction after collision denoted by pr. The kinetic energy of approach of a siRNA towards the RISC molecule with relative velocity U12 is E = m12.U122/2, where m12 = m1.m2/(m1+ m2) = the reduced mass, m1 = mass (in gm) of a siRNA molecule and m2 = mass (in gm) of a RISC molecule. We also assume that as E increases above EAct, the number of collisions that result in reaction also increases linearly. Thus the probability for a reaction to occur, pr, is given by:
and hence, the joint probability, p, for collision and reaction is given by:
Until now we were working with a fixed relative velocity U12 for the siRNA molecules. The velocity distribution of the macromolecules inside a cell capturing the effects of the different forces as obtained from Molecular Dynamic Simulation is generally found to be comparable to the Maxwell-Boltzmann distribution [18,19]. Hence, to approximate the relative velocity of the siRNA molecules, we use the Maxwell-Boltzmann distribution of molecular velocities for a species of mass m given by:
where kb = Boltzmann's constant = 1.381 × 10-23 kg-m2/s2/K/molecule and T is the absolute temperature at which the reaction occurs. Replacing m with the reduced mass m12 of the molecules, we get,
The term on the left hand side of the above equation denotes the fraction of siRNA molecules with relative velocities between U and (U+dU). Summing up the collisions for the siRNA molecules for all velocities we get the probability of reaction, p, as a function of temperature only as follows:
Now, recalling E = m12.U122/2, i.e., dE = m12U12dU and substituting into the expression for f(U, T)dU, we get:
Thus we get:
To consider a certain concentration of siRNA molecules, we assume that n1 number of siRNA molecules are present in the cell. This will increase the probability of a successful reaction by a factor of n1, and hence we have:
We discretize the temporal reaction process as a Bernoulli trial process. Next we compute the average time taken to complete the reaction with this probability. Let us assume that the molecule composition does not change during the reaction time. This is valid due to the very short time for reaction compared to the time taken for a potential change in the reaction environment for the associated molecules. The molecules try to react through repeated collisions. If the first collision fails to produce a reaction, they collide again after Δt time units and so on. We can interpret p as the probability of a successful reaction in time Δt. Thus the average time of reaction, Tavg, can be approximated by summing up the times taken for a successful reaction by the first collision, or that by the second collision and so on. Thus we get:
Similarly, the corresponding second moment, T2ndmoment, can be formalized by
This gives us the first and second moments of the reaction time, which is considered to be a random variable. The binding rate for this reaction can be easily estimated as 1/Tavg although the reaction is essentially stochastic. Note that, the expression for the binding rate from our model is exactly the same as that derived by chemical kinetic modeling of reaction rates. However, our model allows us to study the inherent stochasticity of the reaction considered.
Estimating the reaction rate constant for the siRNA molecules
The above model was used to study the binding rate of siRNA molecules. Note that the model requires an estimate of the activation energy for the siRNA-RISC complex formation reaction. However, it is very difficult to experimentally measure the activation energy required for a reaction. The activation energy is generally derived from the heat of reaction (-ΔH0) by using the Polyani equation . But it is equally difficult to measure the heat of reaction specifically for those occurring inside the cell.
We observed that the siRNA-RISC reaction occurs when the siRNA duplex breaks down into two single-stranded molecules one of which enters the reaction (the antisense strand), while the other (the sense strand) disintegrates. The antisense strand enters into a docking reaction with the RISC complex. Thus, we can rewrite the siRNA-RISC complex formation reaction in the following form:
Note that we treat the double stranded siRNA as a complex between the sense (siRNAsense) and antisense (siRNAantisense) strands. Referring to the reaction model discussed above, the two colliding macromolecules will be the siRNAsense siRNAantisense and RISC complexes. This also motivates the use of the Polyani equation for measuring the activation energy, which particularly works well for the following family of reactions:
Spectometric measurements on the thermodynamics of double-helix formation reported in Refs.  and  reveal the contribution of the 3' end in increasing the stability of the helix. Ref.  also relates the sequence dependence to the Gibb's free energy change to study the energetics of dangling ends and terminal base pairs in ribonucleic acid.
Hence, to measure the heat of reaction, we can use the heat of reaction required to dissociate the siRNA duplex. Note that, the second part of the reaction being essentially a docking reaction does not require any activation energy . However, instead of using the actual heat of reaction (measured in terms of the change in enthalpy for the reaction), we plan to use the Gibb's free energy change (ΔG370) that measures the "useful" or process-initiating work. The idea is that only the Gibb's free energy (which is a subset of the enthalpy or heat of reaction) is required for the overall reaction to occur, as the remaining heat will be generated by the docking part of the reaction. The docking reaction as mentioned before does not require any activation energy as it does not require changes in the chemical structure (i.e., no reaction is required), and in fact most binding events involve non-bonding interactions. As a consequence, such complexes may be formed by crossing small energy barriers that are within the range of thermal fluctuations (kT range at physiological temperature), and would not necessarily require a large kinetic energy during the collision. However, the electrostatic and solvation effects are the driving force for the formation of macromolecular complexes, and hence the docking part of the reaction. Hence, the activation energy, EAct will have the following components:
where, Edissociation is the energy required for the dissociation of the double-stranded siRNAs while Eelectrostatic and Esolvationcorrespond to the electrostatic and solvation energy requirements for the docking reaction.
The Polyani equation for exothermic reactions is used for estimating Edissociation from the Gibb's free energy measurements as follows:
A recent study on siRNA duplex stability  gives the thermodynamic energy requirements for 6-nt long siRNA duplex stability. We used the Gibb's free energy change measured by them for different types of residues present in the 2-nt overhangs to generate the binding rates predicted by our model. An 8-fold difference in the reaction rate was observed suggesting a purine residue followed by a purine/pyrimidine residue in the overhang results in an 8-fold slowdown of the reaction rate (the first 20 data points in Figure 4) compared to a pyrimidine residue followed by a purine/pyrimidine residue which supports the findings in Ref. . The other parameters from the model (e.g., mass and radii of siRNAs and RISC enzymes) were arbitrarily assumed as we are interested in the relative difference in the binding rates for the different siRNAs instead of the actual values. Also, it should be noted that the mass and radii of the different siRNAs in Ref.  will be comparable and will not affect the binding rate appreciably.
Figure 4. Binding rate constant plot.
Ref.  also allows us to find out the distribution of the reaction time for the siRNA molecules reported from the Gibb's free energy estimates reported in Table 1 in Ref. . In Figure 5, we plot the standard deviation-to-mean ratio of the reaction time from our model using the above-mentioned activation energy estimate. We observe that the standard deviation and mean times for reaction are the same, and hence the reaction time follows an exponential distribution. Thus, we can use the standard Gillespie stochastic simulator for studying the RNAi system where each reaction is considered stochastic, with the time for reaction completion following an exponential distribution.
Table 1. Simple RNAi system equations
Figure 5. The mean and standard deviation of reaction time are fairly equal and hence the reaction time is assumed to follow an exponential distribution.
With the above findings, we estimated the Gibb's free energy change for dissociating the 21-nt siRNA duplexes of the following two types: 1) 21-nt siRNAs with 2-nt overhangs and 2) 21-nt blunt-ended siRNAs. The corresponding ΔG370 estimates are -21.2 kcal/mol and -20.4 kcal/mol respectively. The experimental setup is explained in the Appendix. This results in a 2-fold difference in the rate constant using these two types of molecules with the blunt-ended ones having a lower rate constant than their 2-nt overhang counterparts (i.e. Tavgblunt-end/Tavg2-nt = 2, where Tavgblunt-endand Tavg2-nt are the average reaction times for blunt-ended and 2-nt siRNAs respectively). Again we have assumed that the mass and radii of these two types of siRNAs are comparable and only the activation energy parameter from our reaction model makes the rate constant different. Also, note that the cell volume parameter cancels out as we only need to compute the ratio (Tavgblunt-end/Tavg2-nt) and not the actual rate constants. Similarly, we have assumed that Eelectrostatic and Esolvation will be the comparable for the two types of siRNAs and will cancel out as we compute Tavgblunt-end/Tavg2-nt. Note that Eelectrostatic can be computed using the siRNA 3-d structures using standard software like Delphi  but will be very similar for the siRNAs that we have considered (which only differ by 2-nt at the 3' end). We next use this predicted difference in rate constants to compare the potency effects of these two types of molecules in an entire RNAi system simulation. Note that, because of the unavailability of some of the parameters required in the model (that we have assumed to be comparable for the two types of siRNAs), we have used the rate constant reported in Ref  for 2-nt siRNAs directly in the simulation model. To get the rate constant for blunt-ended siRNAs, we halved the rate constant as predicted by our reaction model.
siRNA duplex stability measurements
Blunt-end and 2-nt 3' overhang siRNA's were reconstituted in buffer at 20 μM. Stock solutions were then diluted in buffer to 10 μM for data collection and determination of the thermodynamic parameters for siRNA unfolding. 120 μL of solution was added to a micro-volume UV cuvette with a pathlength of 1.0 cm. The cuvette was placed in an AVIV UV/Vis spectrophotometer with a Peltier temperature controller. A reference cuvette filled with buffer solution was placed in the reference beam. The temperature was lowered to 10°C and the absorbance at 280 nm was auto zeroed. Temperature was then increased at 1°C/min and the change in absorbance recorded until the temperature reached 95°C. Absorbance data were exported to Meltwin 3.5 software  for data analysis. Data were fit to the self-complimentary model available in Meltwin 3.5.
Stochastic simulation of the RNAi system
We primarily employed Dizzy  and the SimBiology toolkit from Matlab for the stochastic simulation of the RNAi system. Let us consider a simple RNAi system from Table 1. mRNA is transcribed according to reaction 4 and degraded according to reaction 6. dsRNA is synthesized from mRNA by RDR and is cleaved into 10 siRNAs according to reactions 3 and 2 respectively. siRNA can associate with mRNA according to reaction 1 to cleave the mRNA and produce garbage RNA (gRNA). For simplicity, we do not implement the formation of RISC explicitly in our model; instead, the siRNA-mRNA complex is directly degraded into aberrant gRNA. Thus, we assume that the RISC enzymes are not a rate limiting component of the RNAi system. Reactions 5 and 7 describe the degradation of siRNAs and aberrant garbage pieces, respectively. This simple RNAi system is motivated from  from where we get an estimate of the basic rate constants. We did not consider any nonspecific effects following the observation in Ref.  that siRNA-mediated inhibition of gene expression is generally independent of nonspecific interference pathways triggered by the dsRNAs. Our goal in this study is to present a systems biology framework for studying the potency effects of the two types of siRNA based on the difference in their rate constants. Thus, the rate constant of reaction 1 is 0.008 for siRNAs with 2-nt overhangs while it will be approximately 0.004 for the blunt-ended molecules.
Table 2 presents the initial number of molecules considered in the cell for the simulation. These values are just for illustrative purposes but serve well to compare the two types of siRNAs. In Figures 6 and 7 we report the change in concentration of the dsRNA and mRNA molecules with time (in hours) for the 2-nt overhangs and blunt-ended molecules respectively. Note that the dsRNA molecules get used up quite quickly (in less than 3 hours) in both cases depending on the siRNA production rate from dsRNAs and dsRNA degradation. However, for blunt-ended siRNAs, the mRNA levels go down to about 250 molecules (in about 3 hours) compared to about 150 molecules when the 2-nt overhang siRNAs are used. This is obvious as the siRNA-RISC complex formation reaction rate for blunt ended molecules is about half of that for 2-nt overhang siRNAs. Another point worth noting is that the mRNA levels return to control levels or untreated levels after a longer period when blunt-ended molecules are used. This can be attributed to the fact that the rate constant for the siRNA-RISC complex formation reaction being low for blunt-ended molecules, the blunt-ended siRNAs persist in the system for a longer time than the 2-nt 3' overhang counterparts. Figure 8 illustrates this phenomenon, where the mRNA levels using blunt-ended siRNAs seem to lag behind that observed for the 2-nt overhang species. The plot was generated with the initial concentration of dsRNA molecules set to 100 for illustrative purposes.
Table 2. Initial conditions
Figure 6. Concentration change using 2-nt 3'siRNAs.
Figure 7. Concentration change blunt-ended siRNAs.
Figure 8. mRNA levels saturate after longer time while blunt-ended siRNAs are used.
Next we consider the transcription/translation machinery to study the potential effects of stochastic mRNA production on the RNAi system . The additional reactions (7–18) were added to consider the transcription/translation processes (presented in Table 3) and the initial conditions are reported in Table 4. A single gene is transcribed into mRNA by RNA polymerase (RNAP). The process is initiated with the binding of RNAP to the promoter, usually near the beginning of the transcribed sequence. Expression of most genes are regulated at the level of transcription and more specifically during the initiation of transcription, that is, before the first phosphodiester bond is formed.
Table 3. RNAi system with transcription/translation machinery
Table 4. Initial Conditions for the combined RNAi system
Figures 9 and 10 report the mRNA and dsRNA concentration changes with time for this combined RNAi system. We observe similar effects (both quantitative and qualitative) as in the simple system without considering the transcription/translation reactions. Note that the simple system considered a certain mRNA production rate and does not capture the uneven mRNA production observed recently . Because of stochastic gene expression in eukaryotic cells, mRNA and hence protein production is erratic and occurs in bursts. The average intervals between bursts is longer in eukaryotes compared to prokaryotes . Thus, the stochasticity in gene expression does not have any observable effects in the RNAi system. One plausible reason behind this is the fact that the stochasticity is observed at a much lower time scale (generally in seconds) while we in this simulation, are interested in studying the mRNA concentrations over a longer time scale (hours). Initial findings show about 1.5 times more blunt-ended molecules will be required to keep the mRNA at a certain low level. Figure 11 plots the mRNA levels for different types of siRNAs in order to quantify the gene silencing period. With 10,000 dsRNA molecules for the two types of siRNAs, the silenced period is larger for 2-nt overhang siRNAs. Using 15,000 dsRNAs of blunt-ended type we can achieve the same amount of silencing as the 2-nt overhang siRNAs.
Figure 9. Combined RNAi system plot using 2-nt 3'siRNAs.
Figure 10. Combined RNAi system plot using blunt-ended siRNAs.
Figure 11. Quantification of gene silencing period with different types of siRNAs.
Next consider the effects of the RDR mechanism of replenishing the dsRNAs in the cell. Using the simple RNAi system (of Table 1), as the transcription/translation machinery failed to show any significant changes in the system behavior. Table 5 lists the additional reactions to incorporate the RDR mechanism. Reaction 8 describes the unprimed amplification – the synthesis of dsRNA from aberrant garbage RNA by RDR; and reaction 9 describes primed amplification – the synthesis of dsRNA primed by the presence of a siRNA on mRNA. We consider the pathway with the amplification terms at both low and high amplification rates (Figures 12 and 13 respectively). Thus sustained gene silencing can only occur through the RNAi pathway if there is a way of replenishing the dsRNAs in the cell. This can occur through the RDR mechanism (with high amplification as illustrated in Figure 13), or through artificial dsRNA injection into the cell (Figure 14) at certain time periods. There is no conclusive evidence of the presence of the RDR mechanism in mammalian cells, and hence the artifical dsRNA injection will play a very important role in any RNAi based therapeutic tool.
Table 5. Additional reactions from Table 1 to consider the RDR mechanism
Figure 12. 2-nt 3' and blunt-ended siRNAs with RDR show the mRNA/dsRNA concentration changes with time at low amplification through RDR for the two types of siRNAs considered. We do find an expected increase in the potency effects for both types of siRNAs (compared to Figures 6 and 7).
Figure 13. 2-nt 3' and blunt-ended siRNAs with RDR. Figure 13 shows the effects of high amplification through RDR and we can observe sustained gene silencing. The actual mRNA levels will however depend on the amplification rates and initial number of dsRNA molecules inserted into the cell.
Figure 14. 2-nt 3'and blunt-ended siRNAs with artificial dsRNA injection. Fiureg 14 show the effects of artificial dsRNA injection (for the simple RNAi system of Table 1) for the two types of siRNAs with 200 dsRNA molecules inserted every hour into the cell.
There is a slight difference in the mRNA levels with the 2-nt overhang siRNAs fairing a little better than their blunt-ended counterparts. However, we can get a relative quantification of the number of dsRNA molecules required to be inserted into the cell for two types of siRNAs considered to keep the mRNAs down at the same level. This will require more real-life initial concentration values for the other components of the RNAi system and also a model to predict the average number of dsRNA molecules that actually enter a single cell depending on the dosage amount and intervals. Our model provides a simulation framework that will help us study the gene silencing duration using siRNAs in the future.
Handling delayed reactions
We have computed the time taken to complete the reaction. The underlying assumption is that reactant collisions occur with some probability and once a collision of sufficient energy occurs, a reaction takes place instantaneously. Hence, we assume that there is no time delay to form an activated complex. If there is some time delay associated with initiation and completion of the reaction, the probability evolution becomes more complicated . Our reaction model cannot directly handle such delayed reactions, which would require comprehensive modeling of the delayed states. However, the reaction model in this paper was primarily developed for the siRNA-RISC complex formation reaction, which does not require the handling of delayed reactions. However, this involves an implicit approximation that the reactant molecules are not available for other reactions once it has entered the present reaction event. But because both the original reaction event time and the delay can be random variables incorporating the probability of successfully completing the delayed reaction event, this approximation should be small. Further analysis is required to study the effect of this approximation.
Maxwell-Boltzmann distribution of molecular velocities
The Maxwell-Boltzmann distribution gives a good estimate of molecular velocities where we have spatial homogeneity and is widely used in practice. Molecular dynamic (MD) simulation measurements during protein reactions show that the velocity distribution of proteins in the cytoplasm closely matched the Maxwell-Boltzmann distribution [18,19]. However, its application in our model may not give the predicted results for cases that violate the assumption of uniform distribution in a volume. Ideally the velocity distribution should incorporate the properties of the reaction space (nucleus/membrane/cytoplasm for reactions occurring in the membrane, nucleus or cytoplasm) and the effect on velocity distribution due to its space shape and irregularities. We plan to explore the possibility to improve this velocity distribution by considering the other biological factors that can influence the velocity of the reacting molecules.
Activation energy threshold
The activation energy has been measured for many reactions and we need an estimate of this parameter to be able to predict the nature of the reaction time. We used the Polyani equation to compute Edissociation for the reaction by measuring the Gibb's free energy change to dissociate a siRNA duplex. As discussed before, this approximation can give us a relative difference in the reaction rates for the two types of siRNAs used in our study, but will fail to give us a direct quantification of the actual reaction rates (which will also require estimates of Eelectrostatic, Esolvation, mass, radii etc. for the siRNAs). We are exploring ways to measure the actual heat of reaction for the entire siRNA-RISC complex formation reaction and the other parameters to make the model predictions more accurate.
We did not consider the reverse reaction conditions in our model because we assumed the siRNA-RISC complex formation reaction as non-reversible. However, the Gillespie stochastic simulation framework allows us to incorporate reversible reactions if we know the forward and backward reaction rates . We plan to make the simulation framework presented in this paper more comprehensive once we have experimental measure of these reaction rates.
In addition, there is increasing evidence of sub-compartmental (i.e., intra-compartmental localization) in cells, so local neighborhoods of reactions will have higher apparent concentrations than simply the number of molecules divided by the size of the compartment. However, this would require more in depth modeling of the different molecular concentrations inside the cell that reduces the scalability of the stochastic simulation framework. Indeed, the Gillespie simulator fails to address this issue as it requires different rate constants for specific neighborhoods of the reaction type. Nevertheless, our model can be easily extended to incorporate such reaction neighborhoods by limiting the movement of the reactant molecules inside a reaction space while computing the probability. However, the applicability of the Maxwell-Boltzmann velocity distribution in the neighborhood requires further research.
Role of GW/P-bodies in RNA silencing
A related concern with this simulation model is the recent finding that RISC and other enzymes involved in RNA degradation tend to localize to discrete cytoplasmic foci known as P-bodies [32,33]. Moreover, siRNAs were also observed to localize to specific cytoplasmic compartments in the periphery of the nucleus in granular-like structures . This compartmentalization will have an effect on the reaction model that we have developed based on the entire cell volume. However, it should be noted that the reaction model is only used to find the relative difference in the reaction rates between the two types of siRNAs and not to predict the overall reaction rate. We have used the reaction rates reported in Ref  in our simulation model and changed it accordingly for blunt ended siRNAs using the reaction model. Also, finding the relative difference in the reaction rates would mean the cancellation of the parameter V, and hence is independent of the effects of this compartmentalization. It can be argued that the simulation model however requires additional reactions to address the effect of GW/P-bodies (e.g., by adding additional delay for siRNAs to reach the cytoplasmic foci and additional reactions for forming the complex with GW/P-bodies). This will require additional research to estimate the kinetic parameters for these reactions. In this paper, however, we can neglect the effect of GW/P-body formations because the reaction rate mentioned in Ref  is for the overall siRNA-mRNA complex formation reaction and should already incorporate the effects of this compartmentalization.
27-mer RNA duplexes show higher potency than 21-mer siRNAs
Many researchers today employ synthetic 21 mer RNA duplexes as their RNAi reagents, which mimic the natural siRNAs that result from Dicer processing of dsRNAs. An alternative approach is to use synthetic RNA duplexes that are greater than 21 mer length which are substrates for Dicer. These duplexes are typically 27-nt long and are processed by Dicer into 21 mer siRNAs. It has been reported that synthetic Dicer-substrate RNAs can have significantly increased potency (~100-fold) when compared with 21 mer duplexes [35,36]. The 100-fold increase in potency of the 27 mer is due to a combination of Dicer cleavage resulting in "a better 21 mer" (10-fold increase) plus some other effect that required use of the intact 27 mer (another 10-fold increase).
More recent works targeted additional sites in other genes, to find examples where the longer RNAs had greater potency, the same potency, and lower potency than 21 mers at the same site. This wide variation in performance was primarily attributed to the differences in dicing patterns: sometimes Dicer processing resulted in a "good" 21 mer while other times it resulted in a "bad" 21 mer (Ref ). Small shifts in sequence can have a large effect on siRNA potency. Combination of asymmetric 2-base 3'-overhang with 3'-DNA residues on the blunt end result in a duplex form which directs dicing to predictably yield a single primary cleavage product. Using this strategy, strand targeting experiments in Ref  show that Dicer processing confers functional polarity within the RNAi pathway.
In this work, we have not studied the effects on potency while using synthetic 27 mers instead of the dsRNAs. This will require the ratios of the "good" and "bad" 21 mers that are diced from the 27 mers entering the cell as well as a similar thermodynamic stability study of these 21 mers to study their effects on the rate constant estimates.
The RNAi system simulation framework developed here presents some qualitative and quantitative results on the RNAi pathway. The simulation also explores the potency effects of the two types of siRNAs considered. We primarily focused on identifying the most important components of the system that play a role in identifying the potency effects. This framework can be extended to predict very important features of the RNAi pathway including the duration of gene silencing for the different siRNA molecules. This will require models to study the effective rate of dsRNA entry into the cell depending on the amount of dosage and dosage intervals. We plan to study these important aspects of the RNAi system, once we build more fidelity on this simulation framework.
The proposed model computes the reaction time for the siRNA-RISC complex formation reaction as a stochastic variable that appropriately reflects the cell environment. The concept of the model is to transform the reaction process from a continuous deterministic process to a discrete random one. The model allows the transformation of biological reactions to the stochastic domain and makes it suitable for a stochastic simulation of the RNAi system. The average reaction time estimated from this method is exactly the same as the reaction rate estimates of kinetic modeling. In addition, we are able to estimate the first two moments of the reaction time to capture the stochastic nature of the reaction. The reaction model was used to study the difference in binding rate for the 2-nt 3' overhang siRNAs and the blunt-ended siRNAs, and we found that the reaction rate is predicted to double when 2-nt 3' overhang siRNAs are used. We next built an RNAi stochastic system simulation using the Gillespie simulator to study the overall potency effects of the two types of siRNAs. Initial findings suggest that about 1.5 times more blunt-ended siRNAs are required to keep down the mRNA at the same level as using 2-nt 3' overhang siRNAs. The additional blunt-end siRNAs may be needed because the siRNA-RISC complex formation reaction is not the only part of the RNAi pathway that is affected by the presence of the 2-nt overhangs in the siRNA molecules. The reaction in which the bound RISC complex cleaves the mRNA might be affected more due to the difference in the siRNA structures. Further research is required to study this aspect of the RNAi system.
The stochastic simulation framework presented here allows us to predict some quantitative and qualitative aspects of the RNAi system for the two different types of siRNAs used in our study. More importantly, it allows us to build a quantitative framework for studying the length of the gene silencing period and number of siRNA molecules (of both types) required for the same for a cost-potency study. This will require a model to estimate the average number of dsRNA molecules actually entering the cell depending on the dosage interval and amount. The initial predictions from our work are promising, as we plan to build a computational tool for studying the therapeutic effects of the RNAi pathway.
The authors declare that they have no competing interests.
All authors designed, analyzed, implemented and tested the proposed models. Each author contributed equally in writing the paper. All authors read and approved the final manuscript.
This article has been published as part of BMC Genomics Volume 10 Supplement 1, 2009: The 2008 International Conference on Bioinformatics & Computational Biology (BIOCOMP'08). The full contents of the supplement are available online at http://www.biomedcentral.com/1471-2164/10?issue=S1.
Nat Biotechnol 2004, 22:326-330.
Epub 2004 Feb 1.PubMed Abstract | Publisher Full Text
Nucleic Acids Res 2004, 32:691-699.
Print 2004.PubMed Abstract | Publisher Full Text | PubMed Central Full Text
Nucleic Acids Res 2005, 33(13):4140-56.
Print 2005.PubMed Abstract | Publisher Full Text | PubMed Central Full Text
Fogler H, Gurmen M: [http://www.engin.umich.edu/~cre/03chap/html/collision/] webcite
Biochemistry 1987, 26:4554-4558. PubMed Abstract
Proc Natl Acad Sci USA 2001, 98:9742-9747.
Epub 2001 Jul 31.PubMed Abstract | Publisher Full Text | PubMed Central Full Text
Biochimie 2007, 89(10):1245-51.
Epub 2007 May 17.PubMed Abstract | Publisher Full Text
Nature 2004, 431(7006):343-9.
Review.PubMed Abstract | Publisher Full Text
Nat Biotechnol 2005, 23(2):222-6.
Epub 2004 Dec 26.PubMed Abstract | Publisher Full Text