An in vitro model of tumor heterogeneity resolves genetic, epigenetic, and stochastic sources of cell state variability.

Tumor heterogeneity is a primary cause of treatment failure and acquired resistance in cancer patients. Even in cancers driven by a single mutated oncogene, variability in response to targeted therapies is well known. The existence of additional genomic alterations among tumor cells can only partially explain this variability. As such, nongenetic factors are increasingly seen as critical contributors to tumor relapse and acquired resistance in cancer. Here, we show that both genetic and nongenetic factors contribute to targeted drug response variability in an experimental model of tumor heterogeneity. We observe significant variability to epidermal growth factor receptor (EGFR) inhibition among and within multiple versions and clonal sublines of PC9, a commonly used EGFR mutant nonsmall cell lung cancer (NSCLC) cell line. We resolve genetic, epigenetic, and stochastic components of this variability using a theoretical framework in which distinct genetic states give rise to multiple epigenetic "basins of attraction," across which cells can transition driven by stochastic noise. Using mutational impact analysis, single-cell differential gene expression, and correlations among Gene Ontology (GO) terms to connect genomics to transcriptomics, we establish a baseline for genetic differences driving drug response variability among PC9 cell line versions. Applying the same approach to clonal sublines, we conclude that drug response variability in all but one of the sublines is due to epigenetic differences; in the other, it is due to genetic alterations. Finally, using a clonal drug response assay together with stochastic simulations, we attribute subclonal drug response variability within sublines to stochastic cell fate decisions and confirm that one subline likely contains genetic resistance mutations that emerged in the absence of drug treatment.

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.

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.

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.

A Nonquiescent "Idling" Population State in Drug-Treated, BRAF-Mutated Melanoma.

Targeted therapy is an effective standard of care in BRAF-mutated malignant melanoma. However, the duration of tumor remission varies unpredictably among patients, and relapse is almost inevitable. Here, we examine the responses of several BRAF-mutated melanoma cell lines (including isogenic subclones) to BRAF inhibitors. We observe complex response dynamics across cell lines, with short-term responses (<100 h) varying from cell line to cell line. In the long term, however, we observe equilibration of all drug-treated populations into a nonquiescent state characterized by a balanced rate of death and division, which we term the "idling" state, and to our knowledge, this state has not been previously reported. Using mathematical modeling, we propose that the observed population-level dynamics are the result of cells transitioning between basins of attraction within a drug-modified phenotypic landscape. Each basin is associated with a drug-induced proliferation rate, a recently introduced metric of an antiproliferative drug effect. The idling population state represents a new dynamic equilibrium in which cells are distributed across the landscape such that the population achieves zero net growth. By fitting our model to experimental drug-response data, we infer the phenotypic landscapes of all considered melanoma cell lines and provide a unifying view of how BRAF-mutated melanomas respond to BRAF inhibition. We hypothesize that the residual disease observed in patients after targeted therapy is composed of a significant number of idling cells. Thus, defining molecular determinants of the phenotypic landscape that idling populations occupy may lead to "targeted landscaping" therapies based on rational modification of the landscape to favor basins with greater drug susceptibility.

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.

BioNetGen 2.2: advances in rule-based modeling.

: BioNetGen is an open-source software package for rule-based modeling of complex biochemical systems. Version 2.2 of the software introduces numerous new features for both model specification and simulation. Here, we report on these additions, discussing how they facilitate the construction, simulation and analysis of larger and more complex models than previously possible. AVAILABILITY AND IMPLEMENTATION: Stable BioNetGen releases (Linux, Mac OS/X and Windows), with documentation, are available at http://bionetgen.org Source code is available at http://github.com/RuleWorld/bionetgen CONTACT: bionetgen.help@gmail.comSupplementary information: Supplementary data are available at Bioinformatics online.

Modeling for (physical) biologists: an introduction to the rule-based approach.

Models that capture the chemical kinetics of cellular regulatory networks can be specified in terms of rules for biomolecular interactions. A rule defines a generalized reaction, meaning a reaction that permits multiple reactants, each capable of participating in a characteristic transformation and each possessing certain, specified properties, which may be local, such as the state of a particular site or domain of a protein. In other words, a rule defines a transformation and the properties that reactants must possess to participate in the transformation. A rule also provides a rate law. A rule-based approach to modeling enables consideration of mechanistic details at the level of functional sites of biomolecules and provides a facile and visual means for constructing computational models, which can be analyzed to study how system-level behaviors emerge from component interactions.

Exact hybrid particle/population simulation of rule-based models of biochemical systems.

