Comparing different ODE modelling approaches for gene regulatory networks

Oscillations in delayed positive feedback systems

Positive feedback loops exist in many biological circuits important for organismal function. In this work, we investigate how temporal delay affects the dynamics of two canonical positive feedback models. We consider models of a genetic toggle switch and a one-way switch with delay added to the feedback terms. We show that long-lasting transient oscillations exist in both models under general conditions and that the duration depends strongly on the magnitude of the delay and initial conditions. We then show the existence of long-lasting oscillations in specific biological examples: the Cdc2-Cyclin B/Wee1 system and a genetic regulatory network. Our results challenge fundamental assumptions underlying oscillatory behavior in biological systems. While generally delayed negative feedback systems are canonical in generating oscillations, we show that delayed positive feedback systems are a mechanism for generating oscillations as well.

Dynamics and bifurcations in genetic circuits with fibration symmetries

Circuit building blocks of gene regulatory networks (GRN) have been identified through the fibration symmetries of the underlying biological graph. Here, we analyse analytically six of these circuits that occur as functional and synchronous building blocks in these networks. Of these, the lock-on, toggle switch, Smolen oscillator, feed-forward fibre and Fibonacci fibre circuits occur in living organisms, notably Escherichia coli; the sixth, the repressilator, is a synthetic GRN. We consider synchronous steady states determined by a fibration symmetry (or balanced colouring) and determine analytic conditions for local bifurcation from such states, which can in principle be either steady-state or Hopf bifurcations. We identify conditions that characterize the first bifurcation, the only one that can be stable near the bifurcation point. We model the state of each gene in terms of two variables: mRNA and protein concentration. We consider all possible 'admissible' models-those compatible with the network structure-and then specialize these general results to simple models based on Hill functions and linear degradation. The results systematically classify using graph symmetries the complexity and dynamics of these circuits, which are relevant to understand the functionality of natural and synthetic cells.

Structural and practical identifiability analysis in bioengineering: a beginner's guide

Advancements in digital technology have brought modelling to the forefront in many disciplines from healthcare to architecture. Mathematical models, often represented using parametrised sets of ordinary differential equations, can be used to characterise different processes. To infer possible estimates for the unknown parameters, these models are usually calibrated using associated experimental data. Structural and practical identifiability analyses are a key component that should be assessed prior to parameter estimation. This is because identifiability analyses can provide insights as to whether or not a parameter can take on single, multiple, or even infinitely or countably many values which will ultimately have an impact on the reliability of the parameter estimates. Also, identifiability analyses can help to determine whether the data collected are sufficient or of good enough quality to truly estimate the parameters or if more data or even reparameterization of the model is necessary to proceed with the parameter estimation process. Thus, such analyses also provide an important role in terms of model design (structural identifiability analysis) and the collection of experimental data (practical identifiability analysis). Despite the popularity of using data to estimate the values of unknown parameters, structural and practical identifiability analyses of these models are often overlooked. Possible reasons for non-consideration of application of such analyses may be lack of awareness, accessibility, and usability issues, especially for more complicated models and methods of analysis. The aim of this study is to introduce and perform both structural and practical identifiability analyses in an accessible and informative manner via application to well established and commonly accepted bioengineering models. This will help to improve awareness of the importance of this stage of the modelling process and provide bioengineering researchers with an understanding of how to utilise the insights gained from such analyses in future model development.

A review on gene regulatory network reconstruction algorithms based on single cell RNA sequencing

BACKGROUND

Understanding gene regulatory networks (GRNs) is essential for unraveling the molecular mechanisms governing cellular behavior. With the advent of high-throughput transcriptome measurement technology, researchers have aimed to reverse engineer the biological systems, extracting gene regulatory rules from their outputs, which represented by gene expression data. Bulk RNA sequencing, a widely used method for measuring gene expression, has been employed for GRN reconstruction. However, it falls short in capturing dynamic changes in gene expression at the level of individual cells since it averages gene expression across mixed cell populations.

OBJECTIVE

In this review, we provide an overview of 15 GRN reconstruction tools and discuss their respective strengths and limitations, particularly in the context of single cell RNA sequencing (scRNA-seq).

METHODS

Recent advancements in scRNA-seq break new ground of GRN reconstruction. They offer snapshots of the individual cell transcriptomes and capturing dynamic changes. We emphasize how these technological breakthroughs have enhanced GRN reconstruction.

CONCLUSION

GRN reconstructors can be classified based on their requirement for cellular trajectory, which represents a dynamical cellular process including differentiation, aging, or disease progression. Benchmarking studies support the superiority of GRN reconstructors that do not require trajectory analysis in identifying regulator-target relationships. However, methods equipped with trajectory analysis demonstrate better performance in identifying key regulatory factors. In conclusion, researchers should select a suitable GRN reconstructor based on their specific research objectives.

Gene regulatory network reconstruction: harnessing the power of single-cell multi-omic data

Inferring gene regulatory networks (GRNs) is a fundamental challenge in biology that aims to unravel the complex relationships between genes and their regulators. Deciphering these networks plays a critical role in understanding the underlying regulatory crosstalk that drives many cellular processes and diseases. Recent advances in sequencing technology have led to the development of state-of-the-art GRN inference methods that exploit matched single-cell multi-omic data. By employing diverse mathematical and statistical methodologies, these methods aim to reconstruct more comprehensive and precise gene regulatory networks. In this review, we give a brief overview on the statistical and methodological foundations commonly used in GRN inference methods. We then compare and contrast the latest state-of-the-art GRN inference methods for single-cell matched multi-omics data, and discuss their assumptions, limitations and opportunities. Finally, we discuss the challenges and future directions that hold promise for further advancements in this rapidly developing field.

From time-series transcriptomics to gene regulatory networks: A review on inference methods

Inference of gene regulatory networks has been an active area of research for around 20 years, leading to the development of sophisticated inference algorithms based on a variety of assumptions and approaches. With the ever increasing demand for more accurate and powerful models, the inference problem remains of broad scientific interest. The abstract representation of biological systems through gene regulatory networks represents a powerful method to study such systems, encoding different amounts and types of information. In this review, we summarize the different types of inference algorithms specifically based on time-series transcriptomics, giving an overview of the main applications of gene regulatory networks in computational biology. This review is intended to give an updated reference of regulatory networks inference tools to biologists and researchers new to the topic and guide them in selecting the appropriate inference method that best fits their questions, aims, and experimental data.

