Quantification of short term signaling by the epidermal growth factor receptor

Leveraging public AI tools to explore systems biology resources in mathematical modeling

Predictive mathematical modeling is an essential part of systems biology and is interconnected with information management. Systems biology information is often stored in specialized formats to facilitate data storage and analysis. These formats are not designed for easy human readability and thus require specialized software to visualize and interpret results. Therefore, comprehending modeling and underlying networks and pathways is contingent on mastering systems biology tools, which is particularly challenging for users with no or little background in data science or system biology. To address this challenge, we investigated the usage of public Artificial Intelligence (AI) tools in exploring systems biology resources in mathematical modeling. We tested public AI's understanding of mathematics in models, related systems biology data, and the complexity of model structures. Our approach can enhance the accessibility of systems biology for non-system biologists and help them understand systems biology without a deep learning curve.

Mathematical basis and toolchain for hierarchical optimization of biochemical networks

Biological signalling systems are complex, and efforts to build mechanistic models must confront a huge parameter space, indirect and sparse data, and frequently encounter multiscale and multiphysics phenomena. We present HOSS, a framework for Hierarchical Optimization of Systems Simulations, to address such problems. HOSS operates by breaking down extensive systems models into individual pathway blocks organized in a nested hierarchy. At the first level, dependencies are solely on signalling inputs, and subsequent levels rely only on the preceding ones. We demonstrate that each independent pathway in every level can be efficiently optimized. Once optimized, its parameters are held constant while the pathway serves as input for succeeding levels. We develop an algorithmic approach to identify the necessary nested hierarchies for the application of HOSS in any given biochemical network. Furthermore, we devise two parallelizable variants that generate numerous model instances using stochastic scrambling of parameters during initial and intermediate stages of optimization. Our results indicate that these variants produce superior models and offer an estimate of solution degeneracy. Additionally, we showcase the effectiveness of the optimization methods for both abstracted, event-based simulations and ODE-based models.

Computational Modeling and Analysis of the TGF-β-induced ERK and SMAD Pathways

TGF-β, an important cytokine that plays a key role in many diseases regulates a wide array of cellular and physiologic processes via several TGF-β-driven signaling cascades, including the SMAD and non-SMAD-driven pathways. However, the detailed mechanisms by which TGF-β induces such diverse responses remain poorly understood. In particular, compared to the SMAD-dependent pathway, SMAD-independent pathways such as the ERK/MAPK pathway, which is critical in cancer progression, are less characterized. Here, we develop an integrated mechanistic model of the TGF-β-triggered ERK activation pathway and its crosstalk with the SMAD pathway, an analysis of which demonstrates how SMAD dynamics can be significantly modulated and regulated by the ERK pathway. In particular, SMAD-mediated transcription can be altered and delayed due to expedited phosphorylation of the linker of SMAD by TGF-β-activated ERK; and enhanced ERK activity, but attenuated SMAD activity, can be achieved simultaneously by fast turnover of TGF-β receptors via lipid-rafts. Also, in silico mutations of the TGF-β pathways reveal that the dynamic characteristics of both SMAD and ERK signaling may change significantly during cancer development. Specifically, normal cells may exhibit enhanced and sustained SMAD signaling with transient ERK activation, whereas cancerous cells may produce elevated and prolonged ERK signaling with enervated SMAD activation. These distinctive differences between normal and cancerous signaling behavior provide clues concerning, and potential explanations for, the seemingly contradictory roles played by TGF-β during cancer progression. We demonstrate how crosstalk among various branch pathways of TGF-β can influence overall cellular behavior. Based on model analysis, we hypothesize that aberrant molecular alterations drive changes in the intensity and duration of SMAD and ERK signaling during cancer progression and ultimately lead to an imbalance between the SMAD and ERK pathways in favor of tumor promotion. Thus, to treat cancer patients with a genetic signature of oncogenic Ras effectively may require at least a combination therapy to restore both the expression of TGF-β receptors and the GTPase activity of Ras.

Cell-specific models reveal conformation-specific RAF inhibitor combinations that synergistically inhibit ERK signaling in pancreatic cancer cells

Pancreatic ductal adenocarcinoma (PDAC) presents significant challenges for targeted clinical interventions due to prevalent KRAS mutations, rendering PDAC resistant to RAF and MEK inhibitors (RAFi and MEKi). In addition, responses to targeted therapies vary between patients. Here, we explored the differential sensitivities of PDAC cell lines to RAFi and MEKi and developed an isogenic pair comprising the most sensitive and resistant PDAC cells. To simulate patient- or tumor-specific variations, we constructed cell-line-specific mechanistic models based on protein expression profiling and differential properties of KRAS mutants. These models predicted synergy between two RAFi with different conformation specificity (type I½ and type II RAFi) in inhibiting phospho-ERK (ppERK) and reducing PDAC cell viability. This synergy was experimentally validated across all four studied PDAC cell lines. Our findings underscore the need for combination approaches to inhibit the ERK pathway in PDAC.

Mathematical Modeling and Inference of Epidermal Growth Factor-Induced Mitogen-Activated Protein Kinase Cell Signaling Pathways

The mitogen-activated protein kinase (MAPK) pathway is an important intracellular signaling cascade that plays a key role in various cellular processes. Understanding the regulatory mechanisms of this pathway is essential for developing effective interventions and targeted therapies for related diseases. Recent advances in single-cell proteomic technologies have provided unprecedented opportunities to investigate the heterogeneity and noise within complex, multi-signaling networks across diverse cells and cell types. Mathematical modeling has become a powerful interdisciplinary tool that bridges mathematics and experimental biology, providing valuable insights into these intricate cellular processes. In addition, statistical methods have been developed to infer pathway topologies and estimate unknown parameters within dynamic models. This review presents a comprehensive analysis of how mathematical modeling of the MAPK pathway deepens our understanding of its regulatory mechanisms, enhances the prediction of system behavior, and informs experimental research, with a particular focus on recent advances in modeling and inference using single-cell proteomic data.

Cell State Transition Models Stratify Breast Cancer Cell Phenotypes and Reveal New Therapeutic Targets

Understanding signaling patterns of transformation and controlling cell phenotypes is a challenge of current biology. Here we applied a cell State Transition Assessment and Regulation (cSTAR) approach to a perturbation dataset of single cell phosphoproteomic patterns of multiple breast cancer (BC) and normal breast tissue-derived cell lines. Following a separation of luminal, basal, and normal cell states, we identified signaling nodes within core control networks, delineated causal connections, and determined the primary drivers underlying oncogenic transformation and transitions across distinct BC subtypes. Whereas cell lines within the same BC subtype have different mutational and expression profiles, the architecture of the core network was similar for all luminal BC cells, and mTOR was a main oncogenic driver. In contrast, core networks of basal BC were heterogeneous and segregated into roughly four major subclasses with distinct oncogenic and BC subtype drivers. Likewise, normal breast tissue cells were separated into two different subclasses. Based on the data and quantified network topologies, we derived mechanistic cSTAR models that serve as digital cell twins and allow the deliberate control of cell movements within a Waddington landscape across different cell states. These cSTAR models suggested strategies of normalizing phosphorylation networks of BC cell lines using small molecule inhibitors.