Detailed modeling and simulation of biochemical systems is complicated by the problem of combinatorial complexity, an explosion in the number of species and reactions due to myriad protein-protein interactions and post-translational modifications. Rule-based modeling overcomes this problem by representing molecules as structured objects and encoding their interactions as pattern-based rules. This greatly simplifies the process of model specification, avoiding the tedious and error prone task of manually enumerating all species and reactions that can potentially exist in a system. From a simulation perspective, rule-based models can be expanded algorithmically into fully-enumerated reaction networks and simulated using a variety of network-based simulation methods, such as ordinary differential equations or Gillespie's algorithm, provided that the network is not exceedingly large. Alternatively, rule-based models can be simulated directly using particle-based kinetic Monte Carlo methods. This "network-free" approach produces exact stochastic trajectories with a computational cost that is independent of network size. However, memory and run time costs increase with the number of particles, limiting the size of system that can be feasibly simulated. Here, we present a hybrid particle/population simulation method that combines the best attributes of both the network-based and network-free approaches. The method takes as input a rule-based model and a user-specified subset of species to treat as population variables rather than as particles. The model is then transformed by a process of "partial network expansion" into a dynamically equivalent form that can be simulated using a population-adapted network-free simulator. The transformation method has been implemented within the open-source rule-based modeling platform BioNetGen, and resulting hybrid models can be simulated using the particle-based simulator NFsim. Performance tests show that significant memory savings can be achieved using the new approach and a monetary cost analysis provides a practical measure of its utility.

MOSBIE: a tool for comparison and analysis of rule-based biochemical models.

BACKGROUND: Mechanistic models that describe the dynamical behaviors of biochemical systems are common in computational systems biology, especially in the realm of cellular signaling. The development of families of such models, either by a single research group or by different groups working within the same area, presents significant challenges that range from identifying structural similarities and differences between models to understanding how these differences affect system dynamics. RESULTS: We present the development and features of an interactive model exploration system, MOSBIE, which provides utilities for identifying similarities and differences between models within a family. Models are clustered using a custom similarity metric, and a visual interface is provided that allows a researcher to interactively compare the structures of pairs of models as well as view simulation results. CONCLUSIONS: We illustrate the usefulness of MOSBIE via two case studies in the cell signaling domain. We also present feedback provided by domain experts and discuss the benefits, as well as the limitations, of the approach.

Concepts of multi-level dynamical modelling: understanding mechanisms of squamous cell carcinoma development in Fanconi anemia.

Fanconi anemia (FA) is a rare disease (incidence of 1:300,000) primarily based on the inheritance of pathogenic variants in genes of the FA/BRCA (breast cancer) pathway. These variants ultimately reduce the functionality of different proteins involved in the repair of DNA interstrand crosslinks and DNA double-strand breaks. At birth, individuals with FA might present with typical malformations, particularly radial axis and renal malformations, as well as other physical abnormalities like skin pigmentation anomalies. During the first decade of life, FA mostly causes bone marrow failure due to reduced capacity and loss of the hematopoietic stem and progenitor cells. This often makes hematopoietic stem cell transplantation necessary, but this therapy increases the already intrinsic risk of developing squamous cell carcinoma (SCC) in early adult age. Due to the underlying genetic defect in FA, classical chemo-radiation-based treatment protocols cannot be applied. Therefore, detecting and treating the multi-step tumorigenesis process of SCC in an early stage, or even its progenitors, is the best option for prolonging the life of adult FA individuals. However, the small number of FA individuals makes classical evidence-based medicine approaches based on results from randomized clinical trials impossible. As an alternative, we introduce here the concept of multi-level dynamical modelling using large, longitudinally collected genome, proteome- and transcriptome-wide data sets from a small number of FA individuals. This mechanistic modelling approach is based on the "hallmarks of cancer in FA", which we derive from our unique database of the clinical history of over 750 FA individuals. Multi-omic data from healthy and diseased tissue samples of FA individuals are to be used for training constituent models of a multi-level tumorigenesis model, which will then be used to make experimentally testable predictions. In this way, mechanistic models facilitate not only a descriptive but also a functional understanding of SCC in FA. This approach will provide the basis for detecting signatures of SCCs at early stages and their precursors so they can be efficiently treated or even prevented, leading to a better prognosis and quality of life for the FA individual.

Involvement of Epithelial-Mesenchymal Transition Genes in Small Cell Lung Cancer Phenotypic Plasticity.

Small cell lung cancer (SCLC) is an aggressive cancer recalcitrant to treatment, arising predominantly from epithelial pulmonary neuroendocrine (NE) cells. Intratumor heterogeneity plays critical roles in SCLC disease progression, metastasis, and treatment resistance. At least five transcriptional SCLC NE and non-NE cell subtypes were recently defined by gene expression signatures. Transition from NE to non-NE cell states and cooperation between subtypes within a tumor likely contribute to SCLC progression by mechanisms of adaptation to perturbations. Therefore, gene regulatory programs distinguishing SCLC subtypes or promoting transitions are of great interest. Here, we systematically analyze the relationship between SCLC NE/non-NE transition and epithelial to mesenchymal transition (EMT)-a well-studied cellular process contributing to cancer invasiveness and resistance-using multiple transcriptome datasets from SCLC mouse tumor models, human cancer cell lines, and tumor samples. The NE SCLC-A2 subtype maps to the epithelial state. In contrast, SCLC-A and SCLC-N (NE) map to a partial mesenchymal state (M1) that is distinct from the non-NE, partial mesenchymal state (M2). The correspondence between SCLC subtypes and the EMT program paves the way for further work to understand gene regulatory mechanisms of SCLC tumor plasticity with applicability to other cancer types.