Analysis of tumor-immune functional responses in a mathematical model of neoantigen cancer vaccines

Cancer neoantigen vaccines have emerged as a promising approach to stimulating the immune system to fight cancer. We propose a simple model including key elements of cancer-immune interactions and conduct a phase plane analysis to understand the immunological mechanisms of cancer neoantigen vaccines. Analytical results are obtained for two widely used functional forms that represent the killing rate of tumor cells by immune cells: the law of mass action (LMA) and the dePillis-Radunskaya Law (LPR). Using the LMA, our results reveal that a slowly growing tumor can escape the immune surveillance and that there is a unique periodic solution. The LPR offers richer dynamics, in which tumor elimination and uncontrolled tumor growth are both present. We show that tumor elimination requires sufficient number of initial activated T cells in relationship to the malignant cells, which lends support to using the neoantigen cancer vaccine as an adjuvant therapy after the primary tumor is surgically removed or treated using radiotherapy. We also derive a sufficient condition for uncontrolled tumor growth under the assumption of the LPR. The juxtaposition of analyses with these two different choices for the killing rate function highlights their importance on model behavior and biological implications, by which we hope to spur further theoretical and experimental work to understand mechanisms underlying different functional forms for the killing rate.

Numerical Stability and Performance of Semi-Explicit and Semi-Implicit Predictor–Corrector Methods

Semi-implicit multistep methods are an efficient tool for solving large-scale ODE systems. This recently emerged technique is based on modified Adams–Bashforth–Moulton (ABM) methods. In this paper, we introduce new semi-explicit and semi-implicit predictor–corrector methods based on the backward differentiation formula and Adams–Bashforth methods. We provide a thorough study of the numerical stability and performance of new methods and compare their stability with semi-explicit and semi-implicit Adams–Bashforth–Moulton methods and their performance with conventional linear multistep methods: Adams–Bashforth, Adams–Moulton, and BDF. The numerical stability of the investigated methods was assessed by plotting stability regions and their performances were assessed by plotting error versus CPU time plots. The mathematical developments leading to the increase in numerical stability and performance are carefully reported. The obtained results show the potential superiority of semi-explicit and semi-implicit methods over conventional linear multistep algorithms.

Identification of condition-specific regulatory mechanisms in normal and cancerous human lung tissue

Background

Lung cancer is the leading cause of cancer death in both men and women. The most common lung cancer subtype is non-small cell lung carcinoma (NSCLC) comprising about 85% of all cases. NSCLC can be further divided into three subtypes: adenocarcinoma (LUAD), squamous cell carcinoma (LUSC), and large cell lung carcinoma. Specific genetic mutations and epigenetic aberrations play an important role in the developmental transition to a specific tumor subtype. The elucidation of normal lung versus lung tumor gene expression patterns and regulatory targets yields biomarker systems that discriminate lung phenotypes (i.e., biomarkers) and provide a foundation for the discovery of normal and aberrant gene regulatory mechanisms.

Results

We built condition-specific gene co-expression networks (csGCNs) for normal lung, LUAD, and LUSC conditions. Then, we integrated normal lung tissue-specific gene regulatory networks (tsGRNs) to elucidate control-target biomarker systems for normal and cancerous lung tissue. We characterized co-expressed gene edges, possibly under common regulatory control, for relevance in lung cancer.

Conclusions

Our approach demonstrates the ability to elucidate csGCN:tsGRN merged biomarker systems based on gene expression correlation and regulation. The biomarker systems we describe can be used to classify and further describe lung specimens. Our approach is generalizable and can be used to discover and interpret complex gene expression patterns for any condition or species.

Supplementary Information

The online version contains available at 10.1186/s12864-022-08591-9.

Inference on the structure of gene regulatory networks

In this paper, we conduct theoretical analyses on inferring the structure of gene regulatory networks. Depending on the experimental method and data type, the inference problem is classified into 20 different scenarios. For each scenario, we discuss the problem that with enough data, under what assumptions, what can be inferred about the structure. For scenarios that have been covered in the literature, we provide a brief review. For scenarios that have not been covered in literature, if the structure can be inferred, we propose new mathematical inference methods and evaluate them on simulated data. Otherwise, we prove that the structure cannot be inferred.

The Interdependency and Co-Regulation of the Vitamin D and Cholesterol Metabolism

Vitamin D and cholesterol metabolism overlap significantly in the pathways that contribute to their biosynthesis. However, our understanding of their independent and co-regulation is limited. Cardiovascular disease is the leading cause of death globally and atherosclerosis, the pathology associated with elevated cholesterol, is the leading cause of cardiovascular disease. It is therefore important to understand vitamin D metabolism as a contributory factor. From the literature, we compile evidence of how these systems interact, relating the understanding of the molecular mechanisms involved to the results from observational studies. We also present the first systems biology pathway map of the joint cholesterol and vitamin D metabolisms made available using the Systems Biology Graphical Notation (SBGN) Markup Language (SBGNML). It is shown that the relationship between vitamin D supplementation, total cholesterol, and LDL-C status, and between latitude, vitamin D, and cholesterol status are consistent with our knowledge of molecular mechanisms. We also highlight the results that cannot be explained with our current knowledge of molecular mechanisms: (i) vitamin D supplementation mitigates the side-effects of statin therapy; (ii) statin therapy does not impact upon vitamin D status; and critically (iii) vitamin D supplementation does not improve cardiovascular outcomes, despite improving cardiovascular risk factors. For (iii), we present a hypothesis, based on observations in the literature, that describes how vitamin D regulates the balance between cellular and plasma cholesterol. Answering these questions will create significant opportunities for advancement in our understanding of cardiovascular health.

Sliding homoclinic bifurcations in a Lorenz-type system: Analytic proofs

Non-smooth systems can generate dynamics and bifurcations that are drastically different from their smooth counterparts. In this paper, we study such homoclinic bifurcations in a piecewise-smooth analytically tractable Lorenz-type system that was recently introduced by Belykh et al. [Chaos 29, 103108 (2019)]. Through a rigorous analysis, we demonstrate that the emergence of sliding motions leads to novel bifurcation scenarios in which bifurcations of unstable homoclinic orbits of a saddle can yield stable limit cycles. These bifurcations are in sharp contrast with their smooth analogs that can generate only unstable (saddle) dynamics. We construct a Poincaré return map that accounts for the presence of sliding motions, thereby rigorously characterizing sliding homoclinic bifurcations that destroy a chaotic Lorenz-type attractor. In particular, we derive an explicit scaling factor for period-doubling bifurcations associated with sliding multi-loop homoclinic orbits and the formation of a quasi-attractor. Our analytical results lay the foundation for the development of non-classical global bifurcation theory in non-smooth flow systems.

