Fast integration-based prediction bands for ordinary differential equation models.

MOTIVATION: To gain a deeper understanding of biological processes and their relevance in disease, mathematical models are built upon experimental data. Uncertainty in the data leads to uncertainties of the model's parameters and in turn to uncertainties of predictions. Mechanistic dynamic models of biochemical networks are frequently based on nonlinear differential equation systems and feature a large number of parameters, sparse observations of the model components and lack of information in the available data. Due to the curse of dimensionality, classical and sampling approaches propagating parameter uncertainties to predictions are hardly feasible and insufficient. However, for experimental design and to discriminate between competing models, prediction and confidence bands are essential. To circumvent the hurdles of the former methods, an approach to calculate a profile likelihood on arbitrary observations for a specific time point has been introduced, which provides accurate confidence and prediction intervals for nonlinear models and is computationally feasible for high-dimensional models. RESULTS: In this article, reliable and smooth point-wise prediction and confidence bands to assess the model's uncertainty on the whole time-course are achieved via explicit integration with elaborate correction mechanisms. The corresponding system of ordinary differential equations is derived and tested on three established models for cellular signalling. An efficiency analysis is performed to illustrate the computational benefit compared with repeated profile likelihood calculations at multiple time points. AVAILABILITY AND IMPLEMENTATION: The integration framework and the examples used in this article are provided with the software package Data2Dynamics, which is based on MATLAB and freely available at http://www.data2dynamics.org CONTACT: helge.hass@fdm.uni-freiburg.de SUPPLEMENTARY INFORMATION: Supplementary data are available at Bioinformatics online.

Population Pharmacokinetic Modeling of Adavosertib (AZD1775) in Patients with Solid Tumors.

Adavosertib (AZD1775) is a potent small-molecule inhibitor of Wee1 kinase. This analysis utilized pharmacokinetic data from 8 Phase I/II studies of adavosertib to characterize the population pharmacokinetics of adavosertib in patients (n = 538) with solid tumors and evaluate the impact of covariates on exposure. A nonlinear mixed-effects modeling approach was employed to estimate population and individual parameters from the clinical trial data. The model for time dependency of apparent clearance (CL) was developed in a stepwise manner and the final model validated by visual predictive checks (VPCs). Using an adavosertib dose of 300 mg once daily on a 5 days on/2 days off dosing schedule given 2 weeks out of a 3-week cycle, simulation analyses evaluated the impact of covariates on the following exposure metrics at steady state: maximum concentration during a 21-day cycle, area under the curve (AUC) during a 21-day cycle, AUC during the second week of a treatment cycle, and AUC on day 12 of a treatment cycle. The final model was a linear 2-compartment model with lag time into the dosing compartment and first-order absorption into the central compartment, time-varying CL, and random effects on all model parameters. VPCs and steady-state observations confirmed that the final model satisfied all the requirements for reliable simulation of randomly sampled Phase I and II populations with different covariate characteristics. Simulation-based analyses revealed that body weight, renal impairment status, and race were key factors determining the variability of drug-exposure metrics.

Division of labor by dual feedback regulators controls JAK2/STAT5 signaling over broad ligand range.

Cellular signal transduction is governed by multiple feedback mechanisms to elicit robust cellular decisions. The specific contributions of individual feedback regulators, however, remain unclear. Based on extensive time-resolved data sets in primary erythroid progenitor cells, we established a dynamic pathway model to dissect the roles of the two transcriptional negative feedback regulators of the suppressor of cytokine signaling (SOCS) family, CIS and SOCS3, in JAK2/STAT5 signaling. Facilitated by the model, we calculated the STAT5 response for experimentally unobservable Epo concentrations and provide a quantitative link between cell survival and the integrated response of STAT5 in the nucleus. Model predictions show that the two feedbacks CIS and SOCS3 are most effective at different ligand concentration ranges due to their distinct inhibitory mechanisms. This divided function of dual feedback regulation enables control of STAT5 responses for Epo concentrations that can vary 1000-fold in vivo. Our modeling approach reveals dose-dependent feedback control as key property to regulate STAT5-mediated survival decisions over a broad range of ligand concentrations.

