Estimate hidden dynamic profiles of siRNA effect on apoptosis
 Takanori Ueda^{1},
 Daisuke Tominaga^{2}Email author,
 Noriko Araki^{1} and
 Tomohiro Yoshikawa^{1}
DOI: 10.1186/147121051497
© Ueda et al.; licensee BioMed Central Ltd. 2013
Received: 20 September 2012
Accepted: 4 March 2013
Published: 15 March 2013
Abstract
Background
For the representation of RNA interference (RNAi) dynamics, several mathematical models based on systems of ordinary differential equations (ODEs) have been proposed. These models consist of equations for each molecule that are involved in RNAi phenomena. Therefore, many realvalue parameters must be optimized to identify the models. They also have many ‘hidden variables’, which cannot be observed directly through experimentation. Calculation of the values of the hidden variables is generally very difficult, if not impossible in some special cases. Identification of the ODE models is also quite difficult.
Results
We show that the simplified logistic LotkaVolterra model, a wellestablished ODE model for biological and biochemical phenomena, can represent RNAi dynamics as a predatorprey system. Although a hidden variable exists in the model, its values can be determined and made visible as dynamic profiles of RNAdecomposing effects of siRNAs. Correlation analysis shows that the model parameters correlate highly with the total effect of the siRNA.
Conclusions
The results suggest that analyses using our model are useful to estimate dynamic profiles of siRNA effects on apoptosis and to score siRNA by its effects on apoptosis, namely ‘phenotypic scoring’.
Keywords
siRNA RNA interference Preypredator model Ordinary differential equation Parameter estimation Genetic algorithmBackground
In a body of living eukaryotic cells, RNA interference (RNAi) is a phenomenon caused by small interfering RNA (siRNA, a doublestrand RNA consisting of 2123 base pairs) in which there is decomposition of singlestrand RNA that has a sequence that is complementary with the siRNA [1, 2]. Artificially introducing siRNA into cells is widely used in experiments to suppress gene expression or to interrupt gene regulatory networks. Therefore, RNAi is a very useful and popular technique for approaching the molecular mechanisms of life.
Eukaryotic cells can incorporate external RNA. The incorporated doublestrand RNA is then separated into two single strands that combine with protein molecules to form an RNAinduced silencing complex (RISC). The RISC combines with an RNA molecule in a cell that has a complementary sequence with the siRNA and decomposes it into fragments. When this target RNA is messenger RNA (mRNA), then either expression of the protein that is translated from the mRNA is suppressed or malfunctioning protein molecules will be produced.
The molecular mechanisms of RNAi are complex and are not completely revealed. In addition, quantitative measurements of the amount of incorporated siRNA and its effectiveness or strength of RNA decomposition are very difficult, especially in a time series. The authors believe that quantitative mathematical models can be applied to address this problem.
Several mathematical models have been proposed to represent the mechanisms and dynamics of RNAi [36]. These models are systems of linear ordinary differential equations (ODEs). Each equation in the system represents kinetics of a chemical reaction that constitutes the RNAi mechanism.
Scales of previously proposed models
Equations  Parameters  

