Model certainty in cellular network-driven processes with missing data.

Mathematical models are often used to explore network-driven cellular processes from a systems perspective. However, a dearth of quantitative data suitable for model calibration leads to models with parameter unidentifiability and questionable predictive power. Here we introduce a combined Bayesian and Machine Learning Measurement Model approach to explore how quantitative and non-quantitative data constrain models of apoptosis execution within a missing data context. We find model prediction accuracy and certainty strongly depend on rigorous data-driven formulations of the measurement, and the size and make-up of the datasets. For instance, two orders of magnitude more ordinal (e.g., immunoblot) data are necessary to achieve accuracy comparable to quantitative (e.g., fluorescence) data for calibration of an apoptosis execution model. Notably, ordinal and nominal (e.g., cell fate observations) non-quantitative data synergize to reduce model uncertainty and improve accuracy. Finally, we demonstrate the potential of a data-driven Measurement Model approach to identify model features that could lead to informative experimental measurements and improve model predictive power.

Unified tumor growth mechanisms from multimodel inference and dataset integration.

Mechanistic models of biological processes can explain observed phenomena and predict responses to a perturbation. A mathematical model is typically constructed using expert knowledge and informal reasoning to generate a mechanistic explanation for a given observation. Although this approach works well for simple systems with abundant data and well-established principles, quantitative biology is often faced with a dearth of both data and knowledge about a process, thus making it challenging to identify and validate all possible mechanistic hypothesis underlying a system behavior. To overcome these limitations, we introduce a Bayesian multimodel inference (Bayes-MMI) methodology, which quantifies how mechanistic hypotheses can explain a given experimental datasets, and concurrently, how each dataset informs a given model hypothesis, thus enabling hypothesis space exploration in the context of available data. We demonstrate this approach to probe standing questions about heterogeneity, lineage plasticity, and cell-cell interactions in tumor growth mechanisms of small cell lung cancer (SCLC). We integrate three datasets that each formulated different explanations for tumor growth mechanisms in SCLC, apply Bayes-MMI and find that the data supports model predictions for tumor evolution promoted by high lineage plasticity, rather than through expanding rare stem-like populations. In addition, the models predict that in the presence of cells associated with the SCLC-N or SCLC-A2 subtypes, the transition from the SCLC-A subtype to the SCLC-Y subtype through an intermediate is decelerated. Together, these predictions provide a testable hypothesis for observed juxtaposed results in SCLC growth and a mechanistic interpretation for tumor treatment resistance.

A biochemical necroptosis model explains cell-type-specific responses to cell death cues.

Necroptosis is a form of regulated cell death associated with degenerative disorders, autoimmune and inflammatory diseases, and cancer. To better understand the biochemical mechanisms regulating necroptosis, we constructed a detailed computational model of tumor necrosis factor-induced necroptosis based on known molecular interactions from the literature. Intracellular protein levels, used as model inputs, were quantified using label-free mass spectrometry, and the model was calibrated using Bayesian parameter inference to experimental protein time course data from a well-established necroptosis-executing cell line. The calibrated model reproduced the dynamics of phosphorylated mixed lineage kinase domain-like protein, an established necroptosis reporter. A subsequent dynamical systems analysis identified four distinct modes of necroptosis signal execution, distinguished by rate constant values and the roles of the RIP1 deubiquitinating enzymes A20 and CYLD. In one case, A20 and CYLD both contribute to RIP1 deubiquitination, in another RIP1 deubiquitination is driven exclusively by CYLD, and in two modes either A20 or CYLD acts as the driver with the other enzyme, counterintuitively, inhibiting necroptosis. We also performed sensitivity analyses of initial protein concentrations and rate constants to identify potential targets for modulating necroptosis sensitivity within each mode. We conclude by associating numerous contrasting and, in some cases, counterintuitive experimental results reported in the literature with one or more of the model-predicted modes of necroptosis execution. In all, we demonstrate that a consensus pathway model of tumor necrosis factor-induced necroptosis can provide insights into unresolved controversies regarding the molecular mechanisms driving necroptosis execution in numerous cell types under different experimental conditions.