BlotIt-Optimal alignment of Western blot and qPCR experiments.

Biological systems are frequently analyzed by means of mechanistic mathematical models. In order to infer model parameters and provide a useful model that can be employed for systems understanding and hypothesis testing, the model is often calibrated on quantitative, time-resolved data. To do so, it is typically important to compare experimental measurements over broad time ranges and various experimental conditions, e.g. perturbations of the biological system. However, most of the established experimental techniques such as Western blot, or quantitative real-time polymerase chain reaction only provide measurements on a relative scale, since different sample volumes, experimental adjustments or varying development times of a gel lead to systematic shifts in the data. In turn, the number of measurements corresponding to the same scale enabling comparability is limited. Here, we present a new flexible method to align measurement data that obeys different scaling factors and compare it to existing normalization approaches. We propose an alignment model to estimate these scaling factors and provide the possibility to adapt this model depending on the measurement technique of interest. In addition, an error model can be specified to adequately weight the different data points and obtain scaling-model based confidence intervals of the finally scaled data points. Our approach is applicable to all sorts of relative measurements and does not need a particular experimental condition that has been measured over all available scales. An implementation of the method is provided with the R package blotIt including refined ways of visualization.

Efficient simulation of clinical target response surfaces.

Simulation of combination therapies is challenging due to computational complexity. Either a simple model is used to simulate the response for many combinations of concentration to generate a response surface but parameter variability and uncertainty are neglected and the concentrations are constant-the link to the doses to be administered is difficult to make-or a population pharmacokinetic/pharmacodynamic model is used to predict the response to combination therapy in a clinical trial taking into account the time-varying concentration profile, interindividual variability (IIV), and parameter uncertainty but simulations are limited to only a few selected doses. We devised new algorithms to efficiently search for the combination doses that achieve a predefined efficacy target while taking into account the IIV and parameter uncertainty. The result of this method is a response surface of confidence levels, indicating for all dose combinations the likelihood of reaching the specified efficacy target. We highlight the importance to simulate across a population rather than focus on an individual. Finally, we provide examples of potential applications, such as informing experimental design.

Dynamic modeling of Nrf2 pathway activation in liver cells after toxicant exposure.

Cells are exposed to oxidative stress and reactive metabolites every day. The Nrf2 signaling pathway responds to oxidative stress by upregulation of antioxidants like glutathione (GSH) to compensate the stress insult and re-establish homeostasis. Although mechanisms describing the interaction between the key pathway constituents Nrf2, Keap1 and p62 are widely reviewed and discussed in literature, quantitative dynamic models bringing together these mechanisms with time-resolved data are limited. Here, we present an ordinary differential equation (ODE) based dynamic model to describe the dynamic response of Nrf2, Keap1, Srxn1 and GSH to oxidative stress caused by the soft-electrophile diethyl maleate (DEM). The time-resolved data obtained by single-cell confocal microscopy of green fluorescent protein (GFP) reporters and qPCR of the Nrf2 pathway components complemented with siRNA knock down experiments, is accurately described by the calibrated mathematical model. We show that the quantitative model can describe the activation of the Nrf2 pathway by compounds with a different mechanism of activation, including drugs which are known for their ability to cause drug induced liver-injury (DILI) i.e., diclofenac (DCF) and omeprazole (OMZ). Finally, we show that our model can reveal differences in the processes leading to altered activation dynamics amongst DILI inducing drugs.

Local Riemannian geometry of model manifolds and its implications for practical parameter identifiability.

