Addressing barriers in comprehensiveness, accessibility, reusability, interoperability and reproducibility of computational models in systems biology.

Computational models are often employed in systems biology to study the dynamic behaviours of complex systems. With the rise in the number of computational models, finding ways to improve the reusability of these models and their ability to reproduce virtual experiments becomes critical. Correct and effective model annotation in community-supported and standardised formats is necessary for this improvement. Here, we present recent efforts toward a common framework for annotated, accessible, reproducible and interoperable computational models in biology, and discuss key challenges of the field.

Setting the basis of best practices and standards for curation and annotation of logical models in biology-highlights of the [BC]2 2019 CoLoMoTo/SysMod Workshop.

The fast accumulation of biological data calls for their integration, analysis and exploitation through more systematic approaches. The generation of novel, relevant hypotheses from this enormous quantity of data remains challenging. Logical models have long been used to answer a variety of questions regarding the dynamical behaviours of regulatory networks. As the number of published logical models increases, there is a pressing need for systematic model annotation, referencing and curation in community-supported and standardised formats. This article summarises the key topics and future directions of a meeting entitled 'Annotation and curation of computational models in biology', organised as part of the 2019 [BC]2 conference. The purpose of the meeting was to develop and drive forward a plan towards the standardised annotation of logical models, review and connect various ongoing projects of experts from different communities involved in the modelling and annotation of molecular biological entities, interactions, pathways and models. This article defines a roadmap towards the annotation and curation of logical models, including milestones for best practices and minimum standard requirements.

Metabolic reprogramming in Rheumatoid Arthritis Synovial Fibroblasts: A hybrid modeling approach.

Rheumatoid Arthritis (RA) is an autoimmune disease characterized by a highly invasive pannus formation consisting mainly of Synovial Fibroblasts (RASFs). This pannus leads to cartilage, bone, and soft tissue destruction in the affected joint. RASFs' activation is associated with metabolic alterations resulting from dysregulation of extracellular signals' transduction and gene regulation. Deciphering the intricate mechanisms at the origin of this metabolic reprogramming may provide significant insight into RASFs' involvement in RA's pathogenesis and offer new therapeutic strategies. Qualitative and quantitative dynamic modeling can address some of these features, but hybrid models represent a real asset in their ability to span multiple layers of biological machinery. This work presents the first hybrid RASF model: the combination of a cell-specific qualitative regulatory network with a global metabolic network. The automated framework for hybrid modeling exploits the regulatory network's trap-spaces as additional constraints on the metabolic network. Subsequent flux balance analysis allows assessment of RASFs' regulatory outcomes' impact on their metabolic flux distribution. The hybrid RASF model reproduces the experimentally observed metabolic reprogramming induced by signaling and gene regulation in RASFs. Simulations also enable further hypotheses on the potential reverse Warburg effect in RA. RASFs may undergo metabolic reprogramming to turn into "metabolic factories", producing high levels of energy-rich fuels and nutrients for neighboring demanding cells through the crucial role of HIF1.

Model learning to identify systemic regulators of the peripheral circadian clock.

MOTIVATION: Personalized medicine aims at providing patient-tailored therapeutics based on multi-type data toward improved treatment outcomes. Chronotherapy that consists in adapting drug administration to the patient's circadian rhythms may be improved by such approach. Recent clinical studies demonstrated large variability in patients' circadian coordination and optimal drug timing. Consequently, new eHealth platforms allow the monitoring of circadian biomarkers in individual patients through wearable technologies (rest-activity, body temperature), blood or salivary samples (melatonin, cortisol) and daily questionnaires (food intake, symptoms). A current clinical challenge involves designing a methodology predicting from circadian biomarkers the patient peripheral circadian clocks and associated optimal drug timing. The mammalian circadian timing system being largely conserved between mouse and humans yet with phase opposition, the study was developed using available mouse datasets. RESULTS: We investigated at the molecular scale the influence of systemic regulators (e.g. temperature, hormones) on peripheral clocks, through a model learning approach involving systems biology models based on ordinary differential equations. Using as prior knowledge our existing circadian clock model, we derived an approximation for the action of systemic regulators on the expression of three core-clock genes: Bmal1, Per2 and Rev-Erbα. These time profiles were then fitted with a population of models, based on linear regression. Best models involved a modulation of either Bmal1 or Per2 transcription most likely by temperature or nutrient exposure cycles. This agreed with biological knowledge on temperature-dependent control of Per2 transcription. The strengths of systemic regulations were found to be significantly different according to mouse sex and genetic background. AVAILABILITY AND IMPLEMENTATION: https://gitlab.inria.fr/julmarti/model-learning-mb21eccb. SUPPLEMENTARY INFORMATION: Supplementary data are available at Bioinformatics online.

