Mathematical Analysis of Chemical Reaction Systems

PCA-based synthetic sensitivity coefficients for chemical reaction network in cancer

Chemical reaction networks are powerful tools for modeling cell signaling and its disruptions in diseases like cancer. Realistic chemical reaction networks involve hundreds of proteins and reactions, resulting in a model depending on a consistently large number of kinetic parameters. Since finely calibrating all the parameters would require an unrealistic amount of data, proper sensitivity analysis is required to identify a subset of parameters for which fine tuning is needed and thus provide a fundamental tool for the qualitative analysis of the network. We present a multidisciplinary approach for computing a set of synthetic sensitivity indices. These indices rank the kinetic parameters, based on the impact that errors in their values would have on the protein concentration profile at equilibrium. Our tests on a chemical reaction network devised for colorectal cells demonstrate the effectiveness of the considered sensitivity indices in different scenarios including in-silico drug dosage and novel therapeutic target discovery. The Matlab code for computing the synthetic sensitivity indices and the data concerning the network for colorectal cells are available at https://github.com/theMIDAgroup/CRN_sensitivity.

A mathematically rigorous algorithm to define, compute and assess relevance of the probable dissociation constants in characterizing a biochemical network

Metabolism results from enzymatic- and non-enzymatic interactions of several molecules, is easily parameterized with the dissociation constant and occurs via biochemical networks. The dissociation constant is an empirically determined parameter and cannot be used directly to investigate in silico models of biochemical networks. Here, we develop and present an algorithm to define, compute and assess the relevance of the probable dissociation constant for every reaction of a biochemical network. The reactants and reactions of this network are modelled by a stoichiometry number matrix. The algorithm computes the null space and then serially generates subspaces by combinatorially summing the spanning vectors that are non-trivial and unique. This is done until the terms of each row either monotonically diverge or form an alternating sequence whose terms can be partitioned into subsets with almost the same number of oppositely signed terms. For a selected null space-generated subspace the algorithm utilizes several statistical and mathematical descriptors to select and bin terms from each row into distinct outcome-specific subsets. The terms of each subset are summed, mapped to the real-valued open interval [Formula: see text] and used to populate a reaction-specific outcome vector. The p1-norm for this vector is then the probable dissociation constant for this reaction. These steps are continued until every reaction of a modelled network is unambiguously annotated. The assertions presented are complemented by computational studies of a biochemical network for aerobic glycolysis. The fundamental premise of this work is that every row of a null space-generated subspace is a valid reaction and can therefore, be modelled as a reaction-specific sequence vector with a dimension that corresponds to the cardinality of the subspace after excluding all trivial- and redundant-vectors. A major finding of this study is that the row-wise sum or the sum of the terms contained in each reaction-specific sequence vector is mapped unambiguously to a positive real number. This means that the probable dissociation constants, for all reactions, can be directly computed from the stoichiometry number matrix and are suitable indicators of outcome for every reaction of the modelled biochemical network. Additionally, we find that the unambiguous annotation for a biochemical network will require a minimum number of iterations and will determine computational complexity.

Oscillations in three-reaction quadratic mass-action systems

It is known that rank-two bimolecular mass-action systems do not admit limit cycles. With a view to understanding which small mass-action systems admit oscillation, in this paper we study rank-two networks with bimolecular source complexes but allow target complexes with higher molecularities. As our goal is to find oscillatory networks of minimal size, we focus on networks with three reactions, the minimum number that is required for oscillation. However, some of our intermediate results are valid in greater generality. One key finding is that an isolated periodic orbit cannot occur in a three-reaction, trimolecular, mass-action system with bimolecular sources. In fact, we characterize all networks in this class that admit a periodic orbit; in every case, all nearby orbits are periodic too. Apart from the well-known Lotka and Ivanova reactions, we identify another network in this class that admits a center. This new network exhibits a vertical Andronov-Hopf bifurcation. Furthermore, we characterize all two-species, three-reaction, bimolecular-sourced networks that admit an Andronov-Hopf bifurcation with mass-action kinetics. These include two families of networks that admit a supercritical Andronov-Hopf bifurcation and hence a stable limit cycle. These networks necessarily have a target complex with a molecularity of at least four, and it turns out that there are exactly four such networks that are tetramolecular.