Extending BioMASS to construct mathematical models from external knowledge

Motivation

Mechanistic modeling based on ordinary differential equations has led to numerous findings in systems biology by integrating prior knowledge and experimental data. However, the manual curation of knowledge necessary when constructing models poses a bottleneck. As the speed of knowledge accumulation continues to grow, there is a demand for a scalable means of constructing executable models.

Results

We previously introduced BioMASS-an open-source, Python-based framework-to construct, simulate, and analyze mechanistic models of signaling networks. With one of its features, Text2Model, BioMASS allows users to define models in a natural language-like format, thereby facilitating the construction of large-scale models. We demonstrate that Text2Model can serve as a tool for integrating external knowledge for mathematical modeling by generating Text2Model files from a pathway database or through the use of a large language model, and simulating its dynamics through BioMASS. Our findings reveal the tool's capabilities to encourage exploration from prior knowledge and pave the way for a fully data-driven approach to constructing mathematical models.

Availability and implementation

The code and documentation for BioMASS are available at https://github.com/biomass-dev/biomass and https://biomass-core.readthedocs.io, respectively. The code used in this article are available at https://github.com/okadalabipr/text2model-from-knowledge.

Index analysis: An approach to understand signal transduction with application to the EGFR signalling pathway

In systems biology and pharmacology, large-scale kinetic models are used to study the dynamic response of a system to a specific input or stimulus. While in many applications, a deeper understanding of the input-response behaviour is highly desirable, it is often hindered by the large number of molecular species and the complexity of the interactions. An approach that identifies key molecular species for a given input-response relationship and characterises dynamic properties of states is therefore highly desirable. We introduce the concept of index analysis; it is based on different time- and state-dependent quantities (indices) to identify important dynamic characteristics of molecular species. All indices are defined for a specific pair of input and response variables as well as for a specific magnitude of the input. In application to a large-scale kinetic model of the EGFR signalling cascade, we identified different phases of signal transduction, the peculiar role of Phosphatase3 during signal activation and Ras recycling during signal onset. In addition, we discuss the challenges and pitfalls of interpreting the relevance of molecular species based on knock-out simulation studies, and provide an alternative view on conflicting results on the importance of parallel EGFR downstream pathways. Beyond the applications in model interpretation, index analysis is envisioned to be a valuable tool in model reduction.

A guide to ERK dynamics, part 1: mechanisms and models

Extracellular signal-regulated kinase (ERK) has long been studied as a key driver of both essential cellular processes and disease. A persistent question has been how this single pathway is able to direct multiple cell behaviors, including growth, proliferation, and death. Modern biosensor studies have revealed that the temporal pattern of ERK activity is highly variable and heterogeneous, and critically, that these dynamic differences modulate cell fate. This two-part review discusses the current understanding of dynamic activity in the ERK pathway, how it regulates cellular decisions, and how these cell fates lead to tissue regulation and pathology. In part 1, we cover the optogenetic and live-cell imaging technologies that first revealed the dynamic nature of ERK, as well as current challenges in biosensor data analysis. We also discuss advances in mathematical models for the mechanisms of ERK dynamics, including receptor-level regulation, negative feedback, cooperativity, and paracrine signaling. While hurdles still remain, it is clear that higher temporal and spatial resolution provide mechanistic insights into pathway circuitry. Exciting new algorithms and advanced computational tools enable quantitative measurements of single-cell ERK activation, which in turn inform better models of pathway behavior. However, the fact that current models still cannot fully recapitulate the diversity of ERK responses calls for a deeper understanding of network structure and signal transduction in general.

Multiscale modeling of drug resistance in glioblastoma with gene mutations and angiogenesis

Drug resistance is a prominent impediment to the efficacy of targeted therapies across various cancer types, including glioblastoma (GBM). However, comprehending the intricate intracellular and extracellular mechanisms underlying drug resistance remains elusive. Empirical investigations have elucidated that genetic aberrations, such as gene mutations, along with microenvironmental adaptation, notably angiogenesis, act as pivotal drivers of tumor progression and drug resistance. Nonetheless, mathematical models frequently compartmentalize these factors in isolation. In this study, we present a multiscale agent-based model of GBM, encompassing cellular dynamics, intricate signaling pathways, gene mutations, angiogenesis, and therapeutic interventions. This integrative framework facilitates an exploration of the interplay between genetic mutations and the vascular microenvironment in shaping the dynamic evolution of tumors during treatment with tyrosine kinase inhibitor. Our simulations unveil that mutations influencing the migration and proliferation of tumor cells expedite the emergence of phenotype heterogeneity, thereby exacerbating tumor invasion under both treated and untreated conditions. Moreover, angiogenesis proximate to the tumor fosters a protumoral milieu, augmenting mutation-induced drug resistance by increasing the survival rate of tumor cells. Collectively, our findings underscore the dual roles of intrinsic genetic mutations and extrinsic microenvironmental adaptations in steering tumor growth and drug resistance. Finally, we substantiate our model predictions concerning the impact of gene mutations and angiogenesis on the responsiveness of targeted therapies by integrating single-cell RNA-seq, spatial transcriptomics, bulk RNA-seq, and clinical data from GBM patients. The multidimensional approach enhances our understanding of the complexities governing drug resistance in glioma and offers insights into potential therapeutic strategies.

PI3K signaling through a biochemical systems lens

Following 3 decades of extensive research into PI3K signaling, it is now evidently clear that the underlying network does not equate to a simple ON/OFF switch. This is best illustrated by the multifaceted nature of the many diseases associated with aberrant PI3K signaling, including common cancers, metabolic disease, and rare developmental disorders. However, we are still far from a complete understanding of the fundamental control principles that govern the numerous phenotypic outputs that are elicited by activation of this well-characterized biochemical signaling network, downstream of an equally diverse set of extrinsic inputs. At its core, this is a question on the role of PI3K signaling in cellular information processing and decision making. Here, we review the determinants of accurate encoding and decoding of growth factor signals and discuss outstanding questions in the PI3K signal relay network. We emphasize the importance of quantitative biochemistry, in close integration with advances in single-cell time-resolved signaling measurements and mathematical modeling.

Quantitative modeling of EGF receptor ligand discrimination via internalization proofreading

The epidermal growth factor receptor (EGFR) is a central regulator of cell physiology that is stimulated by multiple distinct ligands. Although ligands bind to EGFR while the receptor is exposed on the plasma membrane, EGFR incorporation into endosomes following receptor internalization is an important aspect of EGFR signaling, with EGFR internalization behavior dependent upon the type of ligand bound. We develop quantitative modeling for EGFR recruitment to and internalization from clathrin domains, focusing on how internalization competes with ligand unbinding from EGFR. We develop two model versions: a kinetic model with EGFR behavior described as transitions between discrete states and a spatial model with EGFR diffusion to circular clathrin domains. We find that a combination of spatial and kinetic proofreading leads to enhanced EGFR internalization ratios in comparison to unbinding differences between ligand types. Various stages of the EGFR internalization process, including recruitment to and internalization from clathrin domains, modulate the internalization differences between receptors bound to different ligands. Our results indicate that following ligand binding, EGFR may encounter multiple clathrin domains before successful recruitment and internalization. The quantitative modeling we have developed describes competition between EGFR internalization and ligand unbinding and the resulting proofreading.

