Comparative study of parameter sensitivity analyses of the TCR-activated Erk-MAPK signalling pathway

Mathematical models of TCR initial triggering

T cell receptors (TCRs) play crucial roles in regulating T cell response by rapidly and accurately recognizing foreign and non-self antigens. The process involves multiple molecules and regulatory mechanisms, forming a complex network to achieve effective antigen recognition. Mathematical modeling techniques can help unravel the intricate network of TCR signaling and identify key regulators that govern it. In this review, we introduce and briefly discuss relevant mathematical models of TCR initial triggering, with a focus on kinetic proofreading (KPR) models with different modified structures. We compare the topology structures, biological hypotheses, parameter choices, and simulation performance of each model, and summarize the advantages and limitations of them. Further studies on TCR modeling design, aiming for an optimized balance of specificity and sensitivity, are expected to contribute to the development of new therapeutic strategies.

Machine learning ranking of plausible (un)explored synergistic gene combinations using sensitivity indices of time series measurements of Wnt signaling pathway

Combinations of genes or proteins work in synergy at different times and durations in a signaling pathway. However, which combinations are prevalent at a particular time point or duration is mostly not known. Sensitivity analysis plays a major role in computing the strength of the influence of involved factors in any phenomena under investigation. When applied to expression profiles of various intra/extracellular factors that work in a signaling pathway, the variance- and density-based analysis yields a range of sensitivity indices for individual and various combinations of factors. These combinations denote the higher order interactions among the involved factors, which might be of interest. In this work, after estimating the individual effects of factors for a higher order combination, the individual indices are considered as discriminative features. Exploiting the analogy of prioritizing webpages using ranking algorithms, for a particular order, a full set of combinations of genes can be prioritized based on these features using a powerful support vector ranking algorithm. Recording the changing rankings of the combinations over time points and durations reveals which higher order combinations influence the pathway and when and where an intervention might be necessary to affect the pathway. Integration, innovation, and insight Combinations of genes or proteins work in synergy at different times and durations in a signaling pathway. However, which combinations are prevalent at a particular time point or duration is mostly not known. This work develops a search engine that reveals ground-breaking results in the form of higher order (un)explored/(un)tested combinations (as biological hypotheses), based on sensitivity indices. These indices capture the strength of influence of factors (here genes/proteins) that affect a signaling pathway. Recording the changing rankings of these combinations over time points and durations reveals how higher order combinations behave within the pathway. Significance The manuscript develops a search engine that reveals ground-breaking results in the form of higher order (un)explored/(un)tested combinations of genes/proteins (as biological hypotheses), based on sensitivity indices that capture the strength of influence of factors (here genes/proteins) that affect the Wnt signaling pathway. The pipeline uses kernel-based sensitivity indices to capture the influence of the factors in a pathway and employs powerful support vector ranking algorithm. Because of the above point, biologists/oncologists will be able to narrow down their search to particular combinations that are ranked and, if a synergistic functioning is confirmed, will be able to study the mechanism between the components of a combination, in the Wnt pathway. The search engine design is not only limited to one dataset and a range of combinations of genes/proteins. The framework can be applied/modified to all problems where one is interested in searching for particular combinations of factors involved in a particular phenomena. Recording the changing rankings of the combinations over time points and durations reveals how higher order interactions behave within the pathway and when and where an intervention might be necessary to influence the pathway, for therapeutic purpose. It reveals the various unexplored FZD-WNT combinations that have been untested till now in the Wnt pathway.

The impact of a power law-induced memory effect on the SARS-CoV-2 transmission

It is well established that COVID-19 incidence data follows some power law growth pattern. Therefore, it is natural to believe that the COVID-19 transmission process follows some power law. However, we found no existing model on COVID-19 with a power law effect only in the disease transmission process. Inevitably, it is not clear how this power law effect in disease transmission can influence multiple COVID-19 waves in a location. In this context, we developed a completely new COVID-19 model where a force of infection function in disease transmission follows some power law. Furthermore, different realistic epidemiological scenarios like imperfect social distancing among home-quarantined individuals, disease awareness, vaccination, treatment, and possible reinfection of the recovered population are also considered in the model. Applying some recent techniques, we showed that the proposed system converted to a COVID-19 model with fractional order disease transmission, where order of the fractional derivative ( α ) in the force of infection function represents the memory effect in disease transmission. We studied some mathematical properties of this newly formulated model and determined the basic reproduction number (R 0 ). Furthermore, we estimated several epidemiological parameters of the newly developed fractional order model (including memory index α ) by fitting the model to the daily reported COVID-19 cases from Russia, South Africa, UK, and USA, respectively, for the time period March 01, 2020, till December 01, 2021. Variance-based Sobol's global sensitivity analysis technique is used to measure the effect of different important model parameters (including α ) on the number of COVID-19 waves in a location (W C ). Our findings suggest that α along with the average transmission rate of the undetected (symptomatic and asymptomatic) cases in the community (β 1 ) are mainly influencing multiple COVID-19 waves in those four locations. Numerically, we identified the regions in the parameter space of α andβ 1 for which multiple COVID-19 waves are occurring in those four locations. Furthermore, our findings suggested that increasing memory effect in disease transmission ( α → 0) may decrease the possibility of multiple COVID-19 waves and as well as reduce the severity of disease transmission in those four locations. Based on all the results, we try to identify a few non-pharmaceutical control strategies that may reduce the risk of further SARS-CoV-2 waves in Russia, South Africa, UK, and USA, respectively.

Global sensitivity analysis of a single-cell HBV model for viral dynamics in the liver

The predictive accuracy of mathematical models representing anything ranging from the meteorological to the biological system profoundly depends on the quality of model parameters derived from experimental data. Hence, robust sensitivity analysis (SA) of these critical model parameters aids in sifting the influential from the negligible out of typically vast parameter regimes, thus illuminating key components of the system under study. We here move beyond traditional local sensitivity analysis to the adoption of global SA techniques. Partial rank correlation coefficient (PRCC) based on Latin hypercube sampling is compared with the variance-based Sobol method. We selected for this SA investigation an infection model for the hepatitis-B virus (HBV) that describes infection dynamics and clearance of HBV in the liver [Murray & Goyal, 2015]. The model tracks viral particles such as the tenacious and nearly ineradicable covalently closed circular DNA (cccDNA) embedded in infected nuclei and an HBV protein known as p36. Our application of these SA methods to the HBV model illuminates, especially over time, the quantitative relationships between cccDNA synthesis rate and p36 synthesis and export. Our results reinforce previous observations that the viral protein, p36, is by far the most influential factor for cccDNA replication. Moreover, both methods are capable of finding crucial parameters of the model. Though the Sobol method is independent of model structure (e.g., linearity and monotonicity) and well suited for SA, our results ensure that LHS-PRCC suffices for SA of a non-linear model if it is monotonic.

A quantitative systems pharmacology model of hyporesponsiveness to erythropoietin in rats

Recombinant human erythropoietin (rHuEPO) is effective in managing chronic kidney disease and chemotherapy-induced anemia. However, hyporesponsiveness to rHuEPO treatment was reported in about 10% of the patients. A decreased response in rats receiving a single or multiple doses of rHuEPO was also observed. In this study, we aimed to develop a quantitative systems pharmacology (QSP) model to examine hyporesponsiveness to rHuEPO in rats. Pharmacokinetic (PK) and pharmacodynamic (PD) data after a single intravenous dose of rHuEPO (100 IU/kg) was obtained from a previous study (Yan et al. in Pharm Res, 30:1026-1036, 2013) including rHuEPO plasma concentrations, erythroid precursors counts in femur bone marrow and spleen, reticulocytes (RETs), red blood cells (RBCs), and hemoglobin (HGB) in circulation. Parameter values were obtained from literature or calibrated with experimental data. Global sensitivity analysis and model-based simulations were performed to assess parameter sensitivity and hyporesponsiveness. The final QSP model adequately characterizes time courses of rHuEPO PK and nine PD endpoints in both control and treatment groups simultaneously. The model indicates that negative feedback regulation, neocytolysis, and depletion of erythroid precursors are major factors leading to hyporesponsiveness to rHuEPO treatment in rats.

