Comparing 2-nt 3' overhangs against blunt-ended siRNAs: a systems biology based study
© Ghosh et al. 2009
Published: 7 July 2009
Skip to main content
© Ghosh et al. 2009
Published: 7 July 2009
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.
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.
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.
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).
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 = π r12 2 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.
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.
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.
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.
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.
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.
Simple RNAi system equations
mRNA + siRNA -> gRNA
dsRNA -> 10*siRNA
mRNA -> dsRNA
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 ΔG37 0 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. Tavg blunt-end/Tavg 2-nt = 2, where Tavg blunt-endand Tavg 2-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 (Tavg blunt-end/Tavg 2-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 Tavg blunt-end/Tavg 2-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.
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.
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.
Number of molecules
RNAi system with transcription/translation machinery
mRNA + siRNA -> gRNA
dsRNA -> 10*siRNA
mRNA -> dsRNA
O + R -> O_R
Regulatory molecule R binds to the operator region O to form the bound complex O_R
O_R -> O + R
O_R dissociates into free R and O
P + RNAP -> P_RNAPC
RNAP binds to promoter region P forming closed complex P_RNAPC
P_RNAPC -> P + RNAP
P_RNAPC dissociates into free RNAP and P
P_RNAPC -> P_RNAPO
Isomerization of closed to open complex (P_RNAPO)
P_RNAPO -> TrRNAP + mRNA + P
RNAP clears promoter region and mRNA chain synthesis starts. TrRNAP denotes transcribing RNA polymerase.
TrRNAP -> RNAP
RNAP completes transcription and is released from DNA.
mRNA + Ribosome -> RibRBS
Ribosome binds to mRNA forming bound complex RibRBS
RibRBS -> mRNA+ Ribosome
Ribosome dissociation from RibRBS
RibRBS -> EIRib + mRNA
Ribosome EIRib initiates translation of mRNA chain
EIRib -> protein
Protein synthesis by transcribing ribosome
Protein product degradation
Initial Conditions for the combined RNAi system
Number of molecules
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.
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.
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.
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.
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.
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.
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.
This article is published under license to BioMed Central Ltd. This is an open access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/2.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.