Realistic simulation of time-course measurements in systems biology

In systems biology, the analysis of complex nonlinear systems faces many methodological challenges. For the evaluation and comparison of the performances of novel and competing computational methods, one major bottleneck is the availability of realistic test problems. We present an approach for performing realistic simulation studies for analyses of time course data as they are typically measured in systems biology. Since the design of experiments in practice depends on the process of interest, our approach considers the size and the dynamics of the mathematical model which is intended to be used for the simulation study. To this end, we used 19 published systems biology models with experimental data and evaluated the relationship between model features (e.g., the size and the dynamics) and features of the measurements such as the number and type of observed quantities, the number and the selection of measurement times, and the magnitude of measurement errors. Based on these typical relationships, our novel approach enables suggestions of realistic simulation study designs in the systems biology context and the realistic generation of simulated data for any dynamic model. The approach is demonstrated on three models in detail and its performance is validated on nine models by comparing ODE integration, parameter optimization, and parameter identifiability. The presented approach enables more realistic and less biased benchmark studies and thereby constitutes an important tool for the development of novel methods for dynamic modeling.

Mechanistic model of MAPK signaling reveals how allostery and rewiring contribute to drug resistance

BRAF is prototypical of oncogenes that can be targeted therapeutically and the treatment of BRAF^V600E melanomas with RAF and MEK inhibitors results in rapid tumor regression. However, drug-induced rewiring generates a drug adapted state thought to be involved in acquired resistance and disease recurrence. In this article, we study mechanisms of adaptive rewiring in BRAF^V600E melanoma cells using an energy-based implementation of ordinary differential equation (ODE) modeling in combination with proteomic, transcriptomic and imaging data. We develop a method for causal tracing of ODE models and identify two parallel MAPK reaction channels that are differentially sensitive to RAF and MEK inhibitors due to differences in protein oligomerization and drug binding. We describe how these channels, and timescale separation between immediate-early signaling and transcriptional feedback, create a state in which the RAS-regulated MAPK channel can be activated by growth factors under conditions in which the BRAF^V600E -driven channel is fully inhibited. Further development of the approaches in this article is expected to yield a unified model of adaptive drug resistance in melanoma.

Comparison of rule- and ordinary differential equation-based dynamic model of DARPP-32 signalling network

Dynamic modelling has considerably improved our understanding of complex molecular mechanisms. Ordinary differential equations (ODEs) are the most detailed and popular approach to modelling the dynamics of molecular systems. However, their application in signalling networks, characterised by multi-state molecular complexes, can be prohibitive. Contemporary modelling methods, such as rule- based (RB) modelling, have addressed these issues. The advantages of RB modelling over ODEs have been presented and discussed in numerous reviews. In this study, we conduct a direct comparison of the time courses of a molecular system founded on the same reaction network but encoded in the two frameworks. To make such a comparison, a set of reactions that underlie an ODE model was manually encoded in the Kappa language, one of the RB implementations. A comparison of the models was performed at the level of model specification and dynamics, acquired through model simulations. In line with previous reports, we confirm that the Kappa model recapitulates the general dynamics of its ODE counterpart with minor differences. These occur when molecules have multiple sites binding the same interactor. Furthermore, activation of these molecules in the RB model is slower than in the ODE one. As reported for other molecular systems, we find that, also for the DARPP-32 reaction network, the RB representation offers a more expressive and flexible syntax that facilitates access to fine details of the model, easing model reuse. In parallel with these analyses, we report a refactored model of the DARPP-32 interaction network that can serve as a canvas for the development of more complex dynamic models to study this important molecular system.

Utility of constraints reflecting system stability on analyses for biological models

Simulating complex biological models consisting of multiple ordinary differential equations can aid in the prediction of the pharmacological/biological responses; however, they are often hampered by the availability of reliable kinetic parameters. In the present study, we aimed to discover the properties of behaviors without determining an optimal combination of kinetic parameter values (parameter set). The key idea was to collect as many parameter sets as possible. Given that many systems are biologically stable and resilient (BSR), we focused on the dynamics around the steady state and formulated objective functions for BSR by partial linear approximation of the focused region. Using the objective functions and modified global cluster Newton method, we developed an algorithm for a thorough exploration of the allowable parameter space for biological systems (TEAPS). We first applied TEAPS to the NF-κB signaling model. This system shows a damped oscillation after stimulation and seems to fit the BSR constraint. By applying TEAPS, we found several directions in parameter space which stringently determines the BSR property. In such directions, the experimentally fitted parameter values were included in the range of the obtained parameter sets. The arachidonic acid metabolic pathway model was used as a model related to pharmacological responses. The pharmacological effects of nonsteroidal anti-inflammatory drugs were simulated using the parameter sets obtained by TEAPS. The structural properties of the system were partly extracted by analyzing the distribution of the obtained parameter sets. In addition, the simulations showed inter-drug differences in prostacyclin to thromboxane A2 ratio such that aspirin treatment tends to increase the ratio, while rofecoxib treatment tends to decrease it. These trends are comparable to the clinical observations. These results on real biological models suggest that the parameter sets satisfying the BSR condition can help in finding biologically plausible parameter sets and understanding the properties of biological systems.

We propose a new method to analyze the properties of biological dynamic models, which we named TEAPS (Thorough Exploration of Allowable Parameter Space). TEAPS can thoroughly determine combinations of parameter values for ordinary differential equations with which an initial state in a certain range converges to a particular fixed point. This stable and resilient behavior is a characteristic shared with many biological systems, including metabolic systems and intracellular signaling systems. Therefore, this thorough search outlined the possible parameter space as biological systems for target models, which helps to understand the system constraints when the target systems behave dynamically. The obtained parameter space can be used as an initial space for parameter tuning. For models that include a large number of parameters, the parameter space to be searched in the parameter tuning process is too large; therefore, narrowing down the space by TEAPS potentially contributes to the analysis of the dynamics of complicated biological models. Thus, our approach can partly overcome the current problem in parameter tuning and can advance the computational dynamic analyses of biological systems.

Control of cell state transitions

Understanding cell state transitions and purposefully controlling them is a longstanding challenge in biology. Here we present cell state transition assessment and regulation (cSTAR), an approach for mapping cell states, modelling transitions between them and predicting targeted interventions to convert cell fate decisions. cSTAR uses omics data as input, classifies cell states, and develops a workflow that transforms the input data into mechanistic models that identify a core signalling network, which controls cell fate transitions by influencing whole-cell networks. By integrating signalling and phenotypic data, cSTAR models how cells manoeuvre in Waddington's landscape1 and make decisions about which cell fate to adopt. Notably, cSTAR devises interventions to control the movement of cells in Waddington's landscape. Testing cSTAR in a cellular model of differentiation and proliferation shows a high correlation between quantitative predictions and experimental data. Applying cSTAR to different types of perturbation and omics datasets, including single-cell data, demonstrates its flexibility and scalability and provides new biological insights. The ability of cSTAR to identify targeted perturbations that interconvert cell fates will enable designer approaches for manipulating cellular development pathways and mechanistically underpinned therapeutic interventions.