Modeling Temporal Biomarkers With Semiparametric Nonlinear Dynamical Systems

Dynamical systems based on differential equations are useful for modeling the temporal evolution of biomarkers. These systems can characterize the temporal patterns of biomarkers and inform the detection of interactions among biomarkers. Existing statistical methods for dynamical systems mostly target single time-course data based on a linear model or generalized additive model. Hence, they cannot adequately capture the complex interactions among biomarkers; neither can they take into account the heterogeneity between systems or subjects. in this work, we propose a semiparametric dynamical system based on multi-index models for multiple subjects time-course data. Our model accounts for between-subject heterogeneity by introducing system-level or subject-level covariates to dynamic systems, and it allows for nonlinear relationship and interaction between the combined biomarkers and the temporal rate of each biomarker. For estimation and inference, we consider a two-step procedure based on integral equations from the proposed model. We propose an algorithm that iterates between the estimation of the link function through splines and the estimation of index parameters and that allows for regularization to achieve sparsity. We prove model identifiability and derive the asymptotic properties of the estimated model parameters. A benefit of our approach is to pool information from multiple subjects to identify the interaction among biomarkers. We apply the method to analyze electroencephalogram (EEG) data for patients affected by alcohol dependence. The results reveal new insight on patients' brain activities and demonstrate differential interaction patterns in patients compared to health control subjects.

Single-cell network biology for resolving cellular heterogeneity in human diseases

Understanding cellular heterogeneity is the holy grail of biology and medicine. Cells harboring identical genomes show a wide variety of behaviors in multicellular organisms. Genetic circuits underlying cell-type identities will facilitate the understanding of the regulatory programs for differentiation and maintenance of distinct cellular states. Such a cell-type-specific gene network can be inferred from coregulatory patterns across individual cells. Conventional methods of transcriptome profiling using tissue samples provide only average signals of diverse cell types. Therefore, reconstructing gene regulatory networks for a particular cell type is not feasible with tissue-based transcriptome data. Recently, single-cell omics technology has emerged and enabled the capture of the transcriptomic landscape of every individual cell. Although single-cell gene expression studies have already opened up new avenues, network biology using single-cell transcriptome data will further accelerate our understanding of cellular heterogeneity. In this review, we provide an overview of single-cell network biology and summarize recent progress in method development for network inference from single-cell RNA sequencing (scRNA-seq) data. Then, we describe how cell-type-specific gene networks can be utilized to study regulatory programs specific to disease-associated cell types and cellular states. Moreover, with scRNA data, modeling personal or patient-specific gene networks is feasible. Therefore, we also introduce potential applications of single-cell network biology for precision medicine. We envision a rapid paradigm shift toward single-cell network analysis for systems biology in the near future.

Gene regulatory networks reconstructed from single-cell RNA sequencing datasets are allowing researchers to better understand the molecular circuits and cell states that contribute to complex human disease. Junha Cha and Insuk Lee from Yonsei University in Seoul, South Korea, review the concept of ‘single-cell network biology’, which involves using computational algorithms on genetic expression data from thousands of cells to infer functional interactions in various biological contexts. This systems biology approach to analyzing the profiles of messenger RNA in single cells is helping researchers discover new signaling pathways that could serve as disease biomarkers or therapeutic targets. In the future, patient-specific models of personal gene networks could explain why certain genetic variants affect disease risk. This research could also eventually lead to new types of individualized medical treatments.

Inference of gene regulatory networks based on nonlinear ordinary differential equations

Motivation

Gene regulatory networks (GRNs) capture the regulatory interactions between genes, resulting from the fundamental biological process of transcription and translation. In some cases, the topology of GRNs is not known, and has to be inferred from gene expression data. Most of the existing GRNs reconstruction algorithms are either applied to time-series data or steady-state data. Although time-series data include more information about the system dynamics, steady-state data imply stability of the underlying regulatory networks.

Results

In this article, we propose a method for inferring GRNs from time-series and steady-state data jointly. We make use of a non-linear ordinary differential equations framework to model dynamic gene regulation and an importance measurement strategy to infer all putative regulatory links efficiently. The proposed method is evaluated extensively on the artificial DREAM4 dataset and two real gene expression datasets of yeast and Escherichia coli. Based on public benchmark datasets, the proposed method outperforms other popular inference algorithms in terms of overall score. By comparing the performance on the datasets with different scales, the results show that our method still keeps good robustness and accuracy at a low computational complexity.

Availability and implementation

The proposed method is written in the Python language, and is available at: https://github.com/lab319/GRNs_nonlinear_ODEs

Supplementary information

Supplementary data are available at Bioinformatics online.

Reprogramming multistable monotone systems with application to cell fate control

Multistability is a key property of dynamical systems modeling cellular regulatory networks implicated in cell fate decisions, where, different stable steady states usually represent distinct cell phenotypes. Monotone network motifs are highly represented in these regulatory networks. In this paper, we leverage the properties of monotone dynamical systems to provide theoretical results that guide the selection of inputs that trigger a transition, i.e., reprogram the network, to a desired stable steady state. We first show that monotone dynamical systems with bounded trajectories admit a minimum and a maximum stable steady state. Then, we provide input choices that are guaranteed to reprogram the system to these extreme steady states. For intermediate states, we provide an input space that is guaranteed to contain an input that reprograms the system to the desired state. We then provide implementation guidelines for finite-time procedures that search this space for such an input, along with rules to prune parts of the space during search. We demonstrate these results on simulations of two recurrent regulatory network motifs: self-activation within mutual antagonism and self-activation within mutual cooperation. Our results depend uniquely on the structure of the network and are independent of specific parameter values.

A Multi-Laplacian Prior and Augmented Lagrangian Approach to the Exploratory Analysis of Time-Varying Gene and Transcriptional Regulatory Networks for Gene Microarray Data

