Generating Effective Models and Parameters for RNA Genetic Circuits

RNA genetic circuitry is emerging as a powerful tool to control gene expression. However, little work has been done to create a theoretical foundation for RNA circuit design. A prerequisite to this is a quantitative modeling framework that accurately describes the dynamics of RNA circuits. In this work, we develop an ordinary differential equation model of transcriptional RNA genetic circuitry, using an RNA cascade as a test case. We show that parameter sensitivity analysis can be used to design a set of four simple experiments that can be performed in parallel using rapid cell-free transcription-translation (TX-TL) reactions to determine the 13 parameters of the model. The resulting model accurately recapitulates the dynamic behavior of the cascade, and can be easily extended to predict the function of new cascade variants that utilize new elements with limited additional characterization experiments. Interestingly, we show that inconsistencies between model predictions and experiments led to the model-guided discovery of a previously unknown maturation step required for RNA regulator function. We also determine circuit parameters in two different batches of TX-TL, and show that batch-to-batch variation can be attributed to differences in parameters that are directly related to the concentrations of core gene expression machinery. We anticipate the RNA circuit models developed here will inform the creation of computer aided genetic circuit design tools that can incorporate the growing number of RNA regulators, and that the parametrization method will find use in determining functional parameters of a broad array of natural and synthetic regulatory systems

Reduced order complement model training.

<p>An ensemble of model parameters was estimated using multiobjective optimization from dynamic C3a and C5a measurements with and without zymosan. The model was trained using C3a and C5a measurements in the absence of zymosan (<b>A</b>–<b>B</b>) or in the presence of 1 mg/ml zymosan (<b>C</b>–<b>D</b>). The solid black lines show the simulated mean value of C3a or C5a for the ensemble, while the dark shaded region denotes the 99% confidence interval of mean. The light shaded region denotes the 99% confidence interval of the simulated C3a and C5a concentration. The experimental training data (points) was taken from Morad et al [<a href="http://www.plosone.org/article/info:doi/10.1371/journal.pone.0187373#pone.0187373.ref035" target="_blank">35</a>]. All initial conditions not specified by the experimental condition were assumed to be at zero or their physiological serum levels unless otherwise noted (<a href="http://www.plosone.org/article/info:doi/10.1371/journal.pone.0187373#pone.0187373.s001" target="_blank">S1 Table</a>).</p

Reduced order complement model predictions.

<p>Simulations of C3a and C5a generated using 0.1 mg/ml, 0.01 mg/ml, and 0.001 mg/ml zymosan were compared with the corresponding experimental measurements. The solid black lines show the simulated mean value of C3a or C5a for the ensemble, while the dark shaded region denotes the 99% confidence interval of mean. The light shaded region denotes the 99% confidence interval of the simulated C3a and C5a concentration. The experimental validation data (points) was taken from Morad et al [<a href="http://www.plosone.org/article/info:doi/10.1371/journal.pone.0187373#pone.0187373.ref035" target="_blank">35</a>]. All initial conditions not specified by the experimental condition were assumed to be at zero or their physiological serum levels unless otherwise noted (<a href="http://www.plosone.org/article/info:doi/10.1371/journal.pone.0187373#pone.0187373.s001" target="_blank">S1 Table</a>).</p

Reduced order modeling and analysis of the human complement system

<div><p>Complement is an important pathway in innate immunity, inflammation, and many disease processes. However, despite its importance, there are few validated mathematical models of complement activation. In this study, we developed an ensemble of experimentally validated reduced order complement models. We combined ordinary differential equations with logical rules to produce a compact yet predictive model of complement activation. The model, which described the lectin and alternative pathways, was an order of magnitude smaller than comparable models in the literature. We estimated an ensemble of model parameters from <i>in vitro</i> dynamic measurements of the C3a and C5a complement proteins. Subsequently, we validated the model on unseen C3a and C5a measurements not used for model training. Despite its small size, the model was surprisingly predictive. Global sensitivity and robustness analysis suggested complement was robust to any single therapeutic intervention. Only the simultaneous knockdown of both C3 and C5 consistently reduced C3a and C5a formation from all pathways. Taken together, we developed a validated mathematical model of complement activation that was computationally inexpensive, and could easily be incorporated into pre-existing or new pharmacokinetic models of immune system function. The model described experimental data, and predicted the need for multiple points of therapeutic intervention to fully disrupt complement activation.</p></div

Simulated VEGF-A and TGF−<i>β</i>1/2 exposure promoted phenotype heterogeneity.

