Fitting Bayesian Stochastic Differential Equation Models with Mixed Effects through a Filtering Approach.

Recent advances in technology contribute to a fast-growing number of studies utilizing intensive longitudinal data, and call for more flexible methods to address the demands that come with them. One issue that arises from collecting longitudinal data from multiple units in time is nested data, where the variability observed in such data is a mixture of within-unit changes and between-unit differences. This article aims to provide a model-fitting approach that simultaneously models the within-unit changes with differential equation models and accounts for between-unit differences with mixed effects. This approach combines a variant of the Kalman filter, the continuous-discrete extended Kalman filter (CDEKF), and the Markov chain Monte Carlo method often employed in the Bayesian framework through the platform Stan. At the same time, it utilizes Stan's functionality of numerical solvers for the implementation of CDEKF. For an empirical illustration, we applied this method in the context of differential equation models to an empirical dataset to explore the physiological dynamics and co-regulation between couples.

Identifiability analysis for models of the translation kinetics after mRNA transfection.

Mechanistic models are a powerful tool to gain insights into biological processes. The parameters of such models, e.g. kinetic rate constants, usually cannot be measured directly but need to be inferred from experimental data. In this article, we study dynamical models of the translation kinetics after mRNA transfection and analyze their parameter identifiability. That is, whether parameters can be uniquely determined from perfect or realistic data in theory and practice. Previous studies have considered ordinary differential equation (ODE) models of the process, and here we formulate a stochastic differential equation (SDE) model. For both model types, we consider structural identifiability based on the model equations and practical identifiability based on simulated as well as experimental data and find that the SDE model provides better parameter identifiability than the ODE model. Moreover, our analysis shows that even for those parameters of the ODE model that are considered to be identifiable, the obtained estimates are sometimes unreliable. Overall, our study clearly demonstrates the relevance of considering different modeling approaches and that stochastic models can provide more reliable and informative results.

A Method of Correcting Estimation Failure in Latent Differential Equations with Comparisons to Kalman Filtering.

Studies have used the latent differential equation (LDE) model to estimate the parameters of damped oscillation in various phenomena, but it has been shown that correct, non-zero parameter estimates are only obtained when the latent series exhibits little or no process noise. Consequently, LDEs are limited to modeling deterministic processes with measurement error rather than those with random behavior in the true latent state. The reasons for these limitations are considered, and a piecewise deterministic approximation (PDA) algorithm is proposed to treat process noise outliers as functional discontinuities and obtain correct estimates of the damping parameter. Comprehensive, random-effects simulations were used to compare results with those obtained using a state-space model (SSM) based on the Kalman filter. The LDE with the PDA algorithm (LDEPDA) successfully recovered the simulated damping parameter under a variety of conditions when process noise was present in the latent state. The LDEPDA had greater precision and accuracy than the SSM when estimating parameters from data with sparse jump discontinuities, but worse performance for diffusion processes overall. All three methods were applied to a sample of postural sway data. The basic LDE estimated zero damping, while the LDEPDA and SSM estimated moderate to high damping. The SSM estimated the smallest standard errors for both frequency and damping parameter estimates.

CRA toolbox: software package for conditional robustness analysis of cancer systems biology models in MATLAB.

BACKGROUND: In cancer research, robustness of a complex biochemical network is one of the most relevant properties to investigate for the development of novel targeted therapies. In cancer systems biology, biological networks are typically modeled through Ordinary Differential Equation (ODE) models. Hence, robustness analysis consists in quantifying how much the temporal behavior of a specific node is influenced by the perturbation of model parameters. The Conditional Robustness Algorithm (CRA) is a valuable methodology to perform robustness analysis on a selected output variable, representative of the proliferation activity of cancer disease. RESULTS: Here we introduce our new freely downloadable software, the CRA Toolbox. The CRA Toolbox is an Object-Oriented MATLAB package which implements the features of CRA for ODE models. It offers the users the ability to import a mathematical model in Systems Biology Markup Language (SBML), to perturb the model parameter space and to choose the reference node for the robustness analysis. The CRA Toolbox allows the users to visualize and save all the generated results through a user-friendly Graphical User Interface (GUI). The CRA Toolbox has a modular and flexible architecture since it is designed according to some engineering design patterns. This tool has been successfully applied in three nonlinear ODE models: the Prostate-specific Pten-/- mouse model, the Pulse Generator Network and the EGFR-IGF1R pathway. CONCLUSIONS: The CRA Toolbox for MATLAB is an open-source tool implementing the CRA to perform conditional robustness analysis. With its unique set of functions, the CRA Toolbox is a remarkable software for the topological study of biological networks. The source and example code and the corresponding documentation are freely available at the web site: http://gitlab.ict4life.com/SysBiOThe/CRA-Matlab .