Automated inference of Boolean models from molecular interaction maps using CaSQ.

MOTIVATION: Molecular interaction maps have emerged as a meaningful way of representing biological mechanisms in a comprehensive and systematic manner. However, their static nature provides limited insights to the emerging behaviour of the described biological system under different conditions. Computational modelling provides the means to study dynamic properties through in silico simulations and perturbations. We aim to bridge the gap between static and dynamic representations of biological systems with CaSQ, a software tool that infers Boolean rules based on the topology and semantics of molecular interaction maps built with CellDesigner. RESULTS: We developed CaSQ by defining conversion rules and logical formulas for inferred Boolean models according to the topology and the annotations of the starting molecular interaction maps. We used CaSQ to produce executable files of existing molecular maps that differ in size, complexity and the use of Systems Biology Graphical Notation (SBGN) standards. We also compared, where possible, the manually built logical models corresponding to a molecular map to the ones inferred by CaSQ. The tool is able to process large and complex maps built with CellDesigner (either following SBGN standards or not) and produce Boolean models in a standard output format, Systems Biology Marked Up Language-qualitative (SBML-qual), that can be further analyzed using popular modelling tools. References, annotations and layout of the CellDesigner molecular map are retained in the obtained model, facilitating interoperability and model reusability. AVAILABILITY AND IMPLEMENTATION: The present tool is available online: https://lifeware.inria.fr/∼soliman/post/casq/ and distributed as a Python package under the GNU GPLv3 license. The code can be accessed here: https://gitlab.inria.fr/soliman/casq. SUPPLEMENTARY INFORMATION: Supplementary data are available at Bioinformatics online.

Graphical requirements for multistationarity in reaction networks and their verification in BioModels.

Thomas' necessary conditions for the existence of multiple steady states in gene networks have been proved by Soulé with high generality for dynamical systems defined by differential equations. When applied to (protein) reaction networks however, those conditions do not provide information since they are trivially satisfied as soon as there is a bimolecular or a reversible reaction. Refined graphical requirements have been proposed to deal with such cases. In this paper, we present for the first time a graph rewriting algorithm for checking the refined conditions given by Soliman, and evaluate its practical performance by applying it systematically to the curated branch of the BioModels repository. This algorithm analyzes all reaction networks (of size up to 430 species) in less than 0.05 second per network, and permits to conclude to the absence of multistationarity in 160 networks over 506. The short computation times obtained in this graphical approach are in sharp contrast to the Jacobian-based symbolic computation approach. We also discuss the case of one extra graphical condition by arc rewiring that allows us to conclude on 20 more networks of this benchmark but with a high computational cost. Finally, we study with some details the case of phosphorylation cycles and MAPK signalling models which show the importance of modelling the intermediate complexations with the enzymes in order to correctly analyze the multistationarity capabilities of such biochemical reaction networks.

MetaLo: metabolic analysis of Logical models extracted from molecular interaction maps.