Microbench: automated metadata management for systems biology benchmarking and reproducibility in Python.

MOTIVATION: Computational systems biology analyses typically make use of multiple software and their dependencies, which are often run across heterogeneous compute environments. This can introduce differences in performance and reproducibility. Capturing metadata (e.g. package versions, GPU model) currently requires repetitious code and is difficult to store centrally for analysis. Even where virtual environments and containers are used, updates over time mean that versioning metadata should still be captured within analysis pipelines to guarantee reproducibility. RESULTS: Microbench is a simple and extensible Python package to automate metadata capture to a file or Redis database. Captured metadata can include execution time, software package versions, environment variables, hardware information, Python version and more, with plugins. We present three case studies demonstrating Microbench usage to benchmark code execution and examine environment metadata for reproducibility purposes. AVAILABILITY AND IMPLEMENTATION: Install from the Python Package Index using pip install microbench. Source code is available from https://github.com/alubbock/microbench. SUPPLEMENTARY INFORMATION: Supplementary data are available at Bioinformatics online.

Thunor: visualization and analysis of high-throughput dose-response datasets.

High-throughput cell proliferation assays to quantify drug-response are becoming increasingly common and powerful with the emergence of improved automation and multi-time point analysis methods. However, pipelines for analysis of these datasets that provide reproducible, efficient, and interactive visualization and interpretation are sorely lacking. To address this need, we introduce Thunor, an open-source software platform to manage, analyze, and visualize large, dose-dependent cell proliferation datasets. Thunor supports both end-point and time-based proliferation assays as input. It provides a simple, user-friendly interface with interactive plots and publication-quality images of cell proliferation time courses, dose-response curves, and derived dose-response metrics, e.g. IC50, including across datasets or grouped by tags. Tags are categorical labels for cell lines and drugs, used for aggregation, visualization and statistical analysis, e.g. cell line mutation or drug class/target pathway. A graphical plate map tool is included to facilitate plate annotation with cell lines, drugs and concentrations upon data upload. Datasets can be shared with other users via point-and-click access control. We demonstrate the utility of Thunor to examine and gain insight from two large drug response datasets: a large, publicly available cell viability database and an in-house, high-throughput proliferation rate dataset. Thunor is available from www.thunor.net.

Unsupervised logic-based mechanism inference for network-driven biological processes.

Modern analytical techniques enable researchers to collect data about cellular states, before and after perturbations. These states can be characterized using analytical techniques, but the inference of regulatory interactions that explain and predict changes in these states remains a challenge. Here we present a generalizable, unsupervised approach to generate parameter-free, logic-based models of cellular processes, described by multiple discrete states. Our algorithm employs a Hamming-distance based approach to formulate, test, and identify optimized logic rules that link two states. Our approach comprises two steps. First, a model with no prior knowledge except for the mapping between initial and attractor states is built. We then employ biological constraints to improve model fidelity. Our algorithm automatically recovers the relevant dynamics for the explored models and recapitulates key aspects of the biochemical species concentration dynamics in the original model. We present the advantages and limitations of our work and discuss how our approach could be used to infer logic-based mechanisms of signaling, gene-regulatory, or other input-output processes describable by the Boolean formalism.

Arrestin-3 scaffolding of the JNK3 cascade suggests a mechanism for signal amplification.