Follicular competition in cows: the selection of dominant follicles as a synergistic effect.

The reproductive cycle of mono-ovulatory species such as cows or humans is known to show two or more waves of follicular growth and decline between two successive ovulations. Within each wave, there is one dominant follicle escorted by subordinate follicles of varying number. Under the surge of the luteinizing hormone a growing dominant follicle ovulates. Rarely the number of ovulating follicles exceeds one. In the biological literature, the change of hormonal concentrations and individually varying numbers of follicular receptors are made responsible for the selection of exactly one dominant follicle, yet a clear cause has not been identified. In this paper, we suggest a synergistic explanation based on competition, formulated by a parsimoniously defined system of ordinary differential equations (ODEs) that quantifies the time evolution of multiple follicles and their competitive interaction during one wave. Not discriminating between follicles, growth and decline are given by fixed rates. Competition is introduced via a growth-suppressing term, equally supported by all follicles. We prove that the number of dominant follicles is determined exclusively by the ratio of follicular growth and competition. This number turns out to be independent of the number of subordinate follicles. The asymptotic behavior of the corresponding dynamical system is investigated rigorously, where we demonstrate that the [Formula: see text]-limit set only contains fixed points. When also including follicular decline, our ODEs perfectly resemble ultrasound data of bovine follicles. Implications for the involved but not explicitly modeled hormones are discussed.

Dynamic Regulation Responding to an External Stimulus: A Differential Equation Model.

This study examines the dynamic regulation process responding to an external stimulus. The damped oscillator model has been used to describe this process. However, the model does not allow a nonzero steady state, even though the oscillations may continue and do not necessarily damp toward zero. This study introduces the driven damped oscillator model which has an additional parameter to identify different patterns of the steady state. Three methods, generalized local linear approximation, continuous time structural equation modeling, and analytic solutions of differential equations are provided to estimate model parameters. A simulation study indicates that parameters in the driven damped oscillator model are well recovered. The model is then illustrated using a data set on the daily reports of sales after a sale promotion. Potential applications and possible expansions of this model are also discussed.

Reconstruction of disease transmission rates: Applications to measles, dengue, and influenza.

Transmission rates are key in understanding the spread of infectious diseases. Using the framework of compartmental models, we introduce a simple method to reconstruct time series of transmission rates directly from incidence or disease-related mortality data. The reconstruction employs differential equations, which model the time evolution of infective stages and strains. Being sensitive to initial values, the method produces asymptotically correct solutions. The computations are fast, with time complexity being quadratic. We apply the reconstruction to data of measles (England and Wales, 1948-1967), dengue (Thailand, 1982-1999), and influenza (U.S., 1910-1927). The Measles example offers comparison with earlier work. Here we re-investigate reporting corrections, include and exclude demographic information. The dengue example deals with the failure of vector-control measures in reducing dengue hemorrhagic fever (DHF) in Thailand. Two competing mechanisms have been held responsible: strain interaction and demographic transitions. Our reconstruction reveals that both explanations are possible, showing that the increase in DHF cases is consistent with decreasing transmission rates resulting from reduced vector counts. The flu example focuses on the 1918/1919 pandemic, examining the transmission rate evolution for an invading strain. Our analysis indicates that the pandemic strain could have circulated in the population for many months before the pandemic was initiated by an event of highly increased transmission.

Model integration of circadian- and sleep-wake-driven contributions to rhythmic gene expression reveals distinct regulatory principles.

Analyses of gene-expression dynamics in research on circadian rhythms and sleep homeostasis often describe these two processes using separate models. Rhythmically expressed genes are, however, likely to be influenced by both processes. We implemented a driven, damped harmonic oscillator model to estimate the contribution of circadian- and sleep-wake-driven influences on gene expression. The model reliably captured a wide range of dynamics in cortex, liver, and blood transcriptomes taken from mice and humans under various experimental conditions. Sleep-wake-driven factors outweighed circadian factors in driving gene expression in the cortex, whereas the opposite was observed in the liver and blood. Because of tissue- and gene-specific responses, sleep deprivation led to a long-lasting intra- and inter-tissue desynchronization. The model showed that recovery sleep contributed to these long-lasting changes. The results demonstrate that the analyses of the daily rhythms in gene expression must take the complex interactions between circadian and sleep-wake influences into account. A record of this paper's transparent peer review process is included in the supplemental information.

