Kinetic Modeling to Accelerate the Development of Nucleic Acid Formulations.

A critical hurdle in the clinical translation of nucleic acid drugs is the inefficiency in testing formulations for therapeutic potential. Specifically, the ability to quantitatively predict gene expression is lacking when transitioning between cell culture and animal studies. We address this challenge by developing a mathematical framework that can reliably predict short-interfering RNA (siRNA)-mediated gene silencing with as few as one experimental data point as an input, evaluate the efficacies of existing formulations in an expeditious manner, and ultimately guide the design of nanocarriers with optimized performances. The model herein consisted of only essential rate-limiting steps and parameters with easily characterizable values of the RNA interference process, enabling the easy identification of which parameters play dominant roles in determining the potencies of siRNA formulations. Predictions from our framework were in close agreement with in vitro and in vivo experimental results across a retrospective analysis using multiple published data sets. Notably, our findings suggested that siRNA dilution was the primary determinant of gene-silencing kinetics. Our framework shed light on the fact that this dilution rate is governed by different parameters, i.e., cell dilution (in vitro) versus clearance from target tissue (in vivo), highlighting a key reason why in vitro experiments do not always predict in vivo outcomes. Moreover, although our current effort focuses on siRNA, we anticipate that the framework can be modified and applied to other nucleic acids, such as mRNA, that rely on similar biological processes.

Completion of the human genome project and advances in drug delivery technologies have led to skyrocketing interest in nucleic acid therapies over the past decade.[1−3] For example, with an estimated market size of $14.6 billion by 2030, gene therapeutics development is being pursued by numerous biopharmaceutical companies, and there are thousands of ongoing gene therapy clinical trials.[4,5] However, even though nucleic acids are now recognized as the third major class of drugs behind small molecules and antibodies,[6] the translation to practice of nucleic acid based therapies remains limited. Only 20 nucleic acid drugs have been approved by the Food and Drug Administration for commercial use as of May 2021, including two mRNA vaccines for COVID-19.[7,8] One of the most prominent bottlenecks that hampers the development of nucleic acid therapeutics is the inability to formulate and test these drugs in an expedient manner.[9] The determination of appropriate dosing schedules necessary for therapeutically relevant outcomes is a critical yet tedious process requiring extensive preclinical experimentation and a battery of clinical studies.[10−12] To this end, we report on a simple mathematical model that can predict siRNA-mediated gene silencing built from knowledge of fundamental biological processes, i.e., the central dogma of molecular biology. We demonstrate the utility of this streamlined modeling approach for the optimization of dosing schedules needing only limited experimental data.

The above-mentioned mathematical framework is vital because, in many cases, therapeutics are tested in vivo only if the desired effects are first achieved in vitro.[13−15] However, the application of nucleic acid therapeutics still remains a resource-intensive process as a result of the lack of predictive capacity when transitioning from in vitro studies to more complex preclinical animal studies and ultimately clinical application.[16−18] Even when accounting for the fractional biodistribution of nucleic acid therapeutics into different organs/tissues, in vitro studies often fail to accurately predict in vivo potency.[18,19] This gap in actionable understanding can be attributed to a couple of factors. One, differences between key kinetic parameters that affect nucleic acid concentrations in vitro versus in vivo are not considered in most studies, although the associated concentration changes directly correlate to alterations in therapeutic efficacy. Notably, the clearance of nucleic acids at the target tissue in vivo results in dilution rates that are different from the analogous rates applicable to in vitro studies.[20] Two, the collection of comprehensive data sets from cell culture is largely impractical due to the need to test a wide, and sometimes unknown, range of conditions. For example, the time points selected for most experimental analyses are arbitrary, and these time points often do not capture the extent and duration of gene silencing in different application systems, e.g., different target genes, vectors, cell types.[21−23] Mathematical models can address these limitations by predicting gene expression levels over broad time scales in various systems and provide a necessary bridge to link crucial information regarding the potency of nucleic acid therapeutics.

Although a few kinetic models have been developed for the RNA interference (RNAi) process,[24−26] to the authors’ knowledge, none of these prior models have been adopted in a widespread fashion. This lack of comprehensive application is due in part to the need to specify numerous parameters that are neither readily available nor easy to characterize experimentally across different systems.[24−33] For example, previous kinetic models, while robust, required the input of up to 30 parameters, including cell membrane binding rates, endosomal escape rates, and the formation of various RNA-induced silencing complex (RISC) intermediates,[25,29,32] which may be challenging to characterize experimentally. As a result, existing mathematical models have not been widely adopted or used to conduct a retrospective analysis of experimental systems reported in literature. Such studies are essential to validate the reliability and applicability of the model frameworks to predict gene expression behavior spanning multiple genetic targets, cell types, and experimental conditions. Analyses across diverse data sets also can provide insights into the primary parameters governing nucleic acid efficacy.