Scaffold proteins tether and orient components of a signaling cascade to facilitate signaling. Although much is known about how scaffolds colocalize signaling proteins, it is unclear whether scaffolds promote signal amplification. Here, we used arrestin-3, a scaffold of the ASK1-MKK4/7-JNK3 cascade, as a model to understand signal amplification by a scaffold protein. We found that arrestin-3 exhibited >15-fold higher affinity for inactive JNK3 than for active JNK3, and this change involved a shift in the binding site following JNK3 activation. We used systems biochemistry modeling and Bayesian inference to evaluate how the activation of upstream kinases contributed to JNK3 phosphorylation. Our combined experimental and computational approach suggested that the catalytic phosphorylation rate of JNK3 at Thr-221 by MKK7 is two orders of magnitude faster than the corresponding phosphorylation of Tyr-223 by MKK4 with or without arrestin-3. Finally, we showed that the release of activated JNK3 was critical for signal amplification. Collectively, our data suggest a "conveyor belt" mechanism for signal amplification by scaffold proteins. This mechanism informs on a long-standing mystery for how few upstream kinase molecules activate numerous downstream kinases to amplify signaling.

Systems-level network modeling of Small Cell Lung Cancer subtypes identifies master regulators and destabilizers.

Adopting a systems approach, we devise a general workflow to define actionable subtypes in human cancers. Applied to small cell lung cancer (SCLC), the workflow identifies four subtypes based on global gene expression patterns and ontologies. Three correspond to known subtypes (SCLC-A, SCLC-N, and SCLC-Y), while the fourth is a previously undescribed ASCL1+ neuroendocrine variant (NEv2, or SCLC-A2). Tumor deconvolution with subtype gene signatures shows that all of the subtypes are detectable in varying proportions in human and mouse tumors. To understand how multiple stable subtypes can arise within a tumor, we infer a network of transcription factors and develop BooleaBayes, a minimally-constrained Boolean rule-fitting approach. In silico perturbations of the network identify master regulators and destabilizers of its attractors. Specific to NEv2, BooleaBayes predicts ELF3 and NR0B1 as master regulators of the subtype, and TCF3 as a master destabilizer. Since the four subtypes exhibit differential drug sensitivity, with NEv2 consistently least sensitive, these findings may lead to actionable therapeutic strategies that consider SCLC intratumoral heterogeneity. Our systems-level approach should generalize to other cancer types.

Cardiolipin-Dependent Properties of Model Mitochondrial Membranes from Molecular Simulations.

Cardiolipin is an anionic lipid found in the mitochondrial membranes of eukaryotes ranging from unicellular microorganisms to metazoans. This unique lipid contributes to various mitochondrial functions, including metabolism, mitochondrial membrane fusion and/or fission dynamics, and apoptosis. However, differences in cardiolipin content between the two mitochondrial membranes, as well as dynamic fluctuations in cardiolipin content in response to stimuli and cellular signaling events, raise questions about how cardiolipin concentration affects mitochondrial membrane structure and dynamics. Although cardiolipin's structural and dynamic roles have been extensively studied in binary mixtures with other phospholipids, the biophysical properties of cardiolipin in higher number lipid mixtures are still not well resolved. Here, we used molecular dynamics simulations to investigate the cardiolipin-dependent properties of ternary lipid bilayer systems that mimic the major components of mitochondrial membranes. We found that changes to cardiolipin concentration only resulted in minor changes to bilayer structural features but that the lipid diffusion was significantly affected by those alterations. We also found that cardiolipin position along the bilayer surfaces correlated to negative curvature deflections, consistent with the induction of negative curvature stress in the membrane monolayers. This work contributes to a foundational understanding of the role of cardiolipin in altering the properties in ternary lipid mixtures composed of the major mitochondrial phospholipids, providing much-needed insights to help understand how cardiolipin concentration modulates the biophysical properties of mitochondrial membranes.

An unbiased metric of antiproliferative drug effect in vitro.

In vitro cell proliferation assays are widely used in pharmacology, molecular biology, and drug discovery. Using theoretical modeling and experimentation, we show that current metrics of antiproliferative small molecule effect suffer from time-dependent bias, leading to inaccurate assessments of parameters such as drug potency and efficacy. We propose the drug-induced proliferation (DIP) rate, the slope of the line on a plot of cell population doublings versus time, as an alternative, time-independent metric.

PyDREAM: high-dimensional parameter inference for biological models in python.