Fundamental model  4  12 
Considering cell cycle  12  27 
Considering viral effect  17  14 
Selftargeting siRNA  4  8 
Predatorprey model (our model)  2  4 
From the perspective of numerical optimization, scales of these models, i.e. quantities of parameters, are not sufficiently small to identify models for actual amounts of available observed data. In addition, these models include hidden variables that are extremely difficult to observe. Hidden variables increase the degrees of freedom of the models, so quite large amounts of observation data are necessary to identify these models (to make the degrees of freedom of the model zero). Because experimental observation costs money, in many cases the time series data that are necessary to identify the models are insufficient.
The molecular mechanisms of RNAi are the main target of molecular biology today, and new knowledge about that has been growing. Previously established models did not consider such newly found mechanisms. Instead, abstract mathematical models need not change the model formulae.
A mathematical model for a system that has unknown mechanisms, such as for the RNAi phenomenon, can be expected to fit the observed data with fewer degrees of freedom (fewer number of parameters) in the sense of Ockham’s razor rather than those that fit better comparably but with more parameters.
RNAi is a phenomenon by which siRNA degrades target RNA molecules. This relation can be considered like that of prey and predator in a natural food chain even if the siRNA is artificially applied.
This model is called the logistic LotkaVolterra (LLV) model [8].
where x is the number of living cells, and y stands for the strength of siRNA that causes apoptosis (cell death). We will call this ODE system the Ueda model. The variable y is not observable, although x is observable using several experimental techniques. The hidden variable y is readily calculated using numerical integration when the observation data of x are given and the values of all five parameters (a, b, c, d, and e in Equation 3) and the initial value of the hidden variable (y_{0}) are given. Therefore, the trialanderror procedure, or heuristic optimization techniques, can find these parameter values by fitting x to the given time series data.
Here we show that the bestfitted models can clarify the dynamic profile of the invisible hidden variable (y in Equation 3) that implies the strength of the siRNA. We also show that its parameter values are distributed according to the strength of siRNA, and a parameter in the model that highly correlates with the total effect of siRNA. Model parameters are determined based on experimentally observed data. Our method and model proposed in this paper is for evaluation of siRNA, and is not for prediction. The model parameters quantitatively represent the strength of an siRNA, and can be interpreted as a kind of score value. This is the first approach to score siRNA by its effectiveness.
Methods
Negative control experiment
Introduced siRNA to cause apoptosis on HeLa cells
siRNA name  Product name 

ACD  QIAGEN:1027299 (AllStars Hs Cell Death Control siRNA) 
KIF  QIAGEN:SI03019793 (Hs_KIF11_8) 
PLK  QIAGEN:SI02223837 (Hs_PLK1_6) 
VHP  QIAGEN:1027273 (Very High Potency Hs_CDC2 siRNA) 
HP  QIAGEN:1027274 (High Potency Hs_CDC2 siRNA) 
MP  QIAGEN:1027275 (Moderate Potency Hs_CDC2 siRNA) 
Neg  QIAGEN:1027310 (Negative Control) 
The times of observations were at 2.4, 5.4, 8.0, 20.7, 29.0, 44.9, 52.7, 68.7, 76.8, and 92.7 hr after the start of cultivation (ten time points in total). The siRNA introduction started simultaneously with cultivation. The measured values were the intensity of fluorescence in the ATP assay, not the actual numbers of cells. The mean value of four repeated experiments was used as the data for each point of the sampling time.
The parameter estimation in this negative control experiment was optimization of three real parameters (a, b, and x_{0}) from ten realnumber data points at the designated times. The real coded genetic algorithm was used to find the optimal values of the parameters (explained below).
Estimated model parameters
siRNA  a  b  x _{0}  c  d  y _{0} 

Neg        
ACD  2.694e7  8.087e2  2.947e4  
KIF  1.784e7  1.506e1  8.382e3  
PLK  3.949e2  1.445e8  6.512e+5  1.337e7  9.045e2  2.837e3 
VHP  2.498e8  3.091e7  3.458e4  
HP  2.243e8  1.118e7  6.760e4  
MP  7.258e9  6.509e8  2.305e3 
Positive control experiment
Next, we introduced siRNA into HeLa cells to cause apoptosis and observed the changes of the cell population in time. Experimental conditions were the same as those of the negative control observation except that effective siRNA was induced into the cells. The observed cell population changes and fitted curves are shown in Figure 2B through 2G.
The values of the parameters a and b and x_{0} are the same as those estimated at the negative control experiment. Therefore, the number of estimated values for the model of the positive control experiments is three (c, d, and y_{0}), and the numerical difficulties of this optimization are the same as in the case of the negative control experiment. The values of the estimated parameters, c, d, and y_{0}, are shown in Table 3 for each siRNA.
Numerical optimization
Parameter optimization of this case study was accomplished by finding the parameter values that made the curve of the model best fit the given data. The fitness of each model was calculated as the reciprocal of the total sum of the squared relative error between the given data and the model curve, which was calculated using numerical integration. We used the RungeKutta method (4th order) as the method of numerical integration.
The optimization task was undertaken to maximize the fitness of the model. We used the real coded genetic algorithm [11], which introduces UNDX [12] as the crossover operation and MGG [13] as the selection operation. The estimated optimal parameter values are shown in Table 3.
Results and discussion
Estimation of models
The effects of the siRNAs are not defined exactly in common, however, and they are considered ideally as reflecting the rate of the degradation of the target RNA molecules. Assuming that it is reflected in the cell population in our case study, the area under the siRNA strength curve, which is calculated using numerical integration, can be interpreted as the actual total of the siRNA effect. The differences between the curves for cell populations with and without siRNA (positive and negative controls) are also regarded as the siRNA effect.
Correlation between measurements calculated from the estimated dynamics of siRNA strength
Diff.  siRNA  Height  c  d  PC1  

