Identifiability and observability analysis for experimental design in nonlinear dynamical models

In silico model development and optimization of in vitro lung cell population growth

Tissue engineering predominantly relies on trial and error in vitro and ex vivo experiments to develop protocols and bioreactors to generate functional tissues. As an alternative, in silico methods have the potential to significantly reduce the timelines and costs of experimental programs for tissue engineering. In this paper, we propose a methodology to formulate, select, calibrate, and test mathematical models to predict cell population growth as a function of the biochemical environment and to design optimal experimental protocols for model inference of in silico model parameters. We systematically combine methods from the experimental design, mathematical statistics, and optimization literature to develop unique and explainable mathematical models for cell population dynamics. The proposed methodology is applied to the development of this first published model for a population of the airway-relevant bronchio-alveolar epithelial (BEAS-2B) cell line as a function of the concentration of metabolic-related biochemical substrates. The resulting model is a system of ordinary differential equations that predict the temporal dynamics of BEAS-2B cell populations as a function of the initial seeded cell population and the glucose, oxygen, and lactate concentrations in the growth media, using seven parameters rigorously inferred from optimally designed in vitro experiments.

Evolving Improved Sampling Protocols for Dose-Response Modelling Using Genetic Algorithms with a Profile-Likelihood Metric

Practical limitations of quality and quantity of data can limit the precision of parameter identification in mathematical models. Model-based experimental design approaches have been developed to minimise parameter uncertainty, but the majority of these approaches have relied on first-order approximations of model sensitivity at a local point in parameter space. Practical identifiability approaches such as profile-likelihood have shown potential for quantifying parameter uncertainty beyond linear approximations. This research presents a genetic algorithm approach to optimise sample timing across various parameterisations of a demonstrative PK-PD model with the goal of aiding experimental design. The optimisation relies on a chosen metric of parameter uncertainty that is based on the profile-likelihood method. Additionally, the approach considers cases where multiple parameter scenarios may require simultaneous optimisation. The genetic algorithm approach was able to locate near-optimal sampling protocols for a wide range of sample number (n = 3-20), and it reduced the parameter variance metric by 33-37% on average. The profile-likelihood metric also correlated well with an existing Monte Carlo-based metric (with a worst-case r > 0.89), while reducing computational cost by an order of magnitude. The combination of the new profile-likelihood metric and the genetic algorithm demonstrate the feasibility of considering the nonlinear nature of models in optimal experimental design at a reasonable computational cost. The outputs of such a process could allow for experimenters to either improve parameter certainty given a fixed number of samples, or reduce sample quantity while retaining the same level of parameter certainty.

Bayesian Regression Facilitates Quantitative Modeling of Cell Metabolism

This paper presents Maud, a command-line application that implements Bayesian statistical inference for kinetic models of biochemical metabolic reaction networks. Maud takes into account quantitative information from omics experiments and background knowledge as well as structural information about kinetic mechanisms, regulatory interactions, and enzyme knockouts. Our paper reviews the existing options in this area, presents a case study illustrating how Maud can be used to analyze a metabolic network, and explains the biological, statistical, and computational design decisions underpinning Maud.

Minimally sufficient experimental design using identifiability analysis

Mathematical models are increasingly being developed and calibrated in tandem with data collection, empowering scientists to intervene in real time based on quantitative model predictions. Well-designed experiments can help augment the predictive power of a mathematical model but the question of when to collect data to maximize its utility for a model is non-trivial. Here we define data as model-informative if it results in a unique parametrization, assessed through the lens of practical identifiability. The framework we propose identifies an optimal experimental design (how much data to collect and when to collect it) that ensures parameter identifiability (permitting confidence in model predictions), while minimizing experimental time and costs. We demonstrate the power of the method by applying it to a modified version of a classic site-of-action pharmacokinetic/pharmacodynamic model that describes distribution of a drug into the tumor microenvironment (TME), where its efficacy is dependent on the level of target occupancy in the TME. In this context, we identify a minimal set of time points when data needs to be collected that robustly ensures practical identifiability of model parameters. The proposed methodology can be applied broadly to any mathematical model, allowing for the identification of a minimally sufficient experimental design that collects the most informative data.

Likelihood-ratio test statistic for the finite-sample case in nonlinear ordinary differential equation models

Likelihood ratios are frequently utilized as basis for statistical tests, for model selection criteria and for assessing parameter and prediction uncertainties, e.g. using the profile likelihood. However, translating these likelihood ratios into p-values or confidence intervals requires the exact form of the test statistic's distribution. The lack of knowledge about this distribution for nonlinear ordinary differential equation (ODE) models requires an approximation which assumes the so-called asymptotic setting, i.e. a sufficiently large amount of data. Since the amount of data from quantitative molecular biology is typically limited in applications, this finite-sample case regularly occurs for mechanistic models of dynamical systems, e.g. biochemical reaction networks or infectious disease models. Thus, it is unclear whether the standard approach of using statistical thresholds derived for the asymptotic large-sample setting in realistic applications results in valid conclusions. In this study, empirical likelihood ratios for parameters from 19 published nonlinear ODE benchmark models are investigated using a resampling approach for the original data designs. Their distributions are compared to the asymptotic approximation and statistical thresholds are checked for conservativeness. It turns out, that corrections of the likelihood ratios in such finite-sample applications are required in order to avoid anti-conservative results.

Delay-differential SEIR modeling for improved modelling of infection dynamics

SEIR (Susceptible-Exposed-Infected-Recovered) approach is a classic modeling method that is frequently used to study infectious diseases. However, in the vast majority of such models transitions from one population group to another are described using the mass-action law. That causes inability to reproduce observable dynamics of an infection such as the incubation period or progression of the disease's symptoms. In this paper, we propose a new approach to simulate the epidemic dynamics based on a system of differential equations with time delays and instant transitions to approximate durations of transition processes more correctly and make model parameters more clear. The suggested approach can be applied not only to Covid-19 but also to the study of other infectious diseases. We utilized it in the development of the delay-based model of the COVID-19 pandemic in Germany and France. The model takes into account testing of different population groups, symptoms progression from mild to critical, vaccination, duration of protective immunity and new virus strains. The stringency index was used as a generalized characteristic of the non-pharmaceutical government interventions in corresponding countries to contain the virus spread. The parameter identifiability analysis demonstrated that the presented modeling approach enables to significantly reduce the number of parameters and make them more identifiable. Both models are publicly available.

