Efficient probabilistic model personalization integrating uncertainty on data and parameters: Application to Eikonal-Diffusion models in cardiac electrophysiology

Biophysical models are increasingly used for medical applications at the organ scale. However, model predictions are rarely associated with a confidence measure although there are important sources of uncertainty in computational physiology methods. For instance, the sparsity and noise of the clinical data used to adjust the model parameters (personalization), and the difficulty in modeling accurately soft tissue physiology. The recent theoretical progresses in stochastic models make their use computationally tractable, but there is still a challenge in estimating patient-specific parameters with such models. In this work we propose an efficient Bayesian inference method for model personalization using polynomial chaos and compressed sensing. This method makes Bayesian inference feasible in real 3D modeling problems. We demonstrate our method on cardiac electrophysiology. We first present validation results on synthetic data, then we apply the proposed method to clinical data. We demonstrate how this can help in quantifying the impact of the data characteristics on the personalization (and thus prediction) results. Described method can be beneficial for the clinical use of personalized models as it explicitly takes into account the uncertainties on the data and the model parameters while still enabling simulations that can be used to optimize treatment. Such uncertainty handling can be pivotal for the proper use of modeling as a clinical tool, because there is a crucial requirement to know the confidence one can have in personalized models.     

Solvating atomic level fine-grained proteins in supra-molecular level coarse-grained water for molecular dynamics simulations

Simulation of the dynamics of a protein in aqueous solution using an atomic model for both the protein and the many water molecules is still computationally extremely demanding considering the time scale of protein motions. The use of supra-atomic or supra-molecular coarse-grained (CG) models may enhance the computational efficiency, but inevitably at the cost of reduced accuracy. Coarse-graining solvent degrees of freedom is likely to yield a favourable balance between reduced accuracy and enhanced computational speed. Here, the use of a supra-molecular coarse-grained water model that largely preserves the thermodynamic and dielectric properties of atomic level fine-grained (FG) water in molecular dynamics simulations of an atomic model for four proteins is investigated. The results of using an FG, a CG, an implicit, or a vacuum solvent environment of the four proteins are compared, and for hen egg-white lysozyme a comparison to NMR data is made. The mixed-grained simulations do not show large differences compared to the FG atomic level simulations, apart from an increased tendency to form hydrogen bonds between long side chains, which is due to the reduced ability of the supra-molecular CG beads that represent five FG water molecules to make solvent-protein hydrogen bonds. But, the mixed-grained simulations are at least an order of magnitude faster than the atomic level ones.     

Predicting metabolic biomarkers of human inborn errors of metabolism

Early diagnosis of inborn errors of metabolism is commonly performed through biofluid metabolomics, which detects specific metabolic biomarkers whose concentration is altered due to genomic mutations. The identification of new biomarkers is of major importance to biomedical research and is usually performed through data mining of metabolomic data. After the recent publication of the genome‐scale network model of human metabolism, we present a novel computational approach for systematically predicting metabolic biomarkers in stochiometric metabolic models. Applying the method to predict biomarkers for disruptions of red‐blood cell metabolism demonstrates a marked correlation with altered metabolic concentrations inferred through kinetic model simulations. Applying the method to the genome‐scale human model reveals a set of 233 metabolites whose concentration is predicted to be either elevated or reduced as a result of 176 possible dysfunctional enzymes. The method's predictions are shown to significantly correlate with known disease biomarkers and to predict many novel potential biomarkers. Using this method to prioritize metabolite measurement experiments to identify new biomarkers can provide an order of a 10‐fold increase in biomarker detection performance.     

metaGEM: reconstruction of genome scale metabolic models directly from metagenomes