Accuracy assessment methods for physiological model selection toward evaluation of closed-loop controlled medical devices

Physiological closed-loop controlled (PCLC) medical devices are complex systems integrating one or more medical devices with a patient’s physiology through closed-loop control algorithms; introducing many failure modes and parameters that impact performance. These control algorithms should be tested through safety and efficacy trials to compare their performance to the standard of care and determine whether there is sufficient evidence of safety for their use in real care setting. With this aim, credible mathematical models have been constructed and used throughout the development and evaluation phases of a PCLC medical device to support the engineering design and improve safety aspects. Uncertainties about the fidelity of these models and ambiguities about the choice of measures for modeling performance need to be addressed before a reliable PCLC evaluation can be achieved. This research develops tools for evaluating the accuracy of physiological models and establishes fundamental measures for predictive capability assessment across different physiological models. As a case study, we built a refined physiological model of blood volume (BV) response by expanding an original model we developed in our prior work. Using experimental data collected from 16 sheep undergoing hemorrhage and fluid resuscitation, first, we compared the calibration performance of the two candidate physiological models, i.e., original and refined, using root-mean-squared error (RMSE), Akiake information criterion (AIC), and a new multi-dimensional approach utilizing normalized features extracted from the fitting error. Compared to the original model, the refined model demonstrated a significant improvement in calibration performance in terms of RMSE (9%, P = 0.03) and multi-dimensional measure (48%, P = 0.02), while a comparable AIC between the two models verified that the enhanced calibration performance in the refined model is not due to data over-fitting. Second, we compared the physiological predictive capability of the two models under three different scenarios: prediction of subject-specific steady-state BV response, subject-specific transient BV response to hemorrhage perturbation, and leave-one-out inter-subject BV response. Results indicated enhanced accuracy and predictive capability for the refined physiological model with significantly larger proportion of measurements that were within the prediction envelope in the transient and leave-one-out prediction scenarios (P < 0.02). All together, this study helps to identify and merge new methods for credibility assessment and physiological model selection, leading to a more efficient process for PCLC medical device evaluation.

A simple model of immune and muscle cell crosstalk during muscle regeneration

Muscle injury during aging predisposes skeletal muscles to increased damage due to reduced regenerative capacity. Some of the common causes of muscle injury are strains, while other causes are more complex muscle myopathies and other illnesses, and even excessive exercise can lead to muscle damage. We develop a new mathematical model based on ordinary differential equations of muscle regeneration. It includes the interactions between the immune system, healthy and damaged myonuclei as well as satellite cells. Our new mathematical model expands beyond previous ones by accounting for 21 specific parameters, including those parameters that deal with the interactions between the damaged and dead myonuclei, the immune system, and the satellite cells. An important assumption of our model is the replacement of only damaged parts of the muscle fibers and the dead myonuclei. We conduce systematic sensitivity analysis to determine which parameters have larger effects on the model and therefore are more influential for the muscle regeneration process. We propose additional validation for these parameters. We further demonstrate that these simulations are species-, muscle-, and age-dependent. In addition, the knowledge of these parameters and their interactions, may suggest targeting or selecting these interactions for treatments that accelerate the muscle regeneration process.

The prevention and control of tuberculosis: an analysis based on a tuberculosis dynamic model derived from the cases of Americans

Background

Tuberculosis (TB), a preventable and curable disease, is claimed as the second largest number of fatalities, and there are 9,025 cases reported in the United States in 2018. Many researchers have done a lot of research and achieved remarkable results, but TB is still a severe problem for human beings. The study is a further exploration of the prevention and control of tuberculosis.

Methods

In the paper, we propose a new dynamic model to study the transmission dynamics of TB, and then use global differential evolution and local sequential quadratic programming (DESQP) optimization algorithm to estimate parameters of the model. Finally, we use Latin hypercube sampling (LHS) and partial rank correlation coefficients (PRCC) to analyze the influence of parameters on the basic reproduction number (\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$\mathcal R_{0}$\end{document}R0) and the total infectious (including the diagnosed, undiagnosed and incomplete treatment infectious), respectively.

Results

According to the research, the basic reproduction number is computed as 2.3597 from 1984 to 2018, which means TB is also an epidemic in the US. The diagnosed rate is 0.6082, which means the undiagnosed will be diagnosed after 1.6442 years. The diagnosed will recover after an average of 1.9912 years. Moreover, some diagnosed will end the treatment after 1.7550 years for some reason. From the study, it’s shown that 2.40% of the recovered will be reactivated, and 13.88% of the newborn will be vaccinated. However, the immune system will be lost after about 19.6078 years.

Conclusion

Through the results of this study, we give some suggestions to help prevent and control the TB epidemic in the United States, such as prolonging the protection period of the vaccine by developing new and more effective vaccines to prevent TB; using the Chemoprophylaxis for incubation patients to prevent their conversion into active TB; raising people’s awareness of the prevention and control of TB and treatment after illness; isolating the infected to reduce the spread of TB. According to the latest report in the announcement that came at the first WHO Global Ministerial Conference on Ending tuberculosis in the Sustainable Development Era, we predict that it is challenging to control TB by 2030.

Coupled feedback regulation of nuclear factor of activated T-cells (NFAT) modulates activation-induced cell death of T cells

A properly functioning immune system is vital for an organism’s wellbeing. Immune tolerance is a critical feature of the immune system that allows immune cells to mount effective responses against exogenous pathogens such as viruses and bacteria, while preventing attack to self-tissues. Activation-induced cell death (AICD) in T lymphocytes, in which repeated stimulations of the T-cell receptor (TCR) lead to activation and then apoptosis of T cells, is a major mechanism for T cell homeostasis and helps maintain peripheral immune tolerance. Defects in AICD can lead to development of autoimmune diseases. Despite its importance, the regulatory mechanisms that underlie AICD remain poorly understood, particularly at an integrative network level. Here, we develop a dynamic multi-pathway model of the integrated TCR signalling network and perform model-based analysis to characterize the network-level properties of AICD. Model simulation and analysis show that amplified activation of the transcriptional factor NFAT in response to repeated TCR stimulations, a phenomenon central to AICD, is tightly modulated by a coupled positive-negative feedback mechanism. NFAT amplification is predominantly enabled by a positive feedback self-regulated by NFAT, while opposed by a NFAT-induced negative feedback via Carabin. Furthermore, model analysis predicts an optimal therapeutic window for drugs that help minimize proliferation while maximize AICD of T cells. Overall, our study provides a comprehensive mathematical model of TCR signalling and model-based analysis offers new network-level insights into the regulation of activation-induced cell death in T cells.

Vector Preference Annihilates Backward Bifurcation and Reduces Endemicity

We propose and analyze a mathematical model of a vector-borne disease that includes vector feeding preference for carrier hosts and intrinsic incubation in hosts. Analysis of the model reveals the following novel results. We show theoretically and numerically that vector feeding preference for carrier hosts plays an important role for the existence of both the endemic equilibria and backward bifurcation when the basic reproduction number [Formula: see text] is less than one. Moreover, by increasing the vector feeding preference value, backward bifurcation is eliminated and endemic equilibria for hosts and vectors are diminished. Therefore, the vector protects itself and this benefits the host. As an example of these phenomena, we present a case of Andean cutaneous leishmaniasis in Peru. We use parameter values from previous studies, primarily from Peru to introduce bifurcation diagrams and compute global sensitivity of [Formula: see text] in order to quantify and understand the effects of the important parameters of our model. Global sensitivity analysis via partial rank correlation coefficient shows that [Formula: see text] is highly sensitive to both sandflies feeding preference and mortality rate of sandflies.