Molecular interaction maps (MIMs) are static graphical representations depicting complex biochemical networks that can be formalized using one of the Systems Biology Graphical Notation languages. Regardless of their extensive coverage of various biological processes, they are limited in terms of dynamic insights. However, MIMs can serve as templates for developing dynamic computational models. We present MetaLo, an open-source Python package that enables the coupling of Boolean models inferred from process description MIMs with generic core metabolic networks. MetaLo provides a framework to study the impact of signaling cascades, gene regulation processes, and metabolic flux distribution of central energy production pathways. MetaLo computes the Boolean model's asynchronous asymptotic behavior, through the identification of trap-spaces, and extracts metabolic constraints to contextualize the generic metabolic network. MetaLo is able to handle large-scale Boolean models and genome-scale metabolic models without requiring kinetic information or manual tuning. The framework behind MetaLo enables in depth analysis of the regulatory model, and may allow tackling a lack of omics data in poorly addressed biological fields to contextualize generic metabolic networks along with improper automatic reconstructions of cell- and/or disease-specific metabolic networks. MetaLo is available at https://pypi.org/project/metalo/ under the terms of the GNU General Public License v3.

A stronger necessary condition for the multistationarity of chemical reaction networks.

Biochemical reaction networks grow bigger and bigger, fed by the high-throughput data provided by biologists and bred in open repositories of models allowing merging and evolution. Nevertheless, since the available data is still very far from permitting the identification of the increasing number of kinetic parameters of such models, the necessity of structural analyses for describing the dynamics of chemical networks appears stronger every day. Using the structural information, notably from the stoichiometric matrix, of a biochemical reaction system, we state a more strict version of the famous Thomas' necessary condition for multistationarity. In particular, the obvious cases where Thomas' condition was trivially satisfied, mutual inhibition due to a multimolecular reaction and mutual activation due to a reversible reaction, can now easily be ruled out. This more strict condition shall not be seen as some version of Thomas' circuit functionality for the continuous case but rather as related and complementary to the whole domain of the structural analysis of (bio)chemical reaction systems, as pioneered by the chemical reaction network theory.

A large-scale Boolean model of the rheumatoid arthritis fibroblast-like synoviocytes predicts drug synergies in the arthritic joint.

Rheumatoid arthritis (RA) is a complex autoimmune disease with an unknown aetiology. However, rheumatoid arthritis fibroblast-like synoviocytes (RA-FLS) play a significant role in initiating and perpetuating destructive joint inflammation by expressing immuno-modulating cytokines, adhesion molecules, and matrix remodelling enzymes. In addition, RA-FLS are primary drivers of inflammation, displaying high proliferative rates and an apoptosis-resistant phenotype. Thus, RA-FLS-directed therapies could become a complementary approach to immune-directed therapies by predicting the optimal conditions that would favour RA-FLS apoptosis, limit inflammation, slow the proliferation rate and minimise bone erosion and cartilage destruction. In this paper, we present a large-scale Boolean model for RA-FLS that consists of five submodels focusing on apoptosis, cell proliferation, matrix degradation, bone erosion and inflammation. The five-phenotype-specific submodels can be simulated independently or as a global model. In silico simulations and perturbations reproduced the expected biological behaviour of the system under defined initial conditions and input values. The model was then used to mimic the effect of mono or combined therapeutic treatments and predict novel targets and drug candidates through drug repurposing analysis.

Hybrid computational modeling highlights reverse warburg effect in breast cancer-associated fibroblasts.

Cancer-associated fibroblasts (CAFs) are amongst the key players of the tumor microenvironment (TME) and are involved in cancer initiation, progression, and resistance to therapy. They exhibit aggressive phenotypes affecting extracellular matrix remodeling, angiogenesis, immune system modulation, tumor growth, and proliferation. CAFs phenotypic changes appear to be associated with metabolic alterations, notably a reverse Warburg effect that may drive fibroblasts transformation. However, its precise molecular mechanisms and regulatory drivers are still under investigation. Deciphering the reverse Warburg effect in breast CAFs may contribute to a better understanding of the interplay between TME and tumor cells, leading to new treatment strategies. In this regard, dynamic modeling approaches able to span multiple biological layers are essential to capture the emergent properties of various biological entities when complex and intertwined pathways are involved. This work presents the first hybrid large-scale computational model for breast CAFs covering major cellular signaling, gene regulation, and metabolic processes. It was generated by combining a cell- and disease-specific asynchronous Boolean model with a generic core metabolic network leveraging both data-driven and manual curation approaches. This model reproduces the experimentally observed reverse Warburg effect in breast CAFs and further identifies Hypoxia-Inducible Factor 1 (HIF-1) as its key molecular driver. Targeting HIF-1 as part of a TME-centered therapeutic strategy may prove beneficial in the treatment of breast cancer by addressing the reverse Warburg effect. Such findings in CAFs, in light of our previously published results in rheumatoid arthritis synovial fibroblasts, point to a common HIF-1-driven metabolic reprogramming of fibroblasts in breast cancer and rheumatoid arthritis.