ReDirection: an R-package to compute the probable dissociation constant for every reaction of a user-defined biochemical network

Biochemical networks integrate enzyme-mediated substrate conversions with non-enzymatic complex formation and disassembly to accomplish complex biochemical and physiological functions. The choice of parameters and constraints used in most of these studies is numerically motivated and network-specific. Although sound in theory, the outcomes that result depart significantly from the intracellular milieu and are less likely to retain relevance in a clinical setting. There is a need for a computational tool which is biochemically relevant, mathematically rigorous, and unbiased, and can ascribe functionality to and generate potentially testable hypotheses for a user-defined biochemical network. Here, we present "ReDirection," an R-package which computes the probable dissociation constant for every reaction of a biochemical network directly from a null space-generated subspace of the stoichiometry number matrix of the modeled network. "ReDirection" delineates this subspace by excluding all trivial and redundant or duplicate occurrences of non-trivial vectors, combinatorially summing the vectors that remain and verifying that the upper or lower bounds of the sequence of terms formed by each row of this subspace belong to the open real-valued intervals - ∞ , - 1 or 1 , ∞ or whether the number of terms that are differently signed are almost equal. "ReDirection" iterates these steps until these bounds are consistent and unambiguous for all reactions of the modeled biochemical network. Thereafter, "ReDirection" filters the terms from each row of this subspace, bins them to outcome-specific subsets, sums and maps this to an outcome-specific reaction vector, and computes the p1-norm, which is the probable dissociation constant for a reaction. "ReDirection" works on first principles, does not discriminate between enzymatic and non-enzymatic reactions, offers a biochemically relevant and mathematically rigorous environment to explore user-defined biochemical networks under baseline and perturbed conditions, and can be used to address empirically intractable biochemical problems. The utility and relevance of "ReDirection" are highlighted by numerical studies on stoichiometric number models of biochemical networks of galactose metabolism and heme and cholesterol biosynthesis. "ReDirection" is freely available and accessible from the comprehensive R archive network (CRAN) with the URL (https://cran.r-project.org/package=ReDirection).

On The Existence of Nonunique Equilibrium States

Chemists "know" that the equilibrium position in a homogeneous solution is unique and that the one equilibrium is reached, irrespective of how the solution was prepared. This concept has been verified based on the law of mass action applied to concentrations rather than activities. In this contribution, we explore the situation for a very basic equilibrium system, using established approximations for the activity coefficients. We will demonstrate that under the correct application of the law of mass action, bistable equilibria positions are possible, concomitant, of course, with bifurcation in the kinetics of reaching either of the stable equilibria.

Chemical Systems with Limit Cycles

The dynamics of a chemical reaction network (CRN) is often modeled under the assumption of mass action kinetics by a system of ordinary differential equations (ODEs) with polynomial right-hand sides that describe the time evolution of concentrations of chemical species involved. Given an arbitrarily large integer [Formula: see text], we show that there exists a CRN such that its ODE model has at least K stable limit cycles. Such a CRN can be constructed with reactions of at most second-order provided that the number of chemical species grows linearly with K. Bounds on the minimal number of chemical species and the minimal number of chemical reactions are presented for CRNs with K stable limit cycles and at most second order or seventh-order kinetics. We also show that CRNs with only two chemical species can have K stable limit cycles, when the order of chemical reactions grows linearly with K.

Intrinsic catalytic properties of histone H3 lysine-9 methyltransferases preserve monomethylation levels under low S-adenosylmethionine

S-adenosylmethionine (SAM) is the methyl donor for site-specific methylation reactions on histone proteins, imparting key epigenetic information. During SAM-depleted conditions that can arise from dietary methionine restriction, lysine di- and tri-methylation are reduced while sites such as Histone-3 lysine-9 (H3K9) are actively maintained, allowing cells to restore higher-state methylation upon metabolic recovery. Here, we investigated if the intrinsic catalytic properties of H3K9 histone methyltransferases (HMTs) contribute to this epigenetic persistence. We employed systematic kinetic analyses and substrate binding assays using four recombinant H3K9 HMTs (i.e., EHMT1, EHMT2, SUV39H1, and SUV39H2). At both high and low (i.e., sub-saturating) SAM, all HMTs displayed the highest catalytic efficiency (kcat/KM) for monomethylation compared to di- and trimethylation on H3 peptide substrates. The favored monomethylation reaction was also reflected in kcat values, apart from SUV39H2 which displayed a similar kcat regardless of substrate methylation state. Using differentially methylated nucleosomes as substrates, kinetic analyses of EHMT1 and EHMT2 revealed similar catalytic preferences. Orthogonal binding assays revealed only small differences in substrate affinity across methylation states, suggesting that catalytic steps dictate the monomethylation preferences of EHMT1, EHMT2, and SUV39H1. To link in vitro catalytic rates with nuclear methylation dynamics, we built a mathematical model incorporating measured kinetic parameters and a time course of mass spectrometry-based H3K9 methylation measurements following cellular SAM depletion. The model revealed that the intrinsic kinetic constants of the catalytic domains could recapitulate in vivo observations. Together, these results suggest catalytic discrimination by H3K9 HMTs maintains nuclear H3K9me1, ensuring epigenetic persistence after metabolic stress.

Action functional gradient descent algorithm for estimating escape paths in stochastic chemical reaction networks

We first derive the Hamilton-Jacobi theory underlying continuous-time Markov processes, and then we use the construction to develop a variational algorithm for estimating escape (least improbable or first passage) paths for a generic stochastic chemical reaction network that exhibits multiple fixed points. The design of our algorithm is such that it is independent of the underlying dimensionality of the system, the discretization control parameters are updated toward the continuum limit, and there is an easy-to-calculate measure for the correctness of its solution. We consider several applications of the algorithm and verify them against computationally expensive means such as the shooting method and stochastic simulation. While we employ theoretical techniques from mathematical physics, numerical optimization and chemical reaction network theory, we hope that our work finds practical applications with an inter-disciplinary audience including chemists, biologists, optimal control theorists and game theorists.

Homeostasis and injectivity: a reaction network perspective

Homeostasis represents the idea that a feature may remain invariant despite changes in some external parameters. We establish a connection between homeostasis and injectivity for reaction network models. In particular, we show that a reaction network cannot exhibit homeostasis if a modified version of the network (which we call homeostasis-associated network) is injective. We provide examples of reaction networks which can or cannot exhibit homeostasis by analyzing the injectivity of their homeostasis-associated networks.

Foundations of static and dynamic absolute concentration robustness

Absolute Concentration Robustness (ACR) was introduced by Shinar and Feinberg (Science 327:1389-1391, 2010) as robustness of equilibrium species concentration in a mass action dynamical system. Their aim was to devise a mathematical condition that will ensure robustness in the function of the biological system being modeled. The robustness of function rests on what we refer to as empirical robustness-the concentration of a species remains unvarying, when measured in the long run, across arbitrary initial conditions. Even simple examples show that the ACR notion introduced in Shinar and Feinberg (Science 327:1389-1391, 2010) (here referred to as static ACR) is neither necessary nor sufficient for empirical robustness. To make a stronger connection with empirical robustness, we define dynamic ACR, a property related to long-term, global dynamics, rather than only to equilibrium behavior. We discuss general dynamical systems with dynamic ACR properties as well as parametrized families of dynamical systems related to reaction networks. We find necessary and sufficient conditions for dynamic ACR in complex balanced reaction networks, a class of networks that is central to the theory of reaction networks.

Multistationarity in Cyclic Sequestration-Transmutation Networks

We consider a natural class of reaction networks which consist of reactions where either two species can inactivate each other (i.e., sequestration), or some species can be transformed into another (i.e., transmutation), in a way that gives rise to a feedback cycle. We completely characterize the capacity of multistationarity of these networks. This is especially interesting because such networks provide simple examples of "atoms of multistationarity", i.e., minimal networks that can give rise to multiple positive steady states.

Uniqueness of weakly reversible and deficiency zero realizations of dynamical systems

A reaction network together with a choice of rate constants uniquely gives rise to a system of differential equations, according to the law of mass-action kinetics. On the other hand, different networks can generate the same dynamical system under mass-action kinetics. Therefore, the problem of identifying "the" underlying network of a dynamical system is not well-posed, in general. Here we show that the problem of identifying an underlying weakly reversible deficiency zero network is well-posed, in the sense that the solution is unique whenever it exists. This can be very useful in applications because from the perspective of both dynamics and network structure, a weakly reversible deficiency zero (WR₀) realization is the simplest possible one. Moreover, while mass-action systems can exhibit practically any dynamical behavior, including multistability, oscillations, and chaos, WR₀ systems are remarkably stable for any choice of rate constants: they have a unique positive steady state within each invariant polyhedron, and cannot give rise to oscillations or chaotic dynamics. We also prove that both of our hypotheses (i.e., weak reversibility and deficiency zero) are necessary for uniqueness.

Computational quantification of global effects induced by mutations and drugs in signaling networks of colorectal cancer cells

Colorectal cancer (CRC) is one of the most deadly and commonly diagnosed tumors worldwide. Several genes are involved in its development and progression. The most frequent mutations concern APC, KRAS, SMAD4, and TP53 genes, suggesting that CRC relies on the concomitant alteration of the related pathways. However, with classic molecular approaches, it is not easy to simultaneously analyze the interconnections between these pathways. To overcome this limitation, recently these pathways have been included in a huge chemical reaction network (CRN) describing how information sensed from the environment by growth factors is processed by healthy colorectal cells. Starting from this CRN, we propose a computational model which simulates the effects induced by single or multiple concurrent mutations on the global signaling network. The model has been tested in three scenarios. First, we have quantified the changes induced on the concentration of the proteins of the network by a mutation in APC, KRAS, SMAD4, or TP53. Second, we have computed the changes in the concentration of p53 induced by up to two concurrent mutations affecting proteins upstreams in the network. Third, we have considered a mutated cell affected by a gain of function of KRAS, and we have simulated the action of Dabrafenib, showing that the proposed model can be used to determine the most effective amount of drug to be delivered to the cell. In general, the proposed approach displays several advantages, in that it allows to quantify the alteration in the concentration of the proteins resulting from a single or multiple given mutations. Moreover, simulations of the global signaling network of CRC may be used to identify new therapeutic targets, or to disclose unexpected interactions between the involved pathways.

Delay stability of reaction systems

Delay differential equations are used as a model when the effect of past states has to be taken into account. In this work we consider delay models of chemical reaction networks with mass action kinetics. We obtain a sufficient condition for absolute delay stability of equilibrium concentrations, i.e., local asymptotic stability independent of the delay parameters. Several interesting examples on sequestration networks with delays are presented.

Metabolite-mediated modelling of microbial community dynamics captures emergent behaviour more effectively than species-species modelling

Personalized models of the gut microbiome are valuable for disease prevention and treatment. For this, one requires a mathematical model that predicts microbial community composition and the emergent behaviour of microbial communities. We seek a modelling strategy that can capture emergent behaviour when built from sets of universal individual interactions. Our investigation reveals that species-metabolite interaction (SMI) modelling is better able to capture emergent behaviour in community composition dynamics than direct species-species modelling. Using publicly available data, we examine the ability of species-species models and species-metabolite models to predict trio growth experiments from the outcomes of pair growth experiments. We compare quadratic species-species interaction models and quadratic SMI models and conclude that only species-metabolite models have the necessary complexity to explain a wide variety of interdependent growth outcomes. We also show that general species-species interaction models cannot match the patterns observed in community growth dynamics, whereas species-metabolite models can. We conclude that species-metabolite modelling will be important in the development of accurate, clinically useful models of microbial communities.

Strategies for Coexistence in Molecular Communication

Some of the most ambitious applications of molecular communications are expected to lie in nanomedicine and advanced manufacturing. In these domains, the molecular communication system is surrounded by a range of biochemical processes, some of which may be sensitive to chemical species used for communication. Under these conditions, the biological system and the molecular communication system impact each other. As such, the problem of coexistence arises, where both the reliability of the molecular communication system and the function of the biological system must be ensured. In this paper, we study this problem with a focus on interactions with biological systems equipped with chemosensing mechanisms, which arises in a large class of biological systems. We motivate the problem by considering chemosensing mechanisms arising in bacteria chemotaxis, a ubiquitous and well-understood class of biological systems. We, then, propose strategies for a molecular communication system to minimize the disruption of biological system equipped with a chemosensing mechanism. This is achieved by exploiting tools from the theory of chemical reaction networks. To investigate the capabilities of our strategies, we obtain fundamental information theoretic limits by establishing a new connection with the problem of covert communications.