Summary: Biological models contain many parameters whose values are difficult to measure directly via experimentation and therefore require calibration against experimental data. Markov chain Monte Carlo (MCMC) methods are suitable to estimate multivariate posterior model parameter distributions, but these methods may exhibit slow or premature convergence in high-dimensional search spaces. Here, we present PyDREAM, a Python implementation of the (Multiple-Try) Differential Evolution Adaptive Metropolis [DREAM(ZS)] algorithm developed by Vrugt and ter Braak (2008) and Laloy and Vrugt (2012). PyDREAM achieves excellent performance for complex, parameter-rich models and takes full advantage of distributed computing resources, facilitating parameter inference and uncertainty estimation of CPU-intensive biological models. Availability and implementation: PyDREAM is freely available under the GNU GPLv3 license from the Lopez lab GitHub repository at http://github.com/LoLab-VU/PyDREAM. Contact: c.lopez@vanderbilt.edu. Supplementary information: Supplementary data are available at Bioinformatics online.

GPU-powered model analysis with PySB/cupSODA.

SUMMARY: A major barrier to the practical utilization of large, complex models of biochemical systems is the lack of open-source computational tools to evaluate model behaviors over high-dimensional parameter spaces. This is due to the high computational expense of performing thousands to millions of model simulations required for statistical analysis. To address this need, we have implemented a user-friendly interface between cupSODA, a GPU-powered kinetic simulator, and PySB, a Python-based modeling and simulation framework. For three example models of varying size, we show that for large numbers of simulations PySB/cupSODA achieves order-of-magnitude speedups relative to a CPU-based ordinary differential equation integrator. AVAILABILITY AND IMPLEMENTATION: The PySB/cupSODA interface has been integrated into the PySB modeling framework (version 1.4.0), which can be installed from the Python Package Index (PyPI) using a Python package manager such as pip. cupSODA source code and precompiled binaries (Linux, Mac OS/X, Windows) are available at github.com/aresio/cupSODA (requires an Nvidia GPU; developer.nvidia.com/cuda-gpus). Additional information about PySB is available at pysb.org. CONTACT: paolo.cazzaniga@unibg.it or c.lopez@vanderbilt.edu. SUPPLEMENTARY INFORMATION: Supplementary data are available at Bioinformatics online.

Activated Oncogenic Pathway Modifies Iron Network in Breast Epithelial Cells: A Dynamic Modeling Perspective.

Dysregulation of iron metabolism in cancer is well documented and it has been suggested that there is interdependence between excess iron and increased cancer incidence and progression. In an effort to better understand the linkages between iron metabolism and breast cancer, a predictive mathematical model of an expanded iron homeostasis pathway was constructed that includes species involved in iron utilization, oxidative stress response and oncogenic pathways. The model leads to three predictions. The first is that overexpression of iron regulatory protein 2 (IRP2) recapitulates many aspects of the alterations in free iron and iron-related proteins in cancer cells without affecting the oxidative stress response or the oncogenic pathways included in the model. This prediction was validated by experimentation. The second prediction is that iron-related proteins are dramatically affected by mitochondrial ferritin overexpression. This prediction was validated by results in the pertinent literature not used for model construction. The third prediction is that oncogenic Ras pathways contribute to altered iron homeostasis in cancer cells. This prediction was validated by a combination of simulation experiments of Ras overexpression and catalase knockout in conjunction with the literature. The model successfully captures key aspects of iron metabolism in breast cancer cells and provides a framework upon which more detailed models can be built.

Competition and allostery govern substrate selectivity of cyclooxygenase-2.