When non-linear models are fitted to experimental data, parameter estimates can be poorly constrained albeit being identifiable in principle. This means that along certain paths in parameter space, the log-likelihood does not exceed a given statistical threshold but remains bounded. This situation, denoted as practical non-identifiability, can be detected by Monte Carlo sampling or by systematic scanning using the profile likelihood method. In contrast, any method based on a Taylor expansion of the log-likelihood around the optimum, e.g., parameter uncertainty estimation by the Fisher Information Matrix, reveals no information about the boundedness at all. In this work, we present a geometric approach, approximating the original log-likelihood by geodesic coordinates of the model manifold. The Christoffel Symbols in the geodesic equation are fixed to those obtained from second order model sensitivities at the optimum. Based on three exemplary non-linear models we show that the information about the log-likelihood bounds and flat parameter directions can already be contained in this local information. Whereas the unbounded case represented by the Fisher Information Matrix is embedded in the geometric framework as vanishing Christoffel Symbols, non-vanishing constant Christoffel Symbols prove to define prototype non-linear models featuring boundedness and flat parameter directions of the log-likelihood. Finally, we investigate if those models could allow to approximate and replace computationally expensive objective functions originating from non-linear models by a surrogate objective function in parameter estimation problems.

Mathematical modeling of drug-induced receptor internalization in the HER2-positive SKBR3 breast cancer cell-line.

About 20% of breast cancer tumors over-express the HER2 receptor. Trastuzumab, an approved drug to treat this type of breast cancer, is a monoclonal antibody directly binding at the HER2 receptor and ultimately inhibiting cancer cell growth. The goal of our study was to understand the early impact of trastuzumab on HER2 internalization and recycling in the HER2-overexpressing breast cancer cell line SKBR3. To this end, fluorescence microscopy, monitoring the amount of HER2 expression in the plasma membrane, was combined with mathematical modeling to derive the flux of HER2 receptors from and to the membrane. We constructed a dynamic multi-compartment model based on ordinary differential equations. To account for cancer cell heterogeneity, a first, dynamic model was expanded to a second model including two distinct cell phenotypes, with implications for different conformational states of HER2, i.e. monomeric or homodimeric. Our mathematical model shows that the hypothesis of fast constitutive HER2 recycling back to the plasma membrane does not match the experimental data. It conclusively describes the experimental observation that trastuzumab induces sustained receptor internalization in cells with membrane ruffles. It is also concluded that for rare, non-ruffled (flat) cells, HER2 internalization occurs three orders of magnitude slower than for the bulk, ruffled cell population.

Quantitative Systems Pharmacology: An Exemplar Model-Building Workflow With Applications in Cardiovascular, Metabolic, and Oncology Drug Development.

Quantitative systems pharmacology (QSP), a mechanistically oriented form of drug and disease modeling, seeks to address a diverse set of problems in the discovery and development of therapies. These problems bring a considerable amount of variability and uncertainty inherent in the nonclinical and clinical data. Likewise, the available modeling techniques and related software tools are manifold. Appropriately, the development, qualification, application, and impact of QSP models have been similarly varied. In this review, we describe the progressive maturation of a QSP modeling workflow: a necessary step for the efficient, reproducible development and qualification of QSP models, which themselves are highly iterative and evolutive. Furthermore, we describe three applications of QSP to impact drug development; one supporting new indications for an approved antidiabetic clinical asset through mechanistic hypothesis generation, one highlighting efficacy and safety differentiation within the sodium-glucose cotransporter-2 inhibitor drug class, and one enabling rational selection of immuno-oncology drug combinations.

A Dynamic Mathematical Model of Bile Acid Clearance in HepaRG Cells.