Clostridium butyricum population balance model: Predicting dynamic metabolic flux distributions using an objective function related to extracellular glycerol content

BACKGROUND

Extensive experimentation has been conducted to increment 1,3-propanediol (PDO) production using Clostridium butyricum cultures in glycerol, but computational predictions are limited. Previously, we reconstructed the genome-scale metabolic (GSM) model iCbu641, the first such model of a PDO-producing Clostridium strain, which was validated at steady state using flux balance analysis (FBA). However, the prediction ability of FBA is limited for batch and fed-batch cultures, which are the most often employed industrial processes.

RESULTS

We used the iCbu641 GSM model to develop a dynamic flux balance analysis (DFBA) approach to predict the PDO production of the Colombian strain Clostridium sp IBUN 158B. First, we compared the predictions of the dynamic optimization approach (DOA), static optimization approach (SOA), and direct approach (DA). We found no differences between approaches, but the DOA simulation duration was nearly 5000 times that of the SOA and DA simulations. Experimental results at glycerol limitation and glycerol excess allowed for validating dynamic predictions of growth, glycerol consumption, and PDO formation. These results indicated a 4.4% error in PDO prediction and therefore validated the previously proposed objective functions. We performed two global sensitivity analyses, finding that the kinetic input parameters of glycerol uptake flux had the most significant effect on PDO predictions. The other input parameters evaluated during global sensitivity analysis were biomass composition (precursors and macromolecules), death constants, and the kinetic parameters of acetic acid secretion flux. These last input parameters, all obtained from other Clostridium butyricum cultures, were used to develop a population balance model (PBM). Finally, we simulated fed-batch cultures, predicting a final PDO production near to 66 g/L, almost three times the PDO predicted in the best batch culture.

CONCLUSIONS

We developed and validated a dynamic approach to predict PDO production using the iCbu641 GSM model and the previously proposed objective functions. This validated approach was used to propose a population model and then an increment in predictions of PDO production through fed-batch cultures. Therefore, this dynamic model could predict different scenarios, including its integration into downstream processes to predict technical-economic feasibilities and reducing the time and costs associated with experimentation.

Robust Control Design of Heart Rate Response during Treadmill Exercise under Parametric Uncertainty

Rehabilitation (Rehab) exercise can benefit cardiac patients as it can promote the recovery and improve the heart wellness. However, heart failure (HF) patients can only take mild exercise, since excessive exercise may lead to fatal events. It is important to control the exercise intensity at a desired level to maximize exercise benefit. Heart Rate (HR) is an essential factor for measuring exercise intensity. Mathematical models of HR can be used to study exercise physiology. However, HR models involve model uncertainty, resulting from model calibration or variability in patients. It is important to quantify the effect of uncertainty on HR prediction for optimizing exercise intensity, such as treadmill speed. A probabilistic model-based control design is presented in this work to obtain an optimal treadmill speed for Rehab exercise in the presence of uncertainty. To obtain a computationally tractable formulation, the generalized polynomial chaos (gPC) theory is used to propagate uncertainty via a model to HR predictions, and predict slow-acting responses such as peripheral local metabolism that can be used to evaluate exercise outcome for individual patients. The speed control of treadmill is formulated as an optimization problem that can maximize the exercise outcome, while minimizing the slow-acting effects. The effectiveness of the proposed control design was experimentally verified with simulations, showing potentials in the exercise control of individual patients.

Sequential Parameter Estimation for Mammalian Cell Model Based on In Silico Design of Experiments

Due to the complicated metabolism of mammalian cells, the corresponding dynamic mathematical models usually consist of large sets of differential and algebraic equations with a large number of parameters to be estimated. On the other hand, the measured data for estimating the model parameters are limited. Consequently, the parameter estimates may converge to a local minimum far from the optimal ones, especially when the initial guesses of the parameter values are poor. The methodology presented in this paper provides a systematic way for estimating parameters sequentially that generates better initial guesses for parameter estimation and improves the accuracy of the obtained metabolic model. The model parameters are first classified into four subsets of decreasing importance, based on the sensitivity of the model’s predictions on the parameters’ assumed values. The parameters in the most sensitive subset, typically a small fraction of the total, are estimated first. When estimating the remaining parameters with next most sensitive subset, the subsets of parameters with higher sensitivities are estimated again using their previously obtained optimal values as the initial guesses. The power of this sequential estimation approach is illustrated through a case study on the estimation of parameters in a dynamic model of CHO cell metabolism in fed-batch culture. We show that the sequential parameter estimation approach improves model accuracy and that using limited data to estimate low-sensitivity parameters can worsen model performance.

Hilbert-Schmidt and Sobol sensitivity indices for static and time series Wnt signaling measurements in colorectal cancer - part A

Background

Ever since the accidental discovery of Wingless [Sharma R.P., Drosophila information service, 1973, 50, p 134], research in the field of Wnt signaling pathway has taken significant strides in wet lab experiments and various cancer clinical trials, augmented by recent developments in advanced computational modeling of the pathway. Information rich gene expression profiles reveal various aspects of the signaling pathway and help in studying different issues simultaneously. Hitherto, not many computational studies exist which incorporate the simultaneous study of these issues.

Results

This manuscript ∙ explores the strength of contributing factors in the signaling pathway, ∙ analyzes the existing causal relations among the inter/extracellular factors effecting the pathway based on prior biological knowledge and ∙ investigates the deviations in fold changes in the recently found prevalence of psychophysical laws working in the pathway. To achieve this goal, local and global sensitivity analysis is conducted on the (non)linear responses between the factors obtained from static and time series expression profiles using the density (Hilbert-Schmidt Information Criterion) and variance (Sobol) based sensitivity indices.

Conclusion

The results show the advantage of using density based indices over variance based indices mainly due to the former’s employment of distance measures & the kernel trick via Reproducing kernel Hilbert space (RKHS) that capture nonlinear relations among various intra/extracellular factors of the pathway in a higher dimensional space. In time series data, using these indices it is now possible to observe where in time, which factors get influenced & contribute to the pathway, as changes in concentration of the other factors are made. This synergy of prior biological knowledge, sensitivity analysis & representations in higher dimensional spaces can facilitate in time based administration of target therapeutic drugs & reveal hidden biological information within colorectal cancer samples.

Global sensitivity analysis for developing biological models: Application to K+ channel model in mouse ventricular myocytes

Mathematical models of cardiac myocytes are highly nonlinear and involve a large number of model parameters. The parameters are estimated using experimental data, which are often corrupted by noise and uncertainty. Such uncertainty can be propagated onto model parameters during model calibration, which further affects model reliability and credibility. In order to improve model accuracy, it is important to quantify and reduce the uncertainty in model response resulting from parametric uncertainty. Sensitivity analysis is a key technique to investigate the significance of parametric uncertainty and its effect on model responses. This can identify and rank most sensitive parameters, and evaluate the effect of uncertainty on model outputs. In this work, a global sensitivity analysis is developed to determine the significance of parametric uncertainty on model responses using Sobol indices. This method is applied to nonlinear K+ channel models of mouse ventricular myocytes to demonstrate the efficacy of the developed algorithm.

An Integrative multi-lineage model of variation in leukopoiesis and acute myelogenous leukemia

Background

Acute myelogenous leukemia (AML) progresses uniquely in each patient. However, patients are typically treated with the same types of chemotherapy, despite biological differences that lead to differential responses to treatment.