This paper proposes a novel multi-Laplacian prior (MLP) and augmented Lagrangian method (ALM) approach for gene interactions and putative transcription factors (TFs) identification from time-course gene microarray data. It employs a non-linear time-varying auto-regressive (N-TVAR) model and the Maximum-A-Posteriori-Probability method for incorporating the multi-Laplacian prior and the continuity constraint. The MLP allows connections to/from a gene to be better preserved for putative TF identification in non-stationarity gene regulatory network as compared with conventional L₁-based penalties. Moreover, the ALM allows the resultant non-smooth L₁-based penalties to be decoupled from the remaining smooth terms, so that the former and latter can be efficiently solved using a low-complexity proximity operator and smooth optimization technique, respectively. Synthetic and real time-course gene microarray datasets are tested to evaluate the performance of the proposed method. Experimental results show that the proposed method gives better accuracy and higher computational speed than our previous work using smoothed approximation. Moreover, its performance, without the use of ChIP-chip data, is found to be highly comparable with other state-of-the-art methods integrating both ChIP-chip and gene microarray data. It suggests that the proposed method may serve as a useful exploratory tool for putative TF identification with reduced experimental cost.

An in-silico approach for discovery of microRNA-TF regulation of DISC1 interactome mediating neuronal migration

Neuronal migration constitutes an important step in corticogenesis; dysregulation of the molecular mechanisms mediating this crucial step in neurodevelopment may result in various neuropsychiatric disorders. By curating experimental data from published literature, we identified eight functional modules involving Disrupted-in-schizophrenia 1 (DISC1) and its interacting proteins that regulate neuronal migration. We then identified miRNAs and transcription factors (TFs) that form functional feedback loops and regulate gene expression of the DISC1 interactome. Using this curated data, we conducted in-silico modeling of the DISC1 interactome involved in neuronal migration and identified the proteins that either facilitate or inhibit neuronal migrational processes. We also studied the effect of perturbation of miRNAs and TFs in feedback loops on the DISC1 interactome. From these analyses, we discovered that STAT3, TCF3, and TAL1 (through feedback loop with miRNAs) play a critical role in the transcriptional control of DISC1 interactome thereby regulating neuronal migration. To the best of our knowledge, regulation of the DISC1 interactome mediating neuronal migration by these TFs has not been previously reported. These potentially important TFs can serve as targets for undertaking validation studies, which in turn can reveal the molecular processes that cause neuronal migration defects underlying neurodevelopmental disorders. This underscores the importance of the use of in-silico techniques in aiding the discovery of mechanistic evidence governing important molecular and cellular processes. The present work is one such step towards the discovery of regulatory factors of the DISC1 interactome that mediates neuronal migration.

WASABI: a dynamic iterative framework for gene regulatory network inference

BACKGROUND

Inference of gene regulatory networks from gene expression data has been a long-standing and notoriously difficult task in systems biology. Recently, single-cell transcriptomic data have been massively used for gene regulatory network inference, with both successes and limitations.

RESULTS

In the present work we propose an iterative algorithm called WASABI, dedicated to inferring a causal dynamical network from time-stamped single-cell data, which tackles some of the limitations associated with current approaches. We first introduce the concept of waves, which posits that the information provided by an external stimulus will affect genes one-by-one through a cascade, like waves spreading through a network. This concept allows us to infer the network one gene at a time, after genes have been ordered regarding their time of regulation. We then demonstrate the ability of WASABI to correctly infer small networks, which have been simulated in silico using a mechanistic model consisting of coupled piecewise-deterministic Markov processes for the proper description of gene expression at the single-cell level. We finally apply WASABI on in vitro generated data on an avian model of erythroid differentiation. The structure of the resulting gene regulatory network sheds a new light on the molecular mechanisms controlling this process. In particular, we find no evidence for hub genes and a much more distributed network structure than expected. Interestingly, we find that a majority of genes are under the direct control of the differentiation-inducing stimulus.

CONCLUSIONS

Together, these results demonstrate WASABI versatility and ability to tackle some general gene regulatory networks inference issues. It is our hope that WASABI will prove useful in helping biologists to fully exploit the power of time-stamped single-cell data.

Computational Analysis of Altering Cell Fate

The notion of reprogramming cell fate is a direct challenge to the traditional view in developmental biology that a cell's phenotypic identity is sealed after undergoing differentiation. Direct experimental evidence, beginning with the somatic cell nuclear transfer experiments of the twentieth century and culminating in the more recent breakthroughs in transdifferentiation and induced pluripotent stem cell (iPSC) reprogramming, have rewritten the rules for what is possible with cell fate transformation. Research is ongoing in the manipulation of cell fate for basic research in disease modeling, drug discovery, and clinical therapeutics. In many of these cell fate reprogramming experiments, there is often little known about the genetic and molecular changes accompanying the reprogramming process. However, gene regulatory networks (GRNs) can in some cases be implicated in the switching of phenotypes, providing a starting point for understanding the dynamic changes that accompany a given cell fate reprogramming process. In this chapter, we present a framework for computationally analyzing cell fate changes by mathematically modeling these GRNs. We provide a user guide with several tutorials of a set of techniques from dynamical systems theory that can be used to probe the intrinsic properties of GRNs as well as study their responses to external perturbations.

Extreme learning machines for reverse engineering of gene regulatory networks from expression time series

Motivation

The reconstruction of gene regulatory networks (GRNs) from genes profiles has a growing interest in bioinformatics for understanding the complex regulatory mechanisms in cellular systems. GRNs explicitly represent the cause-effect of regulation among a group of genes and its reconstruction is today a challenging computational problem. Several methods were proposed, but most of them require different input sources to provide an acceptable prediction. Thus, it is a great challenge to reconstruct a GRN only from temporal gene expression data.

Results

Extreme Learning Machine (ELM) is a new supervised neural model that has gained interest in the last years because of its higher learning rate and better performance than existing supervised models in terms of predictive power. This work proposes a novel approach for GRNs reconstruction in which ELMs are used for modeling the relationships between gene expression time series. Artificial datasets generated with the well-known benchmark tool used in DREAM competitions were used. Real datasets were used for validation of this novel proposal with well-known GRNs underlying the time series. The impact of increasing the size of GRNs was analyzed in detail for the compared methods. The results obtained confirm the superiority of the ELM approach against very recent state-of-the-art methods in the same experimental conditions.

Availability and implementation