Drug-target identification in COVID-19 disease mechanisms using computational systems biology approaches.

Introduction: The COVID-19 Disease Map project is a large-scale community effort uniting 277 scientists from 130 Institutions around the globe. We use high-quality, mechanistic content describing SARS-CoV-2-host interactions and develop interoperable bioinformatic pipelines for novel target identification and drug repurposing. Methods: Extensive community work allowed an impressive step forward in building interfaces between Systems Biology tools and platforms. Our framework can link biomolecules from omics data analysis and computational modelling to dysregulated pathways in a cell-, tissue- or patient-specific manner. Drug repurposing using text mining and AI-assisted analysis identified potential drugs, chemicals and microRNAs that could target the identified key factors. Results: Results revealed drugs already tested for anti-COVID-19 efficacy, providing a mechanistic context for their mode of action, and drugs already in clinical trials for treating other diseases, never tested against COVID-19. Discussion: The key advance is that the proposed framework is versatile and expandable, offering a significant upgrade in the arsenal for virus-host interactions and other complex pathologies.

A graphical method for reducing and relating models in systems biology.

MOTIVATION: In Systems Biology, an increasing collection of models of various biological processes is currently developed and made available in publicly accessible repositories, such as biomodels.net for instance, through common exchange formats such as SBML. To date, however, there is no general method to relate different models to each other by abstraction or reduction relationships, and this task is left to the modeler for re-using and coupling models. In mathematical biology, model reduction techniques have been studied for a long time, mainly in the case where a model exhibits different time scales, or different spatial phases, which can be analyzed separately. These techniques are however far too restrictive to be applied on a large scale in systems biology, and do not take into account abstractions other than time or phase decompositions. Our purpose here is to propose a general computational method for relating models together, by considering primarily the structure of the interactions and abstracting from their dynamics in a first step. RESULTS: We present a graph-theoretic formalism with node merge and delete operations, in which model reductions can be studied as graph matching problems. From this setting, we derive an algorithm for deciding whether there exists a reduction from one model to another, and evaluate it on the computation of the reduction relations between all SBML models of the biomodels.net repository. In particular, in the case of the numerous models of MAPK signalling, and of the circadian clock, biologically meaningful mappings between models of each class are automatically inferred from the structure of the interactions. We conclude on the generality of our graphical method, on its limits with respect to the representation of the structure of the interactions in SBML, and on some perspectives for dealing with the dynamics. AVAILABILITY: The algorithms described in this article are implemented in the open-source software modeling platform BIOCHAM available at http://contraintes.inria.fr/biocham The models used in the experiments are available from http://www.biomodels.net/.

A general computational method for robustness analysis with applications to synthetic gene networks.

MOTIVATION: Robustness is the capacity of a system to maintain a function in the face of perturbations. It is essential for the correct functioning of natural and engineered biological systems. Robustness is generally defined in an ad hoc, problem-dependent manner, thus hampering the fruitful development of a theory of biological robustness, recently advocated by Kitano. RESULTS: In this article, we propose a general definition of robustness that applies to any biological function expressible in temporal logic LTL (linear temporal logic), and to broad model classes and perturbation types. Moreover, we propose a computational approach and an implementation in BIOCHAM 2.8 for the automated estimation of the robustness of a given behavior with respect to a given set of perturbations. The applicability and biological relevance of our approach is demonstrated by testing and improving the robustness of the timed behavior of a synthetic transcriptional cascade that could be used as a biological timer for synthetic biology applications. AVAILABILITY: Version 2.8 of BIOCHAM and the transcriptional cascade model are available at http://contraintes.inria.fr/BIOCHAM/.

Dynamics of the interlocked positive feedback loops explaining the robust epigenetic switching in Candida albicans.