A dynamic model based on ordinary differential equations that describes uptake, basolateral and canalicular export of taurocholic acid (TCA) in human HepaRG cells is presented. The highly reproducible inter-assay experimental data were used to reliably estimate model parameters. Primary human hepatocytes were similarly evaluated to establish a mathematical model, but with notably higher inter-assay differences in TCA clearance and bile canaliculi dynamics. By use of the HepaRG cell line, the simultaneous TCA clearance associated to basolateral uptake, canalicular and sinusoidal efflux, was predicted. The mathematical model accurately reproduced the dose-dependent inhibition of TCA clearance in the presence and absence of the prototypical cholestatic drugs cyclosporine A (CsA) and chlorpromazine. Rapid inhibition of TCA clearance and recovery were found to be major characteristics of CsA. Conversely, the action of chlorpromazine was described by slow onset of inhibition relative to inhibition of TCA clearance by CsA. The established mathematical model, validated by the use of these 2 prototypical cholestatic drugs and the integration of bile canalicular dynamics, provides an important development for the further study of human hepatobiliary function, through simultaneous temporal and vectorial membrane transport of bile acids in drug-induced cholestasis.

Model-based identification of TNFα-induced IKKß-mediated and IκBα-mediated regulation of NFκB signal transduction as a tool to quantify the impact of drug-induced liver injury compounds.

Drug-induced liver injury (DILI) has become a major problem for patients and for clinicians, academics and the pharmaceutical industry. To date, existing hepatotoxicity test systems are only poorly predictive and the underlying mechanisms are still unclear. One of the factors known to amplify hepatotoxicity is the tumor necrosis factor alpha (TNFα), especially due to its synergy with commonly used drugs such as diclofenac. However, the exact mechanism of how diclofenac in combination with TNFα induces liver injury remains elusive. Here, we combined time-resolved immunoblotting and live-cell imaging data of HepG2 cells and primary human hepatocytes (PHH) with dynamic pathway modeling using ordinary differential equations (ODEs) to describe the complex structure of TNFα-induced NFκB signal transduction and integrated the perturbations of the pathway caused by diclofenac. The resulting mathematical model was used to systematically identify parameters affected by diclofenac. These analyses showed that more than one regulatory module of TNFα-induced NFκB signal transduction is affected by diclofenac, suggesting that hepatotoxicity is the integrated consequence of multiple changes in hepatocytes and that multiple factors define toxicity thresholds. Applying our mathematical modeling approach to other DILI-causing compounds representing different putative DILI mechanism classes enabled us to quantify their impact on pathway activation, highlighting the potential of the dynamic pathway model as a quantitative tool for the analysis of DILI compounds.

Customized Steady-State Constraints for Parameter Estimation in Non-Linear Ordinary Differential Equation Models.

Ordinary differential equation models have become a wide-spread approach to analyze dynamical systems and understand underlying mechanisms. Model parameters are often unknown and have to be estimated from experimental data, e.g., by maximum-likelihood estimation. In particular, models of biological systems contain a large number of parameters. To reduce the dimensionality of the parameter space, steady-state information is incorporated in the parameter estimation process. For non-linear models, analytical steady-state calculation typically leads to higher-order polynomial equations for which no closed-form solutions can be obtained. This can be circumvented by solving the steady-state equations for kinetic parameters, which results in a linear equation system with comparatively simple solutions. At the same time multiplicity of steady-state solutions is avoided, which otherwise is problematic for optimization. When solved for kinetic parameters, however, steady-state constraints tend to become negative for particular model specifications, thus, generating new types of optimization problems. Here, we present an algorithm based on graph theory that derives non-negative, analytical steady-state expressions by stepwise removal of cyclic dependencies between dynamical variables. The algorithm avoids multiple steady-state solutions by construction. We show that our method is applicable to most common classes of biochemical reaction networks containing inhibition terms, mass-action and Hill-type kinetic equations. Comparing the performance of parameter estimation for different analytical and numerical methods of incorporating steady-state information, we show that our approach is especially well-tailored to guarantee a high success rate of optimization.

Driving the Model to Its Limit: Profile Likelihood Based Model Reduction.