Metagenomic analyses of microbial communities have revealed a large degree of interspecies and intraspecies genetic diversity through the reconstruction of metagenome assembled genomes (MAGs). Yet, metabolic modeling efforts mainly rely on reference genomes as the starting point for reconstruction and simulation of genome scale metabolic models (GEMs), neglecting the immense intra- and inter-species diversity present in microbial communities. Here, we present metaGEM (https://github.com/franciscozorrilla/metaGEM), an end-to-end pipeline enabling metabolic modeling of multi-species communities directly from metagenomes. The pipeline automates all steps from the extraction of context-specific prokaryotic GEMs from MAGs to community level flux balance analysis (FBA) simulations. To demonstrate the capabilities of metaGEM, we analyzed 483 samples spanning lab culture, human gut, plant-associated, soil, and ocean metagenomes, reconstructing over 14,000 GEMs. We show that GEMs reconstructed from metagenomes have fully represented metabolism comparable to isolated genomes. We demonstrate that metagenomic GEMs capture intraspecies metabolic diversity and identify potential differences in the progression of type 2 diabetes at the level of gut bacterial metabolic exchanges. Overall, metaGEM enables FBA-ready metabolic model reconstruction directly from metagenomes, provides a resource of metabolic models, and showcases community-level modeling of microbiomes associated with disease conditions allowing generation of mechanistic hypotheses.     

Reconstruction of Danio rerio Metabolic Model Accounting for Subcellular Compartmentalisation

Plant and microbial metabolic engineering is commonly used in the production of functional foods and quality trait improvement. Computational model-based approaches have been used in this important endeavour. However, to date, fish metabolic models have only been scarcely and partially developed, in marked contrast to their prominent success in metabolic engineering. In this study we present the reconstruction of fully compartmentalised models of the Danio rerio (zebrafish) on a global scale. This reconstruction involves extraction of known biochemical reactions in D. rerio for both primary and secondary metabolism and the implementation of methods for determining subcellular localisation and assignment of enzymes. The reconstructed model (ZebraGEM) is amenable for constraint-based modelling analysis, and accounts for 4,988 genes coding for 2,406 gene-associated reactions and only 418 non-gene-associated reactions. A set of computational validations (i.e., simulations of known metabolic functionalities and experimental data) strongly testifies to the predictive ability of the model. Overall, the reconstructed model is expected to lay down the foundations for computational-based rational design of fish metabolic engineering in aquaculture.     

Bayesian calibration of a stochastic,multiscale agent-based model for predicting in vitro tumor growth

Hybrid multiscale agent-based models (ABMs) are unique in their ability to simulate individual cell interactions and microenvironmental dynamics. Unfortunately, the high computational cost of modeling individual cells, the inherent stochasticity of cell dynamics, and numerous model parameters are fundamental limitations of applying such models to predict tumor dynamics. To overcome these challenges, we have developed a coarse-grained two-scale ABM (cgABM) with a reduced parameter space that allows for an accurate and efficient calibration using a set of time-resolved microscopy measurements of cancer cells grown with different initial conditions. The multiscale model consists of a reaction-diffusion type model capturing the spatio-temporal evolution of glucose and growth factors in the tumor microenvironment (at tissue scale), coupled with a lattice-free ABM to simulate individual cell dynamics (at cellular scale). The experimental data consists of BT474 human breast carcinoma cells initialized with different glucose concentrations and tumor cell confluences. The confluence of live and dead cells was measured every three hours over four days. Given this model, we perform a time-dependent global sensitivity analysis to identify the relative importance of the model parameters. The subsequent cgABM is calibrated within a Bayesian framework to the experimental data to estimate model parameters, which are then used to predict the temporal evolution of the living and dead cell populations. To this end, a moment-based Bayesian inference is proposed to account for the stochasticity of the cgABM while quantifying uncertainties due to limited temporal observational data. The cgABM reduces the computational time of ABM simulations by 93% to 97% while staying within a 3% difference in prediction compared to ABM. Additionally, the cgABM can reliably predict the temporal evolution of breast cancer cells observed by the microscopy data with an average error and standard deviation for live and dead cells being 7.61±2.01 and 5.78±1.13, respectively.     

In Silico Investigations of the Anti-Catabolic Effects of Pamidronate and Denosumab on Multiple Myeloma-Induced Bone Disease

It is unclear whether the new anti-catabolic agent denosumab represents a viable alternative to the widely used anti-catabolic agent pamidronate in the treatment of Multiple Myeloma (MM)-induced bone disease. This lack of clarity primarily stems from the lack of sufficient clinical investigations, which are costly and time consuming. However, in silico investigations require less time and expense, suggesting that they may be a useful complement to traditional clinical investigations. In this paper, we aim to (i) develop integrated computational models that are suitable for investigating the effects of pamidronate and denosumab on MM-induced bone disease and (ii) evaluate the responses to pamidronate and denosumab treatments using these integrated models. To achieve these goals, pharmacokinetic models of pamidronate and denosumab are first developed and then calibrated and validated using different clinical datasets. Next, the integrated computational models are developed by incorporating the simulated transient concentrations of pamidronate and denosumab and simulations of their actions on the MM-bone compartment into the previously proposed MM-bone model. These integrated models are further calibrated and validated by different clinical datasets so that they are suitable to be applied to investigate the responses to the pamidronate and denosumab treatments. Finally, these responses are evaluated by quantifying the bone volume, bone turnover, and MM-cell density. This evaluation identifies four denosumab regimes that potentially produce an overall improved bone-related response compared with the recommended pamidronate regime. This in silico investigation supports the idea that denosumab represents an appropriate alternative to pamidronate in the treatment of MM-induced bone disease.     

Numerical approaches for the rapid analysis of prophylactic efficacy against HIV with arbitrary drug-dosing schemes

Pre-exposure prophylaxis (PrEP) is an important pillar to prevent HIV transmission. Because of experimental and clinical shortcomings, mathematical models that integrate pharmacological, viral- and host factors are frequently used to quantify clinical efficacy of PrEP. Stochastic simulations of these models provides sample statistics from which the clinical efficacy is approximated. However, many stochastic simulations are needed to reduce the associated sampling error. To remedy the shortcomings of stochastic simulation, we developed a numerical method that allows predicting the efficacy of arbitrary prophylactic regimen directly from a viral dynamics model, without sampling. We apply the method to various hypothetical dolutegravir (DTG) prophylaxis scenarios. The approach is verified against state-of-the-art stochastic simulation. While the method is more accurate than stochastic simulation, it is superior in terms of computational performance. For example, a continuous 6-month prophylactic profile is computed within a few seconds on a laptop computer. The method’s computational performance, therefore, substantially expands the horizon of feasible analysis in the context of PrEP, and possibly other applications.     

Minimal Conductance-Based Model of Auditory Coincidence Detector Neurons

Sound localization is a fundamental sensory function of a wide variety of animals. The interaural time difference (ITD), an important cue for sound localization, is computed in the auditory brainstem. In our previous modeling study, we introduced a two-compartment Hodgkin-Huxley type model to investigate how cellular and synaptic specializations may contribute to precise ITD computation of the barn owl''s auditory coincidence detector neuron. Although our model successfully reproduced fundamental physiological properties observed in vivo, it was unsuitable for mathematical analyses and large scale simulations because of a number of nonlinear variables. In the present study, we reduce our former model into three types of conductance-based integrate-and-fire (IF) models. We test their electrophysiological properties using data from published in vivo and in vitro studies. Their robustness to parameter changes and computational efficiencies are also examined. Our numerical results suggest that the single-compartment active IF model is superior to other reduced models in terms of physiological reproducibility and computational performance. This model will allow future theoretical studies that use more rigorous mathematical analysis and network simulations.     

Experimental verification of a computational technique for determining ground reactions in human bipedal stance

We have developed a three-dimensional (3D) biomechanical model of human standing that enables us to study the mechanisms of posture and balance simultaneously in various directions in space. Since the two feet are on the ground, the system defines a kinematically closed-chain which has redundancy problems that cannot be resolved using the laws of mechanics alone. We have developed a computational (optimization) technique that avoids the problems with the closed-chain formulation thus giving users of such models the ability to make predictions of joint moments, and potentially, muscle activations using more sophisticated musculoskeletal models. This paper describes the experimental verification of the computational technique that is used to estimate the ground reaction vector acting on an unconstrained foot while the other foot is attached to the ground, thus allowing human bipedal standing to be analyzed as an open-chain system. The computational approach was verified in terms of its ability to predict lower extremity joint moments derived from inverse dynamic simulations performed on data acquired from four able-bodied volunteers standing in various postures on force platforms. Sensitivity analyses performed with model simulations indicated which ground reaction force (GRF) and center of pressure (COP) components were most critical for providing better estimates of the joint moments. Overall, the joint moments predicted by the optimization approach are strongly correlated with the joint moments computed using the experimentally measured GRF and COP (0.78 < or = r(2) < or = 0.99,median,0.96) with a best-fit that was not statistically different from a straight line with unity slope (experimental=computational results) for postures of the four subjects examined. These results indicate that this model-based technique can be relied upon to predict reasonable and consistent estimates of the joint moments using the predicted GRF and COP for most standing postures.     

A computational approach for studying antibody–antigen interactions without prior structural information: the anti‐testosterone binding antibody as a case study

Given the increasing exploitation of antibodies in different contexts such as molecular diagnostics and therapeutics, it would be beneficial to unravel the atomistic level properties of antibody‐antigen complexes with the help of computational modeling. Thus, here we have studied the feasibility of computational tools to gather atomic scale information regarding the antibody‐antigen complexes solely starting from an amino acid sequence. First, we constructed a homology model for the anti‐testosterone binding antibody based on the knowledge based classification of complementary determining regions (CDRs) and implicit solvent molecular dynamics simulations. To further examine whether the generated homology model is suitable for studying antibody‐antigen interactions, docking calculations were carried out followed by binding free‐energy simulations. Our results indicate that with the antibody modeling approach presented here it is possible to construct accurate homology models for antibodies which correctly describes the antibody‐antigen interactions, and produces absolute binding free‐energies that are comparable with experimental values. In addition, our simulations suggest that the conformations of complementary determining regions (CDRs) may considerably change from the X‐ray configuration upon solvation. In conclusion, here we have introduced an antibody modeling workflow that can be used in studying the interactions between antibody and antigen solely based on an amino acid sequence, which in turn provides novel opportunities to tune the properties of antibodies in different applications. Proteins 2017; 85:322–331. © 2016 Wiley Periodicals, Inc.     

A Novel Statistical Model to Estimate Host Genetic Effects Affecting Disease Transmission

There is increasing recognition that genetic diversity can affect the spread of diseases, potentially affecting plant and livestock disease control as well as the emergence of human disease outbreaks. Nevertheless, even though computational tools can guide the control of infectious diseases, few epidemiological models can simultaneously accommodate the inherent individual heterogeneity in multiple infectious disease traits influencing disease transmission, such as the frequently modeled propensity to become infected and infectivity, which describes the host ability to transmit the infection to susceptible individuals. Furthermore, current quantitative genetic models fail to fully capture the heritable variation in host infectivity, mainly because they cannot accommodate the nonlinear infection dynamics underlying epidemiological data. We present in this article a novel statistical model and an inference method to estimate genetic parameters associated with both host susceptibility and infectivity. Our methodology combines quantitative genetic models of social interactions with stochastic processes to model the random, nonlinear, and dynamic nature of infections and uses adaptive Bayesian computational techniques to estimate the model parameters. Results using simulated epidemic data show that our model can accurately estimate heritabilities and genetic risks not only of susceptibility but also of infectivity, therefore exploring a trait whose heritable variation is currently ignored in disease genetics and can greatly influence the spread of infectious diseases. Our proposed methodology offers potential impacts in areas such as livestock disease control through selective breeding and also in predicting and controlling the emergence of disease outbreaks in human populations.     

Agent-based dynamic knowledge representation of Pseudomonas aeruginosa virulence activation in the stressed gut: Towards characterizing host-pathogen interactions in gut-derived sepsis

Background

There is a growing realization that alterations in host-pathogen interactions (HPI) can generate disease phenotypes without pathogen invasion. The gut represents a prime region where such HPI can arise and manifest. Under normal conditions intestinal microbial communities maintain a stable, mutually beneficial ecosystem. However, host stress can lead to changes in environmental conditions that shift the nature of the host-microbe dialogue, resulting in escalation of virulence expression, immune activation and ultimately systemic disease. Effective modulation of these dynamics requires the ability to characterize the complexity of the HPI, and dynamic computational modeling can aid in this task. Agent-based modeling is a computational method that is suited to representing spatially diverse, dynamical systems. We propose that dynamic knowledge representation of gut HPI with agent-based modeling will aid in the investigation of the pathogenesis of gut-derived sepsis.

Methodology/Principal Findings

An agent-based model (ABM) of virulence regulation in Pseudomonas aeruginosa was developed by translating bacterial and host cell sense-and-response mechanisms into behavioral rules for computational agents and integrated into a virtual environment representing the host-microbe interface in the gut. The resulting gut milieu ABM (GMABM) was used to: 1) investigate a potential clinically relevant laboratory experimental condition not yet developed - i.e. non-lethal transient segmental intestinal ischemia, 2) examine the sufficiency of existing hypotheses to explain experimental data - i.e. lethality in a model of major surgical insult and stress, and 3) produce behavior to potentially guide future experimental design - i.e. suggested sample points for a potential laboratory model of non-lethal transient intestinal ischemia. Furthermore, hypotheses were generated to explain certain discrepancies between the behaviors of the GMABM and biological experiments, and new investigatory avenues proposed to test those hypotheses.