Cyclooxygenase-2 (COX-2) oxygenates arachidonic acid (AA) and its ester analog, 2-arachidonoylglycerol (2-AG), to prostaglandins (PGs) and prostaglandin glyceryl esters (PG-Gs), respectively. Although the efficiency of oxygenation of these substrates by COX-2 in vitro is similar, cellular biosynthesis of PGs far exceeds that of PG-Gs. Evidence that the COX enzymes are functional heterodimers suggests that competitive interaction of AA and 2-AG at the allosteric site of COX-2 might result in differential regulation of the oxygenation of the two substrates when both are present. Modulation of AA levels in RAW264.7 macrophages uncovered an inverse correlation between cellular AA levels and PG-G biosynthesis. In vitro kinetic analysis using purified protein demonstrated that the inhibition of 2-AG oxygenation by high concentrations of AA far exceeded the inhibition of AA oxygenation by high concentrations of 2-AG. An unbiased systems-based mechanistic model of the kinetic data revealed that binding of AA or 2-AG at the allosteric site of COX-2 results in a decreased catalytic efficiency of the enzyme toward 2-AG, whereas 2-AG binding at the allosteric site increases COX-2's efficiency toward AA. The results suggest that substrates interact with COX-2 via multiple potential complexes involving binding to both the catalytic and allosteric sites. Competition between AA and 2-AG for these sites, combined with differential allosteric modulation, gives rise to a complex interplay between the substrates, leading to preferential oxygenation of AA.

Integrated, High-Throughput, Multiomics Platform Enables Data-Driven Construction of Cellular Responses and Reveals Global Drug Mechanisms of Action.

An understanding of how cells respond to perturbation is essential for biological applications; however, most approaches for profiling cellular response are limited in scope to pre-established targets. Global analysis of molecular mechanism will advance our understanding of the complex networks constituting cellular perturbation and lead to advancements in areas, such as infectious disease pathogenesis, developmental biology, pathophysiology, pharmacology, and toxicology. We have developed a high-throughput multiomics platform for comprehensive, de novo characterization of cellular mechanisms of action. Platform validation using cisplatin as a test compound demonstrates quantification of over 10 000 unique, significant molecular changes in less than 30 days. These data provide excellent coverage of known cisplatin-induced molecular changes and previously unrecognized insights into cisplatin resistance. This proof-of-principle study demonstrates the value of this platform as a resource to understand complex cellular responses in a high-throughput manner.

Potentiation of antifungal effect of a mixture of two antifungal fractions obtained from Baccharis glutinosa and Jacquinia macrocarpa plants.

The aim of the present work was to evaluate the effect of mixtures of antifungal fractions extracted from Baccharis glutinosa and Jacquinia macrocarpa plants on the development of the filamentous fungi Aspergillus flavus and Fusarium verticillioides. The minimal inhibitory concentration that inhibited 50% of growth (MIC50) of each plant antifungal fraction was determined from the percentage radial growth inhibition of both fungi. Binomial mixtures made with both plant fractions were used at their MIC50 to determine the Fractional Inhibitory Concentration index (FIC index) for each fungus in order to evaluate their synergistic effect. Each synergistic mixture was analyzed in their effect on spore germination, spore size, spore viability, mitotic divisions, hyphal diameter and length, and number of septa per hypha. Some antifungal mixtures, even at low concentrations, showed higher antifungal effect than those of the individual antifungal fraction. The FIC indices of mixtures that showed the highest antifungal activity against A. flavus and F. verticillioides were 0.5272 and 0.4577, respectively, indicating a synergistic effect against both fungi. Only 12% and 8% of the spores of A. flavus and F. verticillioides, respectively, treated with the synergistic mixtures, were able to germinate, although their viability was not affected. An increase in the number of septa per hypha of both fungi was observed. The results indicated that the synergistic mixtures strongly affected the fungal growth even at lower concentrations than those of the individual plant fractions.

Programming biological models in Python using PySB.