Herein, we describe the development and testing of a simple kinetic model that leveraged existing mathematical frameworks to predict siRNA-mediated gene silencing.[25,34−38] Our model consisted of major rate limiting steps in the RNAi process, such as the dilution of siRNA, and it included seven independent kinetic parameters, six of which have known values or can be measured experimentally and one unknown parameter that can be fitted from experimental data. The mathematical model was employed to analyze published data obtained from a series of prior studies that investigated siRNA-mediated gene silencing,[21,25,38−44,100] and the level and duration of knockdown predicted by the model were in good agreement with both in vitro and in vivo experimental results. Notably, a comparison of model predictions generated by fitting one data point versus multiple data points were in close agreement and accurately predicted experimental data, confirming that our mathematical framework can be implemented even when experimental data were limited to a single data point. The alignment between the single-data-point fit and multiple-data-point fit held up across retrospective analyses of multiple data sets reporting a variety of experimental scenarios, further demonstrating the robust predictive capacity of our simple framework. Model predictions also suggested that when one correctly accounts for the dominant kinetic parameters that dictate changes in siRNA concentration in vitro versus in vivo, in vivo gene-silencing kinetics can be predicted from in vitro studies. As one example, a key factor highlighted in this study was that the dilution of siRNA in vitro depended solely on cell doubling rates, whereas the dilution in vivo was governed by the clearance rate of siRNA from the target tissue. Additionally, the ability to predict gene silencing levels at prolonged time scales made it possible to come to more informed conclusions with respect to the efficacy of siRNA treatments. For instance, the model could enable a user to draw more comprehensive conclusions regarding the potency of in vitro formulations, particularly when targeting proteins with longer half-lives (t1/2) that did not show significant silencing at earlier time points. We expect the model to provide guidance on the dosing regimens of existing siRNA formulations, and more importantly, to inform the design of next-generation siRNA vehicles to achieve specific kinetic profiles. The impact of, for example, the size and surface charge of nanocarriers on fractional accumulation in tissue can be estimated from experimental results,[45,46] and this parameter value can be updated within the mathematical framework to predict the silencing profiles of various siRNA formulations. Furthermore, while our current model focuses on siRNA-mediated gene silencing, we anticipate that the framework can be readily modified for use with other nucleic acid therapeutics, such as mRNA- and microRNA-based interventions, that rely on similar biological pathways.

Results and Discussion

For in vitro studies, the effectiveness of an siRNA treatment is determined by characterizing changes in mRNA and/or protein expression levels.[47−49] Time points ranging from less than 24 h to multiple days after treatment are typically selected without consideration of cell type, target mRNA half-life, target protein half-life, or siRNA dose.[43,50,51] To investigate whether these arbitrary time points accurately capture the potency of siRNA treatments in various systems, a mathematical model was designed to predict protein expression levels over an extended time frame. The model was used to conduct a retrospective analysis on the kinetics of siRNA-mediated gene silencing reported in existing literature by comparing model predictions of gene expression levels to those in published results.[21,42,43,52−55]

Model Validation

To validate the mathematical model, predicted gene expression levels were compared to experimental data reported by Paganin-Gioanni et al.[21] This study was chosen for the validation process for two main reasons. First, the study characterized protein expression levels at various time points, whereas many studies report data at only one time point.[52−54] A comparison of multiple data points to model results was needed to confirm that the mathematical framework captured not only the extent but also the kinetics, i.e., duration, of gene silencing in a reliable manner. Second, the target gene (green fluorescence protein (GFP)) and cell type (mouse melanoma cells) that the authors used have been characterized thoroughly,[56−58] which facilitates the determination of model parameters such as protein production/degradation rates and siRNA dilution rates.

Though most parameters within the model were determined experimentally or could be calculated from component half-lives, one parameter, cellular internalization rate of siRNA (kint), was difficult to identify.[25,55,59]kint was hard to determine for two main reasons. First, kint was defined within our model as an aggregated parameter encompassing the cellular internalization rate and any necessary post-uptake processing, including endosomal escape, when relevant, and siRNA release from the particle. In other words, kint describes a complex, multistep process. This definition was made because the exact values of the individual rates are not readily available in literature and require meticulous experimental characterization methods to quantify. Second, kint is unique to each system, e.g., delivery vehicle/method and cell type. As a result, the value for kint was fit so that model predictions were in good agreement with experimental data, and this value was subsequently input to compute gene silencing levels as a function of time. For example, data points for GFP expression reported by Paganin-Gioanni et al. were used to fit for a kint of siRNA in B16F10 mouse melanoma cells after electrotransfer (Figure ).[21] Overall, the model captured siRNA-induced changes in protein concentrations in a reliable manner.