The web demo can be found at http://sinc.unl.edu.ar/web-demo/elm-grnnminer/. The source code is available at https://sourceforge.net/projects/sourcesinc/files/elm-grnnminer.

Contact

mrubiolo@santafe-conicet.gov.ar.

Supplementary information

Supplementary data are available at Bioinformatics online.

Finite-time parametric identification for the model representing the metabolic and genetic regulatory effects of sequential aerobic respiration and anaerobic fermentation processes in Escherichia coli

Mathematical modelling applied to biological systems allows for the inferring of changes in the dynamic behaviour of organisms associated with variations in the environment. Models based on ordinary differential equations are most commonly used because of their ability to describe the mechanisms of biological systems such as transcription. The disadvantage of using this approach is that there is a large number of parameters involved and that it is difficult to obtain them experimentally. This study presents an algorithm to obtain a finite-time parameter characterization of the model used to describe changes in the metabolic behaviour of Escherichia coli associated with environmental changes. In this scheme, super-twisting algorithm was proposed to recover the derivative of all the proteins and mRNA of E. coli associated to changes in the concentration of oxygen available in the growth media. The 75 identified parameters in this study maintain the biological coherence of the system and they were estimated with no more than 20% error with respect to the real ones included in the proposed model.

An integrated framework for identification of effective and synergistic anti-cancer drug combinations

Combination drug therapy is considered a better treatment option for various diseases, such as cancer, HIV, hypertension, and infections as compared to targeted drug therapies. Combination or synergism helps to overcome drug resistance, reduction in drug toxicity and dosage. Considering the complexity and heterogeneity among cancer types, drug combination provides promising treatment strategy. Increase in drug combination data raises a challenge for developing a computational approach that can effectively predict drugs synergism. There is a need to model the combination drug screening data to predict new synergistic drug combinations for successful cancer treatment. In such a scenario, machine learning approaches can be used to alleviate the process of drugs synergy prediction. Experimental data from a single-agent or multi-agent drug screens provides feature data for model training. On the contrary, identification of effective drug combination using clinical trials is a time consuming and resource intensive task. This paper attempts to address the aforementioned challenges by developing a computational approach to effectively predict drug synergy. Single-drug efficacy is used for predicting drug synergism. Our approach obviates the need to understand the underlying drug mechanism to predict drug combination synergy. For this purpose, nine machine learning algorithms are trained. It is observed that the Random forest models, in comparison to other models, have shown significant performance. The <mml:math xmlns:mml="http://www.w3.org/1998/Math/MathML"><mml:mi>K</mml:mi></mml:math> -fold cross-validation is performed to evaluate the robustness of the best predictive model. The proposed approach is applied to mutant-BRAF melanoma and further validated using melanoma cell-lines from AstraZeneca-Sanger Drug Combination Prediction DREAM Challenge dataset.

Estimation of Dynamic Systems for Gene Regulatory Networks from Dependent Time-Course Data

Dynamic system consisting of ordinary differential equations (ODEs) is a well-known tool for describing dynamic nature of gene regulatory networks (GRNs), and the dynamic features of GRNs are usually captured through time-course gene expression data. Owing to high-throughput technologies, time-course gene expression data have complex structures such as heteroscedasticity, correlations between genes, and time dependence. Since gene experiments typically yield highly noisy data with small sample size, for a more accurate prediction of the dynamics, the complex structures should be taken into account in ODE models. Hence, this study proposes an ODE model considering such data structures and a fast and stable estimation method for the ODE parameters based on the generalized profiling approach with data smoothing techniques. The proposed method also provides statistical inference for the ODE estimator and it is applied to a zebrafish retina cell network.

Logic of two antagonizing intra-species quorum sensing systems in bacteria

Bacteria release signaling molecules into the surrounding environment and sense them when present in their proximity. Using this strategy, a cell estimates the number of neighbors in its surrounding. Upon sensing a critical number of individuals, bacteria coordinate a number of cellular processes. This density-dependent control of gene expression and physiology is called quorum sensing (QS). Quorum sensing controls a wide variety of functions in bacteria, including those related to motility, growth, virulence etc. Quorum sensing has been widely observed in bacteria while the individuals of the same species or different species compete and cooperate each other. Interestingly, many species possess more than one QS system (intra-species) and these QS systems interact each other to perform quorum sensing. Thus, several logical arrangements can be possible based on the interaction among intra-species QS systems - parallel, series, antagonizing, and agonizing. In this work, we perform simulations to understand the logic of interaction between two antagonizing intra-species QS systems. In such an interaction, one QS system gets fully expressed and the other only gets partially expressed. This is found to be dictated by the interplay between autoinducer's diffusivity and antagonizing strength. In addition, we speculate an important role of the intracellular regulators (eg. LuxR) in maintaining the uniform response among the individual cells from the different localities. We also expect the interplay between the autoinducer's diffusivity and distribution of cells in fine tuning the collective response. Interestingly, in a localized niche with a heterogeneous cell distribution, the cells are expected to perform a global quorum sensing via fully expressed QS system and a local quorum sensing via partially expressed QS system.

Steady-State-Preserving Simulation of Genetic Regulatory Systems

A novel family of exponential Runge-Kutta (expRK) methods are designed incorporating the stable steady-state structure of genetic regulatory systems. A natural and convenient approach to constructing new expRK methods on the base of traditional RK methods is provided. In the numerical integration of the one-gene, two-gene, and p53-mdm2 regulatory systems, the new expRK methods are shown to be more accurate than their prototype RK methods. Moreover, for nonstiff genetic regulatory systems, the expRK methods are more efficient than some traditional exponential RK integrators in the scientific literature.

A system of recurrent neural networks for modularising, parameterising and dynamic analysis of cell signalling networks