In systems biology, one of the major tasks is to tailor model complexity to information content of the data. A useful model should describe the data and produce well-determined parameter estimates and predictions. Too small of a model will not be able to describe the data whereas a model which is too large tends to overfit measurement errors and does not provide precise predictions. Typically, the model is modified and tuned to fit the data, which often results in an oversized model. To restore the balance between model complexity and available measurements, either new data has to be gathered or the model has to be reduced. In this manuscript, we present a data-based method for reducing non-linear models. The profile likelihood is utilised to assess parameter identifiability and designate likely candidates for reduction. Parameter dependencies are analysed along profiles, providing context-dependent suggestions for the type of reduction. We discriminate four distinct scenarios, each associated with a specific model reduction strategy. Iterating the presented procedure eventually results in an identifiable model, which is capable of generating precise and testable predictions. Source code for all toy examples is provided within the freely available, open-source modelling environment Data2Dynamics based on MATLAB available at http://www.data2dynamics.org/, as well as the R packages dMod/cOde available at https://github.com/dkaschek/. Moreover, the concept is generally applicable and can readily be used with any software capable of calculating the profile likelihood.

Higher-order Lie symmetries in identifiability and predictability analysis of dynamic models.

Parameter estimation in ordinary differential equations (ODEs) has manifold applications not only in physics but also in the life sciences. When estimating the ODE parameters from experimentally observed data, the modeler is frequently concerned with the question of parameter identifiability. The source of parameter nonidentifiability is tightly related to Lie group symmetries. In the present work, we establish a direct search algorithm for the determination of admitted Lie group symmetries. We clarify the relationship between admitted symmetries and parameter nonidentifiability. The proposed algorithm is applied to illustrative toy models as well as a data-based ODE model of the NFκB signaling pathway. We find that besides translations and scaling transformations also higher-order transformations play a role. Enabled by the knowledge about the explicit underlying symmetry transformations, we show how models with nonidentifiable parameters can still be employed to make reliable predictions.

RPPanalyzer toolbox: an improved R package for analysis of reverse phase protein array data.

Analysis of large-scale proteomic data sets requires specialized software tools, tailored toward the requirements of individual approaches. Here we introduce an extension of an open-source software solution for analyzing reverse phase protein array (RPPA) data. The R package RPPanalyzer was designed for data preprocessing followed by basic statistical analyses and proteomic data visualization. In this update, we merged relevant data preprocessing steps into a single user-friendly function and included a new method for background noise correction as well as new methods for noise estimation and averaging of replicates to transform data in such a way that they can be used as input for a new time course plotting function. We demonstrate the robustness of our enhanced RPPanalyzer platform by analyzing longitudinal RPPA data of MET receptor signaling upon stimulation with different hepatocyte growth factor concentrations.

Profile likelihood in systems biology.

Inferring knowledge about biological processes by a mathematical description is a major characteristic of Systems Biology. To understand and predict system's behavior the available experimental information is translated into a mathematical model. Since the availability of experimental data is often limited and measurements contain noise, it is essential to appropriately translate experimental uncertainty to model parameters as well as to model predictions. This is especially important in Systems Biology because typically large and complex models are applied and therefore the limited experimental knowledge might yield weakly specified model components. Likelihood profiles have been recently suggested and applied in the Systems Biology for assessing parameter and prediction uncertainty. In this article, the profile likelihood concept is reviewed and the potential of the approach is demonstrated for a model of the erythropoietin (EPO) receptor.

Pre-clustering of the B cell antigen receptor demonstrated by mathematically extended electron microscopy.