Conclusions/Significance

Agent-based modeling can account for the spatio-temporal dynamics of an HPI, and, even when carried out with a relatively high degree of abstraction, can be useful in the investigation of system-level consequences of putative mechanisms operating at the individual agent level. We suggest that an integrated and iterative heuristic relationship between computational modeling and more traditional laboratory and clinical investigations, with a focus on identifying useful and sufficient degrees of abstraction, will enhance the efficiency and translational productivity of biomedical research.     

Studies of folding and misfolding using simplified models

Computer simulations are as vital to our studies of biological systems as experiments. They bridge and rationalize experimental observations, extend the experimental "field of view", which is often limited to a specific time or length scale, and, most importantly, provide novel insights into biological systems, offering hypotheses about yet-to-be uncovered phenomena. These hypotheses spur further experimental discoveries. Simplified molecular models have a special place in the field of computational biology. Branded as less accurate than all-atom protein models, they have offered what all-atom molecular dynamics simulations could not--the resolution of the length and time scales of biological phenomena. Not only have simplified models proven to be accurate in explaining or reproducing several biological phenomena, they have also offered a novel multiscale computational strategy for accessing a broad range of time and length scales upon integration with traditional all-atom simulations. Recent computer simulations of simplified models have shaken or advanced the established understanding of biological phenomena. It was demonstrated that simplified models can be as accurate as traditional molecular dynamics approaches in identifying native conformations of proteins. Their application to protein structure prediction yielded phenomenal accuracy in recapitulating native protein conformations. New studies that utilize the synergy of simplified protein models with all-atom models and experiments yielded novel insights into complex biological processes, such as protein folding, aggregation and the formation of large protein complexes.     