Stochastic model of ERK-mediated progesterone receptor translocation, clustering and transcriptional activity

Progesterone receptor (PR) transcriptional activity is a key factor in the differentiation of the uterine endometrium. By consequence, progestin has been identified as an important treatment modality for endometrial cancer. PR transcriptional activity is controlled by extracellular-signal-regulated kinase (ERK) mediated phosphorylation, downstream of growth factor receptors such as EGFR. However, phosphorylation of PR also targets it for ubiquitination and destruction in the proteasome. Quantitative studies of these opposing roles are much needed toward validation of potential new progestin-based therapeutics. In this work, we propose a spatial stochastic model to study the effects of the opposing roles for PR phosphorylation on the levels of active transcription factor. Our numerical simulations confirm earlier in vitro experiments in endometrial cancer cell lines, identifying clustering as a mechanism that amplifies the ability of progesterone receptors to influence gene transcription. We additionally show the usefulness of a statistical method we developed to quantify and control variations in stochastic simulations in general biochemical systems, assisting modelers in defining minimal but meaningful numbers of simulations while guaranteeing outputs remain within a pre-defined confidence level.

A scalable, open-source implementation of a large-scale mechanistic model for single cell proliferation and death signaling

Mechanistic models of how single cells respond to different perturbations can help integrate disparate big data sets or predict response to varied drug combinations. However, the construction and simulation of such models have proved challenging. Here, we developed a python-based model creation and simulation pipeline that converts a few structured text files into an SBML standard and is high-performance- and cloud-computing ready. We applied this pipeline to our large-scale, mechanistic pan-cancer signaling model (named SPARCED) and demonstrate it by adding an IFNγ pathway submodel. We then investigated whether a putative crosstalk mechanism could be consistent with experimental observations from the LINCS MCF10A Data Cube that IFNγ acts as an anti-proliferative factor. The analyses suggested this observation can be explained by IFNγ-induced SOCS1 sequestering activated EGF receptors. This work forms a foundational recipe for increased mechanistic model-based data integration on a single-cell level, an important building block for clinically-predictive mechanistic models.

From Seeing to Simulating: A Survey of Imaging Techniques and Spatially-Resolved Data for Developing Multiscale Computational Models of Liver Regeneration

Liver regeneration, which leads to the re-establishment of organ mass, follows a specifically organized set of biological processes acting on various time and length scales. Computational models of liver regeneration largely focused on incorporating molecular and signaling detail have been developed by multiple research groups in the recent years. These modeling efforts have supported a synthesis of disparate experimental results at the molecular scale. Incorporation of tissue and organ scale data using noninvasive imaging methods can extend these computational models towards a comprehensive accounting of multiscale dynamics of liver regeneration. For instance, microscopy-based imaging methods provide detailed histological information at the tissue and cellular scales. Noninvasive imaging methods such as ultrasound, computed tomography and magnetic resonance imaging provide morphological and physiological features including volumetric measures over time. In this review, we discuss multiple imaging modalities capable of informing computational models of liver regeneration at the organ-, tissue- and cellular level. Additionally, we discuss available software and algorithms, which aid in the analysis and integration of imaging data into computational models. Such models can be generated or tuned for an individual patient with liver disease. Progress towards integrated multiscale models of liver regeneration can aid in prognostic tool development for treating liver disease.

A semantics, energy-based approach to automate biomodel composition

Hierarchical modelling is essential to achieving complex, large-scale models. However, not all modelling schemes support hierarchical composition, and correctly mapping points of connection between models requires comprehensive knowledge of each model's components and assumptions. To address these challenges in integrating biosimulation models, we propose an approach to automatically and confidently compose biosimulation models. The approach uses bond graphs to combine aspects of physical and thermodynamics-based modelling with biological semantics. We improved on existing approaches by using semantic annotations to automate the recognition of common components. The approach is illustrated by coupling a model of the Ras-MAPK cascade to a model of the upstream activation of EGFR. Through this methodology, we aim to assist researchers and modellers in readily having access to more comprehensive biological systems models.

OncoboxPD: human 51 672 molecular pathways database with tools for activity calculating and visualization

OncoboxPD (Oncobox pathway databank) available at https://open.oncobox.com is the collection of 51 672 uniformly processed human molecular pathways. Superposition of all pathways formed interactome graph of protein–protein interactions and metabolic reactions containing 361 654 interactions and 64 095 molecular participants. Pathways are uniformly classified by biological processes, and each pathway node is algorithmically functionally annotated by specific activator/repressor role. This enables online calculation of statistically supported pathway activation levels (PALs) with the built-in bioinformatic tool using custom RNA/protein expression profiles. Each pathway can be visualized as static or dynamic graph, where vertices are molecules participating in a pathway and edges are interactions or reactions between them. Differentially expressed nodes in a pathway can be visualized in two-color mode with user-defined color scale. For every comparison, OncoboxPD also generates a graph summarizing top up- and downregulated pathways.

Signalling dynamics, cell decisions, and homeostatic control in health and disease

Cell signalling engenders cells with the capability to receive and process information from the intracellular and extracellular environments, trigger and execute biological responses, and communicate with each other. Ultimately, cell signalling is responsible for maintaining homeostasis at the cellular, tissue and systemic level. For this reason, cell signalling is a topic of intense research efforts aimed to elucidate how cells coordinate transitions between states in developing and adult organisms in physiological and pathological conditions. Here, we review current knowledge of how cell signalling operates at multiple spatial and temporal scales, focusing on how single-cell analytical techniques reveal mechanisms underpinning cell-to-cell variability, signalling plasticity, and collective cellular responses.

A synthetic gene circuit for imaging-free detection of signaling pulses

Cells employ intracellular signaling pathways to sense and respond to changes in their external environment. In recent years, live-cell biosensors have revealed complex pulsatile dynamics in many pathways, but studies of these signaling dynamics are limited by the necessity of live-cell imaging at high spatiotemporal resolution. Here, we describe an approach to infer pulsatile signaling dynamics from a single measurement in fixed cells using a pulse-detecting gene circuit. We computationally screened for circuits with the capability to selectively detect signaling pulses, revealing an incoherent feedforward topology that robustly performs this computation. We implemented the motif experimentally for the Erk signaling pathway using a single engineered transcription factor and fluorescent protein reporter. Our "recorder of Erk activity dynamics" (READer) responds sensitively to spontaneous and stimulus-driven Erk pulses. READer circuits open the door to permanently labeling transient, dynamic cell populations to elucidate the mechanistic underpinnings and biological consequences of signaling dynamics.