Results

Here we present a multi-lineage multi-compartment model of the hematopoietic system that captures patient-to-patient variation in both the concentration and rates of change of hematopoietic cell populations. By constraining the model against clinical hematopoietic cell recovery data derived from patients who have received induction chemotherapy, we identified trends for parameters that must be met by the model; for example, the mitosis rates and the probability of self-renewal of progenitor cells are inversely related. Within the data-consistent models, we found 22,796 parameter sets that meet chemotherapy response criteria. Simulations of these parameter sets display diverse dynamics in the cell populations. To identify large trends in these model outputs, we clustered the simulated cell population dynamics using k-means clustering and identified thirteen ‘representative patient’ dynamics. In each of these patient clusters, we simulated AML and found that clusters with the greatest mitotic capacity experience clinical cancer outcomes more likely to lead to shorter survival times. Conversely, other parameters, including lower death rates or mobilization rates, did not correlate with survival times.

Conclusions

Using the multi-lineage model of hematopoiesis, we have identified several key features that determine leukocyte homeostasis, including self-renewal probabilities and mitosis rates, but not mobilization rates. Other influential parameters that regulate AML model behavior are responses to cytokines/growth factors produced in peripheral blood that target the probability of self-renewal of neutrophil progenitors. Finally, our model predicts that the mitosis rate of cancer is the most predictive parameter for survival time, followed closely by parameters that affect the self-renewal of cancer stem cells; most current therapies target mitosis rate, but based on our results, we propose that additional therapeutic targeting of self-renewal of cancer stem cells will lead to even higher survival rates.

Electronic supplementary material

The online version of this article (doi:10.1186/s12918-017-0469-2) contains supplementary material, which is available to authorized users.

Sensitivity analysis of Repast computational ecology models with R/Repast

Computational ecology is an emerging interdisciplinary discipline founded mainly on modeling and simulation methods for studying ecological systems. Among the existing modeling formalisms, the individual‐based modeling is particularly well suited for capturing the complex temporal and spatial dynamics as well as the nonlinearities arising in ecosystems, communities, or populations due to individual variability. In addition, being a bottom‐up approach, it is useful for providing new insights on the local mechanisms which are generating some observed global dynamics. Of course, no conclusions about model results could be taken seriously if they are based on a single model execution and they are not analyzed carefully. Therefore, a sound methodology should always be used for underpinning the interpretation of model results. The sensitivity analysis is a methodology for quantitatively assessing the effect of input uncertainty in the simulation output which should be incorporated compulsorily to every work based on in‐silico experimental setup. In this article, we present R/Repast a GNU R package for running and analyzing Repast Simphony models accompanied by two worked examples on how to perform global sensitivity analysis and how to interpret the results.

The Integrated Role of Wnt/β-Catenin, N-Glycosylation, and E-Cadherin-Mediated Adhesion in Network Dynamics

The cellular network composed of the evolutionarily conserved metabolic pathways of protein N-glycosylation, Wnt/β-catenin signaling pathway, and E-cadherin-mediated cell-cell adhesion plays pivotal roles in determining the balance between cell proliferation and intercellular adhesion during development and in maintaining homeostasis in differentiated tissues. These pathways share a highly conserved regulatory molecule, β-catenin, which functions as both a structural component of E-cadherin junctions and as a co-transcriptional activator of the Wnt/β-catenin signaling pathway, whose target is the N-glycosylation-regulating gene, DPAGT1. Whereas these pathways have been studied independently, little is known about the dynamics of their interaction. Here we present the first numerical model of this network in MDCK cells. Since the network comprises a large number of molecules with varying cell context and time-dependent levels of expression, it can give rise to a wide range of plausible cellular states that are difficult to track. Using known kinetic parameters for individual reactions in the component pathways, we have developed a theoretical framework and gained new insights into cellular regulation of the network. Specifically, we developed a mathematical model to quantify the fold-change in concentration of any molecule included in the mathematical representation of the network in response to a simulated activation of the Wnt/ β-catenin pathway with Wnt3a under different conditions. We quantified the importance of protein N-glycosylation and synthesis of the DPAGT1 encoded enzyme, GPT, in determining the abundance of cytoplasmic β-catenin. We confirmed the role of axin in β-catenin degradation. Finally, our data suggest that cell-cell adhesion is insensitive to E-cadherin recycling in the cell. We validate the model by inhibiting β-catenin-mediated activation of DPAGT1 expression and predicting changes in cytoplasmic β-catenin concentration and stability of E-cadherin junctions in response to DPAGT1 inhibition. We show the impact of pathway dysregulation through measurements of cell migration in scratch-wound assays. Collectively, our results highlight the importance of numerical analyses of cellular networks dynamics to gain insights into physiological processes and potential design of therapeutic strategies to prevent epithelial cell invasion in cancer.

A global sensitivity analysis approach for morphogenesis models

BACKGROUND

Morphogenesis is a developmental process in which cells organize into shapes and patterns. Complex, non-linear and multi-factorial models with images as output are commonly used to study morphogenesis. It is difficult to understand the relation between the uncertainty in the input and the output of such 'black-box' models, giving rise to the need for sensitivity analysis tools. In this paper, we introduce a workflow for a global sensitivity analysis approach to study the impact of single parameters and the interactions between them on the output of morphogenesis models.

RESULTS

To demonstrate the workflow, we used a published, well-studied model of vascular morphogenesis. The parameters of this cellular Potts model (CPM) represent cell properties and behaviors that drive the mechanisms of angiogenic sprouting. The global sensitivity analysis correctly identified the dominant parameters in the model, consistent with previous studies. Additionally, the analysis provided information on the relative impact of single parameters and of interactions between them. This is very relevant because interactions of parameters impede the experimental verification of the predicted effect of single parameters. The parameter interactions, although of low impact, provided also new insights in the mechanisms of in silico sprouting. Finally, the analysis indicated that the model could be reduced by one parameter.

CONCLUSIONS

We propose global sensitivity analysis as an alternative approach to study the mechanisms of morphogenesis. Comparison of the ranking of the impact of the model parameters to knowledge derived from experimental data and from manipulation experiments can help to falsify models and to find the operand mechanisms in morphogenesis. The workflow is applicable to all 'black-box' models, including high-throughput in vitro models in which output measures are affected by a set of experimental perturbations.

Network reconstruction based on proteomic data and prior knowledge of protein connectivity using graph theory

Modeling of signal transduction pathways is instrumental for understanding cells’ function. People have been tackling modeling of signaling pathways in order to accurately represent the signaling events inside cells’ biochemical microenvironment in a way meaningful for scientists in a biological field. In this article, we propose a method to interrogate such pathways in order to produce cell-specific signaling models. We integrate available prior knowledge of protein connectivity, in a form of a Prior Knowledge Network (PKN) with phosphoproteomic data to construct predictive models of the protein connectivity of the interrogated cell type. Several computational methodologies focusing on pathways’ logic modeling using optimization formulations or machine learning algorithms have been published on this front over the past few years. Here, we introduce a light and fast approach that uses a breadth-first traversal of the graph to identify the shortest pathways and score proteins in the PKN, fitting the dependencies extracted from the experimental design. The pathways are then combined through a heuristic formulation to produce a final topology handling inconsistencies between the PKN and the experimental scenarios. Our results show that the algorithm we developed is efficient and accurate for the construction of medium and large scale signaling networks. We demonstrate the applicability of the proposed approach by interrogating a manually curated interaction graph model of EGF/TNFA stimulation against made up experimental data. To avoid the possibility of erroneous predictions, we performed a cross-validation analysis. Finally, we validate that the introduced approach generates predictive topologies, comparable to the ILP formulation. Overall, an efficient approach based on graph theory is presented herein to interrogate protein–protein interaction networks and to provide meaningful biological insights.