The B cell antigen receptor (BCR) plays a crucial role in adaptive immunity, since antigen-induced signaling by the BCR leads to the activation of the B cell and production of antibodies during an immune response. However, the spatial nano-scale organization of the BCR on the cell surface prior to antigen encounter is still controversial. Here, we fixed murine B cells, stained the BCRs on the cell surface with immuno-gold and visualized the distribution of the gold particles by transmission electron microscopy. Approximately 30% of the gold particles were clustered. However the low staining efficiency of 15% precluded a quantitative conclusion concerning the oligomerization state of the BCRs. To overcome this limitation, we used Monte-Carlo simulations to include or to exclude possible distributions of the BCRs. Our combined experimental-modeling approach assuming the lowest number of different BCR sizes to explain the observed gold distribution suggests that 40% of the surface IgD-BCR was present in dimers and 60% formed large laminar clusters of about 18 receptors. In contrast, a transmembrane mutant of the mIgD molecule only formed IgD-BCR dimers. Our approach complements high resolution fluorescence imaging and clearly demonstrates the existence of pre-formed BCR clusters on resting B cells, questioning the classical cross-linking model of BCR activation.

Lessons learned from quantitative dynamical modeling in systems biology.

Due to the high complexity of biological data it is difficult to disentangle cellular processes relying only on intuitive interpretation of measurements. A Systems Biology approach that combines quantitative experimental data with dynamic mathematical modeling promises to yield deeper insights into these processes. Nevertheless, with growing complexity and increasing amount of quantitative experimental data, building realistic and reliable mathematical models can become a challenging task: the quality of experimental data has to be assessed objectively, unknown model parameters need to be estimated from the experimental data, and numerical calculations need to be precise and efficient. Here, we discuss, compare and characterize the performance of computational methods throughout the process of quantitative dynamic modeling using two previously established examples, for which quantitative, dose- and time-resolved experimental data are available. In particular, we present an approach that allows to determine the quality of experimental data in an efficient, objective and automated manner. Using this approach data generated by different measurement techniques and even in single replicates can be reliably used for mathematical modeling. For the estimation of unknown model parameters, the performance of different optimization algorithms was compared systematically. Our results show that deterministic derivative-based optimization employing the sensitivity equations in combination with a multi-start strategy based on latin hypercube sampling outperforms the other methods by orders of magnitude in accuracy and speed. Finally, we investigated transformations that yield a more efficient parameterization of the model and therefore lead to a further enhancement in optimization performance. We provide a freely available open source software package that implements the algorithms and examples compared here.

Semi-automatic determination of cell surface areas used in systems biology.

Quantitative biology requires high precision measurement of cellular parameters such as surface areas or volumes. Here, we have developed an integrated approach in which the data from 3D confocal microscopy and 2D high-resolution transmission electron microscopy were combined. The volumes and diameters of the cells within one population were automatically measured from the confocal data sets. The perimeter of the cell slices was measured in the TEM images using a semi-automated segmentation into background, cytoplasm and nucleus. These data in conjunction with approaches from stereology allowed for an unbiased estimate of surface areas with high accuracy. We have determined the volumes and surface areas of the cells and nuclei of six different immune cell types. In mast cells for example, the resulting cell surface was 3.5 times larger than the theoretical surface assuming the cell was a sphere with the same volume. Thus, our accurate data can now serve as inputs in modeling approaches in systems immunology.

A variational approach to parameter estimation in ordinary differential equations.

BACKGROUND: Ordinary differential equations are widely-used in the field of systems biology and chemical engineering to model chemical reaction networks. Numerous techniques have been developed to estimate parameters like rate constants, initial conditions or steady state concentrations from time-resolved data. In contrast to this countable set of parameters, the estimation of entire courses of network components corresponds to an innumerable set of parameters. RESULTS: The approach presented in this work is able to deal with course estimation for extrinsic system inputs or intrinsic reactants, both not being constrained by the reaction network itself. Our method is based on variational calculus which is carried out analytically to derive an augmented system of differential equations including the unconstrained components as ordinary state variables. Finally, conventional parameter estimation is applied to the augmented system resulting in a combined estimation of courses and parameters. CONCLUSIONS: The combined estimation approach takes the uncertainty in input courses correctly into account. This leads to precise parameter estimates and correct confidence intervals. In particular this implies that small motifs of large reaction networks can be analysed independently of the rest. By the use of variational methods, elements from control theory and statistics are combined allowing for future transfer of methods between the two fields.