The two element mutual activation and inhibitory positive feedback loops are a common motifs that occur in many biological systems in both isolated and interlocked form, as for example, in the cell division cycle and thymus differentiation in eukaryotes. The properties of three element interlocked positive feedback loops that embeds both mutual activation and inhibition are studied in depth for their bistable properties by performing bifurcation and stochastic simulations. Codimension one and two bifurcations reveal important properties like robustness to parameter variations and adaptability under various conditions by its ability to fine tune the threshold to a wide range of values and to maintain a wide bistable regime. Furthermore, we show that in the interlocked circuit, mutual inhibition controls the decision to switch from OFF to ON state, while mutual activation enforces the decision. This view is supported through a concrete biological example Candida albicans, a human fungal pathogen that can exist in two distinctive cell types; one in the default white state and the other in an opaque form. Stochastic switching between these two forms takes place due to the epigenetic alternation induced by the transcriptional regulators in the circuit, albeit without any rearrangement of the nuclear chromosomes. The transcriptional regulators constitute interlocked mutual activation and inhibition feedback circuits that provide adaptable threshold and wide bistable regime. These positive feedback loops are shown to be responsible for robust noise induced transitions without chattering, persistence of particular phenotypes for many generations and selective exhibition of one particular form of phenotype when mutated. Finally, we propose for synthetic biology constructs to use interlocked positive feedback loops instead of two element positive feedback loops because they are better controlled than isolated mutual activation and mutual inhibition feedback circuits.

Influence Networks Compared with Reaction Networks: Semantics, Expressivity and Attractors.

Biochemical reaction networks are one of the most widely used formalisms in systems biology to describe the molecular mechanisms of high-level cell processes. However, modellers also reason with influence diagrams to represent the positive and negative influences between molecular species and may find an influence network useful in the process of building a reaction network. In this paper, we introduce a formalism of influence networks with forces, and equip it with a hierarchy of Boolean, Petri net, stochastic and differential semantics, similarly to reaction networks with rates. We show that the expressive power of influence networks is the same as that of reaction networks under the differential semantics, but weaker under the discrete semantics. Furthermore, the hierarchy of semantics leads us to consider a (positive) Boolean semantics that cannot test the absence of a species, that we compare with the (negative) Boolean semantics with test for absence of a species in gene regulatory networks à la Thomas. We study the monotonicity properties of the positive semantics and derive from them an algorithm to compute attractors in both the positive and negative Boolean semantics. We illustrate our results on models of the literature about the p53/Mdm2 DNA damage repair system, the circadian clock, and the influence of MAPK signaling on cell-fate decision in urinary bladder cancer.

BIOCHAM: an environment for modeling biological systems and formalizing experimental knowledge.

UNLABELLED: BIOCHAM (the BIOCHemical Abstract Machine) is a software environment for modeling biochemical systems. It is based on two aspects: (1) the analysis and simulation of boolean, kinetic and stochastic models and (2) the formalization of biological properties in temporal logic. BIOCHAM provides tools and languages for describing protein networks with a simple and straightforward syntax, and for integrating biological properties into the model. It then becomes possible to analyze, query, verify and maintain the model with respect to those properties. For kinetic models, BIOCHAM can search for appropriate parameter values in order to reproduce a specific behavior observed in experiments and formalized in temporal logic. Coupled with other methods such as bifurcation diagrams, this search assists the modeler/biologist in the modeling process. AVAILABILITY: BIOCHAM (v. 2.5) is a free software available for download, with example models, at http://contraintes.inria.fr/BIOCHAM/.

Model-based investigation of the circadian clock and cell cycle coupling in mouse embryonic fibroblasts: Prediction of RevErb-α up-regulation during mitosis.