A Computational Study of the Effects of Syk Activity on B Cell Receptor Signaling Dynamics

The kinase Syk is intricately involved in early signaling events in B cells and is required for proper response when antigens bind to B cell receptors (BCRs). Experiments using an analog-sensitive version of Syk (Syk-AQL) have better elucidated its role, but have not completely characterized its behavior. We present a computational model for BCR signaling, using dynamical systems, which incorporates both wild-type Syk and Syk-AQL. Following the use of sensitivity analysis to identify significant reaction parameters, we screen for parameter vectors that produced graded responses to BCR stimulation as is observed experimentally. We demonstrate qualitative agreement between the model and dose response data for both mutant and wild-type kinases. Analysis of our model suggests that the level of NF-κB activation, which is reduced in Syk-AQL cells relative to wild-type, is more sensitive to small reductions in kinase activity than Erkp activation, which is essentially unchanged. Since this profile of high Erkp and reduced NF-κB is consistent with anergy, this implies that anergy is particularly sensitive to small changes in catalytic activity. Also, under a range of forward and reverse ligand binding rates, our model of Erkp and NF-κB activation displays a dependence on a power law affinity: the ratio of the forward rate to a non-unit power of the reverse rate. This dependence implies that B cells may respond to certain details of binding and unbinding rates for ligands rather than simple affinity alone.

The sensitivity and significance analysis of parameters in the model of pH regulation on lactic acid production by Lactobacillus bulgaricus

Background

The excessive production of lactic acid by L. bulgaricus during yogurt storage is a phenomenon we are always tried to prevent. The methods used in industry either control the post-acidification inefficiently or kill the probiotics in yogurt. Genetic methods of changing the activity of one enzyme related to lactic acid metabolism make the bacteria short of energy to growth, although they are efficient ways in controlling lactic acid production.

Results

A model of pH-induced promoter regulation on the production of lactic acid by L. bulgaricus was built. The modelled lactic acid metabolism without pH-induced promoter regulation fitted well with wild type L. bulgaricus (R²_LAC = 0.943, R²_LA = 0.942). Both the local sensitivity analysis and Sobol sensitivity analysis indicated parameters T_max, GR, K_LR, S, V₀, V₁ and d_LR were sensitive. In order to guide the future biology experiments, three adjustable parameters, K_LR, V₀ and V₁, were chosen for further simulations. V₀ had little effect on lactic acid production if the pH-induced promoter could be well induced when pH decreased to its threshold. K_LR and V₁ both exhibited great influence on the producing of lactic acid.

Conclusions

The proposed method of introducing a pH-induced promoter to regulate a repressor gene could restrain the synthesis of lactic acid if an appropriate strength of promoter and/or an appropriate strength of ribosome binding sequence (RBS) in lacR gene has been designed.

Identification of crucial parameters in a mathematical multiscale model of glioblastoma growth

Glioblastomas are highly malignant brain tumours. Mathematical models and their analysis provide a tool to support the understanding of the development of these tumours as well as the design of more effective treatment strategies. We have previously developed a multiscale model of glioblastoma progression that covers processes on the cellular and molecular scale. Here, we present a novel nutrient-dependent multiscale sensitivity analysis of this model that helps to identify those reaction parameters of the molecular interaction network that influence the tumour progression on the cellular scale the most. In particular, those parameters are identified that essentially determine tumour expansion and could be therefore used as potential therapy targets. As indicators for the success of a potential therapy target, a deceleration of the tumour expansion and a reduction of the tumour volume are employed. From the results, it can be concluded that no single parameter variation results in a less aggressive tumour. However, it can be shown that a few combined perturbations of two systematically selected parameters cause a slow-down of the tumour expansion velocity accompanied with a decrease of the tumour volume. Those parameters are primarily linked to the reactions that involve the microRNA-451 and the thereof regulated protein MO25.

Multiple model-informed open-loop control of uncertain intracellular signaling dynamics

Computational approaches to tune the activation of intracellular signal transduction pathways both predictably and selectively will enable researchers to explore and interrogate cell biology with unprecedented precision. Techniques to control complex nonlinear systems typically involve the application of control theory to a descriptive mathematical model. For cellular processes, however, measurement assays tend to be too time consuming for real-time feedback control and models offer rough approximations of the biological reality, thus limiting their utility when considered in isolation. We overcome these problems by combining nonlinear model predictive control with a novel adaptive weighting algorithm that blends predictions from multiple models to derive a compromise open-loop control sequence. The proposed strategy uses weight maps to inform the controller of the tendency for models to differ in their ability to accurately reproduce the system dynamics under different experimental perturbations (i.e. control inputs). These maps, which characterize the changing model likelihoods over the admissible control input space, are constructed using preexisting experimental data and used to produce a model-based open-loop control framework. In effect, the proposed method designs a sequence of control inputs that force the signaling dynamics along a predefined temporal response without measurement feedback while mitigating the effects of model uncertainty. We demonstrate this technique on the well-known Erk/MAPK signaling pathway in T cells. In silico assessment demonstrates that this approach successfully reduces target tracking error by 52% or better when compared with single model-based controllers and non-adaptive multiple model-based controllers. In vitro implementation of the proposed approach in Jurkat cells confirms a 63% reduction in tracking error when compared with the best of the single-model controllers. This study provides an experimentally-corroborated control methodology that utilizes the knowledge encoded within multiple mathematical models of intracellular signaling to design control inputs that effectively direct cell behavior in open-loop.

Specifying informative experiment stimulation conditions for resolving dynamical uncertainty in biological systems

A computationally efficient model-based design of experiments (MBDOE) strategy is developed to plan an optimal experiment by specifying the experimental stimulation magnitudes and measurement points. The strategy is extended from previous work which optimized the experimental design over a space of measurable species and time points. We include system inputs (stimulation conditions) in the experiment design search to investigate if the addition of perturbations enhances the ability of the MBDOE method to resolve uncertainties in system dynamics. The MBDOE problem is made computationally tractable by using a sparse-grid approximation of the model output dynamics, pre-specifying the time points at which the input or experimental perturbations can be applied, and creating scenario trees to explore the endogenous uncertainty. Consecutive scenario trees are used to determine the best input magnitudes and select the optimal associated measurement species and time points. We demonstrate the effectiveness of this strategy on a T-Cell Receptor (TCR) signaling pathway model.

Switch of sensitivity dynamics revealed with DyGloSA toolbox for dynamical global sensitivity analysis as an early warning for system's critical transition

Systems with bifurcations may experience abrupt irreversible and often unwanted shifts in their performance, called critical transitions. For many systems like climate, economy, ecosystems it is highly desirable to identify indicators serving as early warnings of such regime shifts. Several statistical measures were recently proposed as early warnings of critical transitions including increased variance, autocorrelation and skewness of experimental or model-generated data. The lack of automatized tool for model-based prediction of critical transitions led to designing DyGloSA - a MATLAB toolbox for dynamical global parameter sensitivity analysis (GPSA) of ordinary differential equations models. We suggest that the switch in dynamics of parameter sensitivities revealed by our toolbox is an early warning that a system is approaching a critical transition. We illustrate the efficiency of our toolbox by analyzing several models with bifurcations and predicting the time periods when systems can still avoid going to a critical transition by manipulating certain parameter values, which is not detectable with the existing SA techniques. DyGloSA is based on the SBToolbox2 and contains functions, which compute dynamically the global sensitivity indices of the system by applying four main GPSA methods: eFAST, Sobol's ANOVA, PRCC and WALS. It includes parallelized versions of the functions enabling significant reduction of the computational time (up to 12 times). DyGloSA is freely available as a set of MATLAB scripts at http://bio.uni.lu/systems_biology/software/dyglosa. It requires installation of MATLAB (versions R2008b or later) and the Systems Biology Toolbox2 available at www.sbtoolbox2.org. DyGloSA can be run on Windows and Linux systems, -32 and -64 bits.