<p>Simulated response to TGF−<i>β</i>1/2 and VEGF-A exposure with and without axis specific inhibitors. Vimentin and E-cadherin abundances (in nM) were used to quantify the shift in population at 48 hrs. (A-C) VEGF-A (50 a.u.) treatment resulted in a population with enhanced epithelial properties. This was contrary to the addition of TGF−<i>β</i>2 (10 a.u.), which shifted the population towards a mesenchymal phenotype. Interestingly, the combined effects of TGF−<i>β</i>2 and VEGFA was found to increase both ecadherin and vimentin levels, creating a heterogeneous population (black arrow), which can also be seen in a minority of untreated cells (gray arrow). (D-F) To isolate the effect of NFAT, we inhibited NFAT de-phosphorylation in combination with VEGFA. This negated the increase in ecadherin expression and shifted the population towards a mesenchymal phenotype. Likewise, combining NFAT inhibition with TGF−<i>β</i> mitigated all VEGF enhanced ecadherin expression. Lastly, combination of TGF−<i>β</i>2, VEGFA, and NFAT inhibition nearly mitigated all effects of VEGFA, shifting the heterogeneous population towards a mesenchymal phenotype. In whole, high levels of phosphorylated-Sp1 correlated with vimentin expression, while NFAT was responsible for maintaining E-cadherin expression in the presence of other factors, although neither were mutually exclusive.</p

Model connectivity recreates the core architecture during EMT.

<p>The EMT network contains 97 nodes (proteins, mRNA, and genes) interconnected by 169 interactions. Central to EMT induction, activation of the MAPK cascade occurs through TGF-<i>β</i>1/2 binding which activates the AP-1/Sp1 transcriptional axis. AP-1/Sp1 drives an autocrine response of TGF-<i>β</i>3, which activates the Smad cascade, leading to phenotypic change. Conversely, VEGF-A binding can stabilize an epithelial phenotype through NFAT activation. Downstream activation of <i>β</i>-catenin signaling due to E-cadherin loss and GSK3 inactivation of <i>β</i>-catenin confinement is critical to the complete activation of the EMT program. The complete list of molecular interactions that comprise the model is given in the supplement.</p

Pairwise sensitivity and clustering of complement model parameters in the presence of 1 mg/ml zymosan.

<p>The response of the complement model was calculated for each parameter combination following a 10% increase in parameter combinations in the presence of 1 mg/ml zymosan. The model parameters were clustered into high (blue), medium (red) and low (green) response clusters based upon the euclidian distance between the perturbed and nominal system state (no perturbation).</p

Simultaneous TGF-<i>β</i>1/2 and VEGF-A treatment induced phenotype heterogeneity and is dependent upon NFAT activity <i>in-vitro</i>.

<p>(A) In MCF10A, treatment with (10ng/ml) TGF−<i>β</i>2 increased Slug and vimentin, while ecadherin expression was inhibited at both the gene and protein level at 48 hrs. Conversely, VEGFA alone increased both NFATc1 and ecadherin gene expression. Simultaneous TGF−<i>β</i>2 (10ng/ml) and VEGFA (50ng/ml) treatment increased Slug, NFATc1, and vimentin expression, while also increasing ecadherin levels via qPCR. (B-C) Immunofluorescence confirmed these results and nuclear co-localization of both phospho-Sp1 and NFAT were found dependent upon TGF−<i>β</i>2 and VEGFA, respectively. (D) To isolate the effect of NFAT, treatment of VEGFA (50ng/ml) and VIVIT (10<i>μ</i>M) reduced ecadherin expression at 48hrs (control-dashed line). Similarly, combined TGF−<i>β</i>2, VEGFA and VIVIT treatment increased Slug and vimentin expression, while inhibiting ecadherin levels via qPCR. (E) These findings were confirmed via immunofluorescence as the VIVIT peptide inhibited ecadherin and nuclear localization of NFATc1 in all three cases. (F) Quantitative flow cytometry also confirmed this trend. Similar experiments in DLD1 followed a similar trend (supplement). Magnification, 40x. Scale bars: 50<i>μ</i>m. C = Control, T = TGF−<i>β</i>2, V = VEGFA, VI = NFAT inhibitor (VIVIT). Asterisks signify statistical differences from each other according to a one-way ANOVA with Tukey’s post hoc (p≺0.05).</p

Analysis of underlying signaling responses.

<p>(A) We examined the distribution of NFATc1 and AP1/SP1 in cells containing the hybrid phenotype (VEGF-A + TGF−<i>β</i>2 case), showing the potential for cells to express both SP1 and NFATc1 in a non exclusive manner. (B) In the absence of YREG1, most of the population failed to consistently to retain a stable epithelial state. (C) We identified a novel regulator of LEF1 called YREG1 that allows Snail/Slug to emulate an inducer by repressing YREG1, which was required to stabilize the untreated population. YREG1 overexpression revealed an enhanced epithelial phenotype, while some inherently transformed cells moved towards a hybrid phenotype.</p

Robustness analysis of the complement model.

<p>Robustness coefficients were calculated for a 50%, 90% and 99% reduction in C3, C5, or C3 and C5 initial conditions. <b>A</b>: Mean robustness index for C3a and C5a generated in the absence of zymosan. <b>B</b>: Mean robustness index for C3a and C5a generated in the presence of 1 mg/ml zymosan. The color describes the degree of reduction of C3a or C5a following the network perturbation. Robustness coefficients were calculated using all parameter sets with Pareto rank less than five (N = 65). Mean robustness values were reported.</p