A modified Ising model of Barabási–Albert network with gene-type spins

The central question of systems biology is to understand how individual components of a biological system such as genes or proteins cooperate in emerging phenotypes resulting in the evolution of diseases. As living cells are open systems in quasi-steady state type equilibrium in continuous exchange with their environment, computational techniques that have been successfully applied in statistical thermodynamics to describe phase transitions may provide new insights to the emerging behavior of biological systems. Here we systematically evaluate the translation of computational techniques from solid-state physics to network models that closely resemble biological networks and develop specific translational rules to tackle problems unique to living systems. We focus on logic models exhibiting only two states in each network node. Motivated by the apparent asymmetry between biological states where an entity exhibits boolean states i.e. is active or inactive, we present an adaptation of symmetric Ising model towards an asymmetric one fitting to living systems here referred to as the modified Ising model with gene-type spins. We analyze phase transitions by Monte Carlo simulations and propose a mean-field solution of a modified Ising model of a network type that closely resembles a real-world network, the Barabási–Albert model of scale-free networks. We show that asymmetric Ising models show similarities to symmetric Ising models with the external field and undergoes a discontinuous phase transition of the first-order and exhibits hysteresis. The simulation setup presented herein can be directly used for any biological network connectivity dataset and is also applicable for other networks that exhibit similar states of activity. The method proposed here is a general statistical method to deal with non-linear large scale models arising in the context of biological systems and is scalable to any network size.

     