What can we learn from global sensitivity analysis of biochemical systems?

Most biological models of intermediate size, and probably all large models, need to cope with the fact that many of their parameter values are unknown. In addition, it may not be possible to identify these values unambiguously on the basis of experimental data. This poses the question how reliable predictions made using such models are. Sensitivity analysis is commonly used to measure the impact of each model parameter on its variables. However, the results of such analyses can be dependent on an exact set of parameter values due to nonlinearity. To mitigate this problem, global sensitivity analysis techniques are used to calculate parameter sensitivities in a wider parameter space. We applied global sensitivity analysis to a selection of five signalling and metabolic models, several of which incorporate experimentally well-determined parameters. Assuming these models represent physiological reality, we explored how the results could change under increasing amounts of parameter uncertainty. Our results show that parameter sensitivities calculated with the physiological parameter values are not necessarily the most frequently observed under random sampling, even in a small interval around the physiological values. Often multimodal distributions were observed. Unsurprisingly, the range of possible sensitivity coefficient values increased with the level of parameter uncertainty, though the amount of parameter uncertainty at which the pattern of control was able to change differed among the models analysed. We suggest that this level of uncertainty can be used as a global measure of model robustness. Finally a comparison of different global sensitivity analysis techniques shows that, if high-throughput computing resources are available, then random sampling may actually be the most suitable technique.

Computational approaches for understanding energy metabolism

There has been a surge of interest in understanding the regulation of metabolic networks involved in disease in recent years. Quantitative models are increasingly being used to interrogate the metabolic pathways that are contained within this complex disease biology. At the core of this effort is the mathematical modeling of central carbon metabolism involving glycolysis and the citric acid cycle (referred to as energy metabolism). Here, we discuss several approaches used to quantitatively model metabolic pathways relating to energy metabolism and discuss their formalisms, successes, and limitations.

Detecting and removing inconsistencies between experimental data and signaling network topologies using integer linear programming on interaction graphs

Cross-referencing experimental data with our current knowledge of signaling network topologies is one central goal of mathematical modeling of cellular signal transduction networks. We present a new methodology for data-driven interrogation and training of signaling networks. While most published methods for signaling network inference operate on Bayesian, Boolean, or ODE models, our approach uses integer linear programming (ILP) on interaction graphs to encode constraints on the qualitative behavior of the nodes. These constraints are posed by the network topology and their formulation as ILP allows us to predict the possible qualitative changes (up, down, no effect) of the activation levels of the nodes for a given stimulus. We provide four basic operations to detect and remove inconsistencies between measurements and predicted behavior: (i) find a topology-consistent explanation for responses of signaling nodes measured in a stimulus-response experiment (if none exists, find the closest explanation); (ii) determine a minimal set of nodes that need to be corrected to make an inconsistent scenario consistent; (iii) determine the optimal subgraph of the given network topology which can best reflect measurements from a set of experimental scenarios; (iv) find possibly missing edges that would improve the consistency of the graph with respect to a set of experimental scenarios the most. We demonstrate the applicability of the proposed approach by interrogating a manually curated interaction graph model of EGFR/ErbB signaling against a library of high-throughput phosphoproteomic data measured in primary hepatocytes. Our methods detect interactions that are likely to be inactive in hepatocytes and provide suggestions for new interactions that, if included, would significantly improve the goodness of fit. Our framework is highly flexible and the underlying model requires only easily accessible biological knowledge. All related algorithms were implemented in a freely available toolbox SigNetTrainer making it an appealing approach for various applications.

Cellular signal transduction is orchestrated by communication networks of signaling proteins commonly depicted on signaling pathway maps. However, each cell type may have distinct variants of signaling pathways, and wiring diagrams are often altered in disease states. The identification of truly active signaling topologies based on experimental data is therefore one key challenge in systems biology of cellular signaling. We present a new framework for training signaling networks based on interaction graphs (IG). In contrast to complex modeling formalisms, IG capture merely the known positive and negative edges between the components. This basic information, however, already sets hard constraints on the possible qualitative behaviors of the nodes when perturbing the network. Our approach uses Integer Linear Programming to encode these constraints and to predict the possible changes (down, neutral, up) of the activation levels of the involved players for a given experiment. Based on this formulation we developed several algorithms for detecting and removing inconsistencies between measurements and network topology. Demonstrated by EGFR/ErbB signaling in hepatocytes, our approach delivers direct conclusions on edges that are likely inactive or missing relative to canonical pathway maps. Such information drives the further elucidation of signaling network topologies under normal and pathological phenotypes.

Analysis and design of a genetic circuit for dynamic metabolic engineering

Recent advances in synthetic biology have equipped us with new tools for bioprocess optimization at the genetic level. Previously, we have presented an integrated in silico design for the dynamic control of gene expression based on a density-sensing unit and a genetic toggle switch. In the present paper, analysis of a serine-producing Escherichia coli mutant shows that an instantaneous ON-OFF switch leads to a maximum theoretical productivity improvement of 29.6% compared to the mutant. To further the design, global sensitivity analysis is applied here to a mathematical model of serine production in E. coli coupled with a genetic circuit. The model of the quorum sensing and the toggle switch involves 13 parameters of which 3 are identified as having a significant effect on serine concentration. Simulations conducted in this reduced parameter space further identified the optimal ranges for these 3 key parameters to achieve productivity values close to the maximum theoretical values. This analysis can now be used to guide the experimental implementation of a dynamic metabolic engineering strategy and reduce the time required to design the genetic circuit components.

Systematic analysis of the mechanisms of virus-triggered type I IFN signaling pathways through mathematical modeling

Based on biological experimental data, we developed a mathematical model of the virus-triggered signaling pathways that lead to induction of type I IFNs and systematically analyzed the mechanisms of the cellular antiviral innate immune responses, including the negative feedback regulation of ISG56 and the positive feedback regulation of IFNs. We found that the time between 5 and 48 hours after viral infection is vital for the control and/or elimination of the virus from the host cells and demonstrated that the ISG56-induced inhibition of MITA activation is stronger than the ISG56-induced inhibition of TBK1 activation. The global parameter sensitivity analysis suggests that the positive feedback regulation of IFNs is very important in the innate antiviral system. Furthermore, the robustness of the innate immune signaling network was demonstrated using a new robustness index. These results can help us understand the mechanisms of the virus-induced innate immune response at a system level and provide instruction for further biological experiments.

Parameter uncertainty in biochemical models described by ordinary differential equations

Improved mechanistic understanding of biochemical networks is one of the driving ambitions of Systems Biology. Computational modeling allows the integration of various sources of experimental data in order to put this conceptual understanding to the test in a quantitative manner. The aim of computational modeling is to obtain both predictive as well as explanatory models for complex phenomena, hereby providing useful approximations of reality with varying levels of detail. As the complexity required to describe different system increases, so does the need for determining how well such predictions can be made. Despite efforts to make tools for uncertainty analysis available to the field, these methods have not yet found widespread use in the field of Systems Biology. Additionally, the suitability of the different methods strongly depends on the problem and system under investigation. This review provides an introduction to some of the techniques available as well as gives an overview of the state-of-the-art methods for parameter uncertainty analysis.

The A, B, Cs of G-protein-coupled receptor pharmacology in assay development for HTS