Mathematical equations are fundamental to modeling biological networks, but as networks get large and revisions frequent, it becomes difficult to manage equations directly or to combine previously developed models. Multiple simultaneous efforts to create graphical standards, rule-based languages, and integrated software workbenches aim to simplify biological modeling but none fully meets the need for transparent, extensible, and reusable models. In this paper we describe PySB, an approach in which models are not only created using programs, they are programs. PySB draws on programmatic modeling concepts from little b and ProMot, the rule-based languages BioNetGen and Kappa and the growing library of Python numerical tools. Central to PySB is a library of macros encoding familiar biochemical actions such as binding, catalysis, and polymerization, making it possible to use a high-level, action-oriented vocabulary to construct detailed models. As Python programs, PySB models leverage tools and practices from the open-source software community, substantially advancing our ability to distribute and manage the work of testing biochemical hypotheses. We illustrate these ideas using new and previously published models of apoptosis.

Signal execution modes emerge in biochemical reaction networks calibrated to experimental data.

Mathematical models of biomolecular networks are commonly used to study cellular processes; however, their usefulness to explain and predict dynamic behaviors is often questioned due to the unclear relationship between parameter uncertainty and network dynamics. In this work, we introduce PyDyNo (Python dynamic analysis of biochemical networks), a non-equilibrium reaction-flux based analysis to identify dominant reaction paths within a biochemical reaction network calibrated to experimental data. We first show, in a simplified apoptosis execution model, that despite the thousands of parameter vectors with equally good fits to experimental data, our framework identifies the dynamic differences between these parameter sets and outputs three dominant execution modes, which exhibit varying sensitivity to perturbations. We then apply our methodology to JAK2/STAT5 network in colony-forming unit-erythroid (CFU-E) cells and provide previously unrecognized mechanistic explanation for the survival responses of CFU-E cell population that would have been impossible to deduce with traditional protein-concentration based analyses.

Data-driven structural analysis of Small Cell Lung Cancer transcription factor network suggests potential subtype regulators and transition pathways.

Small Cell Lung Cancer (SCLC) is an aggressive disease and challenging to treat due to its mixture of transcriptional subtypes and subtype transitions. Transcription factor (TF) networks have been the focus of studies to identify SCLC subtype regulators via systems approaches. Yet, their structures, which can provide clues on subtype drivers and transitions, are barely investigated. Here, we analyze the structure of an SCLC TF network by using graph theory concepts and identify its structurally important components responsible for complex signal processing, called hubs. We show that the hubs of the network are regulators of different SCLC subtypes by analyzing first the unbiased network structure and then integrating RNA-seq data as weights assigned to each interaction. Data-driven analysis emphasizes MYC as a hub, consistent with recent reports. Furthermore, we hypothesize that the pathways connecting functionally distinct hubs may control subtype transitions and test this hypothesis via network simulations on a candidate pathway and observe subtype transition. Overall, structural analyses of complex networks can identify their functionally important components and pathways driving the network dynamics. Such analyses can be an initial step for generating hypotheses and can guide the discovery of target pathways whose perturbation may change the network dynamics phenotypically.

Data-driven structural analysis of small cell lung cancer transcription factor network suggests potential subtype regulators and transition pathways.

Small cell lung cancer (SCLC) is an aggressive disease and challenging to treat due to its mixture of transcriptional subtypes and subtype transitions. Transcription factor (TF) networks have been the focus of studies to identify SCLC subtype regulators via systems approaches. Yet, their structures, which can provide clues on subtype drivers and transitions, are barely investigated. Here, we analyze the structure of an SCLC TF network by using graph theory concepts and identify its structurally important components responsible for complex signal processing, called hubs. We show that the hubs of the network are regulators of different SCLC subtypes by analyzing first the unbiased network structure and then integrating RNA-seq data as weights assigned to each interaction. Data-driven analysis emphasizes MYC as a hub, consistent with recent reports. Furthermore, we hypothesize that the pathways connecting functionally distinct hubs may control subtype transitions and test this hypothesis via network simulations on a candidate pathway and observe subtype transition. Overall, structural analyses of complex networks can identify their functionally important components and pathways driving the network dynamics. Such analyses can be an initial step for generating hypotheses and can guide the discovery of target pathways whose perturbation may change the network dynamics phenotypically.