Experimental observations have put in evidence autonomous self-sustained circadian oscillators in most mammalian cells, and proved the existence of molecular links between the circadian clock and the cell cycle. Some mathematical models have also been built to assess conditions of control of the cell cycle by the circadian clock. However, recent studies in individual NIH3T3 fibroblasts have shown an unexpected acceleration of the circadian clock together with the cell cycle when the culture medium is enriched with growth factors, and the absence of such acceleration in confluent cells. In order to explain these observations, we study a possible entrainment of the circadian clock by the cell cycle through a regulation of clock genes around the mitosis phase. We develop a computational model and a formal specification of the observed behavior to investigate the conditions of entrainment in period and phase. We show that either the selective activation of RevErb-α or the selective inhibition of Bmal1 transcription during the mitosis phase, allow us to fit the experimental data on both period and phase, while a uniform inhibition of transcription during mitosis seems incompatible with the phase data. We conclude on the arguments favoring the RevErb-α up-regulation hypothesis and on some further predictions of the model.

A constraint solving approach to model reduction by tropical equilibration.

Model reduction is a central topic in systems biology and dynamical systems theory, for reducing the complexity of detailed models, finding important parameters, and developing multi-scale models for instance. While singular perturbation theory is a standard mathematical tool to analyze the different time scales of a dynamical system and decompose the system accordingly, tropical methods provide a simple algebraic framework to perform these analyses systematically in polynomial systems. The crux of these methods is in the computation of tropical equilibrations. In this paper we show that constraint-based methods, using reified constraints for expressing the equilibration conditions, make it possible to numerically solve non-linear tropical equilibration problems, out of reach of standard computation methods. We illustrate this approach first with the detailed reduction of a simple biochemical mechanism, the Michaelis-Menten enzymatic reaction model, and second, with large-scale performance figures obtained on the http://biomodels.net repository.

Invariants and other structural properties of biochemical models as a constraint satisfaction problem.

BACKGROUND: We present a way to compute the minimal semi-positive invariants of a Petri net representing a biological reaction system, as resolution of a Constraint Satisfaction Problem. The use of Petri nets to manipulate Systems Biology models and make available a variety of tools is quite old, and recently analyses based on invariant computation for biological models have become more and more frequent, for instance in the context of module decomposition. RESULTS: In our case, this analysis brings both qualitative and quantitative information on the models, in the form of conservation laws, consistency checking, etc. thanks to finite domain constraint programming. It is noticeable that some of the most recent optimizations of standard invariant computation techniques in Petri nets correspond to well-known techniques in constraint solving, like symmetry-breaking. Moreover, we show that the simple and natural encoding proposed is not only efficient but also flexible enough to encompass sub/sur-invariants, siphons/traps, etc., i.e., other Petri net structural properties that lead to supplementary insight on the dynamics of the biochemical system under study. CONCLUSIONS: A simple implementation based on GNU-Prolog's finite domain solver, and including symmetry detection and breaking, was incorporated into the BIOCHAM modelling environment and in the independent tool Nicotine. Some illustrative examples and benchmarks are provided.

A unique transformation from ordinary differential equations to reaction networks.

Many models in Systems Biology are described as a system of Ordinary Differential Equations, which allows for transient, steady-state or bifurcation analysis when kinetic information is available. Complementary structure-related qualitative analysis techniques have become increasingly popular in recent years, like qualitative model checking or pathway analysis (elementary modes, invariants, flux balance analysis, graph-based analyses, chemical organization theory, etc.). They do not rely on kinetic information but require a well-defined structure as stochastic analysis techniques equally do. In this article, we look into the structure inference problem for a model described by a system of Ordinary Differential Equations and provide conditions for the uniqueness of its solution. We describe a method to extract a structured reaction network model, represented as a bipartite multigraph, for example, a continuous Petri net (CPN), from a system of Ordinary Differential Equations (ODEs). A CPN uniquely defines an ODE, and each ODE can be transformed into a CPN. However, it is not obvious under which conditions the transformation of an ODE into a CPN is unique, that is, when a given ODE defines exactly one CPN. We provide biochemically relevant sufficient conditions under which the derived structure is unique and counterexamples showing the necessity of each condition. Our method is implemented and available; we illustrate it on some signal transduction models from the BioModels database. A prototype implementation of the method is made available to modellers at http://contraintes.inria.fr/~soliman/ode2pn.html, and the data mentioned in the "Results" section at http://contraintes.inria.fr/~soliman/ode2pn_data/. Our results yield a new recommendation for the import/export feature of tools supporting the SBML exchange format.