Dynamic modelling predicts lactate and IL-1β as interventional targets in metabolic-inflammation-clock regulatory loop in glioma

In an attempt to understand the role of dysregulated circadian rhythm in glioma, our recent findings highlighted the existence of a feed-forward loop between tumour metabolite lactate, pro-inflammatory cytokine IL-1β and circadian CLOCK. To further elucidate the implication of this complex interplay, we developed a mathematical model that quantitatively describes this lactate dehydrogenase A (LDHA)-IL-1β-CLOCK/BMAL1 circuit and predicts potential therapeutic targets. The model was calibrated on quantitative western blotting data in two glioma cell lines in response to either lactate inhibition or IL-1β stimulation. The calibrated model described the experimental data well and most of the parameters were identifiable, thus the model was predictive. Sensitivity analysis identified IL-1β and LDHA as potential intervention points. Mathematical models described here can be useful to understand the complex interrelationship between metabolism, inflammation and circadian rhythm, and in designing effective therapeutic strategies. Our findings underscore the importance of including the circadian clock when developing pharmacological approaches that target aberrant tumour metabolism and inflammation. Insight box  The complex interplay of metabolism-inflammation-circadian rhythm in tumours is not well understood. Our recent findings provided evidence of a feed-forward loop between tumour metabolite lactate, pro-inflammatory cytokine IL-1β and circadian CLOCK/BMAL1 in glioma. To elucidate the implication of this complex interplay, we developed a mathematical model that quantitatively describes this LDHA-IL-1β-CLOCK/BMAL1 circuit and integrates experimental data to predict potential therapeutic targets. The study employed a multi-start optimization strategy and profile likelihood estimations for parameter estimation and assessing identifiability. The simulations are in reasonable agreement with the experimental data. Sensitivity analysis found LDHA and IL-1β as potential therapeutic points. Mathematical models described here can provide insights to understand the complex interrelationship between metabolism, inflammation and circadian rhythm, and in identifying effective therapeutic targets.

Approaches to Parameter Estimation from Model Neurons and Biological Neurons

Model optimization in neuroscience has focused on inferring intracellular parameters from time series observations of the membrane voltage and calcium concentrations. These parameters constitute the fingerprints of ion channel subtypes and may identify ion channel mutations from observed changes in electrical activity. A central question in neuroscience is whether computational methods may obtain ion channel parameters with sufficient consistency and accuracy to provide new information on the underlying biology. Finding single-valued solutions in particular, remains an outstanding theoretical challenge. This note reviews recent progress in the field. It first covers well-posed problems and describes the conditions that the model and data need to meet to warrant the recovery of all the original parameters—even in the presence of noise. The main challenge is model error, which reflects our lack of knowledge of exact equations. We report on strategies that have been partially successful at inferring the parameters of rodent and songbird neurons, when model error is sufficiently small for accurate predictions to be made irrespective of stimulation.

Body-Worn IMU-Based Human Hip and Knee Kinematics Estimation during Treadmill Walking

Traditionally, inertial measurement unit (IMU)-based human joint angle estimation techniques are evaluated for general human motion where human joints explore all of their degrees of freedom. Pure human walking, in contrast, limits the motion of human joints and may lead to unobservability conditions that confound magnetometer-free IMU-based methods. This work explores the unobservability conditions emergent during human walking and expands upon a previous IMU-based method for the human knee to also estimate human hip angles relative to an assumed vertical datum. The proposed method is evaluated (N=12) in a human subject study and compared against an optical motion capture system. Accuracy of human knee flexion/extension angle (7.87∘ absolute root mean square error (RMSE)), hip flexion/extension angle (3.70∘ relative RMSE), and hip abduction/adduction angle (4.56∘ relative RMSE) during walking are similar to current state-of-the-art self-calibrating IMU methods that use magnetometers. Larger errors of hip internal/external rotation angle (6.27∘ relative RMSE) are driven by IMU heading drift characteristic of magnetometer-free approaches and non-hinge kinematics of the hip during gait, amongst other error sources. One of these sources of error, soft tissue perturbations during gait, is explored further in the context of knee angle estimation and it was observed that the IMU method may overestimate the angle during stance and underestimate the angle during swing. The presented method and results provide a novel combination of observability considerations, heuristic correction methods, and validation techniques to magnetic-blind, kinematic-only IMU-based skeletal pose estimation during human tasks with degenerate kinematics (e.g., straight line walking).

Reconstructing asynchrony for mechanical ventilation using a hysteresis loop virtual patient model

BACKGROUND

Patient-specific lung mechanics during mechanical ventilation (MV) can be identified from measured waveforms of fully ventilated, sedated patients. However, asynchrony due to spontaneous breathing (SB) effort can be common, altering these waveforms and reducing the accuracy of identified, model-based, and patient-specific lung mechanics.

METHODS

Changes in patient-specific lung elastance over a pressure-volume (PV) loop, identified using hysteresis loop analysis (HLA), are used to detect the occurrence of asynchrony and identify its type and pattern. The identified HLA parameters are then combined with a nonlinear mechanics hysteresis loop model (HLM) to extract and reconstruct ventilated waveforms unaffected by asynchronous breaths. Asynchrony magnitude can then be quantified using an energy-dissipation metric, Easyn, comparing PV loop area between model-reconstructed and original, altered asynchronous breathing cycles. Performance is evaluated using both test-lung experimental data with a known ground truth and clinical data from four patients with varying levels of asynchrony.

RESULTS

Root mean square errors for reconstructed PV loops are within 5% for test-lung experimental data, and 10% for over 90% of clinical data. Easyn clearly matches known asynchrony magnitude for experimental data with RMS errors < 4.1%. Clinical data performance shows 57% breaths having Easyn > 50% for Patient 1 and 13% for Patient 2. Patient 3 only presents 20% breaths with Easyn > 10%. Patient 4 has Easyn = 0 for 96% breaths showing accuracy in a case without asynchrony.

CONCLUSIONS

Experimental test-lung validation demonstrates the method's reconstruction accuracy and generality in controlled scenarios. Clinical validation matches direct observations of asynchrony in incidence and quantifies magnitude, including cases without asynchrony, validating its robustness and potential efficacy as a clinical real-time asynchrony monitoring tool.