Origin of diverse phosphorylation patterns in the ERBB system

Intercellular signals induce various cellular responses, including growth, proliferation, and differentiation, via the dynamic processes of signal transduction pathways. For cell fate decisions, ligand-binding induces the phosphorylation of ERBB receptors, which in turn activate downstream molecules. The ERBB family includes four subtypes, which diverged through two gene duplications from a common ancestor. Differences in the expression patterns of the subtypes have been reported between different organs in the human body. However, how these different expression properties influence the diverse phosphorylation levels of ERBB proteins is not well understood. Here we study the origin of the phosphorylation responses by experimental and mathematical analyses. The experimental measurements clarified that the phosphorylation levels heavily depend on the ERBB expression profiles. We developed a mathematical model consisting of the four subtypes as monomers, homodimers, and heterodimers and estimated the rate constants governing the phosphorylation responses from the experimental data. To understand the origin of the diversity, we analyzed the effects of the expression levels and reaction rates of the ERBB subtypes on the diversity. The difference in phosphorylation rates between ERBB subtypes showed a much greater contribution to the diversity than did the dimerization rates. This result implies that divergent evolution in phosphorylation reactions rather than in dimerization reactions after whole genome duplications was essential for increasing the diversity of the phosphorylation responses.

Scalable reaction network modeling with automatic validation of consistency in Event-B

Constructing a large biological model is a difficult, error-prone process. Small errors in writing a part of the model cascade to the system level and their sources are difficult to trace back. In this paper we extend a recent approach based on Event-B, a state-based formal method with refinement as its central ingredient, allowing us to validate for model consistency step-by-step in an automated way. We demonstrate this approach on a model of the heat shock response in eukaryotes and its scalability on a model of the \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mathsf {ErbB}$$\end{document}ErbB signaling pathway. All consistency properties of the model were proved automatically with computer support.

Mini-batch optimization enables training of ODE models on large-scale datasets

Quantitative dynamic models are widely used to study cellular signal processing. A critical step in modelling is the estimation of unknown model parameters from experimental data. As model sizes and datasets are steadily growing, established parameter optimization approaches for mechanistic models become computationally extremely challenging. Mini-batch optimization methods, as employed in deep learning, have better scaling properties. In this work, we adapt, apply, and benchmark mini-batch optimization for ordinary differential equation (ODE) models, thereby establishing a direct link between dynamic modelling and machine learning. On our main application example, a large-scale model of cancer signaling, we benchmark mini-batch optimization against established methods, achieving better optimization results and reducing computation by more than an order of magnitude. We expect that our work will serve as a first step towards mini-batch optimization tailored to ODE models and enable modelling of even larger and more complex systems than what is currently possible.

Relationship Between Dimensionality and Convergence of Optimization Algorithms: A Comparison Between Data-Driven Normalization and Scaling Factor-Based Methods Using PEPSSBI

Ordinary differential equation models are used to represent intracellular signaling pathways in silico, aiding and guiding biological experiments to elucidate intracellular regulation. To construct such quantitative and predictive models of intracellular interactions, unknown parameters need to be estimated. Most of biological data are expressed in relative or arbitrary units, raising the question of how to compare model simulations with data. It has recently been shown that for models with large number of unknown parameters, fitting algorithms using a data-driven normalization of the simulations (DNS) performs best in terms of the convergence time and parameter identifiability. DNS approach compares model simulations and corresponding data both normalized by the same normalization procedure, without requiring additional parameters to be estimated, as necessary for widely used scaling factor-based methods. However, currently there is no parameter estimation software that directly supports DNS. In this chapter, we show how to apply DNS to dynamic models of systems and synthetic biology using PEPSSBI (Parameter Estimation Pipeline for Systems and Synthetic Biology). PEPSSBI is the first software that supports DNS, through algorithmically supported data normalization and objective function construction. PEPSSBI also supports model import using SBML and repeated parameter estimation runs executed in parallel either on a personal computer or a multi-CPU cluster.

Can Systems Biology Advance Clinical Precision Oncology?

Precision oncology is perceived as a way forward to treat individual cancer patients. However, knowing particular cancer mutations is not enough for optimal therapeutic treatment, because cancer genotype-phenotype relationships are nonlinear and dynamic. Systems biology studies the biological processes at the systems' level, using an array of techniques, ranging from statistical methods to network reconstruction and analysis, to mathematical modeling. Its goal is to reconstruct the complex and often counterintuitive dynamic behavior of biological systems and quantitatively predict their responses to environmental perturbations. In this paper, we review the impact of systems biology on precision oncology. We show examples of how the analysis of signal transduction networks allows to dissect resistance to targeted therapies and inform the choice of combinations of targeted drugs based on tumor molecular alterations. Patient-specific biomarkers based on dynamical models of signaling networks can have a greater prognostic value than conventional biomarkers. These examples support systems biology models as valuable tools to advance clinical and translational oncological research.

Multiscale computational modeling of cancer growth using features derived from microCT images

Advances in medical imaging technologies now allow noninvasive image acquisition from individual patients at high spatiotemporal resolutions. A relatively new effort of predictive oncology is to develop a paradigm for forecasting the future status of an individual tumor given initial conditions and an appropriate mathematical model. The objective of this study was to introduce a comprehensive multiscale computational method to predict cancer and microvascular network growth patterns. A rectangular lattice-based model was designed so different evolutionary scenarios could be simulated and for predicting the impact of diffusible factors on tumor morphology and size. Further, the model allows prediction-based simulation of cell and microvascular behavior. Like a single cell, each agent is fully realized within the model and interactions are governed in part by machine learning methods. This multiscale computational model was developed and incorporated input information from in vivo microscale computed tomography (microCT) images acquired from breast cancer-bearing mice. It was found that as the difference between expansion of the cancer cell population and microvascular network increases, cells undergo proliferation and migration with a greater probability compared to other phenotypes. Overall, multiscale computational model agreed with both theoretical expectations and experimental findings (microCT images) not used during model training.

Long noncoding RNA SNHG1 promotes TERT expression by sponging miR-18b-5p in breast cancer

BACKGROUND

Long noncoding RNA (lncRNA) small nucleolar RNA host gene 1 (SNHG1) plays a positive role in the progression of human malignant tumors. However, the molecular mechanism of SNHG1 remains elusive in breast cancer.

RESULTS

LncRNA SNHG1 was upregulated and had a positive relationship with poor prognosis according to bioinformatics analysis in pan-cancer including breast cancer. Silencing SNHG1 inhibited tumorigenesis in breast cancer both in vitro and in vivo. Mechanistically, SNHG1 functioned as a competing endogenous RNA (ceRNA) to promote TERT expression by sponging miR-18b-5p in breast cancer. miR-18b-5p acted as a tumor repressor in breast cancer. Moreover, the combination of SNHG1 knockdown and TERT inhibitor administration showed a synergistic inhibitory effect on breast cancer growth in vivo. Finally, E2F1 as a transcription factor, binding to SNHG1 promoter and enhanced SNHG1 transcription in breast cancer.