Figure 1

Model predictions of GFP expression levels using a fitted internalization rate. Model prediction of GFP expression levels (black line) using a kint fitted from all reported data points as compared to experimentally determined values (blue dots). All data points and error bars were generated using data reported by the original authors.[21]

Accurate Predictions of Gene Expression Behavior Using a Single-Data-Point Fit

After initial validation of the model, we demonstrated the predictive capacity of our simple model even when experimental data were limited to a single data point. For example, we compared predicted Luciferase expression kinetics obtained using either a multiple-data-point fit or a single-data-point fit of experimental data (Figure ).[25] Specifically, an optimized internalization rate was computed from multiple time points or a single time point, and gene expression levels were predicted with these optimized rates as inputs. As shown in Figure , using a single-data-point fit internalization rate produced a prediction for gene silencing efficiency that was very close to the prediction obtained using a multiple-data-point fit internalization rate.

Figure 2

Model predictions obtained from a single-data-point fit versus a multiple-data-point fit. Luciferase expression kinetics modulated by naked siRNA via Oligofectamime predicted by the kinetic model using either a single-data-point fit (red line) or a multiple-data-point fit (black line) compared to experimental data (blue dots). All data points and error bars were generated using data reported by the original authors.[25]

To further validate the robustness and versatility of our model, this comparison was conducted on additional peer-reviewed experimental scenarios (Table ).[21,25,39] When a single-data-point fit approach, e.g., data from 48 h, was implemented, the maximum gene silencing efficiency only differed by 1–3% compared to a multiple-data-point fit approach. The time point of maximum silencing also only differed by 0–2 h. Therefore, we have established that using a single data point to fit for an unknown rate parameter within our kinetic model resulted in reliable predictions that capture important information, i.e., maximum silencing and time point, needed to determine the efficacy of siRNA formulations. With this ability, we anticipate that our model can expedite the testing of various siRNA formulations.

Table 1

Comparisons of Predicted Gene Silencing Efficiencies Using a Single-Data-Point Fit versus Multiple-Data-Point Fit[21,25,39]

example no.	target gene	cell line	delivery vehicle	fit (all vs. 48 h)	max expression	time (h)
1	GFP	B16F10	naked siRNA (electrotransfer)	all	0.44	59
				48 h	0.47	58
2	Luciferase	Neuro2A	Oligofectamine	all	0.39	23
				48 h	0.37	21
3	GFP	HeLa	calcium phosphate nanoparticles	all	0.67	60
				48 h	0.68	60

Due to the simple nature of our model, we expect the mathematical framework to be accessible to a wide variety of users. There is, however, some limitation to the single-data-point fit approach due to the assumptions made throughout the model. More specifically, we suspect that there are limitations to obtaining accurate fits for kint from earlier time points because we are aggregating multiple steps for the internalization process. For reasonable model predictions, users should select a time point consistent with the expected kinetics of the comprehensive internalization process. A suggested time point would be 48 h post-siRNA exposure. This specific time point is recommended because it is a time point longer than most protein half-lives, i.e., a time point that would reflect changes in protein concentrations after siRNA internalization, and common practice is to assess the efficacy of siRNA formulations by determining the extent of gene silencing at 48 h post exposure.[22,23,60,61]

Retrospective Analysis (1): Effect of Target Protein Half-Life on Gene Silencing In Vitro

Protein and mRNA half-lives are important factors that directly impact the duration of siRNA-mediated gene silencing. Proteins with longer half-lives take longer to degrade and will exhibit slower initial responses to siRNA treatments.[25] However, regardless of what the target protein is, common practice is to measure the protein concentration within the cell between 48 and 72 h after siRNA treatment.[21,22,24] For our framework, three separate studies, i.e., Luciferase activity,t1/2 = 2 h; GFP expression, t1/2 = 26 h; glyceraldehyde 3-phosphate dehydrogenase (GAPDH) expression, t1/2 = 35 h, were modeled (Figure )[25,39,40] to determine the effect of protein half-life on the gene-silencing kinetics of proteins with significantly different half-lives. The three studies were chosen because the experiments were carried out in cell lines with similar doubling times, i.e., LNCaP, doubling time = 34 h; HeLa, doubling time = 32 h; H1299, doubling time = 30 h, to avoid variability introduced by siRNA dilution (discussed in the next section). Experimental protein levels were reported at the time points indicated in Figure A–C, and these values were used to validate our mathematical model. Model predictions were in good agreement with the expression levels of both exogenous (Luciferase; GFP) and endogenous (GAPDH) genes, suggesting the versatility in the application of our mathematical approach.