Diff.  1.0  0.99922  0.97268  0.97605  0.68269  0.69267 
siRNA  1.0  0.97307  0.96764  0.66965  0.67984  
height  1.0  0.92129  0.50365  0.51537  
c  1.0  0.79264  0.80059  
d  1.0  0.99989  
PC1  1.0 
Discussion
The parameter estimation problem of the Ueda model in this case study is similar to the least squares method that minimizes the total sum of absolute errors, whereas our optimization minimizes relative errors. The amount of data (10 sampling points) is considered sufficient to estimate the number of parameters (2 parameters and 1 initial value, total is 3).
The profile of the cost function to be minimized (total sum of the relative error between the given data and the time series data calculated using numerical integration of the model) has a highly nonlinear shape because the model cannot be solved numerically in certain regions of the parameter space. In such regions, numerical computational errors such as overflow, underflow, and division by zero prevent calculation of the cost function. Established analytical optimization methods such as the steepest descent method, the conjugate gradient method, arrangements of the NewtonRaphson method, etc. require smooth search regions that do not contain such regions. However, determining those regions before optimization is extremely difficult in most cases. Reportedly, heuristic search methods are effective for optimization of systems of differential equations [14].
Correlation analysis shows that the parameter c is highly correlated with total siRNA effects. This fact suggests that the value of c can be used to score the siRNAs. Larger values of c signify a stronger effect. This score reflects the actual effect without considering the nucleic acid sequence of siRNA that relates its characters. Therefore, this score can be called a ‘phenotypic score’. Further analysis of additional experimental data is needed to prove the reliability of this score.
According to the general interpretation of the LLV model, the parameter b in the Ueda model is understood as the carrying capacity of the environment in which the cells are living. Here, c represents the positive effect of the cell population to the strength of the siRNA. Because intake of external siRNA is more active on more viable cells, c might be interpreted as the siRNAinducing capacity of the cells, or more generally, the cell viability under the cultivating condition.
Each term of the Ueda model can be interpreted as follows: ax represents selfreproduction rate of cells, −b x denotes the carrying capacity, −e x y is degradation rate of cells by Apoptosis caused by siRNA, cxy is siRNA incorporating rate to cell bodies that reflects cell activity or viability, and −d y is the decreasing rate of siRNA effect by decreasing cell viability.
One difficulty of the Ueda model currently is the fitting difficulty. Results demonstrated that multipoint heuristic searches are effective for systems of ODEs [11, 14]. However, this optimization method demands many computational resources. To apply a fast analytic search, the cost function should be modified.
Long simulation of the Ueda model also should be considered. After the end point of the observation the model shows an oscillating profile. Continuous experiments could not match the oscillating profile of the simulation because siRNA molecules in the medium would be exhausted.
The Ueda model is not based on molecular mechanisms. The Ssystem [15] is the same on this point. Applying this canonical form ODE model to RNAi dynamics is expected to be interesting for its capability for application to dynamical system analyses.
Conclusions
Cell population changes by apoptosis that results from introduction of siRNA were observed as quantitative time course datasets. The Ueda model, the simplified logistic LotkaVolterra model, can fit these datasets. The optimal models represent dynamic profiles of RNA decomposing effects of siRNA in apoptosis that cannot be observed directly through experimentation. Parameter estimation using the Ueda model can be done using the real coded genetic algorithm. One estimated parameter correlates highly with the estimated siRNA strength. We think that this parameter might represent the effectiveness of siRNA.
Abbreviations
 LV model:

LotkaVolterra model
 LLV model:

Logistic LotkaVolterra model
 mRNA:

messenger RNA
 ODE:

Ordinary Differential Equation
 PCA:

Principal Component Analysis
 RISC:

RNAInduced Silencing Complex
 RNAi:

RNA interference
 siRNA:

short interfering RNA.
Declarations
Acknowledgements
All works of the manuscript was done by the authors at their affiliations.
Authors’ Affiliations
References
 Fire A, Xu S, Montgomery M, Kostas S, Driver S, Mello C: Potent and specific genetic interference by doublestranded RNA in Caenorhabditis elegans. Nature 1998, 391: 744745. 10.1038/35750View Article
 Elbashir S, Harborth J, Lendeckel W, Yalcin A, Weber K, Tuschl T: Duplexes of 21nucleotide RNAs mediate RNA interference in cultured mammalian cells. Nature 2001, 411: 494498. 10.1038/35078107View ArticlePubMed
 Bergstrom C, McKittrick E, Antia R: Mathematical models of RNA silencing: Unidirectional amplification limits accidental selfdirected reactions. Proc Natl Acad Sci USA 2003, 100: 1151111516. 10.1073/pnas.1931639100PubMed CentralView ArticlePubMed
 Bartlett D, Davis M: Insights into the kinetics of siRNAmediated gene silencing from livecell and liveanimal bioluminescent imaging. Nucleic Acids Res 2006, 34: 322333. 10.1093/nar/gkj439PubMed CentralView ArticlePubMed
 Groenenboom M, Hogeweg P: The dynamics and efficacy of antiviral RNA silencing: A model study. BMC Syst Biol 2008, 2: 28. 10.1186/17520509228PubMed CentralView ArticlePubMed
 Marshall W: Modeling recursive RNA interference. PLoS Comput Biol 2008, 4: e1000183. 10.1371/journal.pcbi.1000183PubMed CentralView ArticlePubMed
 Berryman A: The origins and evolution of predatorprey theory. Ecology 1992, 73: 15301535. 10.2307/1940005View Article
 Gause G: The Struggle for Existence. NY, USA: Hafner Publishing; 1964.
 Yoshikawa T, Uchimura E, Kishi M, Funeriu D, Miyake M, Miyake J: Transfection microarray of human mesenchymal stem cells and onchip siRNA gene knockdown. J Control Release 2004, 96: 227—232.PubMed
 Crouch S, Kozlowski R, Slater K, Fletcher J: The use of ATP bioluminescence as a measure of cell proliferation and cytotoxicity. J Immunol Methods 1993, 160: 8188. 10.1016/00221759(93)90011UView ArticlePubMed
 Ueda T, Koga N, Ono I, Okamoto M: Application of numerical optimization technique based on realcoded genetic algorithm to inverse problem in biochemical systems. GECCO 2002: Proc Genet Evol Comput Conf 2002, 701701.
 Ono I, Kobayashi S: A realcoded genetic algorithm for function optimization using unimodal normal distribution crossover. Proc Seventh Intl Conf Genet Algo 1997, 246253.
 Satoh H, Yamamura I, Kobayashi S: Minimal generation gap model for gas considering both exploration and exploitation. Proc fourth Intl Conf Soft Comp (Iizuka ’96) 1996, 494497.
 Kikuchi S, Tominaga D, Arita M, Takahashi K, Tomita M: Dynamic modeling of genetic networks using genetic algorithm and Ssystem. Bioinformatics 2003, 19: 643—650.View ArticlePubMed
 Voit E: Canonical Nonlinear Modeling: SSystem Approach to Understanding Complexity. NY, USA: Van Nostrand Reinhold; 1991.
Copyright
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.