Quantifying CFSE Label Decay in Flow Cytometry Data.

We developed a series of models for the label decay in cell proliferation assays when the intracellular dye carboxyfluorescein succinimidyl ester (CFSE) is used as a staining agent. Data collected from two healthy patients were used to validate the models and to compare the models with the Akiake Information Criteria. The distinguishing features of multiple decay rates in the data are readily characterized and explained via time dependent decay models such as the logistic and Gompertz models.

Optimally controlling nutrition and propulsion force in a long distance running race.

Introduction: Runners competing in races are looking to optimize their performance. In this paper, a runner's performance in a race, such as a marathon, is formulated as an optimal control problem where the controls are: the nutrition intake throughout the race and the propulsion force of the runner. As nutrition is an integral part of successfully running long distance races, it needs to be included in models of running strategies. Methods: We formulate a system of ordinary differential equations to represent the velocity, fat energy, glycogen energy, and nutrition for a runner competing in a long-distance race. The energy compartments represent the energy sources available in the runner's body. We allocate the energy source from which the runner draws, based on how fast the runner is moving. The food consumed during the race is a source term for the nutrition differential equation. With our model, we are investigating strategies to manage the nutrition and propulsion force in order to minimize the running time in a fixed distance race. This requires the solution of a nontrivial singular control problem. Results: As the goal of an optimal control model is to determine the optimal strategy, comparing our results against real data presents a challenge; however, in comparing our results to the world record for the marathon, our results differed by 0.4%, 31 seconds. Per each additional gel consumed, the runner is able to run 0.5 to 0.7 kilometers further in the same amount of time, resulting in a 7.75% increase in taking five 100 calorie gels vs no nutrition. Discussion: Our results confirm the belief that the most effective way to run a race is to run approximately the same pace the entire race without letting one's energies hit zero, by consuming in-race nutrition. While this model does not take all factors into account, we consider it a building block for future models, considering our novel energy representation, and in-race nutrition.

Optimizing differential equations to fit data and predict outcomes.

Many scientific problems focus on observed patterns of change or on how to design a system to achieve particular dynamics. Those problems often require fitting differential equation models to target trajectories. Fitting such models can be difficult because each evaluation of the fit must calculate the distance between the model and target patterns at numerous points along a trajectory. The gradient of the fit with respect to the model parameters can be challenging to compute. Recent technical advances in automatic differentiation through numerical differential equation solvers potentially change the fitting process into a relatively easy problem, opening up new possibilities to study dynamics. However, application of the new tools to real data may fail to achieve a good fit. This article illustrates how to overcome a variety of common challenges, using the classic ecological data for oscillations in hare and lynx populations. Models include simple ordinary differential equations (ODEs) and neural ordinary differential equations (NODEs), which use artificial neural networks to estimate the derivatives of differential equation systems. Comparing the fits obtained with ODEs versus NODEs, representing small and large parameter spaces, and changing the number of variable dimensions provide insight into the geometry of the observed and model trajectories. To analyze the quality of the models for predicting future observations, a Bayesian-inspired preconditioned stochastic gradient Langevin dynamics (pSGLD) calculation of the posterior distribution of predicted model trajectories clarifies the tendency for various models to underfit or overfit the data. Coupling fitted differential equation systems with pSGLD sampling provides a powerful way to study the properties of optimization surfaces, raising an analogy with mutation-selection dynamics on fitness landscapes.

Optimization of Transcription Factor Genetic Circuits.

Transcription factors (TFs) affect the production of mRNAs. In essence, the TFs form a large computational network that controls many aspects of cellular function. This article introduces a computational method to optimize TF networks. The method extends recent advances in artificial neural network optimization. In a simple example, computational optimization discovers a four-dimensional TF network that maintains a circadian rhythm over many days, successfully buffering strong stochastic perturbations in molecular dynamics and entraining to an external day-night signal that randomly turns on and off at intervals of several days. This work highlights the similar challenges in understanding how computational TF and neural networks gain information and improve performance.

