BACKGROUND
rna interference  is a well characterized regulatory mechanism in eukaryotes  <cit>  as well as a powerful tool for understanding gene function, thanks to the discovery that synthetic small interfering rna oligomers  can efficiently induce rnai in mammalian cells  <cit> .
rnai has also been used extensively as a novel "biological part" to design synthetic biological circuits in synthetic biology  <cit> .
artificial gene silencing has the potential to become a major genetic-based therapeutic tool for viral infections  <cit> , cancer  <cit>  or inherited genetic disorders  <cit> .
despite its widespread experimental application, the best way to quantitatively model rna interference is still under debate.
in systems and synthetic biology, mathematical models are essential to carry out in silico investigations of biological pathways, or novel synthetic circuits.
the aim of this work is to find the most appropriate quantitative mathematical model that can correctly describe the rnai phenomenon in mammalian cells, for varying concentrations of the sirna oligomers.
a schematic representation of the rna interference mechanism is illustrated in figure  <dig>  in step  <dig>  the presence of double stranded rna  elicits a response in the cell mediated by the dicer enzyme, which binds and cleaves the dsrna into fragments of 21- <dig> base pairs, called small interfering rna .
in step  <dig>  sirnas are loaded onto a multiprotein complex called rna induced silencing complex  and then separated into single strands of which one  is discarded and degraded  <cit> , while the guide strand remains within risc and serves as a template in the silencing reaction.
in step  <dig>  the guide strand assembles into a functional sirna-risc complex, which contains the sirna bound to the ago protein  <cit> .
target mrnas are then recognized by watson-crick base pairing  <cit>  and bound by the sirna-risc complex.
finally, in step  <dig>  mrna degradation is induced, target mrna is dissociated from the sirna, and the sirna-risc complex is released to process further mrna targets  <cit> .
here, we will focus on quantitatively modeling step  <dig> to step  <dig> of the rna interference process using, as an experimental tool, synthetic sirna oligomers.
RESULTS
in order to model the effects of rna interference on mrna expression levels at different concentration of sirna oligomers, we carried out in-vivo experiments of rna interference on two mammalian cell-lines stably expressing the egfp protein or the tta protein, respectively.
in the first set of experiments , human embryonic kidney cells stably expressing egfp , were transfected with varying quantities of synthetic sirna oligomers directed against the egfp mrna in the range of  <dig> to  <dig> pmol.
figure  <dig> shows the ratio of egfp mrna levels between treated cells and negative control cells  measured  <dig> hours after transfection.
error-bars represent the standard-error from three biological replicates for each point.
similarly, in the second set of experiments , egfp protein levels were measured via facs analysis  <dig> hours after transfection of hek293-egfp cells with synthetic sirna oligomers directed against the egfp mrna, or with non-specific sirna oligomers as negative controls.
the ratio between egfp protein levels in sirna-treated cells versus negative controls, for the same sirna quantities as in the set i experiments, is reported in figure  <dig>  in the third set of experiments , we tested chinese hamster ovary cells  stably expressing the tetracycline-regulated transactivator tta at a low level, using the same protocol used previously for hek cells.
cells were transfected with varying quantities of synthetic sirna oligomers against the tta mrna in the range of  <dig> to  <dig> pmol.
since the tta is not fluorescent, no facs measurement were performed in this cell-line.
in additional file  <dig> figure 1we shows the ratio of tta mrna levels between cells treated with the silencing tta oligomer and cells treated with the negative control .
rnai modeling
we were interested in formulating a model that can quantitatively describe the effects of varying quantities of sirna oligomers onto the degradation of the target mrna species, and of its corresponding protein product.
a general dynamical model describing transcription of the mrna species, its sirna-mediated degradation, and translation of its protein products, can be described by a system of ordinary differential equations .
let xm, xp and xs be the mrna, protein and sirna concentrations, respectively.
the evolution of their time-dependent concentrations can be described by the following odes:
:  dxmdt=km−dmxm−δ
:  dxpdt=ktxm−dpxp,
the parameter km, represents the transcription rate from the promoter that transcribes the mrna species targeted by the sirna oligomer; dm represents the basal degradation rate of the mrna species.
rnai can be considered as a mechanism that enhances the degradation of the targeted mrna, therefore the function δ  is an extra degradation term that represents the rate at which mrnas are degraded due to rnai.
this function, δ , depends on both the mrna and sirna levels, xm and xs respectively.
the parameter kt is the protein translation rate, whereas dp represents the basal protein degradation rate.
at least four different models have been proposed in the literature, for the rna interference mechanism  <cit> .
each of these models is based on the general approach described by equations  but each has a different functional form for δ .
table  <dig>  lists all the different models studied here with a description of their corresponding parameters.
c3: cleavage and dissociation rate of mrna-sirna*
θ4: michaelis-menten like constant
model 1: the stoichiometric model
one possible way of modeling the effects of rnai on the mrna degradation is to consider a stoichiometric reaction between the sirnas and mrnas.
let sirna* denote the concentration of the sirna-risc complex, namely the fraction of the sirnas that are loaded into risc complexes .
then the stoichiometric reaction between the sirna* and the mrna  can be described as follows:
+←→∅.
namely, the sirna-loaded risc binds to the complementary mrna and then both are degraded.
according to this model, following the law of mass action, we predict that the sirna mediated degradation will be proportional to the product of the concentration of sirna oligomers and the targeted mrna species:
δ=k1xmxs
where parameter k <dig> represents the proportionality constant.
in this modeling approach, the sirna-risc complex is assumed not to be recycled, but the risc needs to be reloaded before it can degrade another mrna molecule .
indeed, this model was suggested for rna interference in prokaryotes  <cit> .
rnai in prokaryotes has important differences with rnai in eukaryotes.
one of these is that the interaction between the sirna and its mrna target is non-catalytic in nature  <cit> .
this is not the case in mammalian cells.
once the mrna is cleaved, the sirna-risc complex is dissociated from it, in an atp-dependent manner  <cit>  and it is free to process further targets  <cit> .
additionally, in some organisms, such as c.elegans, the primary dsrna trigger induces synthesis of secondary sirnas  through the action of rdrp enzymes, strengthening and perpetuating the silencing response  <cit> .
model 2: stoichiometric model with co-operativity
this model is a straightforward extension of model  <dig>  which additionally takes into account the presence of multiple sites on the targeted mrna where the sirna-loaded risc can bind.
model  <dig> can be easily extended to include co-operativity:
++...+︸h2↔mrna−h2sirna*→∅
as before, the rate of rnai-driven degradation, can be easily obtained applying the law of mass action:
δ=k2xsh2xm
where k <dig> is the proportionality constant and h <dig> is the number of sirna binding sites on the targeted mrna species.
this model was suggested in  <cit>  for modeling rnai by mirnas, since it is experimentally proven that multiple sites for the same mirna can boost target mrna repression  <cit> .
this model however suffers from the same limitations as model  <dig>
model 3: enzymatic model
a detailed model of rnai specific for mammalian cells was proposed by malphettes et al  <cit> .
this model accounts for the catalytic nature of rnai in mammalian cells, thus modeling step  <dig>  step  <dig> and step  <dig> of figure  <dig>  including the recycling of the risc complex .
specifically, it assumes that once the cleavage and degradation of the targeted mrna is completed, the sirna-risc complex  dissociates from its mrna target, and it is free to degrade further mrna molecules.
additionally, this model also considers that the targeted mrna may have a multiple number, h <dig>  of sirna binding sites.
the sirna-risc complex  can bind to any site on the mrna to form an intermediate mrna-sirna* complex, which can either accommodate further sirna-risc complexes on any other free binding sites, or cleave the target mrna and dissociate from the cleavage products.
the reaction of the complex formation of the target mrna with the sirna-risc complex is described as follows :
+↔.
the model considers the following reaction between an intermediate mrna-sirna* i- <dig> complex with another sirna-risc complex :
+↔.
the generic cleavage and degradation reaction of the mrna by any interacting sirna-risc complex   +  → cleaved mrna+i.
in  <cit> , the authors show that the rate of rnai-driven mrna degradation, δ , is given by:
δ=c3k3h3c3+k3h3xmxsxm
where k <dig> is the rate of mrna-sirna* complex formation for a single sirna target site  and c <dig> is the cleavage and dissociation rate of mrna-sirna* complex .
this functional form, for a constant xs, is a classic michaelis-mentenu enzymatic reaction, as can be observed by simply defining vm = c3xs and km=c3k3h <dig>  and rewriting eq.
as:
δ=vmxmkm+xm
this is in perfect agreement with the experimental finding on the enzymatic activity of both non-mammalian and mammalian risc on mrna degradation, as reported in  <cit> , where the product of the reaction  was measured at varying the concentration of the substrate  for a constant amount of the risc enzyme.
model  <dig> basically assumes that only the maximal rate vm of the michaelis-menten will be affected when the concentration of the active risc  changes.
in  <cit>  the authors show experimentally that changing the concentration of risc indeed changes the value of vm but also of km.
this is not captured by this model.
note also that when xm ≪ km, i.e.
low amounts of mrna compared to the michaelis-menten constant km, the model can be simplified to: δ=vmkmxm=k3h3xmxs.
this is equivalent to model  <dig>  on the other hand, when the mrna concentration is very high, the model can be simplified to: δ  = vm = c3xs.
namely, the rate of rnai-driven mrna degradation will depend only on the amount of active enzyme, i.e.
the sirna concentration.
the higher the sirna concentration, the higher the degradation rate that can be achieved, without any saturation effect.
model 4: phenomenological model
in  <cit>  the authors proposed a standard hill-kinetic model to describe the post-transcriptional effects of micrornas on the gene expression.
mirnas are processed by the cell to produce dsrnas which then follow the typical rnai pathway, schematically described in figure  <dig> <cit> .
by considering a hill-type enzymatic model with an hill coefficient h ≥  <dig>  the model can be extended to account, either, for multiple binding sites of the sirna on the same target mrna, or, for the cooperativity of protein complexes involved in rnai  <cit> .
differently from the other three model presented above, this model is a phenomenological model not derived from specific biochemical reactions:
δ=d4xsh4θ4h4+xsh4xm
this model, despite being phenomenological has interesting properties.
the kinetic parameters d <dig> and θ <dig> depend on the efficiency of sirna binding to its sites on the target mrna  <cit> : d <dig> represents the maximal degradation rate of the mrna due to rna interference; θ <dig> the concentration of sirna oligomers needed to achieve half of the maximal degradation rate.
the above equation implies that for xs ≪ θ <dig>  the increase in the rnai mediated degradation is linear with xsh <dig> , while it saturates at higher levels of xs, reaching the maximal degradation rate d <dig>  differently from model  <dig>
parameter identification
the four models were fitted to the three mrna and protein experimental datasets , by searching for the parameter values for which the model-generated data best fitted the experimental data, according to a squared error measure.
the results of the fitting procedure for each of the models, together with the optimized values of their parameters, are given in table  <dig> for egfp and in additional file  <dig> table a <dig> for tta.
the relative value of the error  and the relative value of the prediction error  for each model are given, together with the corresponding optimized values of its parameters.
the unit of measurements are reported for the dimensional parameters.
.
the fitting results for the mrna levels  are shown in figure  <dig>  model  <dig> gives a significant smaller error than the other three models , in fact, compared to models  <dig> and  <dig>  the error of model  <dig> is two orders of magnitude smaller.
model  <dig>  gives a better fitting than models  <dig> and  <dig>  however it is worse than model  <dig>  the optimized value for the parameter d <dig> in model  <dig>  is d <dig> =  <dig>  min- <dig>  indicating that the strength of sirna mediated mrna degradation is comparable to the strength of basal mrna degradation .
note also that the parameters found for model  <dig> include a coefficient h <dig> =  <dig> , hence less than unity.
since h <dig> describes the number of sirna binding sites on the targeted mrna, it should be greater than, or equal to,  <dig> in order to have a clear biological interpretation.
however, if we constrain this parameter to be greater or equal to one, then the model optimizes at the value of h <dig> =  <dig>  which makes model  <dig> identical to model  <dig>
we have observed that in all numerical simulations, models  <dig> and  <dig> are almost indistinguishable.
the large optimized value of parameter c <dig> of model  <dig> , the low value of parameter k3h <dig> =  <dig>  × 10- <dig> suggest c <dig> ≫ k3h3xm, and hence the function δ  of model  <dig> can be approximated by k3h3xsxm, which is nothing else than model  <dig> .
therefore, in this parameter space, model  <dig> is almost identical to model  <dig>  similar results were obtained when searching for the parameters which best fitted the measured protein levels , shown in figure  <dig>  model  <dig> is again the one with the smallest error.
model  <dig> and  <dig> are again unable to capture efficiently the experimental data.
model  <dig> is better than models  <dig> and  <dig>  but it is still worse than model  <dig>
when we fitted the models to the third experimental dataset , which was performed on a different cell-line, with both a different target mrna and a different sirna oligomer, model  <dig> still performed better than the others with the smallest error, although model  <dig> was a close match.
model  <dig> and  <dig> were again behaving very similarly and had the largest error.
it should be noted that as it happened with the previous fitting results, also in this case model  <dig> has a hill coefficient smaller than unity .
finally we observed, as shown in additional file 1figure a <dig>  that the experimental error for set iii experiment was larger when compared to the egfp experiment .
this was caused by the relative low expression of the tta in the cho cell lines, when compared to the egfp expression in hek cells, which made measurements more noisy.
assessing the model predictive ability
the models described above have the same number of unknown parameters to be learned , but for model  <dig>  which has one extra parameter.
to be sure that the improved performance of model  <dig> in describing the experimental data was not due to overfitting, we computed for each model and for each experimental dataset, the prediction error, which allows to assess the generalisation performance of the models  <cit> .
we followed a "leave-one-out" cross validation strategy, where for each model and for each dataset, the parameter identification procedure was repeated each time removing one of the experimental points and then predicting the missing point with the identified parameters.
we thus could estimate a prediction error for each model in each experimental dataset.
this value is reported in table 2and additional file 1table a <dig>  model  <dig> is again the one with the smallest prediction error, once again confirming this model superior performance in describing the data.
CONCLUSIONS
our findings show that the simple hill function described by model  <dig> is sufficient to quantitatively describe the effect of rna interference, at the mrna and protein level, in mammalian cells in vitro, for varying concentration of sirna oligomers.
one significant feature of model  <dig> is that it can predict the saturation effect of the rnai process that we observed experimentally.
we considered the possibility that this saturation could be in fact due to the inability of the cell to uptake high concentration of sirna oligomers, however recent experiments  <cit> , prove that uptake of sirna oligomers in cells is linear with the concentration of sirna oligomers transfected, at least in the concentration range we used.
additionally, khan et al in  <cit> , observed upregulation of mrna targets of endogenous micro-rna when transfecting sirna oligomers in mammalian cells.
in order to explain this effect, they suggested a saturation of the risc complex .
it has been demonstrated in  <cit>  that the enzymatic activity of risc can be efficiently modeled in-vitro as a classic michaelis-menten reaction, where the target mrna is the substrate, the sirna-loaded risc is the active enzyme , and the product is the degraded mrna.
this is one feature that model  <dig> does not capture; namely for a fixed amount of sirna-loaded risc , eq.
should approximate the michaelis-menten in eq., instead model  <dig> becomes simply proportional to xm, as model  <dig> and  <dig>  nevertheless, model  <dig> approximates very well the experimental data.
we believe this happens because enzymatic reactions have a typical km much greater that the physiological concentration of their substrate   <cit> , and the same happens for the risc complex  <cit> .
in this condition, the michaelis-menten equation becomes vmxmkm+xm≊vmkmxm and therefore it is linear in xm, as predicted by model  <dig>  when the sirna concentration varies, model  <dig> predicts that the parameters vmkm will change as a function of xs as described by eq..
the three parameters of model  <dig> have a straightforward biological interpretation, and their values can be easily tuned to accommodate for different efficiencies of rnai.
for example, the parameter d <dig> can be used to weigh the degradation due to the rnai compared to the endogeneous mrna degradation, and its strength, i.e.
what is the maximal degradation rate that can be achieved.
θ <dig> quantifies the sirna oligomers concentration needed to achieve half of the maximal degradation of the targeted mrna.
the h <dig> coefficient can accommodate for multiple target sites on the same mrna, or for the cooperativity of the risc complex.
clearly, the rnai process is very complex and no one-to-one relationship can be found between parameters of model  <dig> and rnai biological components.
nevertheless, it has been shown in  <cit>  that between  <dig> and  <dig> sirna oligomers per cell  are sufficient to reach half-maximal mrna target degradation.
model  <dig> predicts that half-maximal degradation is achieved for an amount of sirna oligomers equal to θ <dig>  the value of this parameter when fitting mrna levels  is θ <dig> ≈  <dig>  pmol despite of the different cell-lines and mrna-sirna pairs tested .
this value corresponds to a concentration of  <dig> pm in our experimental setting, hence in good agreement with the previously reported range.
altogether these observations suggest that the quantity θ <dig> could be cell-type, mrna, and sirna indipendent.
it is estimated that the concentration of active risc in a cell is about  <dig> -  <dig> nm  <cit> .
taking into account that the volume of a mammalian cell is in the range 10-13l - 10-12l, then we can estimate that the number of active risc in a cell is in the range  <dig> -  <dig>  the above observations suggest that saturation begins when the number of sirna oligomers in a cell becomes comparable to the number of risc molecules.
we observed that the parameters of model  <dig> estimated when fitting protein levels  are very close to the ones estimated when fitting mrna levels .
namely, the optimized values of d <dig> and h <dig> are very similar for both experimental data.
this is important since these are two independent biological experiments.
this proves the mathematical robustness of model  <dig>  the only parameter changing between the two sets of experiment is θ <dig>  which represents the concentration of sirna oligomers needed to achieve half of the maximal degradation rate .
this is reflected in figure 2and figure  <dig>  where it is clear that saturation is achieved at about  <dig> pmol for the mrna data  and at about  <dig> pmol for the protein data .
this difference may be due to biological variability, or to the simplified model of protein translation dynamics we used .
we also conformed that model  <dig> is cell-line-independent, mrna-independent, and sirna-independent, since it can accurately describe the rna interference process on a different cell-line  expressing a different mrna , silenced by a different sirna oligomer.
interestingly, the difference in model  <dig> parameters, when testing a different mrna-sirna pair , shows that only d <dig>  and h <dig>  change significantly, suggesting that these two parameters can be used to describe changes in sirna-mrna silencing specific strength, whereas θ <dig> may be kept constant.
recently it has been proposed that sirna and microrna efficacy, defined as the percentage decrease in the target mrna level due to the silencing reaction, could be limited due to mrna abundance  <cit> or to mrna degradation rate  <cit> .
model  <dig> predicts that the percentage decrease in target mrna level  simply dividing by km/dm) is indeed sensitive to dm , with a higher degradation rate corresponding to a weaker effect of the silencing reaction, and vice-versa.
this result is in line with the experimental observation described in larsson et al  <cit> .
in addition, according to model  <dig>  the transcription rate km of the target mrna does not have any influence on the silencing reaction efficacy.
the target mrna abundance, in absence of the silencing reaction, is simply obtained from eq.
<dig> as km/dm.
smaller dm will correspond to a higher mrna abundance  therefore a correlation between mrna abundance and sensitivity to mrna can be found  <cit> , but this is only an indirect effect mediated by the degradation rate, at least according to our model.
our conclusion is that sirna-mediated degradation in mammalian cells can always be best represented as an enzymatic reaction described by an hill function, whose parameters have to be tuned to the specific sirna-mrna pair.
the models discussed so far consider the average behavior of a population of cells.
in the case of singe-cell experiments, these models might not be efficient enough due to their deterministic nature and will not be able to capture any stochastic effects.
since rna has a plethora of functional properties and plays many of roles in regulating gene expression, it has been used in a number of different studies as a tool for elucidating gene functions.
in fact with rnai it is possible to selectively knock-down any gene and even modulate its dosage  <cit> .
rna has also been used in the design of therapeutic molecules as well as metabolic reprogramming  <cit> .
the potential uses of this versatile molecule are still very much under study, but their effectiveness depend on many variables such as, the concentration of the silencing reagent, the transfection techniques, the cell type used and the target type selection.
in the present study we biologically validated for the first time a mathematical model  that has a simple mathematical form, amenable to analytical investigations and a small set of parameters with an intuitive physical meaning that can be used both by the computational and the experimental community interested in the analysis and application of rna interference.

 rna interference is a well characterized regulatory mechanism in eukaryotes as well as a powerful tool for understanding gene function, thanks to the discovery that synthetic small interfering rna oligomers can efficiently induce rnai in mammalian cells. rnai has also been used extensively as a novel "biological part" to design synthetic biological circuits in synthetic biology. artificial gene silencing has the potential to become a major genetic-based therapeutic tool for viral infections, cancer or inherited genetic disorders. despite its widespread experimental application, the best way to quantitatively model rna interference is still under debate. in systems and synthetic biology, mathematical models are essential to carry out in silico investigations of biological pathways, or novel synthetic circuits. the aim of this work is to find the most appropriate quantitative mathematical model that can correctly describe the rnai phenomenon in mammalian cells, for varying concentrations of the sirna oligomers. a schematic representation of the rna interference mechanism is illustrated in figure in step the presence of double stranded rna elicits a response in the cell mediated by the dicer enzyme, which binds and cleaves the dsrna into fragments of 21- base pairs, called small interfering rna. in step sirnas are loaded onto a multiprotein complex called rna induced silencing complex and then separated into single strands of which one is discarded and degraded, while the guide strand remains within risc and serves as a template in the silencing reaction. in step the guide strand assembles into a functional sirna-risc complex, which contains the sirna bound to the ago protein. target mrnas are then recognized by watson-crick base pairing and bound by the sirna-risc complex. finally, in step mrna degradation is induced, target mrna is dissociated from the sirna, and the sirna-risc complex is released to process further mrna targets. here, we will focus on quantitatively modeling step to step of the rna interference process using, as an experimental tool, synthetic sirna oligomers.  in order to model the effects of rna interference on mrna expression levels at different concentration of sirna oligomers, we carried out in-vivo experiments of rna interference on two mammalian cell-lines stably expressing the egfp protein or the tta protein, respectively. in the first set of experiments, human embryonic kidney cells stably expressing egfp, were transfected with varying quantities of synthetic sirna oligomers directed against the egfp mrna in the range of to pmol. figure shows the ratio of egfp mrna levels between treated cells and negative control cells measured hours after transfection. error-bars represent the standard-error from three biological replicates for each point. similarly, in the second set of experiments, egfp protein levels were measured via facs analysis hours after transfection of hek293-egfp cells with synthetic sirna oligomers directed against the egfp mrna, or with non-specific sirna oligomers as negative controls. the ratio between egfp protein levels in sirna-treated cells versus negative controls, for the same sirna quantities as in the set i experiments, is reported in figure in the third set of experiments, we tested chinese hamster ovary cells stably expressing the tetracycline-regulated transactivator tta at a low level, using the same protocol used previously for hek cells. cells were transfected with varying quantities of synthetic sirna oligomers against the tta mrna in the range of to pmol. since the tta is not fluorescent, no facs measurement were performed in this cell-line. in additional file figure 1we shows the ratio of tta mrna levels between cells treated with the silencing tta oligomer and cells treated with the negative control. rnai modeling we were interested in formulating a model that can quantitatively describe the effects of varying quantities of sirna oligomers onto the degradation of the target mrna species, and of its corresponding protein product. a general dynamical model describing transcription of the mrna species, its sirna-mediated degradation, and translation of its protein products, can be described by a system of ordinary differential equations. let xm, xp and xs be the mrna, protein and sirna concentrations, respectively. the evolution of their time-dependent concentrations can be described by the following odes::  dxmdt=km−dmxm−δ:  dxpdt=ktxm−dpxp, the parameter km, represents the transcription rate from the promoter that transcribes the mrna species targeted by the sirna oligomer; dm represents the basal degradation rate of the mrna species. rnai can be considered as a mechanism that enhances the degradation of the targeted mrna, therefore the function δ is an extra degradation term that represents the rate at which mrnas are degraded due to rnai. this function, δ, depends on both the mrna and sirna levels, xm and xs respectively. the parameter kt is the protein translation rate, whereas dp represents the basal protein degradation rate. at least four different models have been proposed in the literature, for the rna interference mechanism. each of these models is based on the general approach described by equations but each has a different functional form for δ. table lists all the different models studied here with a description of their corresponding parameters. c3: cleavage and dissociation rate of mrna-sirna* θ4: michaelis-menten like constant model 1: the stoichiometric model one possible way of modeling the effects of rnai on the mrna degradation is to consider a stoichiometric reaction between the sirnas and mrnas. let sirna* denote the concentration of the sirna-risc complex, namely the fraction of the sirnas that are loaded into risc complexes. then the stoichiometric reaction between the sirna* and the mrna can be described as follows: +←→∅. namely, the sirna-loaded risc binds to the complementary mrna and then both are degraded. according to this model, following the law of mass action, we predict that the sirna mediated degradation will be proportional to the product of the concentration of sirna oligomers and the targeted mrna species: δ=k1xmxs where parameter k represents the proportionality constant. in this modeling approach, the sirna-risc complex is assumed not to be recycled, but the risc needs to be reloaded before it can degrade another mrna molecule. indeed, this model was suggested for rna interference in prokaryotes. rnai in prokaryotes has important differences with rnai in eukaryotes. one of these is that the interaction between the sirna and its mrna target is non-catalytic in nature. this is not the case in mammalian cells. once the mrna is cleaved, the sirna-risc complex is dissociated from it, in an atp-dependent manner and it is free to process further targets. additionally, in some organisms, such as c.elegans, the primary dsrna trigger induces synthesis of secondary sirnas through the action of rdrp enzymes, strengthening and perpetuating the silencing response. model 2: stoichiometric model with co-operativity this model is a straightforward extension of model which additionally takes into account the presence of multiple sites on the targeted mrna where the sirna-loaded risc can bind. model can be easily extended to include co-operativity: ++...+︸h2↔mrna−h2sirna*→∅ as before, the rate of rnai-driven degradation, can be easily obtained applying the law of mass action: δ=k2xsh2xm where k is the proportionality constant and h is the number of sirna binding sites on the targeted mrna species. this model was suggested in for modeling rnai by mirnas, since it is experimentally proven that multiple sites for the same mirna can boost target mrna repression. this model however suffers from the same limitations as model model 3: enzymatic model a detailed model of rnai specific for mammalian cells was proposed by malphettes et al. this model accounts for the catalytic nature of rnai in mammalian cells, thus modeling step step and step of figure including the recycling of the risc complex. specifically, it assumes that once the cleavage and degradation of the targeted mrna is completed, the sirna-risc complex dissociates from its mrna target, and it is free to degrade further mrna molecules. additionally, this model also considers that the targeted mrna may have a multiple number, h of sirna binding sites. the sirna-risc complex can bind to any site on the mrna to form an intermediate mrna-sirna* complex, which can either accommodate further sirna-risc complexes on any other free binding sites, or cleave the target mrna and dissociate from the cleavage products. the reaction of the complex formation of the target mrna with the sirna-risc complex is described as follows: +↔. the model considers the following reaction between an intermediate mrna-sirna* i- complex with another sirna-risc complex: +↔. the generic cleavage and degradation reaction of the mrna by any interacting sirna-risc complex  +  → cleaved mrna+i. in, the authors show that the rate of rnai-driven mrna degradation, δ, is given by: δ=c3k3h3c3+k3h3xmxsxm where k is the rate of mrna-sirna* complex formation for a single sirna target site and c is the cleavage and dissociation rate of mrna-sirna* complex. this functional form, for a constant xs, is a classic michaelis-mentenu enzymatic reaction, as can be observed by simply defining vm = c3xs and km=c3k3h and rewriting eq. as: δ=vmxmkm+xm this is in perfect agreement with the experimental finding on the enzymatic activity of both non-mammalian and mammalian risc on mrna degradation, as reported in, where the product of the reaction was measured at varying the concentration of the substrate for a constant amount of the risc enzyme. model basically assumes that only the maximal rate vm of the michaelis-menten will be affected when the concentration of the active risc changes. in the authors show experimentally that changing the concentration of risc indeed changes the value of vm but also of km. this is not captured by this model. note also that when xm ≪ km, i.e. low amounts of mrna compared to the michaelis-menten constant km, the model can be simplified to: δ=vmkmxm=k3h3xmxs. this is equivalent to model on the other hand, when the mrna concentration is very high, the model can be simplified to: δ = vm = c3xs. namely, the rate of rnai-driven mrna degradation will depend only on the amount of active enzyme, i.e. the sirna concentration. the higher the sirna concentration, the higher the degradation rate that can be achieved, without any saturation effect. model 4: phenomenological model in the authors proposed a standard hill-kinetic model to describe the post-transcriptional effects of micrornas on the gene expression. mirnas are processed by the cell to produce dsrnas which then follow the typical rnai pathway, schematically described in figure. by considering a hill-type enzymatic model with an hill coefficient h ≥ the model can be extended to account, either, for multiple binding sites of the sirna on the same target mrna, or, for the cooperativity of protein complexes involved in rnai. differently from the other three model presented above, this model is a phenomenological model not derived from specific biochemical reactions: δ=d4xsh4θ4h4+xsh4xm this model, despite being phenomenological has interesting properties. the kinetic parameters d and θ depend on the efficiency of sirna binding to its sites on the target mrna: d represents the maximal degradation rate of the mrna due to rna interference; θ the concentration of sirna oligomers needed to achieve half of the maximal degradation rate. the above equation implies that for xs ≪ θ the increase in the rnai mediated degradation is linear with xsh, while it saturates at higher levels of xs, reaching the maximal degradation rate d differently from model parameter identification the four models were fitted to the three mrna and protein experimental datasets, by searching for the parameter values for which the model-generated data best fitted the experimental data, according to a squared error measure. the results of the fitting procedure for each of the models, together with the optimized values of their parameters, are given in table for egfp and in additional file table a for tta. the relative value of the error and the relative value of the prediction error for each model are given, together with the corresponding optimized values of its parameters. the unit of measurements are reported for the dimensional parameters.. the fitting results for the mrna levels are shown in figure model gives a significant smaller error than the other three models, in fact, compared to models and the error of model is two orders of magnitude smaller. model gives a better fitting than models and however it is worse than model the optimized value for the parameter d in model is d = min- indicating that the strength of sirna mediated mrna degradation is comparable to the strength of basal mrna degradation. note also that the parameters found for model include a coefficient h =, hence less than unity. since h describes the number of sirna binding sites on the targeted mrna, it should be greater than, or equal to, in order to have a clear biological interpretation. however, if we constrain this parameter to be greater or equal to one, then the model optimizes at the value of h = which makes model identical to model we have observed that in all numerical simulations, models and are almost indistinguishable. the large optimized value of parameter c of model, the low value of parameter k3h = × 10- suggest c ≫ k3h3xm, and hence the function δ of model can be approximated by k3h3xsxm, which is nothing else than model. therefore, in this parameter space, model is almost identical to model similar results were obtained when searching for the parameters which best fitted the measured protein levels, shown in figure model is again the one with the smallest error. model and are again unable to capture efficiently the experimental data. model is better than models and but it is still worse than model when we fitted the models to the third experimental dataset, which was performed on a different cell-line, with both a different target mrna and a different sirna oligomer, model still performed better than the others with the smallest error, although model was a close match. model and were again behaving very similarly and had the largest error. it should be noted that as it happened with the previous fitting results, also in this case model has a hill coefficient smaller than unity. finally we observed, as shown in additional file 1figure a that the experimental error for set iii experiment was larger when compared to the egfp experiment. this was caused by the relative low expression of the tta in the cho cell lines, when compared to the egfp expression in hek cells, which made measurements more noisy. assessing the model predictive ability the models described above have the same number of unknown parameters to be learned, but for model which has one extra parameter. to be sure that the improved performance of model in describing the experimental data was not due to overfitting, we computed for each model and for each experimental dataset, the prediction error, which allows to assess the generalisation performance of the models. we followed a "leave-one-out" cross validation strategy, where for each model and for each dataset, the parameter identification procedure was repeated each time removing one of the experimental points and then predicting the missing point with the identified parameters. we thus could estimate a prediction error for each model in each experimental dataset. this value is reported in table 2and additional file 1table a model is again the one with the smallest prediction error, once again confirming this model superior performance in describing the data.  our findings show that the simple hill function described by model is sufficient to quantitatively describe the effect of rna interference, at the mrna and protein level, in mammalian cells in vitro, for varying concentration of sirna oligomers. one significant feature of model is that it can predict the saturation effect of the rnai process that we observed experimentally. we considered the possibility that this saturation could be in fact due to the inability of the cell to uptake high concentration of sirna oligomers, however recent experiments, prove that uptake of sirna oligomers in cells is linear with the concentration of sirna oligomers transfected, at least in the concentration range we used. additionally, khan et al in, observed upregulation of mrna targets of endogenous micro-rna when transfecting sirna oligomers in mammalian cells. in order to explain this effect, they suggested a saturation of the risc complex. it has been demonstrated in that the enzymatic activity of risc can be efficiently modeled in-vitro as a classic michaelis-menten reaction, where the target mrna is the substrate, the sirna-loaded risc is the active enzyme, and the product is the degraded mrna. this is one feature that model does not capture; namely for a fixed amount of sirna-loaded risc, eq. should approximate the michaelis-menten in eq., instead model becomes simply proportional to xm, as model and nevertheless, model approximates very well the experimental data. we believe this happens because enzymatic reactions have a typical km much greater that the physiological concentration of their substrate, and the same happens for the risc complex. in this condition, the michaelis-menten equation becomes vmxmkm+xm≊vmkmxm and therefore it is linear in xm, as predicted by model when the sirna concentration varies, model predicts that the parameters vmkm will change as a function of xs as described by eq.. the three parameters of model have a straightforward biological interpretation, and their values can be easily tuned to accommodate for different efficiencies of rnai. for example, the parameter d can be used to weigh the degradation due to the rnai compared to the endogeneous mrna degradation, and its strength, i.e. what is the maximal degradation rate that can be achieved. θ quantifies the sirna oligomers concentration needed to achieve half of the maximal degradation of the targeted mrna. the h coefficient can accommodate for multiple target sites on the same mrna, or for the cooperativity of the risc complex. clearly, the rnai process is very complex and no one-to-one relationship can be found between parameters of model and rnai biological components. nevertheless, it has been shown in that between and sirna oligomers per cell are sufficient to reach half-maximal mrna target degradation. model predicts that half-maximal degradation is achieved for an amount of sirna oligomers equal to θ the value of this parameter when fitting mrna levels is θ ≈ pmol despite of the different cell-lines and mrna-sirna pairs tested. this value corresponds to a concentration of pm in our experimental setting, hence in good agreement with the previously reported range. altogether these observations suggest that the quantity θ could be cell-type, mrna, and sirna indipendent. it is estimated that the concentration of active risc in a cell is about - nm. taking into account that the volume of a mammalian cell is in the range 10-13l - 10-12l, then we can estimate that the number of active risc in a cell is in the range - the above observations suggest that saturation begins when the number of sirna oligomers in a cell becomes comparable to the number of risc molecules. we observed that the parameters of model estimated when fitting protein levels are very close to the ones estimated when fitting mrna levels. namely, the optimized values of d and h are very similar for both experimental data. this is important since these are two independent biological experiments. this proves the mathematical robustness of model the only parameter changing between the two sets of experiment is θ which represents the concentration of sirna oligomers needed to achieve half of the maximal degradation rate. this is reflected in figure 2and figure where it is clear that saturation is achieved at about pmol for the mrna data and at about pmol for the protein data. this difference may be due to biological variability, or to the simplified model of protein translation dynamics we used. we also conformed that model is cell-line-independent, mrna-independent, and sirna-independent, since it can accurately describe the rna interference process on a different cell-line expressing a different mrna, silenced by a different sirna oligomer. interestingly, the difference in model parameters, when testing a different mrna-sirna pair, shows that only d and h change significantly, suggesting that these two parameters can be used to describe changes in sirna-mrna silencing specific strength, whereas θ may be kept constant. recently it has been proposed that sirna and microrna efficacy, defined as the percentage decrease in the target mrna level due to the silencing reaction, could be limited due to mrna abundance or to mrna degradation rate. model predicts that the percentage decrease in target mrna level simply dividing by km/dm) is indeed sensitive to dm, with a higher degradation rate corresponding to a weaker effect of the silencing reaction, and vice-versa. this result is in line with the experimental observation described in larsson et al. in addition, according to model the transcription rate km of the target mrna does not have any influence on the silencing reaction efficacy. the target mrna abundance, in absence of the silencing reaction, is simply obtained from eq. as km/dm. smaller dm will correspond to a higher mrna abundance therefore a correlation between mrna abundance and sensitivity to mrna can be found, but this is only an indirect effect mediated by the degradation rate, at least according to our model. our conclusion is that sirna-mediated degradation in mammalian cells can always be best represented as an enzymatic reaction described by an hill function, whose parameters have to be tuned to the specific sirna-mrna pair. the models discussed so far consider the average behavior of a population of cells. in the case of singe-cell experiments, these models might not be efficient enough due to their deterministic nature and will not be able to capture any stochastic effects. since rna has a plethora of functional properties and plays many of roles in regulating gene expression, it has been used in a number of different studies as a tool for elucidating gene functions. in fact with rnai it is possible to selectively knock-down any gene and even modulate its dosage. rna has also been used in the design of therapeutic molecules as well as metabolic reprogramming. the potential uses of this versatile molecule are still very much under study, but their effectiveness depend on many variables such as, the concentration of the silencing reagent, the transfection techniques, the cell type used and the target type selection. in the present study we biologically validated for the first time a mathematical model that has a simple mathematical form, amenable to analytical investigations and a small set of parameters with an intuitive physical meaning that can be used both by the computational and the experimental community interested in the analysis and application of rna interference. 