In this paper, we show how to extend our previously proposed novel continuous time Recurrent Neural Networks (RNN) approach that retains the advantage of continuous dynamics offered by Ordinary Differential Equations (ODE) while enabling parameter estimation through adaptation, to larger signalling networks using a modular approach. Specifically, the signalling network is decomposed into several sub-models based on important temporal events in the network. Each sub-model is represented by the proposed RNN and trained using data generated from the corresponding ODE model. Trained sub-models are assembled into a whole system RNN which is then subjected to systems dynamics and sensitivity analyses. The concept is illustrated by application to G1/S transition in cell cycle using Iwamoto et al. (2008) ODE model. We decomposed the G1/S network into 3 sub-models: (i) E2F transcription factor release; (ii) E2F and CycE positive feedback loop for elevating cyclin levels; and (iii) E2F and CycA negative feedback to degrade E2F. The trained sub-models accurately represented system dynamics and parameters were in good agreement with the ODE model. The whole system RNN however revealed couple of parameters contributing to compounding errors due to feedback and required refinement to sub-model 2. These related to the reversible reaction between CycE/CDK2 and p27, its inhibitor. The revised whole system RNN model very accurately matched dynamics of the ODE system. Local sensitivity analysis of the whole system model further revealed the most dominant influence of the above two parameters in perturbing G1/S transition, giving support to a recent hypothesis that the release of inhibitor p27 from Cyc/CDK complex triggers cell cycle stage transition. To make the model useful in a practical setting, we modified each RNN sub-model with a time relay switch to facilitate larger interval input data (≈20min) (original model used data for 30s or less) and retrained them that produced parameters and protein concentrations similar to the original RNN system. Results thus demonstrated the reliability of the proposed RNN method for modelling relatively large networks by modularisation for practical settings. Advantages of the method are its ability to represent accurate continuous system dynamics and ease of: parameter estimation through training with data from a practical setting, model analysis (40% faster than ODE), fine tuning parameters when more data are available, sub-model extension when new elements and/or interactions come to light and model expansion with addition of sub-models.

The recurrent architecture of tumour initiation, progression and drug sensitivity

Recent studies across multiple tumour types are starting to reveal a recurrent regulatory architecture in which genomic alterations cluster upstream of functional master regulator (MR) proteins, the aberrant activity of which is both necessary and sufficient to maintain tumour cell state. These proteins form small, hyperconnected and autoregulated modules (termed tumour checkpoints) that are increasingly emerging as optimal biomarkers and therapeutic targets. Crucially, as their activity is mostly dysregulated in a post-translational manner, rather than by mutations in their corresponding genes or by differential expression, the identification of MR proteins by conventional methods is challenging. In this Opinion article, we discuss novel methods for the systematic analysis of MR proteins and of the modular regulatory architecture they implement, including their use as a valuable reductionist framework to study the genetic heterogeneity of human disease and to drive key translational applications.

Inference of Gene Regulatory Networks Using Bayesian Nonparametric Regression and Topology Information

Gene regulatory networks (GRNs) play an important role in cellular systems and are important for understanding biological processes. Many algorithms have been developed to infer the GRNs. However, most algorithms only pay attention to the gene expression data but do not consider the topology information in their inference process, while incorporating this information can partially compensate for the lack of reliable expression data. Here we develop a Bayesian group lasso with spike and slab priors to perform gene selection and estimation for nonparametric models. B-spline basis functions are used to capture the nonlinear relationships flexibly and penalties are used to avoid overfitting. Further, we incorporate the topology information into the Bayesian method as a prior. We present the application of our method on DREAM3 and DREAM4 datasets and two real biological datasets. The results show that our method performs better than existing methods and the topology information prior can improve the result.

Estimation of delays in generalized asynchronous Boolean networks

A new generalized asynchronous Boolean network (GABN) model has been proposed in this paper. This continuous-time discrete-state model captures the biological reality of cellular dynamics without compromising the computational efficiency of the Boolean framework. The GABN synthesis procedure is based on the prior knowledge of the logical structure of the regulatory network, and the experimental transcriptional parameters. The novelty of the proposed methodology lies in considering different delays associated with the activation and deactivation of a particular protein (especially the transcription factors). A few illustrative examples of some well-studied network motifs have been provided to explore the scope of using the GABN model for larger networks. The GABN model of the p53-signaling pathway in response to γ-irradiation has also been simulated in the current paper to provide an indirect validation of the proposed schema.

Hydrogen Peroxide Induced Cell Death: The Major Defences Relative Roles and Consequences in E. coli

We recently developed a mathematical model for predicting reactive oxygen species (ROS) concentration and macromolecules oxidation in vivo. We constructed such a model using Escherichia coli as a model organism and a set of ordinary differential equations. In order to evaluate the major defences relative roles against hydrogen peroxide (H₂O₂), we investigated the relative contributions of the various reactions to the dynamic system and searched for approximate analytical solutions for the explicit expression of changes in H₂O₂ internal or external concentrations. Although the key actors in cell defence are enzymes and membrane, a detailed analysis shows that their involvement depends on the H₂O₂ concentration level. Actually, the impact of the membrane upon the H₂O₂ stress felt by the cell is greater when micromolar H₂O₂ is present (9-fold less H₂O₂ in the cell than out of the cell) than when millimolar H₂O₂ is present (about 2-fold less H₂O₂ in the cell than out of the cell). The ratio between maximal external H₂O₂ and internal H₂O₂ concentration also changes, reducing from 8 to 2 while external H₂O₂ concentration increases from micromolar to millimolar. This non-linear behaviour mainly occurs because of the switch in the predominant scavenger from Ahp (Alkyl Hydroperoxide Reductase) to Cat (catalase). The phenomenon changes the internal H₂O₂ maximal concentration, which surprisingly does not depend on cell density. The external H₂O₂ half-life and the cumulative internal H₂O₂ exposure do depend upon cell density. Based on these analyses and in order to introduce a concept of dose response relationship for H₂O₂-induced cell death, we developed the concepts of “maximal internal H₂O₂ concentration” and “cumulative internal H₂O₂ concentration” (e.g. the total amount of H₂O₂). We predict that cumulative internal H₂O₂ concentration is responsible for the H₂O₂-mediated death of bacterial cells.

Stable Gene Regulatory Network Modeling From Steady-State Data

Gene regulatory networks represent an abstract mapping of gene regulations in living cells. They aim to capture dependencies among molecular entities such as transcription factors, proteins and metabolites. In most applications, the regulatory network structure is unknown, and has to be reverse engineered from experimental data consisting of expression levels of the genes usually measured as messenger RNA concentrations in microarray experiments. Steady-state gene expression data are obtained from measurements of the variations in expression activity following the application of small perturbations to equilibrium states in genetic perturbation experiments. In this paper, the least absolute shrinkage and selection operator-vector autoregressive (LASSO-VAR) originally proposed for the analysis of economic time series data is adapted to include a stability constraint for the recovery of a sparse and stable regulatory network that describes data obtained from noisy perturbation experiments. The approach is applied to real experimental data obtained for the SOS pathway in Escherichia coli and the cell cycle pathway for yeast Saccharomyces cerevisiae. Significant features of this method are the ability to recover networks without inputting prior knowledge of the network topology, and the ability to be efficiently applied to large scale networks due to the convex nature of the method.