Qualitative Modeling, Analysis and Control of Synthetic Regulatory Circuits.

Qualitative modeling approaches are promising and still underexploited tools for the analysis and design of synthetic circuits. They can make predictions of circuit behavior in the absence of precise, quantitative information. Moreover, they provide direct insight into the relation between the feedback structure and the dynamical properties of a network. We review qualitative modeling approaches by focusing on two specific formalisms, Boolean networks and piecewise-linear differential equations, and illustrate their application by means of three well-known synthetic circuits. We describe various methods for the analysis of state transition graphs, discrete representations of the network dynamics that are generated in both modeling frameworks. We also briefly present the problem of controlling synthetic circuits, an emerging topic that could profit from the capacity of qualitative modeling approaches to rapidly scan a space of design alternatives.

A one-step method for modelling longitudinal data with differential equations.

Differential equation models are frequently used to describe non-linear trajectories of longitudinal data. This study proposes a new approach to estimate the parameters in differential equation models. Instead of estimating derivatives from the observed data first and then fitting a differential equation to the derivatives, our new approach directly fits the analytic solution of a differential equation to the observed data, and therefore simplifies the procedure and avoids bias from derivative estimations. A simulation study indicates that the analytic solutions of differential equations (ASDE) approach obtains unbiased estimates of parameters and their standard errors. Compared with other approaches that estimate derivatives first, ASDE has smaller standard error, larger statistical power and accurate Type I error. Although ASDE obtains biased estimation when the system has sudden phase change, the bias is not serious and a solution is also provided to solve the phase problem. The ASDE method is illustrated and applied to a two-week study on consumers' shopping behaviour after a sale promotion, and to a set of public data tracking participants' grammatical facial expression in sign language. R codes for ASDE, recommendations for sample size and starting values are provided. Limitations and several possible expansions of ASDE are also discussed.

Guidelines for benchmarking of optimization-based approaches for fitting mathematical models.

Insufficient performance of optimization-based approaches for the fitting of mathematical models is still a major bottleneck in systems biology. In this article, the reasons and methodological challenges are summarized as well as their impact in benchmark studies. Important aspects for achieving an increased level of evidence for benchmark results are discussed. Based on general guidelines for benchmarking in computational biology, a collection of tailored guidelines is presented for performing informative and unbiased benchmarking of optimization-based fitting approaches. Comprehensive benchmark studies based on these recommendations are urgently required for the establishment of a robust and reliable methodology for the systems biology community.

Classification-Based Inference of Dynamical Models of Gene Regulatory Networks.

Cell-fate decisions during development are controlled by densely interconnected gene regulatory networks (GRNs) consisting of many genes. Inferring and predictively modeling these GRNs is crucial for understanding development and other physiological processes. Gene circuits, coupled differential equations that represent gene product synthesis with a switch-like function, provide a biologically realistic framework for modeling the time evolution of gene expression. However, their use has been limited to smaller networks due to the computational expense of inferring model parameters from gene expression data using global non-linear optimization. Here we show that the switch-like nature of gene regulation can be exploited to break the gene circuit inference problem into two simpler optimization problems that are amenable to computationally efficient supervised learning techniques. We present FIGR (Fast Inference of Gene Regulation), a novel classification-based inference approach to determining gene circuit parameters. We demonstrate FIGR's effectiveness on synthetic data generated from random gene circuits of up to 50 genes as well as experimental data from the gap gene system of Drosophila melanogaster, a benchmark for inferring dynamical GRN models. FIGR is faster than global non-linear optimization by a factor of 600 and its computational complexity scales much better with GRN size. On a practical level, FIGR can accurately infer the biologically realistic gap gene network in under a minute on desktop-class hardware instead of requiring hours of parallel computing. We anticipate that FIGR would enable the inference of much larger biologically realistic GRNs than was possible before.

Effective implicit finite-difference method for sensitivity analysis of stiff stochastic discrete biochemical systems.