Simulations of ion current in realistic models of ion channels: the KcsA potassium channel

Realistic studies of ion current in biologic channels present a major challenge for computer simulation approaches. All-atom molecular dynamics simulations involve serious time limitations that prevent their use in direct evaluation of ion current in channels with significant barriers. The alternative use of Brownian dynamics (BD) simulations can provide the current for simplified macroscopic models. However, the time needed for accurate calculations of electrostatic energies can make BD simulations of ion current expensive. The present work develops an approach that overcomes some of the above challenges and allows one to simulate ion currents in models of biologic channels. Our method provides a fast and reliable estimate of the energetics of the system by combining semimacroscopic calculations of the self-energy of each ion and an implicit treatment of the interactions between the ions, as well as the interactions between the ions and the protein-ionizable groups. This treatment involves the use of the semimacroscopic version of the protein dipole Langevin dipole (PDLD/S) model in its linear response approximation (LRA) implementation, which reduces the uncertainties about the value of the protein "dielectric constant." The resulting free energy surface is used to generate the forces for on-the-fly BD simulations of the corresponding ion currents. Our model is examined in a preliminary simulation of the ion current in the KcsA potassium channel. The complete free energy profile for a single ion transport reflects reasonable energetics and captures the effect of the protein-ionized groups. This calculated profile indicates that we are dealing with the channel in its closed state. Reducing the barrier at the gate region allows us to simulate the ion current in a reasonable computational time. Several limiting cases are examined, including those that reproduce the observed current, and the nature of the productive trajectories is considered. The ability to simulate the current in realistic models of ion channels should provide a powerful tool for studies of the biologic function of such systems, including the analysis of the effect of mutations, pH, and electric potentials.     