Real-time luminescence enables continuous drug-response analysis in adherent and suspension cell lines.

The drug-induced proliferation (DIP) rate is a metric of in vitro drug response that avoids inherent biases in commonly used metrics such as 72 h viability. However, DIP rate measurements rely on direct cell counting over time, a laborious task that is subject to numerous challenges, including the need to fluorescently label cells and automatically segment nuclei. Moreover, it is incredibly difficult to directly count cells and accurately measure DIP rates for cell populations in suspension. As an alternative, we use real-time luminescence measurements derived from the cellular activity of NAD(P)H oxidoreductase to efficiently estimate drug response in both adherent and suspension cell populations to a panel of known anticancer agents. For the adherent cell lines, we collect both luminescence reads and direct cell counts over time simultaneously to assess their congruency. Our results demonstrate that the proposed approach significantly speeds up data collection, avoids the need for cellular labels and image segmentation, and opens the door to significant advances in high-throughput screening of anticancer drugs.

Processes in DNA damage response from a whole-cell multi-omics perspective.

Technological advances have made it feasible to collect multi-condition multi-omic time courses of cellular response to perturbation, but the complexity of these datasets impedes discovery due to challenges in data management, analysis, visualization, and interpretation. Here, we report a whole-cell mechanistic analysis of HL-60 cellular response to bendamustine. We integrate both enrichment and network analysis to show the progression of DNA damage and programmed cell death over time in molecular, pathway, and process-level detail using an interactive analysis framework for multi-omics data. Our framework, Mechanism of Action Generator Involving Network analysis (MAGINE), automates network construction and enrichment analysis across multiple samples and platforms, which can be integrated into our annotated gene-set network to combine the strengths of networks and ontology-driven analysis. Taken together, our work demonstrates how multi-omics integration can be used to explore signaling processes at various resolutions and demonstrates multi-pathway involvement beyond the canonical bendamustine mechanism.

Accurate implementation of leaping in space: the spatial partitioned-leaping algorithm.

There is a great need for accurate and efficient computational approaches that can account for both the discrete and stochastic nature of chemical interactions as well as spatial inhomogeneities and diffusion. This is particularly true in biology and nanoscale materials science, where the common assumptions of deterministic dynamics and well-mixed reaction volumes often break down. In this article, we present a spatial version of the partitioned-leaping algorithm, a multiscale accelerated-stochastic simulation approach built upon the tau-leaping framework of Gillespie. We pay special attention to the details of the implementation, particularly as it pertains to the time step calculation procedure. We point out conceptual errors that have been made in this regard in prior implementations of spatial tau-leaping and illustrate the manifestation of these errors through practical examples. Finally, we discuss the fundamental difficulties associated with incorporating efficient exact-stochastic techniques, such as the next-subvolume method, into a spatial leaping framework and suggest possible solutions.

Systems biology markup language (SBML) level 3 package: multistate, multicomponent and multicompartment species, version 1, release 2.

Rule-based modeling is an approach that permits constructing reaction networks based on the specification of rules for molecular interactions and transformations. These rules can encompass details such as the interacting sub-molecular domains and the states and binding status of the involved components. Conceptually, fine-grained spatial information such as locations can also be provided. Through "wildcards" representing component states, entire families of molecule complexes sharing certain properties can be specified as patterns. This can significantly simplify the definition of models involving species with multiple components, multiple states, and multiple compartments. The systems biology markup language (SBML) Level 3 Multi Package Version 1 extends the SBML Level 3 Version 1 core with the "type" concept in the Species and Compartment classes. Therefore, reaction rules may contain species that can be patterns and exist in multiple locations. Multiple software tools such as Simmune and BioNetGen support this standard that thus also becomes a medium for exchanging rule-based models. This document provides the specification for Release 2 of Version 1 of the SBML Level 3 Multi package. No design changes have been made to the description of models between Release 1 and Release 2; changes are restricted to the correction of errata and the addition of clarifications.

Quantifying stochastic effects in biochemical reaction networks using partitioned leaping.

"Leaping" methods show great promise for significantly accelerating stochastic simulations of complex biochemical reaction networks. However, few practical applications of leaping have appeared in the literature to date. Here, we address this issue using the "partitioned leaping algorithm" (PLA) [L. A. Harris and P. Clancy, J. Chem. Phys. 125, 144107 (2006)], a recently introduced multiscale leaping approach. We use the PLA to investigate stochastic effects in two model biochemical reaction networks. The networks that we consider are simple enough so as to be accessible to our intuition but sufficiently complex so as to be generally representative of real biological systems. We demonstrate how the PLA allows us to quantify subtle effects of stochasticity in these systems that would be difficult to ascertain otherwise as well as not-so-subtle behaviors that would strain commonly used "exact" stochastic methods. We also illustrate bottlenecks that can hinder the approach and exemplify and discuss possible strategies for overcoming them. Overall, our aim is to aid and motivate future applications of leaping by providing stark illustrations of the benefits of the method while at the same time elucidating obstacles that are often encountered in practice.