Figure 3

Model predictions of the silencing of different proteins within mammalian cells with similar cell doubling times. Model predictions (lines) and experimental values (dots) for changes in in vitro protein expression levels of (A) Luciferase in LNCaP cells using Oligofectamine, (B) GFP in HeLa cells using calcium phosphate nanoparticles, (C) and GAPDH in H1299 cells using chitosan nanoparticles. (D) Comparison of normalized predictions of gene expression levels for Luciferase (black line), GFP (red line), and GAPDH (blue line). All data points and error bars were generated using data reported by the original authors.[25,39,40]

The predicted time to reach maximum silencing levels increased with protein half-life as expected, with maximum silencing occurring at 1, 2, and 4.5 days post-siRNA treatment for Luciferase, GFP, and GAPDH, respectively. Despite these differences, the expression levels of both GFP and GAPDH often are analyzed experimentally at the same time point, i.e., 24 or 48 h after initial siRNA exposure.[52−54] This choice to look only at one or two time points is based on experimental limitations, and the lack of tools to predict expression kinetics necessitates the selection of generic time points. The selection of generic time points, in turn, frequently results in an incomplete understanding of the effectiveness of an siRNA treatment. For example, Ki et al. examined the effects of siRNA in primary microglial cells (doubling time = 24–48 h) that are difficult to maintain in culture for prolonged periods.[62] As a result, the authors concluded that lower doses of siRNA do not induce sufficient GFP silencing in microglial cells at earlier time points, e.g., 16 h post-siRNA treatment.[63] Considering that the half-life of GFP (t1/2 = 26 h) is relatively long in comparison to other reporter genes such as Luciferase (enzymatic activity t1/2 = 2 h), we modeled GFP expression levels up to longer time points after cells were exposed to siRNA. Informative experimental data reported by Ki et al. were used concurrently with our mathematical model to determine the optimal concentration of siRNA needed to induce significant gene inhibition in primary cells, such as microglial cells.[63]

Retrospective Analysis (2): Effect of Cell Doubling Time on Gene Silencing In Vitro

The impacts of siRNA on gene expression are directly dependent on biological parameters that influence the nucleic acid concentration.[64] Notably, cell division, which results in the dilution of siRNA within the cell, is a parameter that strongly influences siRNA concentration in vitro but is often overlooked in experimental designs when selecting time points to characterize gene expression levels.[25,65] To investigate the significance of cell doubling time on the kinetics of siRNA-mediated gene silencing, we compared model results obtained from two studies that targeted the same protein, i.e., Luciferase, in cell lines of significantly different doubling times, i.e., Neuro2A, doubling time = 19 h; HeLa, doubling time = 32 h (Figure ).[25,100] The Luciferase expression level in Neuro2A cells was predicted to recover back to 100% approximately 4 days after the cells were initially exposed to siRNA. Meanwhile, only 70% of Luciferase expression was predicted 4 days after initial siRNA exposure in slower dividing cells, such as HeLa cells. This trend was in agreement with published work that has shown that gene silencing will differ significantly for fast-dividing, slow-dividing, and nondividing cells.[34,38,65] As an example, Raemdonck et al. noted that in fast-dividing cell lines, in vitro gene silencing usually lasted for only approximately 4–7 days, whereas gene inhibition was sustained for several weeks in nondividing cells.[65] These findings suggest that experiments need to be customized on the basis of the cell type to design effective dosing regimens. The mathematical model presented in this work can guide selection of appropriate time points to measure mRNA and/or protein expression levels. Ultimately, this model can expedite the evaluation process of siRNA treatments by identifying the most essential characterization time points.

Figure 4

Model predictions of the silencing of Luciferase in mammalian cells with different cell doubling times. Model predictions (lines) and experimental values (dots) for changes in in vitro Luciferase expression levels induced by siRNA in (A) Neuro2A cells using Oligofectamine (B) and HeLa cells using multifunctional envelope-type nano devices. (C) Comparison of normalized predictions of Luciferase expression levels in Neuro2A cells (black line) and HeLa cells (red line). All data points and error bars were generated using data reported by the original authors.[25,100]

Retrospective Analysis (3): Accounting for Changes in Cell Proliferation Caused by siRNA-Induced Gene Silencing