Inference of asynchronous Boolean network from biological pathways

Gene regulation is a complex process with multiple levels of interactions. In order to describe this complex dynamical system with tractable parameterization, the choice of the dynamical system model is of paramount importance. The right abstraction of the modeling scheme can reduce the complexity in the inference and intervention design, both computationally and experimentally. This article proposes an asynchronous Boolean network framework to capture the transcriptional regulation as well as the protein-protein interactions in a genetic regulatory system. The inference of asynchronous Boolean network from biological pathways information and experimental evidence are explained using an algorithm. The suitability of this paradigm for the variability of several reaction rates is also discussed. This methodology and model selection open up new research challenges in understanding gene-protein interactive system in a coherent way and can be beneficial for designing effective therapeutic intervention strategy.

Large scale gene regulatory network inference with a multi-level strategy

Transcriptional regulation is a basis of many crucial molecular processes and an accurate inference of the gene regulatory network is a helpful and essential task to understand cell functions and gain insights into biological processes of interest in systems biology.

Hydrogen peroxide induced cell death: One or two modes of action?

Imlay and Linn show that exposure of logarithmically growing Escherichia coli to hydrogen peroxide (H₂O₂) leads to two kinetically distinguishable modes of cell killing. Mode one killing is pronounced near 1 mM concentration of H₂O₂ and is caused by DNA damage, whereas mode-two killing requires higher concentration (>10 mM). The second mode seems to be essentially due to damage to all macromolecules. This phenomenon has also been observed in Fenton in vitro systems with DNA nicking caused by hydroxyl radical (HO•).

To our knowledge, there is currently no mathematical model for predicting mode one killing in vitro or in vivo after H₂O₂ exposure.

We propose a simple model, using Escherichia coli as a model organism and a set of ordinary differential equations. Using this model, we show that available iron and cell density, two factors potentially involved in ROS dynamics, play a major role in the prediction of the experimental results obtained by our team and in previous studies. Indeed the presence of the mode one killing is strongly related to those two parameters.

To our knowledge, mode-one death has not previously been explained. Imlay and Linn (Imlay and Linn, 1986) suggested that perhaps the amount of the toxic species was reduced at high concentrations of H₂O₂ because hydroxyl (or other) radicals might be quenched directly by hydrogen peroxide with the concomitant formation of superoxide anion (a less toxic species). We demonstrate (mathematically and numerically) that free available iron decrease is necessary to explain mode one killing which cannot appear without it and that H₂O₂ quenching or consumption is not responsible for mode-one death.

We are able to follow ROS concentration (particularly responsible for mode one killing) after exposure to H₂O₂. This model therefore allows us to understand two major parameters involved in the presence or not of the first killing mode.

The limiting dynamics of a bistable molecular switch with and without noise

We consider the dynamics of a population of organisms containing two mutually inhibitory gene regulatory networks, that can result in a bistable switch-like behaviour. We completely characterize their local and global dynamics in the absence of any noise, and then go on to consider the effects of either noise coming from bursting (transcription or translation), or Gaussian noise in molecular degradation rates when there is a dominant slow variable in the system. We show analytically how the steady state distribution in the population can range from a single unimodal distribution through a bimodal distribution and give the explicit analytic form for the invariant stationary density which is globally asymptotically stable. Rather remarkably, the behaviour of the stationary density with respect to the parameters characterizing the molecular behaviour of the bistable switch is qualitatively identical in the presence of noise coming from bursting as well as in the presence of Gaussian noise in the degradation rate. This implies that one cannot distinguish between either the dominant source or nature of noise based on the stationary molecular distribution in a population of cells. We finally show that the switch model with bursting but two dominant slow genes has an asymptotically stable stationary density.

Exponentially Fitted Two-Derivative Runge-Kutta Methods for Simulation of Oscillatory Genetic Regulatory Systems

Oscillation is one of the most important phenomena in the chemical reaction systems in living cells. The general purpose simulation algorithms fail to take into account this special character and produce unsatisfying results. In order to enhance the accuracy of the integrator, the second-order derivative is incorporated in the scheme. The oscillatory feature of the solution is captured by the integrators with an exponential fitting property. Three practical exponentially fitted TDRK (EFTDRK) methods are derived. To test the effectiveness of the new EFTDRK methods, the two-gene system with cross-regulation and the circadian oscillation of the period protein in Drosophila are simulated. Each EFTDRK method has the best fitting frequency which minimizes the global error. The numerical results show that the new EFTDRK methods are more accurate and more efficient than their prototype TDRK methods or RK methods of the same order and the traditional exponentially fitted RK method in the literature.

Time-Delayed Models of Gene Regulatory Networks

We discuss different mathematical models of gene regulatory networks as relevant to the onset and development of cancer. After discussion of alternative modelling approaches, we use a paradigmatic two-gene network to focus on the role played by time delays in the dynamics of gene regulatory networks. We contrast the dynamics of the reduced model arising in the limit of fast mRNA dynamics with that of the full model. The review concludes with the discussion of some open problems.

Systematic synergy modeling: understanding drug synergy from a systems biology perspective

Owing to drug synergy effects, drug combinations have become a new trend in combating complex diseases like cancer, HIV and cardiovascular diseases. However, conventional synergy quantification methods often depend on experimental dose–response data which are quite resource-demanding. In addition, these methods are unable to interpret the explicit synergy mechanism. In this review, we give representative examples of how systems biology modeling offers strategies toward better understanding of drug synergy, including the protein-protein interaction (PPI) network-based methods, pathway dynamic simulations, synergy network motif recognitions, integrative drug feature calculations, and “omic”-supported analyses. Although partially successful in drug synergy exploration and interpretation, more efforts should be put on a holistic understanding of drug-disease interactions, considering integrative pharmacology and toxicology factors. With a comprehensive and deep insight into the mechanism of drug synergy, systems biology opens a novel avenue for rational design of effective drug combinations.

Circadian systems biology: When time matters