the findings here presented will be a useful instrument for better understanding rnai biology and as modelling tool in systems and synthetic biology. our model has a simple mathematical form, amenable to analytical investigations and a small set of parameters with an intuitive physical meaning, that makes it a unique and reliable mathematical tool. rna interference is a regulatory cellular process that controls post-transcriptional gene silencing. we tested and validated four different mathematical models of rna interference by quantitatively fitting models' parameters to best capture the in vitro experimental data. we show that a simple hill kinetic model is the most efficient way to model rna interference. we performed in vitro experiments in human and hamster cell lines constitutively expressing respectively egfp protein or tta protein, measuring both mrna levels, by quantitative real-time pcr, and protein levels, by facs analysis, for a large range of concentrations of sirna oligomers. sirnas are then loaded onto the rna-induced silencing multiprotein complex, which uses the sirna antisense strand to specifically recognize mrna species which exhibit a complementary sequence. here, we address this pressing open problem performing in vitro experiments of rnai in mammalian cells and testing and comparing different mathematical models fitting experimental data to in-silico generated data. once the sirna loaded-risc binds the target mrna, the mrna is cleaved and degraded, and the sirna loaded-risc can degrade additional mrna molecules. despite the widespread use of sirnas for gene silencing, and the importance of dosage for its efficiency and to avoid off target effects, none of the numerous mathematical models proposed in literature was validated to quantitatively capture the effects of rnai on the target mrna degradation for different concentrations of sirnas.
 rna interference is a regulatory cellular process that controls post-transcriptional gene silencing. during rnai double-stranded rna induces sequence-specific degradation of homologous mrna via the generation of smaller dsrna oligomers of length between 21-23nt. sirnas are then loaded onto the rna-induced silencing multiprotein complex, which uses the sirna antisense strand to specifically recognize mrna species which exhibit a complementary sequence. once the sirna loaded-risc binds the target mrna, the mrna is cleaved and degraded, and the sirna loaded-risc can degrade additional mrna molecules. despite the widespread use of sirnas for gene silencing, and the importance of dosage for its efficiency and to avoid off target effects, none of the numerous mathematical models proposed in literature was validated to quantitatively capture the effects of rnai on the target mrna degradation for different concentrations of sirnas. here, we address this pressing open problem performing in vitro experiments of rnai in mammalian cells and testing and comparing different mathematical models fitting experimental data to in-silico generated data. we performed in vitro experiments in human and hamster cell lines constitutively expressing respectively egfp protein or tta protein, measuring both mrna levels, by quantitative real-time pcr, and protein levels, by facs analysis, for a large range of concentrations of sirna oligomers.  we tested and validated four different mathematical models of rna interference by quantitatively fitting models' parameters to best capture the in vitro experimental data. we show that a simple hill kinetic model is the most efficient way to model rna interference. our experimental and modeling findings clearly show that the rnai-mediated degradation of mrna is subject to saturation effects.  our model has a simple mathematical form, amenable to analytical investigations and a small set of parameters with an intuitive physical meaning, that makes it a unique and reliable mathematical tool. the findings here presented will be a useful instrument for better understanding rnai biology and as modelling tool in systems and synthetic biology. 