Allometric scaling of wall shear stress from mice to humans: quantification using cine phase-contrast MRI and computational fluid dynamics

Allometric scaling laws relate structure or function between species of vastly different sizes. They have rarely been derived for hemodynamic parameters known to affect the cardiovascular system, e.g., wall shear stress (WSS). This work describes noninvasive methods to quantify and determine a scaling law for WSS. Geometry and blood flow velocities in the infrarenal aorta of mice and rats under isoflurane anesthesia were quantified using two-dimensional magnetic resonance angiography and phase-contrast magnetic resonance imaging at 4.7 tesla. Three-dimensional models constructed from anatomic data were discretized and used for computational fluid dynamic simulations using phase-contrast velocity imaging data as inlet boundary conditions. WSS was calculated along the infrarenal aorta and compared between species to formulate an allometric equation for WSS. Mean WSS along the infrarenal aorta was significantly greater in mice and rats compared with humans (87.6, 70.5, and 4.8 dyn/cm(2), P < 0.01), and a scaling exponent of -0.38 (R(2) = 0.92) was determined. Manipulation of the murine genome has made small animal models standard surrogates for better understanding the healthy and diseased human cardiovascular system. It has therefore become increasingly important to understand how results scale from mouse to human. This noninvasive methodology provides the opportunity to serially quantify changes in WSS during disease progression and/or therapeutic intervention.     

Can we project changes in fish abundance and distribution in response to climate?

Large‐scale and long‐term changes in fish abundance and distribution in response to climate change have been simulated using both statistical and process‐based models. However, national and regional fisheries management requires also shorter term projections on smaller spatial scales, and these need to be validated against fisheries data. A 26‐year time series of fish surveys with high spatial resolution in the North‐East Atlantic provides a unique opportunity to assess the ability of models to correctly simulate the changes in fish distribution and abundance that occurred in response to climate variability and change. We use a dynamic bioclimate envelope model forced by physical–biogeochemical output from eight ocean models to simulate changes in fish abundance and distribution at scales down to a spatial resolution of 0.5°. When comparing with these simulations with annual fish survey data, we found the largest differences at the 0.5° scale. Differences between fishery model runs driven by different biogeochemical models decrease dramatically when results are aggregated to larger scales (e.g. the whole North Sea), to total catches rather than individual species or when the ensemble mean instead of individual simulations are used. Recent improvements in the fidelity of biogeochemical models translate into lower error rates in the fisheries simulations. However, predictions based on different biogeochemical models are often more similar to each other than they are to the survey data, except for some pelagic species. We conclude that model results can be used to guide fisheries management at larger spatial scales, but more caution is needed at smaller scales.