In previous sections, the mathematical framework was implemented primarily to predict the expression levels of reporter genes. Reporter genes, such as GFP and Luciferase, are commonly used to develop and optimize siRNA experiments because changes in protein levels can be characterized easily, in situ, via fluorescence imaging or luminescence intensity measurements. However, unlike reporter genes, the target genes for siRNAs that are being developed for clinical applications alter the phenotype and/or behavior of cells to achieve desired therapeutic effects. For instance, siRNA therapeutics for cancer treatment often are designed to reduce growth factor expression and thereby inhibit tumor growth.[66] Considering that we have demonstrated the importance of biological factors, including cell doubling times, when determining the efficacy of siRNA treatments, we hypothesized that changes in cell growth need to be considered within the mathematical model for more accurate predictions.

To test this hypothesis, we compared model predictions that were generated with either fixed or updated values for cell doubling times to experimental results reported by Greco et al. and Arif et al. (Figure ).[38,41] Greco et al. determined that the knockdown of interleukin 1 beta (IL1β) slowed cellular growth rates by ∼30% after a single dose of siRNA.[38] Thus, the kinetic model was encoded with a slower siRNA dilution rate to account for the longer cell doubling time (AoAF cells; doubling time without change in cell proliferation = 38 h, doubling time with change in cell proliferation = 54 ± 6 h) (Figure A). Model predictions utilizing a longer doubling time were in closer agreement with experimental results versus predictions using the original doubling time.

Figure 5

Model predictions of protein expression levels with and without changes in cell proliferation considered. In vitro model predictions generated with the original (dotted lines) and updated (solid lines) cell doubling times were compared to experimental results (blue dots) for (A) IL1β silencing using photoresponsive polymer nanoparticles and (B) voltage dependent anion channel 1 (VDAC1) silencing using jetPRIME. Shaded area: Standard deviation of model predictions associated with the error in measured cell doubling times. All data points and error bars were generated using data reported by the original authors.[38,41]

In the above analysis, it is worth nothing that Greco et al. characterized the change in cell proliferation at a single time point (72 h),[38] and as a result, gene-silencing kinetics induced by a second dose of siRNA were predicted using a step function for the cell doubling time parameter. To test the validity of assuming a step change in cell doubling time, we evaluated an alternative approach in which doubling time was updated more gradually. For this reason, protein expression levels reported by Arif et al. were modeled because this paper characterized changes in the growth of A549 cells as a function of time.[41] The authors reported a 90% decrease in cell proliferation rate following delivery, and this reduced proliferation rate remained constant for at least 144 h following siRNA exposure.[41] Again, the mathematical model was encoded with this information, and model predictions were in good agreement with experimental results once the slower siRNA dilution rate was input (Figure B). Thus, our mathematical framework can be applied to predict expression levels of proteins that affect cell growth by accounting for changes in doubling time via the model. Furthermore, we have confirmed that attenuation in cell proliferation after a dose of siRNA can be modeled as a step change using experimental results from both our own laboratory (Greco et al.,Figure A) and Arif et al. (Figure B).[38,41] Although model predictions were more accurate when cell doubling times were updated more frequently (Figure B), a single update in doubling time also resulted in reasonable agreement between model and experimental results (Figure A).

Translation of In Vitro Results to In Vivo Animal Studies

One of the most prominent challenges in the clinical application of siRNA treatments is that formulations that are initially promising in vitro do not retain efficacy in animal studies.[18] This lack of correlation of in vitro and in vivo results can stem from not considering differences in biological parameters that directly affect intracellular siRNA concentration. More specifically, a dominant factor for changes in siRNA concentration shifts from the doubling of cells in vitro to clearance of siRNA at the target site in vivo. Here, we demonstrate that when siRNA dilution as a result of clearance from target tissue is accounted for quantitatively, experimental results can be used to predict siRNA-mediated gene-silencing kinetics in murine models.

For the modeling of in vivo studies, additional terms were added to the mathematical framework to account for changes in extracellular siRNA concentration. Biological factors that impact siRNA concentration before cellular uptake include extracellular delivery barriers that are unique to in vivo studies, such as rapid clearance of nanoparticles from the bloodstream and nonspecific accumulation of nanoparticles in tissues.[67] Parameters prior to cellular processing were simplified to two terms: accumulation and clearance from the tissue of interest. One advantage of simplifying parameters was that we aggregated the series of hard-to-measure transport steps into a single accumulation parameter.[68] Another advantage was that the biodistribution of nanoparticles throughout the body of an animal model is measurable.[45,46,69,70] Thus, instead of using a kinetic rate, we decided to fix the accumulation based on fractional biodistribution values characterized in previous studies.[45,46,69,70] We chose to do so because transport steps before accumulation into the target tissue, such as clearance from blood, occur rapidly and can be assumed to be instantaneous in comparison to the rates of later steps.[71] By fixing tissue accumulation, our model focused on capturing slower steps, such as the clearance of siRNA from target tissue and intracellular processing after cellular uptake, that govern the overall kinetics of siRNA-mediated gene silencing. The fractional tissue accumulation was estimated from published literature on the basis of particle size.[45,46,69,70]