GAMES: A Dynamic Model Development Workflow for Rigorous Characterization of Synthetic Genetic Systems

Mathematical modeling is invaluable for advancing understanding and design of synthetic biological systems. However, the model development process is complicated and often unintuitive, requiring iteration on various computational tasks and comparisons with experimental data. Ad hoc model development can pose a barrier to reproduction and critical analysis of the development process itself, reducing the potential impact and inhibiting further model development and collaboration. To help practitioners manage these challenges, we introduce the Generation and Analysis of Models for Exploring Synthetic Systems (GAMES) workflow, which includes both automated and human-in-the-loop processes. We systematically consider the process of developing dynamic models, including model formulation, parameter estimation, parameter identifiability, experimental design, model reduction, model refinement, and model selection. We demonstrate the workflow with a case study on a chemically responsive transcription factor. The generalizable workflow presented in this tutorial can enable biologists to more readily build and analyze models for various applications.

Sensitivity and identifiability analysis of COVID-19 pandemic models

The paper presents the results of sensitivity-based identifiability analysis of the COVID-19 pandemic spread models in the Novosibirsk region using the systems of differential equations and mass balance law. The algorithm is built on the sensitivity matrix analysis using the methods of differential and linear algebra. It allows one to determine the parameters that are the least and most sensitive to data changes to build a regularization for solving an identification problem of the most accurate pandemic spread scenarios in the region. The performed analysis has demonstrated that the virus contagiousness is identifiable from the number of daily confirmed, critical and recovery cases. On the other hand, the predicted proportion of the admitted patients who require a ventilator and the mortality rate are determined much less consistently. It has been shown that building a more realistic forecast requires adding additional information about the process such as the number of daily hospital admissions. In our study, the problems of parameter identification using additional information about the number of daily confirmed, critical and mortality cases in the region were reduced to minimizing the corresponding misfit functions. The minimization problem was solved through the differential evolution method that is widely applied for stochastic global optimization. It has been demonstrated that a more general COVID-19 spread compartmental model consisting of seven ordinary differential equations describes the main trend of the spread and is sensitive to the peaks of confirmed cases but does not qualitatively describe small statistical datasets such as the number of daily critical cases or mortality that can lead to errors in forecasting. A more detailed agent-oriented model has been able to capture statistical data with additional noise to build scenarios of COVID-19 spread in the region.

Assessing the role of initial conditions in the local structural identifiability of large dynamic models

Structural identifiability is a binary property that determines whether or not unique parameter values can, in principle, be estimated from error-free input-output data. The many papers that have been written on this topic collectively stress the importance of this a priori analysis in the model development process. The story however, often ends with a structurally unidentifiable model. This may leave a model developer with no plan of action on how to address this potential issue. We continue this model exploration journey by identifying one of the possible sources of a model's unidentifiability: problematic initial conditions. It is well-known that certain initial values may result in the loss of local structural identifiability. Nevertheless, literature on this topic has been limited to the analysis of small toy models. Here, we present a systematic approach to detect problematic initial conditions of real-world systems biology models, that are usually not small. A model's identifiability can often be reinstated by changing the value of such problematic initial conditions. This provides modellers an option to resolve the "unidentifiable model" problem. Additionally, a good understanding of which initial values should rather be avoided can be very useful during experimental design. We show how our approach works in practice by applying it to five models. First, two small benchmark models are studied to get the reader acquainted with the method. The first one shows the effect of a zero-valued problematic initial condition. The second one illustrates that the approach also yields correct results in the presence of input signals and that problematic initial conditions need not be zero-values. For the remaining three examples, we set out to identify key initial values which may result in the structural unidentifiability. The third and fourth examples involve a systems biology Epo receptor model and a JAK/STAT model, respectively. In the final Pharmacokinetics model, of which its global structural identifiability has only recently been confirmed, we indicate that there are still sets of initial values for which this property does not hold.

Cell-to-cell variability in JAK2/STAT5 pathway components and cytoplasmic volumes defines survival threshold in erythroid progenitor cells

Survival or apoptosis is a binary decision in individual cells. However, at the cell-population level, a graded increase in survival of colony-forming unit-erythroid (CFU-E) cells is observed upon stimulation with erythropoietin (Epo). To identify components of Janus kinase 2/signal transducer and activator of transcription 5 (JAK2/STAT5) signal transduction that contribute to the graded population response, we extended a cell-population-level model calibrated with experimental data to study the behavior in single cells. The single-cell model shows that the high cell-to-cell variability in nuclear phosphorylated STAT5 is caused by variability in the amount of Epo receptor (EpoR):JAK2 complexes and of SHP1, as well as the extent of nuclear import because of the large variance in the cytoplasmic volume of CFU-E cells. 24-118 pSTAT5 molecules in the nucleus for 120 min are sufficient to ensure cell survival. Thus, variability in membrane-associated processes is sufficient to convert a switch-like behavior at the single-cell level to a graded population-level response.

Body-Worn IMU Human Skeletal Pose Estimation Using a Factor Graph-Based Optimization Framework

Traditionally, inertial measurement units- (IMU) based human joint angle estimation requires a priori knowledge about sensor alignment or specific calibration motions. Furthermore, magnetometer measurements can become unreliable indoors. Without magnetometers, however, IMUs lack a heading reference, which leads to unobservability issues. This paper proposes a magnetometer-free estimation method, which provides desirable observability qualities under joint kinematics that sufficiently excite the lower body degrees of freedom. The proposed lower body model expands on the current self-calibrating human-IMU estimation literature and demonstrates a novel knee hinge model, the inclusion of segment length anthropometry, segment cross-leg length discrepancy, and the relationship between the knee axis and femur/tibia segment. The maximum a posteriori problem is formulated as a factor graph and inference is performed via post-hoc, on-manifold global optimization. The method is evaluated (N = 12) for a prescribed human motion profile task. Accuracy of derived knee flexion/extension angle (4.34∘ root mean square error (RMSE)) without magnetometers is similar to current state-of-the-art with magnetometer use. The developed framework can be expanded for modeling additional joints and constraints.

Estimation of neuron parameters from imperfect observations