CONCLUSIONS

Our results provide a comprehensive understanding of the oncogenic mechanism of lncRNA SNHG1 in breast cancer. Importantly, we identified a novel E2F1-SNHG1-miR-18b-5p-TERT axis, which may be a potential therapeutic target for breast cancer. Our results also provided a potential treatment for breast cancer when knockdown SNHG1 and TERT inhibitor administration simultaneously.

Using proteomic and transcriptomic data to assess activation of intracellular molecular pathways

Analysis of molecular pathway activation is the recent instrument that helps to quantize activities of various intracellular signaling, structural, DNA synthesis and repair, and biochemical processes. This may have a deep impact in fundamental research, bioindustry, and medicine. Unlike gene ontology analyses and numerous qualitative methods that can establish whether a pathway is affected in principle, the quantitative approach has the advantage of exactly measuring the extent of a pathway up/downregulation. This results in emergence of a new generation of molecular biomarkers-pathway activation levels, which reflect concentration changes of all measurable pathway components. The input data can be the high-throughput proteomic or transcriptomic profiles, and the output numbers take both positive and negative values and positively reflect overall pathway activation. Due to their nature, the pathway activation levels are more robust biomarkers compared to the individual gene products/protein levels. Here, we review the current knowledge of the quantitative gene expression interrogation methods and their applications for the molecular pathway quantization. We consider enclosed bioinformatic algorithms and their applications for solving real-world problems. Besides a plethora of applications in basic life sciences, the quantitative pathway analysis can improve molecular design and clinical investigations in pharmaceutical industry, can help finding new active biotechnological components and can significantly contribute to the progressive evolution of personalized medicine. In addition to the theoretical principles and concepts, we also propose publicly available software for the use of large-scale protein/RNA expression data to assess the human pathway activation levels.

LncRNA TINCR favors tumorigenesis via STAT3-TINCR-EGFR-feedback loop by recruiting DNMT1 and acting as a competing endogenous RNA in human breast cancer

The long noncoding RNA (lncRNA) TINCR has recently been found to be associated with the progression of human malignancies, but the molecular mechanism of TINCR action remains elusive, particularly in breast cancer. The oncogenic role of TINCR was examined in vitro and in vivo in breast cancer. Next, the interaction between TINCR, DNMT1, and miR-503-5p methylation was explored. Moreover, the mechanism by which TINCR enhances EGFR expression and downstream signaling via an RNA–RNA interaction was comprehensively investigated. Furthermore, upstream transcriptional regulation of TINCR expression by STAT3 was examined by performing chromatin immunoprecipitation. Finally, feedback signaling in the STAT3–TINCR–EGFR downstream cascade was also investigated. TINCR is upregulated in human breast cancer tissues, and TINCR knockdown suppresses tumorigenesis in vitro and in vivo. Mechanistically, TINCR recruits DNMT1 to the miR-503-5p locus promoter, which increases the methylation and suppresses the transcriptional expression of miR-503-5p. Furthermore, TINCR also functions as a competing endogenous RNA to upregulate EGFR expression by sponging miR-503-5p. In addition, TINCR stimulates JAK2–STAT3 signaling downstream from EGFR, and STAT3 reciprocally enhances the transcriptional expression of TINCR. Our findings broaden the current understanding of the diverse manners in which TINCR functions in cancer biology. The newly identified STAT3–TINCR–EGFR-feedback loop could serve as a potential therapeutic target for human cancer.

Impact of explicit area scaling on kinetic models involving multiple compartments

BACKGROUND

Computational modelling of cell biological processes is a frequently used technique to analyse the underlying mechanisms and to generally understand the behaviour of these processes in the context of a pathway, network or even the whole cell. The most common technique in this context is the usage of ordinary differential equations that describe the kinetics of the relevant processes in mechanistic detail. Here, it is usually assumed that the content of the cell is well-stirred and thus homogeneous - which is of course an over-simplification, but often worked in the past. However, many processes happen at membranes and thus not in 3D, but in 2D. The scaling of the rates of these processes poses a special problem, if volumes of compartments are changed. They will typically scale with an area, but not with the volume of the involved compartment. However, commonly, this is neglected when setting up models and/or volume scaling also sometimes automatically happens when using modelling software in the field.

RESULTS

Here, we investigate generic as well as specific, realistic cases to find out, how strong the impact of the wrong scaling is for the outcome of simulations. We show that the importance of correct area scaling depends on the architecture of the reaction site and its changes upon volume alterations and it is hard to foresee, if it has a significant impact or not just by looking at the original model set-up. Moreover, scaled rates might exhibit more or less control over the behaviour of the system and therefore, accordingly, incorrect scaling will have more or less influence.

CONCLUSIONS

Working with multi-compartment reactions requires a careful consideration of the correct scaling of the rates when changing the volumes of the involved compartments. The error following incorrect scaling - often done by scaling with the volume of the respective compartments can lead to significant aberrations of model behaviour.

To be more precise: the role of intracellular trafficking in development and pattern formation

Living cells interpret a variety of signals in different contexts to elucidate functional responses. While the understanding of signalling molecules, their respective receptors and response at the gene transcription level have been relatively well-explored, how exactly does a single cell interpret a plethora of time-varying signals? Furthermore, how their subsequent responses at the single cell level manifest in the larger context of a developing tissue is unknown. At the same time, the biophysics and chemistry of how receptors are trafficked through the complex dynamic transport network between the plasma membrane-endosome-lysosome-Golgi-endoplasmic reticulum are much more well-studied. How the intracellular organisation of the cell and inter-organellar contacts aid in orchestrating trafficking, as well as signal interpretation and modulation by the cells are beginning to be uncovered. In this review, we highlight the significant developments that have strived to integrate endosomal trafficking, signal interpretation in the context of developmental biology and relevant open questions with a few chosen examples. Furthermore, we will discuss the imaging technologies that have been developed in the recent past that have the potential to tremendously accelerate knowledge gain in this direction while shedding light on some of the many challenges.

Kinetics of receptor tyrosine kinase activation define ERK signaling dynamics

In responses to activation of receptor tyrosine kinases (RTKs), crucial cell fate decisions depend on the duration and dynamics of ERK signaling. In PC12 cells, epidermal growth factor (EGF) induces transient ERK activation that leads to cell proliferation, whereas nerve growth factor (NGF) promotes sustained ERK activation and cell differentiation. These differences have typically been assumed to reflect distinct feedback mechanisms in the Raf-MEK-ERK signaling network, with the receptors themselves acting as simple upstream inputs. We failed to confirm the expected differences in feedback type when investigating transient versus sustained signaling downstream of the EGF receptor (EGFR) and NGF receptor (TrkA). Instead, we found that ERK signaling faithfully followed RTK dynamics when receptor signaling was modulated in different ways. EGFR activation kinetics, and consequently ERK signaling dynamics, were switched from transient to sustained when receptor internalization was inhibited with drugs or mutations, or when cells expressed a chimeric receptor likely to have impaired dimerization. In addition, EGFR and ERK signaling both became more sustained when substoichiometric levels of erlotinib were added to reduce duration of EGFR kinase activation. Our results argue that RTK activation kinetics play a crucial role in determining MAP kinase cascade signaling dynamics and cell fate decisions, and that signaling outcome can be modified by activating a given RTK in different ways.

