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 (p_c) 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 (p_r).

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, n₂, 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 [17] to model the chemical kinetics of the reaction. As shown in Figure 2, we model each reactant molecule as a rigid sphere with radii r₁ and r₂ 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 U₁₂. The siRNA molecule moves through the cell cytoplasm to sweep out a collision cross section A = π r₁₂² (as illustrated in Figure 3), where r₁₂ is the collision radius given by: r₁₂ = r₁ + r₂. If the center of a RISC molecule comes within a distance of r₁₂ of the center of a siRNA molecule, they will collide which might result in a successful reaction.

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 = π r₁₂² U₁₂ Δ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 U₁₂ Δt.

Now, the probability of a siRNA molecule being present in the collision volume ΔV is p_siRNA = 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 p_RISC = ΔV.n₂/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:

(1)

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 E_Act, 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 p_r. The kinetic energy of approach of a siRNA towards the RISC molecule with relative velocity U₁₂ is E = m₁₂.U₁₂²/2, where m₁₂ = m₁.m₂/(m₁+ m₂) = the reduced mass, m₁ = mass (in gm) of a siRNA molecule and m₂ = mass (in gm) of a RISC molecule. We also assume that as E increases above E_Act, the number of collisions that result in reaction also increases linearly. Thus the probability for a reaction to occur, p_r, is given by:

(2)

and hence, the joint probability, p, for collision and reaction is given by:

(3)

Until now we were working with a fixed relative velocity U₁₂ 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:

(4)

where k_b = Boltzmann's constant = 1.381 × 10^-23 kg-m²/s²/K/molecule and T is the absolute temperature at which the reaction occurs. Replacing m with the reduced mass m₁₂ of the molecules, we get,

(5)

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:

(6)

Now, recalling E = m₁₂.U₁₂²/2, i.e., dE = m₁₂U₁₂dU and substituting into the expression for f(U, T)dU, we get:

(7)

Thus we get:

(8)

To consider a certain concentration of siRNA molecules, we assume that n₁ number of siRNA molecules are present in the cell. This will increase the probability of a successful reaction by a factor of n₁, and hence we have:

(9)

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, T_avg, 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:

(10)

Similarly, the corresponding second moment, T_2ndmoment, can be formalized by

(11)

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/T_avg 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 (-ΔH⁰) by using the Polyani equation [17]. 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 (siRNA_sense) and antisense (siRNA_antisense) strands. Referring to the reaction model discussed above, the two colliding macromolecules will be the siRNA_sense siRNA_antisense 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. [20] and [21] reveal the contribution of the 3' end in increasing the stability of the helix. Ref. [22] 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 [23]. 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 (ΔG₃₇⁰) 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, E_Act will have the following components:

(12)

where, E_dissociation is the energy required for the dissociation of the double-stranded siRNAs while E_{electrostatic} and E_solvationcorrespond to the electrostatic and solvation energy requirements for the docking reaction.

The Polyani equation for exothermic reactions is used for estimating E_dissociation from the Gibb's free energy measurements as follows:

(13)

A recent study on siRNA duplex stability [24] 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. [24]. 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. [24] will be comparable and will not affect the binding rate appreciably.

Ref. [24] 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. [24]. 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.

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 ΔG₃₇⁰ 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. T_avg^blunt-end/T_avg^2-nt = 2, where T_avg^blunt-endand T_avg^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 (T_avg^blunt-end/T_avg^2-nt) and not the actual rate constants. Similarly, we have assumed that E_{electrostatic} and E_solvation will be the comparable for the two types of siRNAs and will cancel out as we compute T_avg^blunt-end/T_avg^2-nt. Note that E_{electrostatic} can be computed using the siRNA 3-d structures using standard software like Delphi [25] 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 [26] 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 [27] for data analysis. Data were fit to the self-complimentary model available in Meltwin 3.5.

We primarily employed Dizzy [28] 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 [26] from where we get an estimate of the basic rate constants. We did not consider any nonspecific effects following the observation in Ref. [29] 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.

Next we consider the transcription/translation machinery to study the potential effects of stochastic mRNA production on the RNAi system [30]. 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.

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 [16]. 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 [16]. 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.

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.

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 [31]. 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 E_dissociation 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 E_{electrostatic}, E_solvation, 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 [15]. 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 [34]. 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 [26] 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 [26] 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 [14]). 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 [14] 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.