The estimation of parameters controlling the electrical properties of biological neurons is essential to determine their complement of ion channels and to understand the function of biological circuits. By synchronizing conductance models to time series observations of the membrane voltage, one may construct models capable of predicting neuronal dynamics. However, identifying the actual set of parameters of biological ion channels remains a formidable theoretical challenge. Here, we present a regularization method that improves convergence towards this optimal solution when data are noisy and the model is unknown. Our method relies on the existence of an offset in parameter space arising from the interplay between model nonlinearity and experimental error. By tuning this offset, we induce saddle-node bifurcations from sub-optimal to optimal solutions. This regularization method increases the probability of finding the optimal set of parameters from 67% to 94.3%. We also reduce parameter correlations by implementing adaptive sampling and stimulation protocols compatible with parameter identifiability requirements. Our results show that the optimal model parameters may be inferred from imperfect observations provided the conditions of observability and identifiability are fulfilled.

BioUML: an integrated environment for systems biology and collaborative analysis of biomedical data

BioUML (homepage: http://www.biouml.org, main public server: https://ict.biouml.org) is a web-based integrated environment (platform) for systems biology and the analysis of biomedical data generated by omics technologies. The BioUML vision is to provide a computational platform to build virtual cell, virtual physiological human and virtual patient. BioUML spans a comprehensive range of capabilities, including access to biological databases, powerful tools for systems biology (visual modelling, simulation, parameters fitting and analyses), a genome browser, scripting (R, JavaScript) and a workflow engine. Due to integration with the Galaxy platform and R/Bioconductor, BioUML provides powerful possibilities for the analyses of omics data. The plug-in-based architecture allows the user to add new functionalities using plug-ins. To facilitate a user focus on a particular task or database, we have developed several predefined perspectives that display only those web interface elements that are needed for a specific task. To support collaborative work on scientific projects, there is a central authentication and authorization system (https://bio-store.org). The diagram editor enables several remote users to simultaneously edit diagrams.

Mechanistic modeling explains the dsRNA length-dependent activation of the RIG-I mediated immune response

In cell-intrinsic antiviral immunity, cytoplasmic receptors such as retinoic acid-inducible gene I (RIG-I) detect viral double-stranded RNA (dsRNA) and trigger a signaling cascade activating the interferon (IFN) system. This leads to the transcription of hundreds of interferon-stimulated genes (ISGs) with a wide range of antiviral effects. This recognition of dsRNA not only has to be very specific to discriminate foreign from self but also highly sensitive to detect even very low numbers of pathogenic dsRNA molecules. Previous work indicated an influence of the dsRNA length on the binding behavior of RIG-I and its potential to elicit antiviral signaling. However, the molecular mechanisms behind the binding process are still under debate. We compare two hypothesized RIG-I binding mechanisms by translating them into mathematical models and analyzing their potential to describe published experimental data. The models consider the length of the dsRNA as well as known RIG-I binding motifs and describe RIG-I pathway activation after stimulation with dsRNA. We show that internal RIG-I binding sites in addition to cooperative RIG-I oligomerization are essential to describe the experimentally observed RIG-I binding behavior and immune response activation for different dsRNA lengths and concentrations. The combination of RIG-I binding to internal sites on the dsRNA and cooperative oligomerization compensates for a lack of high-affinity binding motifs and triggers a strong antiviral response for long dsRNAs. Model analysis reveals dsRNA length-dependency as a potential mechanism to discriminate between different types of dsRNAs: It allows for sensitive detection of small numbers of long dsRNAs, a typical by-product of viral replication, while ensuring tolerance against non-harming small dsRNAs.

Mathematical Modeling and Forecasting of COVID-19 in Moscow and Novosibirsk Region

We investigate inverse problems of finding unknown parameters of mathematical models SEIR-HCD and SEIR-D of COVID-19 spread with additional information about the number of detected cases, mortality, self-isolation coefficient, and tests performed for the city of Moscow and Novosibirsk region since 23.03.2020. In SEIR-HCD the population is divided into seven groups, and in SEIR-D into five groups with similar characteristics and transition probabilities depending on the specific region of interest. An identifiability analysis of SEIR-HCD is made to reveal the least sensitive unknown parameters as related to the additional information. The parameters are corrected by minimizing some objective functionals which is made by stochastic methods (simulated annealing, differential evolution, and genetic algorithm). Prognostic scenarios for COVID-19 spread in Moscow and in Novosibirsk region are developed, and the applicability of the models is analyzed.

Information sensitivity functions to assess parameter information gain and identifiability of dynamical systems

A new class of functions, called the ‘information sensitivity functions’ (ISFs), which quantify the information gain about the parameters through the measurements/observables of a dynamical system are presented. These functions can be easily computed through classical sensitivity functions alone and are based on Bayesian and information-theoretic approaches. While marginal information gain is quantified by decrease in differential entropy, correlations between arbitrary sets of parameters are assessed through mutual information. For individual parameters, these information gains are also presented as marginal posterior variances, and, to assess the effect of correlations, as conditional variances when other parameters are given. The easy to interpret ISFs can be used to (a) identify time intervals or regions in dynamical system behaviour where information about the parameters is concentrated; (b) assess the effect of measurement noise on the information gain for the parameters; (c) assess whether sufficient information in an experimental protocol (input, measurements and their frequency) is available to identify the parameters; (d) assess correlation in the posterior distribution of the parameters to identify the sets of parameters that are likely to be indistinguishable; and (e) assess identifiability problems for particular sets of parameters.

Derivation and identification of a mechanistic model for a branched enzyme-catalyzed carboligation

The kinetic description of enzyme-catalyzed reactions is a core task in biotechnology and biochemical engineering. In particular, mechanistic kinetic models help from the discovery of the biocatalyst throughout its application. Chemo- or enantioselective enzyme reactions often undergo two alternative pathways for the release of two different products from a central intermediate. For these types of reaction, no explicit reaction equations have been derived so far. To this end, we extend the commonly used Cleland's notation for branched reaction pathways and explicitly derive the rate expressions for two-coupled ordered bi-uni reactions. This mechanism also leads to a ping-pong bi-bi mechanism for a transfer reaction between the two products via the same central intermediate of the reaction system. Using the cross-ligation of benzaldehyde and propanal catalyzed by the thiamine diphosphate-dependent enzyme benzaldehyde lyase from Pseudomonas fluorescens yielding (R)-2-hydroxy-1-phenylbutan-1-one as a case study, we performed model-based experimental analysis to show that such a reaction mechanism can be modeled mechanistically and leads to reasonable kinetic parameters.

Particle filtering of dynamical networks: Highlighting observability issues

In a network of high-dimensionality, it is not feasible to measure every single node. Thus, an important goal is to define the optimal choice of sensor nodes that provides a reliable state reconstruction of the network system state-space. This is an observability problem. In this paper, we propose a particle filtering (PF) framework as a way to assess observability properties of a dynamical network, where each node is composed of an individual dynamical system. The PF framework is applied to two benchmarks, networks of Kuramoto and Rössler oscillators, to investigate how the interplay between dynamics and topology impacts the network observability. Based on the numerical results, we conjecture that, when the network nodal dynamics are heterogeneous, better observability is conveyed for sets of sensor nodes that share some dynamical affinity to its neighbourhood. Moreover, we also investigate how the choice of an internal measured variable of a multidimensional sensor node affects the PF performance. The PF framework effectiveness as an observability measure is compared with a well-consolidated nonlinear observability metric for a small network case and some chaotic system benchmarks.

Parameter Identification For Continuous Fluidized Bed Spray Agglomeration

Agglomeration represents an important particle formation process used in many industries. One particularly attractive process setup is continuous fluidized bed spray agglomeration, which features good mixing as well as high heat and mass transfer on the one hand and constant product throughput with constant quality as well as high flow rates compared to batch mode on the other hand. Particle properties such as agglomerate size or porosity significantly affect overall product properties such as re-hydration behavior and dissolubility. These can be influenced by different operating parameters. In this manuscript, a population balance model for a continuous fluidized bed spray agglomeration is presented and adapted to experimental data. Focus is on the description of the dynamic behavior in continuous operation mode in a certain neighborhood around steady-state. Different kernel candidates are evaluated and it is shown that none of the kernels are able to match the first six minutes with time independent parameters. Afterwards, a good fit can be obtained, where the Brownian and the volume independent kernel models match best with the experimental data. Model fit is improved for identification on a shifted time domain neglecting the initial start-up phase. Here, model identifiability is shown and parameter confidence intervals are computed via parametric bootstrap.

Determining minimal output sets that ensure structural identifiability

The process of inferring parameter values from experimental data can be a cumbersome task. In addition, the collection of experimental data can be time consuming and costly. This paper covers both these issues by addressing the following question: “Which experimental outputs should be measured to ensure that unique model parameters can be calculated?”. Stated formally, we examine the topic of minimal output sets that guarantee a model’s structural identifiability. To that end, we introduce an algorithm that guides a researcher as to which model outputs to measure. Our algorithm consists of an iterative structural identifiability analysis and can determine multiple minimal output sets of a model. This choice in different output sets offers researchers flexibility during experimental design. Our method can determine minimal output sets of large differential equation models within short computational times.

Successive complementary model-based experimental designs for parameter estimation of fed-batch bioreactors

When a dynamic model is used for the description of (fed-)batch bioreactors, it is typical that the model parameters are highly correlated to each other. In this case, it is important to keep the parameter correlation as small as possible to obtain a reliable set of parameter estimates. In this study, we propose an anticorrelation parameter estimation scheme that can be best utilized when a number of different batch experiments are sequentially processed. The scheme iteratively performs parameter estimation and model-based design of experiment (MBDOE) at the beginning and between the batches. The important difference from the existing approaches is that the MBDOE objective is defined according to the system analysis performed a priori, so that each new batch supplements what is lacking from the previous batches combined, in terms of information. The use of the scheme is illustrated on a fed-batch bioreactor model.

A Green-Light-Responsive System for the Control of Transgene Expression in Mammalian and Plant Cells

The ever-increasing complexity of synthetic gene networks and applications of synthetic biology requires precise and orthogonal gene expression systems. Of particular interest are systems responsive to light as they enable the control of gene expression dynamics with unprecedented resolution in space and time. While broadly used in mammalian backgrounds, however, optogenetic approaches in plant cells are still limited due to interference of the activating light with endogenous photoreceptors. Here, we describe the development of the first synthetic light-responsive system for the targeted control of gene expression in mammalian and plant cells that responds to the green range of the light spectrum in which plant photoreceptors have minimal activity. We first engineered a system based on the light-sensitive bacterial transcription factor CarH and its cognate DNA operator sequence CarO from Thermus thermophilus to control gene expression in mammalian cells. The system was functional in various mammalian cell lines, showing high induction (up to 350-fold) along with low leakiness, as well as high reversibility. We quantitatively described the systems characteristics by the development and experimental validation of a mathematical model. Finally, we transferred the system into A. thaliana protoplasts and demonstrated gene repression in response to green light. We expect that this system will provide new opportunities in applications based on synthetic gene networks and will open up perspectives for optogenetic studies in mammalian and plant cells.

Profile likelihood-based analyses of infectious disease models

Ordinary differential equation models are frequently applied to describe the temporal evolution of epidemics. However, ordinary differential equation models are also utilized in other scientific fields. We summarize and transfer state-of-the art approaches from other fields like Systems Biology to infectious disease models. For this purpose, we use a simple SIR model with data from an influenza outbreak at an English boarding school in 1978 and a more complex model of a vector-borne disease with data from the Zika virus outbreak in Colombia in 2015-2016. Besides parameter estimation using a deterministic multistart optimization approach, a multitude of analyses based on the profile likelihood are presented comprising identifiability analysis and model reduction. The analyses were performed using the freely available modeling framework Data2Dynamics (data2dynamics.org) which has been awarded as best performing within the DREAM6 parameter estimation challenge and in the DREAM7 network reconstruction challenge.

A symbolic network-based nonlinear theory for dynamical systems observability

When the state of the whole reaction network can be inferred by just measuring the dynamics of a limited set of nodes the system is said to be fully observable. However, as the number of all possible combinations of measured variables and time derivatives spanning the reconstructed state of the system exponentially increases with its dimension, the observability becomes a computationally prohibitive task. Our approach consists in computing the observability coefficients from a symbolic Jacobian matrix whose elements encode the linear, nonlinear polynomial or rational nature of the interaction among the variables. The novelty we introduce in this paper, required for treating large-dimensional systems, is to identify from the symbolic Jacobian matrix the minimal set of variables (together with their time derivatives) candidate to be measured for completing the state space reconstruction. Then symbolic observability coefficients are computed from the symbolic observability matrix. Our results are in agreement with the analytical computations, evidencing the correctness of our approach. Its application to efficiently exploring the dynamics of real world complex systems such as power grids, socioeconomic networks or biological networks is quite promising.

Windowed Granger causal inference strategy improves discovery of gene regulatory networks

Accurate inference of regulatory networks from experimental data facilitates the rapid characterization and understanding of biological systems. High-throughput technologies can provide a wealth of time-series data to better interrogate the complex regulatory dynamics inherent to organisms, but many network inference strategies do not effectively use temporal information. We address this limitation by introducing Sliding Window Inference for Network Generation (SWING), a generalized framework that incorporates multivariate Granger causality to infer network structure from time-series data. SWING moves beyond existing Granger methods by generating windowed models that simultaneously evaluate multiple upstream regulators at several potential time delays. We demonstrate that SWING elucidates network structure with greater accuracy in both in silico and experimentally validated in vitro systems. We estimate the apparent time delays present in each system and demonstrate that SWING infers time-delayed, gene-gene interactions that are distinct from baseline methods. By providing a temporal framework to infer the underlying directed network topology, SWING generates testable hypotheses for gene-gene influences.

Parameter Estimation in Mathematical Models of Viral Infections Using R

In recent years, mathematical modeling approaches have played a central role in understanding and quantifying mechanisms in different viral infectious diseases. In this approach, biology-based hypotheses are expressed via mathematical relations and then tested based on empirical data. The simulation results can be used to either identify underlying mechanisms and provide predictions of infection outcomes or to evaluate the efficacy of a treatment.Conducting parameter estimation for mathematical models is not an easy task. Here we detail an approach to conduct parameter estimation and to evaluate the results using the free software R. The method is applicable to influenza virus dynamics at different complexity levels, widening experimentalists' capabilities in understanding their data. The parameter estimation approach presented here can be also applied to other viral infections or biological applications.

Estimability and dependency analysis of model parameters based on delay coordinates

In data-driven system identification, values of parameters and not observed variables of a given model of a dynamical system are estimated from measured time series. We address the question of estimability and redundancy of parameters and variables, that is, whether unique results can be expected for the estimates or whether, for example, different combinations of parameter values would provide the same measured output. This question is answered by analyzing the null space of the linearized delay coordinates map. Examples with zero-dimensional, one-dimensional, and two-dimensional null spaces are presented employing the Hindmarsh-Rose model, the Colpitts oscillator, and the Rössler system.

Prediction uncertainty and optimal experimental design for learning dynamical systems

Dynamical systems are frequently used to model biological systems. When these models are fit to data, it is necessary to ascertain the uncertainty in the model fit. Here, we present prediction deviation, a metric of uncertainty that determines the extent to which observed data have constrained the model's predictions. This is accomplished by solving an optimization problem that searches for a pair of models that each provides a good fit for the observed data, yet has maximally different predictions. We develop a method for estimating a priori the impact that additional experiments would have on the prediction deviation, allowing the experimenter to design a set of experiments that would most reduce uncertainty. We use prediction deviation to assess uncertainty in a model of interferon-alpha inhibition of viral infection, and to select a sequence of experiments that reduces this uncertainty. Finally, we prove a theoretical result which shows that prediction deviation provides bounds on the trajectories of the underlying true model. These results show that prediction deviation is a meaningful metric of uncertainty that can be used for optimal experimental design.

Model calibration and uncertainty analysis in signaling networks

For a long time the biggest challenges in modeling cellular signal transduction networks has been the inference of crucial pathway components and the qualitative description of their interactions. As a result of the emergence of powerful high-throughput experiments, it is now possible to measure data of high temporal and spatial resolution and to analyze signaling dynamics quantitatively. In addition, this increase of high-quality data is the basis for a better understanding of model limitations and their influence on the predictive power of models. We review established approaches in signal transduction network modeling with a focus on ordinary differential equation models as well as related developments in model calibration. As central aspects of the calibration process we discuss possibilities of model adaptation based on data-driven parameter optimization and the concomitant objective of reducing model uncertainties.

Numerical simulation of vascular tumour growth under antiangiogenic treatment: addressing the paradigm of single-agent bevacizumab therapy with the use of experimental data

Background

Antiangiogenic agents have been recently added to the oncological armamentarium with bevacizumab probably being the most popular representative in current clinical practice. The elucidation of the mode of action of these agents is a prerequisite for personalized prediction of antiangiogenic treatment response and selection of patients who may benefit from this kind of therapy. To this end, having used as a basis a preexisting continuous vascular tumour growth model which addresses the targeted nature of antiangiogenic treatment, we present a paper characterized by the following three features. First, the integration of a two-compartmental bevacizumab specific pharmacokinetic module into the core of the aforementioned preexisting model. Second, its mathematical modification in order to reproduce the asymptotic behaviour of tumour volume in the theoretical case of a total destruction of tumour neovasculature. Third, the exploitation of a range of published animal datasets pertaining to antitumour efficacy of bevacizumab on various tumour types (breast, lung, head and neck, colon).

Results

Results for both the unperturbed growth and the treatment module reveal qualitative similarities with experimental observations establishing the biologically acceptable behaviour of the model. The dynamics of the untreated tumour has been studied via a parameter analysis, revealing the role of each relevant input parameter to tumour evolution. The combined effect of endogenous proangiogenic and antiangiogenic factors on the angiogenic potential of a tumour is also studied, in order to capture the dynamics of molecular competition between the two key-players of tumoural angiogenesis. The adopted methodology also allows accounting for the newly recognized direct antitumour effect of the specific agent.

Conclusions

Interesting observations have been made, suggesting a potential size-dependent tumour response to different treatment modalities and determining the relative timing of cytotoxic versus antiangiogenic agents administration. Insight into the comparative effectiveness of different antiangiogenic treatment strategies is revealed. The results of a series of in vivo experiments in mice bearing diverse types of tumours (breast, lung, head and neck, colon) and treated with bevacizumab are successfully reproduced, supporting thus the validity of the underlying model.

Reviewers

This article was reviewed by L. Hanin, T. Radivoyevitch and L. Edler.

Exploiting intrinsic fluctuations to identify model parameters

Parameterisation of kinetic models plays a central role in computational systems biology. Besides the lack of experimental data of high enough quality, some of the biggest challenges here are identification issues. Model parameters can be structurally non‐identifiable because of functional relationships. Noise in measured data is usually considered to be a nuisance for parameter estimation. However, it turns out that intrinsic fluctuations in particle numbers can make parameters identifiable that were previously non‐identifiable. The authors present a method to identify model parameters that are structurally non‐identifiable in a deterministic framework. The method takes time course recordings of biochemical systems in steady state or transient state as input. Often a functional relationship between parameters presents itself by a one‐dimensional manifold in parameter space containing parameter sets of optimal goodness. Although the system's behaviour cannot be distinguished on this manifold in a deterministic framework it might be distinguishable in a stochastic modelling framework. Their method exploits this by using an objective function that includes a measure for fluctuations in particle numbers. They show on three example models, immigration‐death, gene expression and Epo‐EpoReceptor interaction, that this resolves the non‐identifiability even in the case of measurement noise with known amplitude. The method is applied to partially observed recordings of biochemical systems with measurement noise. It is simple to implement and it is usually very fast to compute. This optimisation can be realised in a classical or Bayesian fashion.

Learning (from) the errors of a systems biology model

Mathematical modelling is a labour intensive process involving several iterations of testing on real data and manual model modifications. In biology, the domain knowledge guiding model development is in many cases itself incomplete and uncertain. A major problem in this context is that biological systems are open. Missed or unknown external influences as well as erroneous interactions in the model could thus lead to severely misleading results. Here we introduce the dynamic elastic-net, a data driven mathematical method which automatically detects such model errors in ordinary differential equation (ODE) models. We demonstrate for real and simulated data, how the dynamic elastic-net approach can be used to automatically (i) reconstruct the error signal, (ii) identify the target variables of model error, and (iii) reconstruct the true system state even for incomplete or preliminary models. Our work provides a systematic computational method facilitating modelling of open biological systems under uncertain knowledge.

Data2Dynamics: a modeling environment tailored to parameter estimation in dynamical systems

Modeling of dynamical systems using ordinary differential equations is a popular approach in the field of systems biology. Two of the most critical steps in this approach are to construct dynamical models of biochemical reaction networks for large datasets and complex experimental conditions and to perform efficient and reliable parameter estimation for model fitting. We present a modeling environment for MATLAB that pioneers these challenges. The numerically expensive parts of the calculations such as the solving of the differential equations and of the associated sensitivity system are parallelized and automatically compiled into efficient C code. A variety of parameter estimation algorithms as well as frequentist and Bayesian methods for uncertainty analysis have been implemented and used on a range of applications that lead to publications.

AVAILABILITY AND IMPLEMENTATION

The Data2Dynamics modeling environment is MATLAB based, open source and freely available at http://www.data2dynamics.org.

CONTACT

andreas.raue@fdm.uni-freiburg.de

SUPPLEMENTARY INFORMATION

Supplementary data are available at Bioinformatics online.

A Priori Identifiability of Target-Mediated Drug Disposition Models and Approximations

A priori identifiability of mathematical models assures that for a given input/output experiment, the parameter set has one unique solution within a defined space, independent of the experimental design. Many biologic therapeutics exhibit target-mediated drug disposition (TMDD), and use of the full compartmental model describing this system is well documented. In practice, estimation of the full parameter set for TMDD models, given real-world clinical data, is characterized by convergence difficulties and unstable solutions. Still, the formal assessment of the a priori identifiability of these systems has yet to be reported. The exact arithmetic rank (EAR) approach was used to test the a priori identifiability of a TMDD model as well as model approximations. The full TMDD and quasi-equilibrium/rapid binding (QE/RB), quasi-steady state (QSS), and Michaelis-Menten (MM) approximations were fully identifiable, a priori, regardless of whether observations were taken from a single or multiple compartments. The results of these identifiability analyses indicated that the difficulty with TMDD model convergence, a posteriori, lies in the experimental design, not in the mathematical identifiability in the lack of samples from several compartments. Experiments can be tailored to resolve these structurally non-identifiable parameters, notwithstanding practical implementation challenges. This work highlights the importance of identifiability analyses, specifically how they can influence experimental design and selection of the appropriate model structure to describe a dynamic biological system.

Cell-specific cardiac electrophysiology models

The traditional cardiac model-building paradigm involves constructing a composite model using data collected from many cells. Equations are derived for each relevant cellular component (e.g., ion channel, exchanger) independently. After the equations for all components are combined to form the composite model, a subset of parameters is tuned, often arbitrarily and by hand, until the model output matches a target objective, such as an action potential. Unfortunately, such models often fail to accurately simulate behavior that is dynamically dissimilar (e.g., arrhythmia) to the simple target objective to which the model was fit. In this study, we develop a new approach in which data are collected via a series of complex electrophysiology protocols from single cardiac myocytes and then used to tune model parameters via a parallel fitting method known as a genetic algorithm (GA). The dynamical complexity of the electrophysiological data, which can only be fit by an automated method such as a GA, leads to more accurately parameterized models that can simulate rich cardiac dynamics. The feasibility of the method is first validated computationally, after which it is used to develop models of isolated guinea pig ventricular myocytes that simulate the electrophysiological dynamics significantly better than does a standard guinea pig model. In addition to improving model fidelity generally, this approach can be used to generate a cell-specific model. By so doing, the approach may be useful in applications ranging from studying the implications of cell-to-cell variability to the prediction of intersubject differences in response to pharmacological treatment.

Mathematical models of cardiac cell electrophysiology are widely used as predictive and illuminatory tools, but have been developed for decades using a suboptimal process. The models are typically constructed by manual adjustment of parameters to fit simple data and therefore often underperform when used to predict complex behavior such as arrhythmias. We present a novel method of model parameterization using automated optimization and dynamically rich fitting data and then demonstrate that this approach is better at finding the “real” model of a cell. Application of the method to cardiac myocytes leads to cell-specific models, which may enable well-controlled studies of both cellular- and subject-level population heterogeneity in disease propensity and response to therapies.

A unified framework for estimating parameters of kinetic biological models

Background

Utilizing kinetic models of biological systems commonly require computational approaches to estimate parameters, posing a variety of challenges due to their highly non-linear and dynamic nature, which is further complicated by the issue of non-identifiability. We propose a novel parameter estimation framework by combining approaches for solving identifiability with a recently introduced filtering technique that can uniquely estimate parameters where conventional methods fail. This framework first conducts a thorough analysis to identify and classify the non-identifiable parameters and provides a guideline for solving them. If no feasible solution can be found, the framework instead initializes the filtering technique with informed prior to yield a unique solution.

Results

This framework has been applied to uniquely estimate parameter values for the sucrose accumulation model in sugarcane culm tissue and a gene regulatory network. In the first experiment the results show the progression of improvement in reliable and unique parameter estimation through the use of each tool to reduce and remove non-identifiability. The latter experiment illustrates the common situation where no further measurement data is available to solve the non-identifiability. These results show the successful application of the informed prior as well as the ease with which parallel data sources may be utilized without increasing the model complexity.

Conclusion

The proposed unified framework is distinct from other approaches by providing a robust and complete solution which yields reliable and unique parameter estimation even in the face of non-identifiability.

Electronic supplementary material

The online version of this article (doi:10.1186/s12859-015-0500-9) contains supplementary material, which is available to authorized users.

Assessing parameter identifiability for dynamic causal modeling of fMRI data

Deterministic dynamic causal modeling (DCM) for fMRI data is a sophisticated approach to analyse effective connectivity in terms of directed interactions between brain regions of interest. To date it is difficult to know if acquired fMRI data will yield precise estimation of DCM parameters. Focusing on parameter identifiability, an important prerequisite for research questions on directed connectivity, we present an approach inferring if parameters of an envisaged DCM are identifiable based on information from fMRI data. With the freely available “attention to motion” dataset, we investigate identifiability of two DCMs and show how different imaging specifications impact on identifiability. We used the profile likelihood, which has successfully been applied in systems biology, to assess the identifiability of parameters in a DCM with specified scanning parameters. Parameters are identifiable when minima of the profile likelihood as well as finite confidence intervals for the parameters exist. Intermediate epoch duration, shorter TR and longer session duration generally increased the information content in the data and thus improved identifiability. Irrespective of biological factors such as size and location of a region, attention should be paid to densely interconnected regions in a DCM, as those seem to be prone to non-identifiability. Our approach, available in the DCMident toolbox, enables to judge if the parameters of an envisaged DCM are sufficiently determined by underlying data without priors as opposed to primarily reflecting the Bayesian priors in a SPM–DCM. Assessments with the DCMident toolbox prior to a study will lead to improved identifiability of the parameters and thus might prevent suboptimal data acquisition. Thus, the toolbox can be used as a preprocessing step to provide immediate statements on parameter identifiability.

A dynamic growth model of Dunaliella salina: parameter identification and profile likelihood analysis

In this work, a photoautotrophic growth model incorporating light and nutrient effects on growth and pigmentation of Dunaliella salina was formulated. The model equations were taken from literature and modified according to the experimental setup with special emphasis on model reduction. The proposed model has been evaluated with experimental data of D. salina cultivated in a flat-plate photobioreactor under stressed and non-stressed conditions. Simulation results show that the model can represent the experimental data accurately. The identifiability of the model parameters was studied using the profile likelihood method. This analysis revealed that three model parameters are practically non-identifiable. However, some of these non-identifiabilities can be resolved by model reduction and additional measurements. As a conclusion, our results suggest that the proposed model equations result in a predictive growth model for D. salina.

Computational strategies for a system-level understanding of metabolism

Cell metabolism is the biochemical machinery that provides energy and building blocks to sustain life. Understanding its fine regulation is of pivotal relevance in several fields, from metabolic engineering applications to the treatment of metabolic disorders and cancer. Sophisticated computational approaches are needed to unravel the complexity of metabolism. To this aim, a plethora of methods have been developed, yet it is generally hard to identify which computational strategy is most suited for the investigation of a specific aspect of metabolism. This review provides an up-to-date description of the computational methods available for the analysis of metabolic pathways, discussing their main advantages and drawbacks. In particular, attention is devoted to the identification of the appropriate scale and level of accuracy in the reconstruction of metabolic networks, and to the inference of model structure and parameters, especially when dealing with a shortage of experimental measurements. The choice of the proper computational methods to derive in silico data is then addressed, including topological analyses, constraint-based modeling and simulation of the system dynamics. A description of some computational approaches to gain new biological knowledge or to formulate hypotheses is finally provided.

Modeling the inflammatory response in the hypothalamus ensuing heat stroke: iterative cycle of model calibration, identifiability analysis, experimental design and data collection

Heat Stroke (HS) is a life-threatening illness caused by prolonged exposure to heat that causes severe hyperthermia and nervous system abnormalities. The long term consequences of HS are poorly understood and deeper insight is required to find possible treatment strategies. Elevated pro- and anti-inflammatory cytokines during HS recovery suggest to play a major role in the immune response. In this study, we developed a mathematical model to understand the interactions and dynamics of cytokines in the hypothalamus, the main thermoregulatory center in the brain. Uncertainty and identifiability analysis of the calibrated model parameters revealed non-identifiable parameters due to the limited amount of data. To overcome the lack of identifiability of the parameters, an iterative cycle of optimal experimental design, data collection, re-calibration and model reduction was applied and further informative experiments were suggested. Additionally, a new method of approximating the prior distribution of the parameters for Bayesian optimal experimental design based on the profile likelihood is presented.

Cause and cure of sloppiness in ordinary differential equation models

Data-based mathematical modeling of biochemical reaction networks, e.g., by nonlinear ordinary differential equation (ODE) models, has been successfully applied. In this context, parameter estimation and uncertainty analysis is a major task in order to assess the quality of the description of the system by the model. Recently, a broadened eigenvalue spectrum of the Hessian matrix of the objective function covering orders of magnitudes was observed and has been termed as sloppiness. In this work, we investigate the origin of sloppiness from structures in the sensitivity matrix arising from the properties of the model topology and the experimental design. Furthermore, we present strategies using optimal experimental design methods in order to circumvent the sloppiness issue and present nonsloppy designs for a benchmark model.