Periodic propagating waves coordinate RhoGTPase network dynamics at the leading and trailing edges during cell migration

Migrating cells need to coordinate distinct leading and trailing edge dynamics but the underlying mechanisms are unclear. Here, we combine experiments and mathematical modeling to elaborate the minimal autonomous biochemical machinery necessary and sufficient for this dynamic coordination and cell movement. RhoA activates Rac1 via DIA and inhibits Rac1 via ROCK, while Rac1 inhibits RhoA through PAK. Our data suggest that in motile, polarized cells, RhoA–ROCK interactions prevail at the rear, whereas RhoA-DIA interactions dominate at the front where Rac1/Rho oscillations drive protrusions and retractions. At the rear, high RhoA and low Rac1 activities are maintained until a wave of oscillatory GTPase activities from the cell front reaches the rear, inducing transient GTPase oscillations and RhoA activity spikes. After the rear retracts, the initial GTPase pattern resumes. Our findings show how periodic, propagating GTPase waves coordinate distinct GTPase patterns at the leading and trailing edge dynamics in moving cells.

Revealing Dynamic Protein Acetylation across Subcellular Compartments

Protein acetylation is a widespread post-translational modification implicated in many cellular processes. Recent advances in mass spectrometry have enabled the cataloging of thousands of sites throughout the cell; however, identifying regulatory acetylation marks have proven to be a daunting task. Knowledge of the kinetics and stoichiometry of site-specific acetylation is an important factor to uncover function. Here, an improved method of quantifying acetylation stoichiometry was developed and validated, providing a detailed landscape of dynamic acetylation stoichiometry within cellular compartments. The dynamic nature of site-specific acetylation in response to serum stimulation was revealed. In two distinct human cell lines, growth factor stimulation led to site-specific, temporal acetylation changes, revealing diverse kinetic profiles that clustered into several groups. Overlap of dynamic acetylation sites among two different human cell lines suggested similar regulatory control points across major cellular pathways that include splicing, translation, and protein homeostasis. Rapid increases in acetylation on protein translational machinery suggest a positive regulatory role under progrowth conditions. Finally, higher median stoichiometry was observed in cellular compartments where active acetyltransferases are well described. Data sets can be accessed through ProteomExchange via the MassIVE repository (ProteomExchange: PXD014453; MassIVE: MSV000084029).

Regulating cellular cyclic adenosine monophosphate: "Sources," "sinks," and now, "tunable valves"

A number of hormones and growth factors stimulate target cells via the second messenger pathways, which in turn regulate cellular phenotypes. Cyclic adenosine monophosphate (cAMP) is a ubiquitous second messenger that facilitates numerous signal transduction pathways; its production in cells is tightly balanced by ligand-stimulated receptors that activate adenylate cyclases (ACs), that is, "source" and by phosphodiesterases (PDEs) that hydrolyze it, that is, "sinks." Because it regulates various cellular functions, including cell growth and differentiation, gene transcription and protein expression, the cAMP signaling pathway has been exploited for the treatment of numerous human diseases. Reduction in cAMP is achieved by blocking "sources"; however, elevation in cAMP is achieved by either stimulating "source" or blocking "sinks." Here we discuss an alternative paradigm for the regulation of cellular cAMP via GIV/Girdin, the prototypical member of a family of modulators of trimeric GTPases, Guanine nucleotide Exchange Modulators (GEMs). Cells upregulate or downregulate cellular levels of GIV-GEM, which modulates cellular cAMP via spatiotemporal mechanisms distinct from the two most often targeted classes of cAMP modulators, "sources" and "sinks." A network-based compartmental model for the paradigm of GEM-facilitated cAMP signaling has recently revealed that GEMs such as GIV serve much like a "tunable valve" that cells may employ to finetune cellular levels of cAMP. Because dysregulated signaling via GIV and other GEMs has been implicated in multiple disease states, GEMs constitute a hitherto untapped class of targets that could be exploited for modulating aberrant cAMP signaling in disease states. This article is categorized under: Models of Systems Properties and Processes > Mechanistic Models Biological Mechanisms > Cell Signaling.

Multisite EGFR phosphorylation is regulated by adaptor protein abundances and dimer lifetimes

Differential epidermal growth factor receptor (EGFR) phosphorylation is thought to couple receptor activation to distinct signaling pathways. However, the molecular mechanisms responsible for biased signaling are unresolved due to a lack of insight into the phosphorylation patterns of full-length EGFR. We extended a single-molecule pull-down technique previously used to study protein-protein interactions to allow for robust measurement of receptor phosphorylation. We found that EGFR is predominantly phosphorylated at multiple sites, yet phosphorylation at specific tyrosines is variable and only a subset of receptors share phosphorylation at the same site, even with saturating ligand concentrations. We found distinct populations of receptors as soon as 1 min after ligand stimulation, indicating early diversification of function. To understand this heterogeneity, we developed a mathematical model. The model predicted that variations in phosphorylation are dependent on the abundances of signaling partners, while phosphorylation levels are dependent on dimer lifetimes. The predictions were confirmed in studies of cell lines with different expression levels of signaling partners, and in experiments comparing low- and high-affinity ligands and oncogenic EGFR mutants. These results reveal how ligand-regulated receptor dimerization dynamics and adaptor protein concentrations play critical roles in EGFR signaling.

Quantitation of Molecular Pathway Activation Using RNA Sequencing Data

Intracellular molecular pathways (IMPs) control all major events in the living cell. IMPs are considered hotspots in biomedical sciences and thousands of IMPs have been discovered for humans and model organisms. Knowledge of IMPs activation is essential for understanding biological functions and differences between the biological objects at the molecular level. Here we describe the Oncobox system for accurate quantitative scoring activities of up to several thousand molecular pathways based on high throughput molecular data. Although initially designed for gene expression and mainly RNA sequencing data, Oncobox is now also applicable for quantitative proteomics, microRNA and transcription factor binding sites mapping data. The Oncobox system includes modules of gene expression data harmonization, aggregation and comparison and a recursive algorithm for automatic annotation of molecular pathways. The universal rationale of Oncobox enables scoring of signaling, metabolic, cytoskeleton, immunity, DNA repair, and other pathways in a multitude of biological objects. The Oncobox system can be helpful to all those working in the fields of genetics, biochemistry, interactomics, and big data analytics in molecular biomedicine.

Mathematical model of hypoxia and tumor signaling interplay reveals the importance of hypoxia and cell-to-cell variability in tumor growth inhibition

Background

Human tumor is a complex tissue with multiple heterogeneous hypoxic regions and significant cell-to-cell variability. Due to the complexity of the disease, the explanation of why anticancer therapies fail cannot be attributed to intrinsic or acquired drug resistance alone. Furthermore, there are inconsistent reports of hypoxia-induced kinase activities in different cancer cell-lines, where increase, decreases, or no change has been observed. Thus, we asked, why are there widely contrasting results in kinase activity under hypoxia in different cancer cell-lines and how does hypoxia play a role in anti-cancer drug sensitivity?