The circadian clock is a powerful endogenous timing system, which allows organisms to fine-tune their physiology and behaviour to the geophysical time. The interplay of a distinct set of core-clock genes and proteins generates oscillations in expression of output target genes which temporally regulate numerous molecular and cellular processes. The study of the circadian timing at the organismal as well as at the cellular level outlines the field of chronobiology, which has been highly interdisciplinary ever since its origins. The development of high-throughput approaches enables the study of the clock at a systems level. In addition to experimental approaches, computational clock models exist which allow the analysis of rhythmic properties of the clock network. Such mathematical models aid mechanistic understanding and can be used to predict outcomes of distinct perturbations in clock components, thereby generating new hypotheses regarding the putative function of particular clock genes. Perturbations in the circadian timing system are linked to numerous molecular dysfunctions and may result in severe pathologies including cancer. A comprehensive knowledge regarding the mechanistic of the circadian system is crucial to develop new procedures to investigate pathologies associated with a deregulated clock.

In this manuscript we review the combination of experimental methodologies, bioinformatics and theoretical models that have been essential to explore this remarkable timing-system. Such an integrative and interdisciplinary approach may provide new strategies with regard to chronotherapeutic treatment and new insights concerning the restoration of the circadian timing in clock-associated diseases.

A Modelling Framework for Gene Regulatory Networks Including Transcription and Translation

Qualitative models of gene regulatory networks have generally considered transcription factors to regulate directly the expression of other transcription factors, without any intermediate variables. In fact, gene expression always involves transcription, which produces mRNA molecules, followed by translation, which produces protein molecules, which can then act as transcription factors for other genes (in some cases after post-transcriptional modifications). Suppressing these multiple steps implicitly assumes that the qualitative behaviour does not depend on them. Here we explore a class of expanded models that explicitly includes both transcription and translation, keeping track of both mRNA and protein concentrations. We mainly deal with regulation functions that are steep sigmoids or step functions, as is often done in protein-only models. We find that flow cannot be constrained to switching domains, though there can still be asymptotic approach to singular stationary points (fixed points in the vicinity of switching thresholds). This avoids the thorny issue of singular flow, but leads to somewhat more complicated possibilities for flow between threshold crossings. In the infinitely fast limit of either mRNA or protein rates, we find that solutions converge uniformly to solutions of the corresponding protein-only model on arbitrary finite time intervals. This leaves open the possibility that the limit system (with one type of variable infinitely fast) may have different asymptotic behaviour, and indeed, we find an example in which stability of a fixed point in the protein-only model is lost in the expanded model. Our results thus show that including mRNA as a variable may change the behaviour of solutions.

Quantitative and logic modelling of molecular and gene networks

Behaviours of complex biomolecular systems are often irreducible to the elementary properties of their individual components. Explanatory and predictive mathematical models are therefore useful for fully understanding and precisely engineering cellular functions. The development and analyses of these models require their adaptation to the problems that need to be solved and the type and amount of available genetic or molecular data. Quantitative and logic modelling are among the main methods currently used to model molecular and gene networks. Each approach comes with inherent advantages and weaknesses. Recent developments show that hybrid approaches will become essential for further progress in synthetic biology and in the development of virtual organisms.

Piecewise linear and Boolean models of chemical reaction networks

Models of biochemical networks are frequently complex and high-dimensional. Reduction methods that preserve important dynamical properties are therefore essential for their study. Interactions in biochemical networks are frequently modeled using Hill functions ([Formula: see text]). Reduced ODEs and Boolean approximations of such model networks have been studied extensively when the exponent [Formula: see text] is large. However, while the case of small constant [Formula: see text] appears in practice, it is not well understood. We provide a mathematical analysis of this limit and show that a reduction to a set of piecewise linear ODEs and Boolean networks can be mathematically justified. The piecewise linear systems have closed-form solutions that closely track those of the fully nonlinear model. The simpler, Boolean network can be used to study the qualitative behavior of the original system. We justify the reduction using geometric singular perturbation theory and compact convergence, and illustrate the results in network models of a toggle switch and an oscillator.

Quantitative dynamic modelling of the gene regulatory network controlling adipogenesis

Gene regulatory networks (GRNs) coherently coordinate the expressions of genes and control the behaviors of cellular systems. The complexity in modeling a quantitative GRN usually results from inaccurate parameter estimation, which is mostly due to small sample sizes. For better modeling of GRNs, we have designed a small-sample iterative optimization algorithm (SSIO) to quantitatively model GRNs with nonlinear regulatory relationships. The algorithm utilizes gene expression data as the primary input and it can be applied in case of small-sized samples. Using SSIO, we have quantitatively constructed the dynamic models for the GRNs controlling human and mouse adipogenesis. Compared with two other commonly-used methods, SSIO shows better performance with relatively lower residual errors, and it generates rational predictions on the adipocyte responses to external signals and steady-states. Sensitivity analysis further indicates the validity of our method. Several differences are observed between the GRNs of human and mouse adipocyte differentiations, suggesting the differences in regulatory efficiencies of the transcription factors between the two species. In addition, we use SSIO to quantitatively determine the strengths of the regulatory interactions as well as to optimize regulatory models. The results indicate that SSIO facilitates better investigation and understanding of gene regulatory processes.

Delays induce novel stochastic effects in negative feedback gene circuits

Stochastic models of reaction networks are widely used to depict gene expression dynamics. However, stochastic does not necessarily imply accurate, as subtle assumptions can yield erroneous results, masking key discrete effects. For instance, transcription and translation are not instantaneous processes-explicit delays separate their initiation from the appearance of their functional products. However, delays are often ignored in stochastic, single-gene expression models. By consequence, effects such as delay-induced stochastic oscillations at the single-cell level have remained relatively unexplored. Here, we present a systematic study of periodicity and multimodality in a simple gene circuit with negative feedback, analyzing the influence of negative feedback strength and transcriptional/translational delays on expression dynamics. We demonstrate that an oscillatory regime emerges through a Hopf bifurcation in both deterministic and stochastic frameworks. Of importance, a shift in the stochastic Hopf bifurcation evidences inaccuracies of the deterministic bifurcation analysis. Furthermore, noise fluctuations within stochastic oscillations decrease alongside increasing values of transcriptional delays and within a specific range of negative feedback strengths, whereas a strong feedback is associated with oscillations triggered by bursts. Finally, we demonstrate that explicitly accounting for delays increases the number of accessible states in the multimodal regime, and also introduces features typical of excitable systems.