Simulation of cellular processes is achieved through a range of mathematical modelling approaches. Deterministic differential equation models are a commonly used first strategy. However, because many biochemical processes are inherently probabilistic, stochastic models are often called for to capture the random fluctuations observed in these systems. In that context, the Chemical Master Equation (CME) is a widely used stochastic model of biochemical kinetics. Use of these models relies on estimates of kinetic parameters, which are often poorly constrained by experimental observations. Consequently, sensitivity analysis, which quantifies the dependence of systems dynamics on model parameters, is a valuable tool for model analysis and assessment. A number of approaches to sensitivity analysis of biochemical models have been developed. In this study, the authors present a novel method for estimation of sensitivity coefficients for CME models of biochemical reaction systems that span a wide range of time-scales. They make use of finite-difference approximations and adaptive implicit tau-leaping strategies to estimate sensitivities for these stiff models, resulting in significant computational efficiencies in comparison with previously published approaches of similar accuracy, as evidenced by illustrative applications.

Linear models of activation cascades: analytical solutions and coarse-graining of delayed signal transduction.

Cellular signal transduction usually involves activation cascades, the sequential activation of a series of proteins following the reception of an input signal. Here, we study the classic model of weakly activated cascades and obtain analytical solutions for a variety of inputs. We show that in the special but important case of optimal gain cascades (i.e. when the deactivation rates are identical) the downstream output of the cascade can be represented exactly as a lumped nonlinear module containing an incomplete gamma function with real parameters that depend on the rates and length of the cascade, as well as parameters of the input signal. The expressions obtained can be applied to the non-identical case when the deactivation rates are random to capture the variability in the cascade outputs. We also show that cascades can be rearranged so that blocks with similar rates can be lumped and represented through our nonlinear modules. Our results can be used both to represent cascades in computational models of differential equations and to fit data efficiently, by reducing the number of equations and parameters involved. In particular, the length of the cascade appears as a real-valued parameter and can thus be fitted in the same manner as Hill coefficients. Finally, we show how the obtained nonlinear modules can be used instead of delay differential equations to model delays in signal transduction.

Estimating biologically relevant parameters under uncertainty for experimental within-host murine West Nile virus infection.

West Nile virus (WNV) is an emerging pathogen that has decimated bird populations and caused severe outbreaks of viral encephalitis in humans. Currently, little is known about the within-host viral kinetics of WNV during infection. We developed mathematical models to describe viral replication, spread and host immune response in wild-type and immunocompromised mice. Our approach fits a target cell-limited model to viremia data from immunocompromised knockout mice and an adaptive immune response model to data from wild-type mice. Using this approach, we first estimate parameters governing viral production and viral spread in the host using simple models without immune responses. We then use these parameters in a more complex immune response model to characterize the dynamics of the humoral immune response. Despite substantial uncertainty in input parameters, our analysis generates relatively precise estimates of important viral characteristics that are composed of nonlinear combinations of model parameters: we estimate the mean within-host basic reproductive number,R0, to be 2.3 (95% of values in the range 1.7-2.9); the mean infectious virion burst size to be 2.9 plaque-forming units (95% of values in the range 1.7-4.7); and the average number of cells infected per infectious virion to be between 0.3 and 0.99. Our analysis gives mechanistic insights into the dynamics of WNV infection and produces estimates of viral characteristics that are difficult to measure experimentally. These models are a first step towards a quantitative understanding of the timing and effectiveness of the humoral immune response in reducing host viremia and consequently the epidemic spread of WNV.

Theoretically exploring direct and indirect chemical effects across ecological and exposure scenarios using mechanistic fate and effects modelling.

Predicting ecosystem response to chemicals is a complex problem in ecotoxicology and a challenge for risk assessors. The variables potentially influencing chemical fate and exposure define the exposure scenario while the variables determining effects at the ecosystem level define the ecological scenario. In absence of any empirical data, the objective of this paper is to present simulations by a fugacity-based fate model and a differential equation-based ecosystem model to theoretically explore how direct and indirect effects on invertebrate shallow pond communities vary with changing ecological and exposure scenarios. These simulations suggest that direct and indirect effects are larger in mesotrophic systems than in oligotrophic systems. In both trophic states, interaction strength (quantified using grazing rates) was suggested a more important driver for the size and recovery from direct and indirect effects than immigration rate. In general, weak interactions led to smaller direct and indirect effects. For chemicals targeting mesozooplankton only, indirect effects were common in (simple) food-chains but rare in (complex) food-webs. For chemicals directly affecting microzooplankton, the dominant zooplankton group in the modelled community, indirect effects occurred both in food-chains and food-webs. We conclude that the choice of the ecological and exposure scenarios in ecotoxicological modelling efforts needs to be justified because of its influence on the prevalence and magnitude of the predicted effects. Overall, more work needs to be done to empirically test the theoretical expectations formulated here.