G-protein-coupled receptors represent one of the most important areas of research in the pharmaceutical industry, being one of the largest druggable gene families. Recognising this fact, manufacturers have developed a huge variety of homogeneous assay technologies that facilitate the quantification of receptor ligand binding events and their downstream signalling cascades. However, while early emphasis was placed on the most sensitive, high-throughput and cost-effective screening technologies to enable identification of the most lead matter for further development, in recent years emphasis has shifted to a focus on maximising the identification of compounds that are new and developing assays that are more biologically/pharmacologically relevant. Therefore, this review provides an overview of the binding and functional techniques available for high-throughput screening, with particular attention on how assay application and configuration can be maximised to ensure their successful identification of relevant chemical matter and thereby optimising project success.

Non Linear Programming (NLP) formulation for quantitative modeling of protein signal transduction pathways

Modeling of signal transduction pathways plays a major role in understanding cells' function and predicting cellular response. Mathematical formalisms based on a logic formalism are relatively simple but can describe how signals propagate from one protein to the next and have led to the construction of models that simulate the cells response to environmental or other perturbations. Constrained fuzzy logic was recently introduced to train models to cell specific data to result in quantitative pathway models of the specific cellular behavior. There are two major issues in this pathway optimization: i) excessive CPU time requirements and ii) loosely constrained optimization problem due to lack of data with respect to large signaling pathways. Herein, we address both issues: the former by reformulating the pathway optimization as a regular nonlinear optimization problem; and the latter by enhanced algorithms to pre/post-process the signaling network to remove parts that cannot be identified given the experimental conditions. As a case study, we tackle the construction of cell type specific pathways in normal and transformed hepatocytes using medium and large-scale functional phosphoproteomic datasets. The proposed Non Linear Programming (NLP) formulation allows for fast optimization of signaling topologies by combining the versatile nature of logic modeling with state of the art optimization algorithms.

Evaluation of design strategies for time course experiments in genetic networks: case study of the XlnR regulon in Aspergillus niger

One of the challenges in genetic network reconstruction is finding experimental designs that maximize the information content in a data set. In this paper, the information value of mRNA transcription time course experiments was used to compare experimental designs. The study concerns the dynamic response of genes in the XlnR regulon of Aspergillus niger, with the goal to find the best moment in time to administer an extra pulse of inducing D-xylose. Low and high D-xylose pulses were used to perturb the XlnR regulon. Evaluation of the experimental methods was based on simulation of the regulon. Models that govern the regulation of the target genes in this regulon were used for the simulations. Parameter sensitivity analysis, the Fisher Information Matrix (FIM) and the modified E-criterion were used to assess the design performances. The results show that the best time to give a second D-xylose pulse is when the D-xylose concentration from the first pulse has not yet completely faded away. Due to the presence of a repression effect the strength of the second pulse must be optimized, rather than maximized. The results suggest that the modified E-criterion is a better metric than the sum of integrals of absolute sensitivity for comparing alternative designs.

Modeling and experimental analyses reveals signaling plasticity in a bi-modular assembly of CD40 receptor activated kinases

Depending on the strength of signal dose, CD40 receptor (CD40) controls ERK-1/2 and p38MAPK activation. At low signal dose, ERK-1/2 is maximally phosphorylated but p38MAPK is minimally phosphorylated; as the signal dose increases, ERK-1/2 phosphorylation is reduced whereas p38MAPK phosphorylation is reciprocally enhanced. The mechanism of reciprocal activation of these two MAPKs remains un-elucidated. Here, our computational model, coupled to experimental perturbations, shows that the observed reciprocity is a system-level behavior of an assembly of kinases arranged in two modules. Experimental perturbations with kinase inhibitors suggest that a minimum of two trans-modular negative feedback loops are required to reproduce the experimentally observed reciprocity. The bi-modular architecture of the signaling pathways endows the system with an inherent plasticity which is further expressed in the skewing of the CD40-induced productions of IL-10 and IL-12, the respective anti-inflammatory and pro-inflammatory cytokines. Targeting the plasticity of CD40 signaling significantly reduces Leishmania major infection in a susceptible mouse strain. Thus, for the first time, using CD40 signaling as a model, we show how a bi-modular assembly of kinases imposes reciprocity to a receptor signaling. The findings unravel that the signalling plasticity is inherent to a reciprocal system and that the principle can be used for designing a therapy.

Modeling and analysis of biopathways dynamics

Cellular processes are governed and coordinated by a multitude of biopathways. A pathway can be viewed as a complex network of biochemical reactions. The dynamics of this network largely determines the functioning of the pathway. Hence the modeling and analysis of biochemical networks dynamics is an important problem and is an active area of research. Here we review quantitative models of biochemical networks based on ordinary differential equations (ODEs). We mainly focus on the parameter estimation and sensitivity analysis problems and survey the current methods for tackling them. In this context we also highlight a recently developed probabilistic approximation technique using which these two problems can be considerably simplified.

Inferring robust gene networks from expression data by a sensitivity-based incremental evolution method

Background

Reconstructing gene regulatory networks (GRNs) from expression data is one of the most important challenges in systems biology research. Many computational models and methods have been proposed to automate the process of network reconstruction. Inferring robust networks with desired behaviours remains challenging, however. This problem is related to network dynamics but has yet to be investigated using network modeling.

Results

We propose an incremental evolution approach for inferring GRNs that takes network robustness into consideration and can deal with a large number of network parameters. Our approach includes a sensitivity analysis procedure to iteratively select the most influential network parameters, and it uses a swarm intelligence procedure to perform parameter optimization. We have conducted a series of experiments to evaluate the external behaviors and internal robustness of the networks inferred by the proposed approach. The results and analyses have verified the effectiveness of our approach.

Conclusions

Sensitivity analysis is crucial to identifying the most sensitive parameters that govern the network dynamics. It can further be used to derive constraints for network parameters in the network reconstruction process. The experimental results show that the proposed approach can successfully infer robust GRNs with desired system behaviors.

Sensitivity analysis approaches applied to systems biology models

With the rising application of systems biology, sensitivity analysis methods have been widely applied to study the biological systems, including metabolic networks, signalling pathways and genetic circuits. Sensitivity analysis can provide valuable insights about how robust the biological responses are with respect to the changes of biological parameters and which model inputs are the key factors that affect the model outputs. In addition, sensitivity analysis is valuable for guiding experimental analysis, model reduction and parameter estimation. Local and global sensitivity analysis approaches are the two types of sensitivity analysis that are commonly applied in systems biology. Local sensitivity analysis is a classic method that studies the impact of small perturbations on the model outputs. On the other hand, global sensitivity analysis approaches have been applied to understand how the model outputs are affected by large variations of the model input parameters. In this review, the author introduces the basic concepts of sensitivity analysis approaches applied to systems biology models. Moreover, the author discusses the advantages and disadvantages of different sensitivity analysis methods, how to choose a proper sensitivity analysis approach, the available sensitivity analysis tools for systems biology models and the caveats in the interpretation of sensitivity analysis results.

A methodology for global-sensitivity analysis of time-dependent outputs in systems biology modelling

One of the main challenges in the development of mathematical and computational models of biological systems is the precise estimation of parameter values. Understanding the effects of uncertainties in parameter values on model behaviour is crucial to the successful use of these models. Global sensitivity analysis (SA) can be used to quantify the variability in model predictions resulting from the uncertainty in multiple parameters and to shed light on the biological mechanisms driving system behaviour. We present a new methodology for global SA in systems biology which is computationally efficient and can be used to identify the key parameters and their interactions which drive the dynamic behaviour of a complex biological model. The approach combines functional principal component analysis with established global SA techniques. The methodology is applied to a model of the insulin signalling pathway, defects of which are a major cause of type 2 diabetes and a number of key features of the system are identified.

Model-based global sensitivity analysis as applied to identification of anti-cancer drug targets and biomarkers of drug resistance in the ErbB2/3 network

