Algebraic Systems Biology: A Case Study for the Wnt Pathway

Quantitative Pharmacology Methods for Bispecific T Cell Engagers

T Cell Engager (TCE)s are an exciting therapeutic modality in immuno-oncology that acts to bypass antigen presentation and forms a direct link between cancer and immune cells in the Tumor Microenvironment (TME). TCEs are efficacious only when the drug is bound to both immune and cancer cell targets. Therefore, approaches that maximize the formation of the drug-target trimer in the TME are expected to increase the drug's efficacy. In this study, we quantitatively investigate how the concentration of ternary complex and its biodistribution depend on both the targets' specific properties and the design characteristics of the TCE, and specifically on the binding kinetics of the drug to its targets. A simplified mathematical model of drug-target interactions is considered here, with insights from the "three-body" problem applied to the model. Parameter identifiability analysis performed on the model demonstrates that steady state data, which is often available at the early pre-clinical stages, is sufficient to estimate the binding affinity of the TCE molecule to both targets. We used the model to analyze several existing antibodies, both clinically approved and under development, to explore their common kinetic features. The manuscript concludes with an assessment of a full quantitative pharmacology model that accounts for drug disposition into the peripheral compartment.

Algebra, Geometry and Topology of ERK Kinetics

The MEK/ERK signalling pathway is involved in cell division, cell specialisation, survival and cell death (Shaul and Seger in Biochim Biophys Acta (BBA)-Mol Cell Res 1773(8):1213-1226, 2007). Here we study a polynomial dynamical system describing the dynamics of MEK/ERK proposed by Yeung et al. (Curr Biol 2019, https://doi.org/10.1016/j.cub.2019.12.052 ) with their experimental setup, data and known biological information. The experimental dataset is a time-course of ERK measurements in different phosphorylation states following activation of either wild-type MEK or MEK mutations associated with cancer or developmental defects. We demonstrate how methods from computational algebraic geometry, differential algebra, Bayesian statistics and computational algebraic topology can inform the model reduction, identification and parameter inference of MEK variants, respectively. Throughout, we show how this algebraic viewpoint offers a rigorous and systematic analysis of such models.

The linear framework: using graph theory to reveal the algebra and thermodynamics of biomolecular systems

The linear framework uses finite, directed graphs with labelled edges to model biomolecular systems. Graph vertices represent biochemical species or molecular states, edges represent reactions or transitions and labels represent rates. The graph yields a linear dynamics for molecular concentrations or state probabilities, with the graph Laplacian as the operator, and the labels encode the nonlinear interactions between system and environment. The labels can be specified by vertices of other graphs or by conservation laws or, when the environment consists of thermodynamic reservoirs, they may be constants. In the latter case, the graphs correspond to infinitesimal generators of Markov processes. The key advantage of the framework has been that steady states are determined as rational algebraic functions of the labels by the Matrix-Tree theorems of graph theory. When the system is at thermodynamic equilibrium, this prescription recovers equilibrium statistical mechanics but it continues to hold for non-equilibrium steady states. The framework goes beyond other graph-based approaches in treating the graph as a mathematical object, for which general theorems can be formulated that accommodate biomolecular complexity. It has been particularly effective at analysing enzyme-catalysed modification systems and input-output responses.

A blood atlas of COVID-19 defines hallmarks of disease severity and specificity

Treatment of severe COVID-19 is currently limited by clinical heterogeneity and incomplete description of specific immune biomarkers. We present here a comprehensive multi-omic blood atlas for patients with varying COVID-19 severity in an integrated comparison with influenza and sepsis patients versus healthy volunteers. We identify immune signatures and correlates of host response. Hallmarks of disease severity involved cells, their inflammatory mediators and networks, including progenitor cells and specific myeloid and lymphocyte subsets, features of the immune repertoire, acute phase response, metabolism, and coagulation. Persisting immune activation involving AP-1/p38MAPK was a specific feature of COVID-19. The plasma proteome enabled sub-phenotyping into patient clusters, predictive of severity and outcome. Systems-based integrative analyses including tensor and matrix decomposition of all modalities revealed feature groupings linked with severity and specificity compared to influenza and sepsis. Our approach and blood atlas will support future drug development, clinical trial design, and personalized medicine approaches for COVID-19.

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.

The Case for Algebraic Biology: from Research to Education

Though it goes without saying that linear algebra is fundamental to mathematical biology, polynomial algebra is less visible. In this article, we will give a brief tour of four diverse biological problems where multivariate polynomials play a central role-a subfield that is sometimes called algebraic biology. Namely, these topics include biochemical reaction networks, Boolean models of gene regulatory networks, algebraic statistics and genomics, and place fields in neuroscience. After that, we will summarize the history of discrete and algebraic structures in mathematical biology, from their early appearances in the late 1960s to the current day. Finally, we will discuss the role of algebraic biology in the modern classroom and curriculum, including resources in the literature and relevant software. Our goal is to make this article widely accessible, reaching the mathematical biologist who knows no algebra, the algebraist who knows no biology, and especially the interested student who is curious about the synergy between these two seemingly unrelated fields.

Robustness and parameter geography in post-translational modification systems

Biological systems are acknowledged to be robust to perturbations but a rigorous understanding of this has been elusive. In a mathematical model, perturbations often exert their effect through parameters, so sizes and shapes of parametric regions offer an integrated global estimate of robustness. Here, we explore this “parameter geography” for bistability in post-translational modification (PTM) systems. We use the previously developed “linear framework” for timescale separation to describe the steady-states of a two-site PTM system as the solutions of two polynomial equations in two variables, with eight non-dimensional parameters. Importantly, this approach allows us to accommodate enzyme mechanisms of arbitrary complexity beyond the conventional Michaelis-Menten scheme, which unrealistically forbids product rebinding. We further use the numerical algebraic geometry tools Bertini, Paramotopy, and alphaCertified to statistically assess the solutions to these equations at ∼10⁹ parameter points in total. Subject to sampling limitations, we find no bistability when substrate amount is below a threshold relative to enzyme amounts. As substrate increases, the bistable region acquires 8-dimensional volume which increases in an apparently monotonic and sigmoidal manner towards saturation. The region remains connected but not convex, albeit with a high visibility ratio. Surprisingly, the saturating bistable region occupies a much smaller proportion of the sampling domain under mechanistic assumptions more realistic than the Michaelis-Menten scheme. We find that bistability is compromised by product rebinding and that unrealistic assumptions on enzyme mechanisms have obscured its parametric rarity. The apparent monotonic increase in volume of the bistable region remains perplexing because the region itself does not grow monotonically: parameter points can move back and forth between monostability and bistability. We suggest mathematical conjectures and questions arising from these findings. Advances in theory and software now permit insights into parameter geography to be uncovered by high-dimensional, data-centric analysis.

Biological organisms are often said to have robust properties but it is difficult to understand how such robustness arises from molecular interactions. Here, we use a mathematical model to study how the molecular mechanism of protein modification exhibits the property of multiple internal states, which has been suggested to underlie memory and decision making. The robustness of this property is revealed by the size and shape, or “geography,” of the parametric region in which the property holds. We use advances in reducing model complexity and in rapidly solving the underlying equations, to extensively sample parameter points in an 8-dimensional space. We find that under realistic molecular assumptions the size of the region is surprisingly small, suggesting that generating multiple internal states with such a mechanism is much harder than expected. While the shape of the region appears straightforward, we find surprising complexity in how the region grows with increasing amounts of the modified substrate. Our approach uses statistical analysis of data generated from a model, rather than from experiments, but leads to precise mathematical conjectures about parameter geography and biological robustness.

Joining and decomposing reaction networks

In systems and synthetic biology, much research has focused on the behavior and design of single pathways, while, more recently, experimental efforts have focused on how cross-talk (coupling two or more pathways) or inhibiting molecular function (isolating one part of the pathway) affects systems-level behavior. However, the theory for tackling these larger systems in general has lagged behind. Here, we analyze how joining networks (e.g., cross-talk) or decomposing networks (e.g., inhibition or knock-outs) affects three properties that reaction networks may possess-identifiability (recoverability of parameter values from data), steady-state invariants (relationships among species concentrations at steady state, used in model selection), and multistationarity (capacity for multiple steady states, which correspond to multiple cell decisions). Specifically, we prove results that clarify, for a network obtained by joining two smaller networks, how properties of the smaller networks can be inferred from or can imply similar properties of the original network. Our proofs use techniques from computational algebraic geometry, including elimination theory and differential algebra.

Inferring reaction network structure from single-cell, multiplex data, using toric systems theory

The goal of many single-cell studies on eukaryotic cells is to gain insight into the biochemical reactions that control cell fate and state. In this paper we introduce the concept of Effective Stoichiometric Spaces (ESS) to guide the reconstruction of biochemical networks from multiplexed, fixed time-point, single-cell data. In contrast to methods based solely on statistical models of data, the ESS method leverages the power of the geometric theory of toric varieties to begin unraveling the structure of chemical reaction networks (CRN). This application of toric theory enables a data-driven mapping of covariance relationships in single-cell measurements into stoichiometric information, one in which each cell subpopulation has its associated ESS interpreted in terms of CRN theory. In the development of ESS we reframe certain aspects of the theory of CRN to better match data analysis. As an application of our approach we process cytomery- and image-based single-cell datasets and identify differences in cells treated with kinase inhibitors. Our approach is directly applicable to data acquired using readily accessible experimental methods such as Fluorescence Activated Cell Sorting (FACS) and multiplex immunofluorescence.

We introduce a new notion, which we call the effective stoichiometric space (ESS), that elucidates network structure from the covariances of single-cell multiplex data. The ESS approach differs from methods that are based on purely statistical models of data: it allows a completely new and data-driven translation of the theory of toric varieties in geometry and specifically their role in chemical reaction networks (CRN). In the process, we reframe certain aspects of the theory of CRN. As illustrations of our approach, we find stoichiometry in different single-cell datasets, and pinpoint dose-dependence of network perturbations in drug-treated cells.

Multistationarity in the Space of Total Concentrations for Systems that Admit a Monomial Parametrization

We apply tools from real algebraic geometry to the problem of multistationarity of chemical reaction networks. A particular focus is on the case of reaction networks whose steady states admit a monomial parametrization. For such systems, we show that in the space of total concentrations multistationarity is scale invariant: If there is multistationarity for some value of the total concentrations, then there is multistationarity on the entire ray containing this value (possibly for different rate constants)-and vice versa. Moreover, for these networks it is possible to decide about multistationarity independent of the rate constants by formulating semi-algebraic conditions that involve only concentration variables. These conditions can easily be extended to include total concentrations. Hence, quantifier elimination may give new insights into multistationarity regions in the space of total concentrations. To demonstrate this, we show that for the distributive phosphorylation of a protein at two binding sites multistationarity is only possible if the total concentration of the substrate is larger than either the total concentration of the kinase or the total concentration of the phosphatase. This result is enabled by the chamber decomposition of the space of total concentrations from polyhedral geometry. Together with the corresponding sufficiency result of Bihan et al., this yields a characterization of multistationarity up to lower-dimensional regions.

Dynamics of Posttranslational Modification Systems: Recent Progress and Future Directions

Posttranslational modification of proteins is important for signal transduction, and hence significant effort has gone toward understanding how posttranslational modification networks process information. This involves, on the theory side, analyzing the dynamical systems arising from such networks. Which networks are, for instance, bistable? Which networks admit sustained oscillations? Which parameter values enable such behaviors? In this Biophysical Perspective, we highlight recent progress in this area and point out some important future directions. Along the way, we summarize several techniques for analyzing general networks, such as eliminating variables to obtain steady-state parameterizations, and harnessing results on how incorporating intermediates affects dynamics.

A Deficiency-Based Approach to Parametrizing Positive Equilibria of Biochemical Reaction Systems

We present conditions which guarantee a parametrization of the set of positive equilibria of a generalized mass-action system. Our main results state that (1) if the underlying generalized chemical reaction network has an effective deficiency of zero, then the set of positive equilibria coincides with the parametrized set of complex-balanced equilibria and (2) if the network is weakly reversible and has a kinetic deficiency of zero, then the equilibrium set is nonempty and has a positive, typically rational, parametrization. Via the method of network translation, we apply our results to classical mass-action systems studied in the biochemical literature, including the EnvZ–OmpR and shuttled WNT signaling pathways. A parametrization of the set of positive equilibria of a (generalized) mass-action system is often a prerequisite for the study of multistationarity and allows an easy check for the occurrence of absolute concentration robustness, as we demonstrate for the EnvZ–OmpR pathway.

Bond graph modelling of the cardiac action potential: implications for drift and non-unique steady states

Mathematical models of cardiac action potentials have become increasingly important in the study of heart disease and pharmacology, but concerns linger over their robustness during long periods of simulation, in particular due to issues such as model drift and non-unique steady states. Previous studies have linked these to violation of conservation laws, but only explored those issues with respect to charge conservation in specific models. Here, we propose a general and systematic method of identifying conservation laws hidden in models of cardiac electrophysiology by using bond graphs, and develop a bond graph model of the cardiac action potential to study long-term behaviour. Bond graphs provide an explicit energy-based framework for modelling physical systems, which makes them well suited for examining conservation within electrophysiological models. We find that the charge conservation laws derived in previous studies are examples of the more general concept of a 'conserved moiety'. Conserved moieties explain model drift and non-unique steady states, generalizing the results from previous studies. The bond graph approach provides a rigorous method to check for drift and non-unique steady states in a wide range of cardiac action potential models, and can be extended to examine behaviours of other excitable systems.

Evidence Synthesis for Stochastic Epidemic Models

In recent years, the role of epidemic models in informing public health policies has progressively grown. Models have become increasingly realistic and more complex, requiring the use of multiple data sources to estimate all quantities of interest. This review summarises the different types of stochastic epidemic models that use evidence synthesis and highlights current challenges.

Numerical algebraic geometry for model selection and its application to the life sciences

Researchers working with mathematical models are often confronted by the related problems of parameter estimation, model validation and model selection. These are all optimization problems, well known to be challenging due to nonlinearity, non-convexity and multiple local optima. Furthermore, the challenges are compounded when only partial data are available. Here, we consider polynomial models (e.g. mass-action chemical reaction networks at steady state) and describe a framework for their analysis based on optimization using numerical algebraic geometry. Specifically, we use probability-one polynomial homotopy continuation methods to compute all critical points of the objective function, then filter to recover the global optima. Our approach exploits the geometrical structures relating models and data, and we demonstrate its utility on examples from cell signalling, synthetic biology and epidemiology.

Identifying parameter regions for multistationarity

Mathematical modelling has become an established tool for studying the dynamics of biological systems. Current applications range from building models that reproduce quantitative data to identifying systems with predefined qualitative features, such as switching behaviour, bistability or oscillations. Mathematically, the latter question amounts to identifying parameter values associated with a given qualitative feature. We introduce a procedure to partition the parameter space of a parameterized system of ordinary differential equations into regions for which the system has a unique or multiple equilibria. The procedure is based on the computation of the Brouwer degree, and it creates a multivariate polynomial with parameter depending coefficients. The signs of the coefficients determine parameter regions with and without multistationarity. A particular strength of the procedure is the avoidance of numerical analysis and parameter sampling. The procedure consists of a number of steps. Each of these steps might be addressed algorithmically using various computer programs and available software, or manually. We demonstrate our procedure on several models of gene transcription and cell signalling, and show that in many cases we obtain a complete partitioning of the parameter space with respect to multistationarity.

The role of the Hes1 crosstalk hub in Notch-Wnt interactions of the intestinal crypt

The Notch pathway plays a vital role in determining whether cells in the intestinal epithelium adopt a secretory or an absorptive phenotype. Cell fate specification is coordinated via Notch's interaction with the canonical Wnt pathway. Here, we propose a new mathematical model of the Notch and Wnt pathways, in which the Hes1 promoter acts as a hub for pathway crosstalk. Computational simulations of the model can assist in understanding how healthy intestinal tissue is maintained, and predict the likely consequences of biochemical knockouts upon cell fate selection processes. Chemical reaction network theory (CRNT) is a powerful, generalised framework which assesses the capacity of our model for monostability or multistability, by analysing properties of the underlying network structure without recourse to specific parameter values or functional forms for reaction rates. CRNT highlights the role of β-catenin in stabilising the Notch pathway and damping oscillations, demonstrating that Wnt-mediated actions on the Hes1 promoter can induce dynamic transitions in the Notch system, from multistability to monostability. Time-dependent model simulations of cell pairs reveal the stabilising influence of Wnt upon the Notch pathway, in which β-catenin- and Dsh-mediated action on the Hes1 promoter are key in shaping the subcellular dynamics. Where Notch-mediated transcription of Hes1 dominates, there is Notch oscillation and maintenance of fate flexibility; Wnt-mediated transcription of Hes1 favours bistability akin to cell fate selection. Cells could therefore regulate the proportion of Wnt- and Notch-mediated control of the Hes1 promoter to coordinate the timing of cell fate selection as they migrate through the intestinal epithelium and are subject to reduced Wnt stimuli. Furthermore, mutant cells characterised by hyperstimulation of the Wnt pathway may, through coupling with Notch, invert cell fate in neighbouring healthy cells, enabling an aberrant cell to maintain its neighbours in mitotically active states.