Results

We took a modeling approach to address these questions by analyzing the model simulation to explain why hypoxia driven signals can have dissimilar impact on tumor growth and alter the efficacy of anti-cancer drugs. Repeated simulations with varying concentrations of biomolecules followed by decision tree analysis reveal that the highly differential effects among heterogeneous subpopulation of tumor cells could be governed by varying concentrations of just a few key biomolecules. These biomolecules include activated serine/threonine-specific protein kinases (pRAF), mitogen-activated protein kinase kinase (pMEK), protein kinase B (pAkt), or phosphoinositide-4,5-bisphosphate 3-kinase (pPI3K). Additionally, the ratio of activated extracellular signal-regulated kinases (pERK) or pAkt to its respective total was a key factor in determining the sensitivity of pERK or pAkt to hypoxia.

Conclusion

This work offers a mechanistic insight into how hypoxia can affect the efficacy of anti-cancer drug that targets tumor signaling and provides a framework to identify the types of tumor cells that are either sensitive or resistant to anti-cancer therapy.

ModelBricks-modules for reproducible modeling improving model annotation and provenance

Most computational models in biology are built and intended for "single-use"; the lack of appropriate annotation creates models where the assumptions are unknown, and model elements are not uniquely identified. Simply recreating a simulation result from a publication can be daunting; expanding models to new and more complex situations is a herculean task. As a result, new models are almost always created anew, repeating literature searches for kinetic parameters, initial conditions and modeling specifics. It is akin to building a brick house starting with a pile of clay. Here we discuss a concept for building annotated, reusable models, by starting with small well-annotated modules we call ModelBricks. Curated ModelBricks, accessible through an open database, could be used to construct new models that will inherit ModelBricks annotations and thus be easier to understand and reuse. Key features of ModelBricks include reliance on a commonly used standard language (SBML), rule-based specification describing species as a collection of uniquely identifiable molecules, association with model specific numerical parameters, and more common annotations. Physical bricks can vary substantively; likewise, to be useful the structure of ModelBricks must be highly flexible-it should encapsulate mechanisms from single reactions to multiple reactions in a complex process. Ultimately, a modeler would be able to construct large models by using multiple ModelBricks, preserving annotations and provenance of model elements, resulting in a highly annotated model. We envision the library of ModelBricks to rapidly grow from community contributions. Persistent citable references will incentivize model creators to contribute new ModelBricks.

Computational Models Accurately Predict Multi-Cell Biomarker Profiles in Inflammation and Cancer

Individual computational models of single myeloid, lymphoid, epithelial, and cancer cells were created and combined into multi-cell computational models and used to predict the collective chemokine, cytokine, and cellular biomarker profiles often seen in inflamed or cancerous tissues. Predicted chemokine and cytokine output profiles from multi-cell computational models of gingival epithelial keratinocytes (GE KER), dendritic cells (DC), and helper T lymphocytes (HTL) exposed to lipopolysaccharide (LPS) or synthetic triacylated lipopeptide (Pam3CSK4) as well as multi-cell computational models of multiple myeloma (MM) and DC were validated using the observed chemokine and cytokine responses from the same cell type combinations grown in laboratory multi-cell cultures with accuracy. Predicted and observed chemokine and cytokine responses of GE KER + DC + HTL exposed to LPS and Pam3CSK4 matched 75% (15/20, p = 0.02069) and 80% (16/20, P = 0.005909), respectively. Multi-cell computational models became ‘personalized’ when cell line-specific genomic data were included into simulations, again validated with the same cell lines grown in laboratory multi-cell cultures. Here, predicted and observed chemokine and cytokine responses of MM cells lines MM.1S and U266B1 matched 75% (3/4) and MM.1S and U266B1 inhibition of DC marker expression in co-culture matched 100% (6/6). Multi-cell computational models have the potential to identify approaches altering the predicted disease-associated output profiles, particularly as high throughput screening tools for anti-inflammatory or immuno-oncology treatments of inflamed multi-cellular tissues and the tumor microenvironment.

Coupled feedback regulation of nuclear factor of activated T-cells (NFAT) modulates activation-induced cell death of T cells

A properly functioning immune system is vital for an organism’s wellbeing. Immune tolerance is a critical feature of the immune system that allows immune cells to mount effective responses against exogenous pathogens such as viruses and bacteria, while preventing attack to self-tissues. Activation-induced cell death (AICD) in T lymphocytes, in which repeated stimulations of the T-cell receptor (TCR) lead to activation and then apoptosis of T cells, is a major mechanism for T cell homeostasis and helps maintain peripheral immune tolerance. Defects in AICD can lead to development of autoimmune diseases. Despite its importance, the regulatory mechanisms that underlie AICD remain poorly understood, particularly at an integrative network level. Here, we develop a dynamic multi-pathway model of the integrated TCR signalling network and perform model-based analysis to characterize the network-level properties of AICD. Model simulation and analysis show that amplified activation of the transcriptional factor NFAT in response to repeated TCR stimulations, a phenomenon central to AICD, is tightly modulated by a coupled positive-negative feedback mechanism. NFAT amplification is predominantly enabled by a positive feedback self-regulated by NFAT, while opposed by a NFAT-induced negative feedback via Carabin. Furthermore, model analysis predicts an optimal therapeutic window for drugs that help minimize proliferation while maximize AICD of T cells. Overall, our study provides a comprehensive mathematical model of TCR signalling and model-based analysis offers new network-level insights into the regulation of activation-induced cell death in T cells.

Local Riemannian geometry of model manifolds and its implications for practical parameter identifiability

When non-linear models are fitted to experimental data, parameter estimates can be poorly constrained albeit being identifiable in principle. This means that along certain paths in parameter space, the log-likelihood does not exceed a given statistical threshold but remains bounded. This situation, denoted as practical non-identifiability, can be detected by Monte Carlo sampling or by systematic scanning using the profile likelihood method. In contrast, any method based on a Taylor expansion of the log-likelihood around the optimum, e.g., parameter uncertainty estimation by the Fisher Information Matrix, reveals no information about the boundedness at all. In this work, we present a geometric approach, approximating the original log-likelihood by geodesic coordinates of the model manifold. The Christoffel Symbols in the geodesic equation are fixed to those obtained from second order model sensitivities at the optimum. Based on three exemplary non-linear models we show that the information about the log-likelihood bounds and flat parameter directions can already be contained in this local information. Whereas the unbounded case represented by the Fisher Information Matrix is embedded in the geometric framework as vanishing Christoffel Symbols, non-vanishing constant Christoffel Symbols prove to define prototype non-linear models featuring boundedness and flat parameter directions of the log-likelihood. Finally, we investigate if those models could allow to approximate and replace computationally expensive objective functions originating from non-linear models by a surrogate objective function in parameter estimation problems.