The process of predicting in vivo gene silencing from in vitro experimental data is depicted in the schematic shown in Figure . Herein, we modeled in vitro and in vivo studies published by Rozema et al. to demonstrate the ability to accurately predict in vivo gene silencing levels using parameter values informed by the in vitro experiments.[42] Similar to the rationale for choosing results published by Paganin-Gioanni et al. for the validation of the model (discussed in the “Methods and Experimental Section” section),[21] Rozema et al. was chosen because it characterized mRNA expression levels at multiple time points post-siRNA treatment.[42] First, the internalization rate was obtained by fitting in vitro data reported by Rozema et al. to the mathematical model using the same procedure detailed within the “Model Validation” section.[42] Next, accumulation within the target tissue, i.e., liver, was estimated on the basis of the particle size, which was 10 nm.[42] According to the studies of Pérez-Campaña et al., for particles with a size of 10 nm, approximately 5% of the injected dose accumulated in the liver after intravenous administration.[70] This value was input into the in vivo mathematical model, and subsequently, the clearance rate was varied to fit data in the work of Rozema et al.[42] When an appropriate value was used for the clearance rate, the mathematical framework gave an accurate prediction of the apolipoprotein B (apoB) mRNA levels in the livers of mice. In particular, the overall trend predicted by the mathematical framework was in good agreement with the trend determined experimentally (Figure B).[42] This trend includes when maximum silencing is expected to happen, along with the duration of silencing, which are crucial pieces of information needed for the design of efficient dosing regimens.

Figure 6

Schematic of how key parameters were obtained for in vivo modeling and the corresponding results. (A) Internalization rate, which is unique to the delivery system and cell type, was obtained by fitting in vitro experimental data to the mathematical framework. The model prediction resulting from this fitted internalization rate (blue) was compared to experimental data (black). (B) In vivo model predictions (lines) were compared to in vivo experimental values (black dots) of mRNA expression levels in response to a single siRNA dose over time induced by siRNA-polymer conjugates (blue line, clearance rate = 0.03 h–1; green line, clearance rate = 0.02 h–1; red line, clearance rate = 0.01 h–1). All data points and error bars were generated using data reported by the original authors.[42]

To further assess the accuracy and robustness of the in vivo mathematical model in various systems, the same process was repeated using experimental results reported by two additional studies.[43,44] The first study tested the gene silencing efficacy of siRNA–lipid nanoparticle (LNP) complexes in mice using siRNA targeted against the apoB gene.[43] The particle size was estimated to be 100 nm on the basis of the composition of the lipid nanoparticle, which would result in 0.1% of the initial injected dose accumulating in the liver.[72] Again, the trends predicted by the mathematical model were in good agreement with experimental data shown in Figure A. The second study used siRNA LNPs of a slightly smaller size (70 nm) to modulate the expression of a different gene (PTEN) in the livers of mice.[44] Considering the size and composition of the LNP, tissue accumulation was estimated to be around 0.08%.[72] As shown in Figure B, model predictions for PTEN expression levels in vivo were close to the values characterized experimentally.

Figure 7

Comparison of model predictions and experimental data of mRNA expression levels in in vivo studies. Model prediction (black line) and experimental values (blue dots) of mRNA expression levels in C57BL/6 mouse models treated with siRNA-LNP formulations. (A) Changes in apoB mRNA expression levels induced by LNPs (ionizable lipid: CLinDMA, size = 100 nm) predicted by our model are compared to experimental values. (B) Changes in PTEN mRNA expression levels induced by LNPs (ionizable lipid: DLin-KC2-DMA, size = 70 nm) predicted by our model are compared to experimental values. All data points and error bars were generated using data reported by the original authors.[43,44]

This facile implementation of our model to describe multiple in vitro and in vivo systems was possible due to the simplicity of the mathematical framework. The framework only included seven crucial rate-limiting kinetic terms that (1) account for changes in intracellular concentrations of key RNAi components and (2) can be easily determined. For example, although previous studies have shown that RISC has multiple intermediates, we assumed that the concentration of the active RISC,[25,73] which can recognize and cleave mRNA, depends solely on the binding of free siRNA to RISC. This assumption was reasonable because the formation of the siRNA–RISC intermediate has the slowest kinetic rate of the steps in active RISC formation.[73] Furthermore, the value of this particular rate is known.[25] Despite this distillation of multiple kinetic processes into a small set of key parameters, our model enabled users to predict gene expression levels with considerable accuracy.