High levels of variability in cancer-related cellular signalling networks and a lack of parameter identifiability in large-scale network models hamper translation of the results of modelling studies into the process of anti-cancer drug development. Recently global sensitivity analysis (GSA) has been recognised as a useful technique, capable of addressing the uncertainty of the model parameters and generating valid predictions on parametric sensitivities. Here we propose a novel implementation of model-based GSA specially designed to explore how multi-parametric network perturbations affect signal propagation through cancer-related networks. We use area-under-the-curve for time course of changes in phosphorylation of proteins as a characteristic for sensitivity analysis and rank network parameters with regard to their impact on the level of key cancer-related outputs, separating strong inhibitory from stimulatory effects. This allows interpretation of the results in terms which can incorporate the effects of potential anti-cancer drugs on targets and the associated biological markers of cancer. To illustrate the method we applied it to an ErbB signalling network model and explored the sensitivity profile of its key model readout, phosphorylated Akt, in the absence and presence of the ErbB2 inhibitor pertuzumab. The method successfully identified the parameters associated with elevation or suppression of Akt phosphorylation in the ErbB2/3 network. From analysis and comparison of the sensitivity profiles of pAkt in the absence and presence of targeted drugs we derived predictions of drug targets, cancer-related biomarkers and generated hypotheses for combinatorial therapy. Several key predictions have been confirmed in experiments using human ovarian carcinoma cell lines. We also compared GSA-derived predictions with the results of local sensitivity analysis and discuss the applicability of both methods. We propose that the developed GSA procedure can serve as a refining tool in combinatorial anti-cancer drug discovery.

Reducing experimental variability in variance-based sensitivity analysis of biochemical reaction systems

Sensitivity analysis is a valuable task for assessing the effects of biological variability on cellular behavior. Available techniques require knowledge of nominal parameter values, which cannot be determined accurately due to experimental uncertainty typical to problems of systems biology. As a consequence, the practical use of existing sensitivity analysis techniques may be seriously hampered by the effects of unpredictable experimental variability. To address this problem, we propose here a probabilistic approach to sensitivity analysis of biochemical reaction systems that explicitly models experimental variability and effectively reduces the impact of this type of uncertainty on the results. The proposed approach employs a recently introduced variance-based method to sensitivity analysis of biochemical reaction systems [Zhang et al., J. Chem. Phys. 134, 094101 (2009)] and leads to a technique that can be effectively used to accommodate appreciable levels of experimental variability. We discuss three numerical techniques for evaluating the sensitivity indices associated with the new method, which include Monte Carlo estimation, derivative approximation, and dimensionality reduction based on orthonormal Hermite approximation. By employing a computational model of the epidermal growth factor receptor signaling pathway, we demonstrate that the proposed technique can greatly reduce the effect of experimental variability on variance-based sensitivity analysis results. We expect that, in cases of appreciable experimental variability, the new method can lead to substantial improvements over existing sensitivity analysis techniques.

Identification of Critical Molecular Components in a Multiscale Cancer Model Based on the Integration of Monte Carlo, Resampling, and ANOVA

To date, parameters defining biological properties in multiscale disease models are commonly obtained from a variety of sources. It is thus important to examine the influence of parameter perturbations on system behavior, rather than to limit the model to a specific set of parameters. Such sensitivity analysis can be used to investigate how changes in input parameters affect model outputs. However, multiscale cancer models require special attention because they generally take longer to run than does a series of signaling pathway analysis tasks. In this article, we propose a global sensitivity analysis method based on the integration of Monte Carlo, resampling, and analysis of variance. This method provides solutions to (1) how to render the large number of parameter variation combinations computationally manageable, and (2) how to effectively quantify the sampling distribution of the sensitivity index to address the inherent computational intensity issue. We exemplify the feasibility of this method using a two-dimensional molecular-microscopic agent-based model previously developed for simulating non-small cell lung cancer; in this model, an epidermal growth factor (EGF)-induced, EGF receptor-mediated signaling pathway was implemented at the molecular level. Here, the cross-scale effects of molecular parameters on two tumor growth evaluation measures, i.e., tumor volume and expansion rate, at the microscopic level are assessed. Analysis finds that ERK, a downstream molecule of the EGF receptor signaling pathway, has the most important impact on regulating both measures. The potential to apply this method to therapeutic target discovery is discussed.

Identifying therapeutic targets in a combined EGFR-TGFβR signalling cascade using a multiscale agent-based cancer model

Applying a previously developed non-small cell lung cancer model, we assess 'cross-scale' the therapeutic efficacy of targeting a variety of molecular components of the epidermal growth factor receptor (EGFR) signalling pathway. Simulation of therapeutic inhibition and amplification allows for the ranking of the implemented downstream EGFR signalling molecules according to their therapeutic values or indices. Analysis identifies mitogen-activated protein kinase and extracellular signal-regulated kinase as top therapeutic targets for both inhibition and amplification-based treatment regimen but indicates that combined parameter perturbations do not necessarily improve the therapeutic effect of the separate parameter treatments as much as might be expected. Potential future strategies using this in silico model to tailor molecular treatment regimen are discussed.

Experiment design through dynamical characterisation of non-linear systems biology models utilising sparse grids

The sparse grid-based experiment design algorithm sequentially selects an experimental design point to discriminate between hypotheses for given experimental conditions. Sparse grids efficiently screen the global uncertain parameter space to identify acceptable parameter subspaces. Clustering the located acceptable parameter vectors by the similarity of the simulated model trajectories characterises the data-compatible model dynamics. The experiment design algorithm capitalizes on the diversity of the experimentally distinguishable system output dynamics to select the design point that best discerns between competing model-structure and parameter-encoded hypotheses. As opposed to designing the experiments to explicitly reduce uncertainty in the model parameters, this approach selects design points to differentiate between dynamical behaviours. This approach further differs from other experimental design methods in that it simultaneously addresses both parameter- and structural-based uncertainty that is applicable to some ill-posed problems where the number of uncertain parameters exceeds the amount of data, places very few requirements on the model type, available data and a priori parameter estimates, and is performed over the global uncertain parameter space. The experiment design algorithm is demonstrated on a mitogen-activated protein kinase cascade model. The results show that system dynamics are highly uncertain with limited experimental data. Nevertheless, the algorithm requires only three additional experimental data points to simultaneously discriminate between possible model structures and acceptable parameter values. This sparse grid-based experiment design process provides a systematic and computationally efficient exploration over the entire uncertain parameter space of potential model structures to resolve the uncertainty in the non-linear systems biology model dynamics.

Systematic calibration of a cell signaling network model

Background

Mathematical modeling is being applied to increasingly complex biological systems and datasets; however, the process of analyzing and calibrating against experimental data is often challenging and a rate limiting step in model development. To address this problem, we developed a systematic methodology for calibrating quantitative models of dynamic biological processes and illustrate its utility by validating a model of TRAIL (Tumor necrosis factor Related Apoptosis-Inducing Ligand)-induced cell death.

Results

We propose a serial framework integrating analysis and calibration modules and we compare various methods for global sensitivity analysis and global parameter estimation. First, adequacy of the network structure is checked by global sensitivity analysis to changes in concentrations of molecular species, validating that the model can reproduce qualitative features of the system behavior derived from experiments or literature surveys. Second, rate parameters are ranked by importance using gradient-based and variance-based sensitivity indices, and we systematically determine the optimal number of parameters to include in model calibration. Third, deterministic, stochastic and hybrid algorithms for global optimization are applied to estimate the values of the most important parameters by fitting to time series data. We compare the performance of these three optimization algorithms.

Conclusions

Our proposed framework covers the entire process from validating a proto-model to establishing a realistic model for in silico experiments and thereby provides a generalized workflow for the construction of predictive models of complex network systems.