The ability to predict gene-silencing kinetics is beneficial because of the high cost of siRNA molecules and extensive time required to conduct gene silencing studies, which has traditionally limited the systematic exploration of parameter space needed to achieve the most effective siRNA dosing schedule. In vitro studies that are designed to assess the potency of siRNA formulations and guide in vivo dosing schedules often fail to shed light on the gene-silencing behavior of different siRNA treatments when translated to in vivo preclinical studies.[18] As a result, researchers test the efficacies of siRNA treatments directly in in vivo studies to identify and optimize disease-specific dosing schedules that achieve therapeutically relevant silencing profiles.[74] However, this process is extremely costly considering that in vivo experiments require large amounts of expensive materials and resources, e.g., siRNA, animals, and time, and these costs ultimately hamper the testing of promising formulations. For instance, replacement of an oxygen in phosphate with boron (boranophosphate) within the siRNA strand has shown to significantly enhance the stability of siRNA against nucleases.[75] Although this modification is promising, the lack of a method to synthesize the large quantities of boranophosphate-modified siRNA needed to assess the effectiveness of this chemically modified siRNA in vivo restricts clinical application.[76] Such limitations and challenges can be addressed by using our simple and versatile mathematical model to maximize the efficiency of testing siRNA formulations both in vitro and in vivo.

Conclusions

We developed a mathematical framework that can predict siRNA-mediated gene silencing in a reliable manner to efficiently assess the efficacies of siRNA formulations across various experimental setups. Despite the simplicity of the model, predictions were in excellent agreement with both in vitro and in vivo experimental results across a retrospective analysis spanning a wide range of literature. In particular, we highlighted the utility and accessibility of our simple model for making accurate predictions using as few as one experimental data point as an input. Model predictions informed the identification of siRNA dilution as the primary determinant of gene-silencing kinetics. Furthermore, this rate was governed by different biological factors within in vitro versus in vivo systems. Specifically, siRNA dilution depended solely on cell division in vitro and on clearance from target tissue in vivo. This discrepancy highlighted by the model shed light on the traditional lack of in vitro–in vivo translatability for siRNA delivery systems when they are analyzed for potency with generic protocols. With our model, when an appropriate value was identified for the dominant dilution factor, in vitro and in vivo efficacies of siRNA formulations reported across a wide range of literature were predicted with the input of only seven and nine parameter values, respectively. Thus, this work provides a straightforward, versatile tool to articulate the critical aspects of siRNA delivery efficacies and provide guidance on the design of both optimal dosing regimens and the design of nanocarrier systems with desirable performances. Furthermore, while this work focuses on siRNA, we anticipate that our model can be easily extended to other nucleic acids, such as mRNA formulations, that rely on similar gene expression processes.

Methods and Experimental Section

Mathematical Model

The objectives of this study were to develop a simple mathematical model that can predict the gene-silencing kinetics induced by siRNA and ultimately guide in vitro and in vivo experimental designs. We determined that only the most important kinetic parameters were needed within the mathematical framework to get quantitative agreement with experimental data. A sketch of the major components included within the model is shown in Figure . As depicted in Figure , the RNAi process began with the internalization of external siRNA into the cytoplasm of a cell. The ordinary differential equation (ODE) for external siRNA consisted of two rate limiting steps: kint and kclear (eq ).

Figure 8

Simplified schematic of the key RNAi components for siRNA-mediated gene silencing. Schematic of the key RNAi components (siRNA, RISC, mRNA, and protein) and the parameters that a user need to characterize to implement the kinetic model.

Although the internalization process of siRNA into cells consisted of multiple steps, e.g., cellular uptake, endosomal escape, and unpackaging of siRNA, we lumped kint into a single parameter because the rates for the individual steps were not characterized in literature.[68] In the in vitro model, kclear was assumed to be zero, whereas the value of kclear was varied in the in vivo model to get good agreement between predictions and experimental results. Subsequently, the concentration of intracellular siRNA directly depended on the amount of internalized siRNA and the amount of siRNA lost either via degradation (ks,deg) or dilution (ks,dil) as cells divide (eq ).[25,37,38,65]

Intracellular siRNA then interacted with the RNAi machinery, in which cytosolic siRNA molecules associated with free RISC to form active RISC. The active RISC was able to recognize and cleave mRNA molecules that have complementary base pairs to the guide siRNA. Although RISC had multiple intermediates,[73] we assumed that the concentration of active RISC, which can cleave mRNA, was limited by the formation of the siRNA–RISC intermediate that had the slowest kinetic rate.[25] Therefore, the ODE captured changes in the active RISC concentration by accounting for the formation (kRISC) and degradation (kr,deg) of the siRNA–RISC species (eq ). The total amount of free RISC proteins within a cell was estimated to remain constant within the cell.[77]

The active RISC was able to bind to and cleave target mRNA molecules within the cell. Accordingly, the concentration of mRNA was modeled to depend on the cleavage rate (kcleav), production rate (kmRNA), and degradation rate (km,deg) of mRNA (eq ).

Changes in mRNA concentrations available for translation within a cell ultimately resulted in protein attenuation. Thus, the protein concentration was considered to depend directly on the translation of intracellular mRNA molecules (kprot) and the inherent degradation of protein within the cytosol (kp,deg) (eq ).

Model Implementation

The series of ODEs used to model the changes in the intracellular concentrations of key RNAi components was solved using Euler’s method in Python, and we computed time courses of species concentrations. Sample codes of the mathematical framework for both in vitro and in vivo modeling can be found in the Supporting Information. The model comprised of eleven kinetic parameters, three geometric parameters, and six model variables, with four nonzero initial conditions.

Publication Selection and Data Extraction

We extracted experimental data from figures found in various published studies using an online data extraction software, WebPlotDigitzer.[78]

Model Parameterization

Initial Parameter Settings

As a starting point, we first set the parameter values based on information from various sources, including experimental studies and previously established mathematical models.[25,34,37,38,55,59] The values used for model variables and parameters are shown in Tables and 3, respectively. Degradation rates were calculated from component half-lives, and the production rates of mRNA and protein were fit so that steady-state values were achieved in the absence of siRNA. The specific parameter values that were used for the examples presented in this work and the citations for those values can be found in the Supporting Information.

Table 2

Model Variables

variable (name)	description (units)
siRNAex	extracellular siRNA (nM)
siRNAin	intracellular siRNA (nM)
RISC	activated RISC (RISC + siRNA) (nM)
mRNA	intracellular mRNA (nM)
DNA	intracellular DNA (nM)
prot	intracellular protein (nM)

Table 3

Model Parameters

parameter (name)	description	determination	value (units)
kint	internalization rate of external siRNA	fit to each system	(h–1)
kclear	clearance rate of external siRNA	Supporting Information	(h–1)
Vex	transfection volume	specified experimentally	(L)
Vin	intracellular volume	Bartlett et al, Wittrup et al.[25,55]	2 × 10–12 (L)
ks,dil	dilution of siRNA due to cell division	Supporting Information	calculated from cell doubling time (h–1)
ks,deg	intracellular degradation of siRNA	Supporting Information	calculated from siRNA half-life (h–1)
RISCtot	intracellular RISC concentration	Bartlett et al.	3.2 (nM)
kRISC	formation of active RISC	Bartlett et al.[25]	1.2 × 10–4 (nM–1 h–1)
kr,deg	intracellular RISC degradation	Bartlett et al.[25]	7.7 × 10–2 (h–1)
kmRNA	intracellular mRNA production rate	steady-state approximation, see the Supporting Information	calculated from mRNA half-life (h–1)
km,deg	intracellular mRNA degradation rate	Supporting Information	calculated from mRNA half-life (h–1)
kcleav	cleavage rate of mRNA	Materials and Methods	5 × 102 (nM–1 h–1)
kprot	intracellular protein production rate	steady-state approximation, Supporting Information	calculated from protein half-life (h–1)
kp,deg	intracellular protein degradation rate	Supporting Information	calculated from protein half-life (h–1)

Model Fitting

For the in vitro mathematical model, the values for two kinetic rates, kcleav and kint, were either not characterized experimentally or unique to each system. The value for kcleav was obtained by fitting the mathematical model to experimental data (Figure ).[55] The referenced experimental study quantified the amount of cytosolic siRNA and its effect on GFP expression levels. The specific publication was chosen to reduce the number of unknown parameters from two (kint and kcleav) to one (kcleav) within the in vitro mathematical framework and thereby obtain an accurate fit for kcleav. Specifically, kint could be eliminated because the intracellular siRNA concentration after cellular uptake was characterized by Wittrup et al.,[55] leaving kcleav as the only unknown parameter. Model predictions accurately captured GFP silencing kinetics and reflected a rapid cleavage rate (Figure ). The fitted value was used in further analyses assuming that RISC has intrinsic activity.[59]kint was obtained using details described in the “Results and Discussion” section.

Figure 9

Comparison of model predictions and experimental data to obtain kcleav. Model predictions (black) and experimental data (blue) of relative GFP levels in cells with an average cytosolic siRNA concentration of 1.6 nM delivered via Lipofectamine. All data points and error bars were generated using data reported by the original authors.[55]

The in vivo kinetic model required knowledge of an additional unknown parameter, kclear. This parameter depended on both the particle size and surface modifications. kclear was varied so that model predictions were in good agreement with experimental results.