Simplification of Mathematical Models of Chemical Reaction Systems

Universal slow dynamics of chemical reaction networks

Understanding the emergent behavior of chemical reaction networks (CRNs) is a fundamental aspect of biology and its origin from inanimate matter. A closed CRN monotonically tends to thermal equilibrium, but when it is opened to external reservoirs, a range of behaviors is possible, including transition to a new equilibrium state, a nonequilibrium state, or indefinite growth. This study shows that slowly driven CRNs are governed by the conserved quantities of the closed system, which are generally far fewer in number than the species. Considering both deterministic and stochastic dynamics, a universal slow-dynamics equation is derived with singular perturbation methods and is shown to be thermodynamically consistent. The slow dynamics is highly robust against microscopic details of the network, which may be unknown in practical situations. In particular, nonequilibrium states of realistic large CRNs can be sought without knowledge of bulk reaction rates. The framework is successfully tested against a suite of networks of increasing complexity and argued to be relevant in the treatment of open CRNs as chemical machines.

Simultaneously Predicting the Pharmacokinetics of CES1-Metabolized Drugs and Their Metabolites Using Physiologically Based Pharmacokinetic Model in Cirrhosis Subjects

Hepatic carboxylesterase 1 (CES1) metabolizes numerous prodrugs into active ingredients or direct-acting drugs into inactive metabolites. We aimed to develop a semi-physiologically based pharmacokinetic (semi-PBPK) model to simultaneously predict the pharmacokinetics of CES1 substrates and their active metabolites in liver cirrhosis (LC) patients. Six prodrugs (enalapril, benazepril, cilazapril, temocapril, perindopril and oseltamivir) and three direct-acting drugs (flumazenil, pethidine and remimazolam) were selected. Parameters such as organ blood flows, plasma-binding protein concentrations, functional liver volume, hepatic enzymatic activity, glomerular filtration rate (GFR) and gastrointestinal transit rate were integrated into the simulation. The pharmacokinetic profiles of these drugs and their active metabolites were simulated for 1000 virtual individuals. The developed semi-PBPK model, after validation in healthy individuals, was extrapolated to LC patients. Most of the observations fell within the 5th and 95th percentiles of simulations from 1000 virtual patients. The estimated AUC and Cmax were within 0.5-2-fold of the observed values. The sensitivity analysis showed that the decreased plasma exposure of active metabolites due to the decreased CES1 was partly attenuated by the decreased GFR. Conclusion: The developed PBPK model successfully predicted the pharmacokinetics of CES1 substrates and their metabolites in healthy individuals and LC patients, facilitating tailored dosing of CES1 substrates in LC patients.

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

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

Embracing fine-root system complexity in terrestrial ecosystem modeling

Projecting the dynamics and functioning of the biosphere requires a holistic consideration of whole-ecosystem processes. However, biases toward leaf, canopy, and soil modeling since the 1970s have constantly left fine-root systems being rudimentarily treated. As accelerated empirical advances in the last two decades establish clearly functional differentiation conferred by the hierarchical structure of fine-root orders and associations with mycorrhizal fungi, a need emerges to embrace this complexity to bridge the data-model gap in still extremely uncertain models. Here, we propose a three-pool structure comprising transport and absorptive fine roots with mycorrhizal fungi (TAM) to model vertically resolved fine-root systems across organizational and spatial-temporal scales. Emerging from a conceptual shift away from arbitrary homogenization, TAM builds upon theoretical and empirical foundations as an effective and efficient approximation that balances realism and simplicity. A proof-of-concept demonstration of TAM in a big-leaf model both conservatively and radically shows robust impacts of differentiation within fine-root systems on simulating carbon cycling in temperate forests. Theoretical and quantitative support warrants exploiting its rich potentials across ecosystems and models to confront uncertainties and challenges for a predictive understanding of the biosphere. Echoing a broad trend of embracing ecological complexity in integrative ecosystem modeling, TAM may offer a consistent framework where modelers and empiricists can work together toward this grand goal.

A compatibility evaluation between the physiologically based pharmacokinetic (PBPK) model and the compartmental PK model using the lumping method with real cases

Pharmacokinetic (PK) modeling is a useful method for investigating drug absorption, distribution, metabolism, and excretion. The most commonly used mathematical models in PK modeling are the compartment model and physiologically based pharmacokinetic (PBPK) model. Although the theoretical characteristics of each model are well known, there have been few comparative studies of the compatibility of the models. Therefore, we evaluated the compatibility of PBPK and compartment models using the lumping method with 20 model compounds. The PBPK model was theoretically reduced to the lumped model using the principle of grouping tissues and organs that show similar kinetic behaviors. The area under the concentration-time curve (AUC) based on the simulated concentration and PK parameters (drug clearance [CL], central volume of distribution [Vc], peripheral volume of distribution [Vp]) in each model were compared, assuming administration to humans. The AUC and PK parameters in the PBPK model were similar to those in the lumped model within the 2-fold range for 17 of 20 model compounds (85%). In addition, the relationship of the calculated Vd/fu (volume of distribution [Vd], drug-unbound fraction [fu]) and the accuracy of AUC between the lumped model and compartment model confirmed their compatibility. Accordingly, the compatibility between PBPK and compartment models was confirmed by the lumping method. This method can be applied depending on the requirement of compatibility between the two models.

The Development of a PBPK Model for Atomoxetine Using Levels in Plasma, Saliva and Brain Extracellular Fluid in Patients with Normal and Deteriorated Kidney Function

BACKGROUND

Atomoxetine is a treatment for attention-deficit hyperactivity disorder. It inhibits Norepinephrine Transporters (NET) in the brain. Renal impairment can reduce hepatic CYP2D6 activity and atomoxetine elimination which may increase its body exposure. Atomoxetine can be secreted in saliva.

OBJECTIVE

The objective of this work was to test the hypothesis that atomoxetine saliva levels (sATX) can be used to predict ATX brain Extracellular Fluid (bECF) levels and their pharmacological effects in healthy subjects and those with End-Stage Renal Disease (ESRD).

METHODS

The pharmacokinetics of atomoxetine after intravenous administration to rats with chemically induced acute and chronic renal impairments were investigated. A physiologically-based pharmacokinetic (PBPK) model was built and verified in rats using previously published measured atomoxetine levels in plasma and brain tissue. The rat PBPK model was then scaled to humans and verified using published measured atomoxetine levels in plasma, saliva, and bECF.

RESULTS

The rat PBPK model predicted the observed reduced atomoxetine clearance due to renal impairment in rats. The PBPK model predicted atomoxetine exposure in human plasma, sATX and bECF. Additionally, it predicted that ATX bECF levels needed to inhibit NET are achieved at 80 mg dose. In ESRD patients, the developed PBPK model predicted that the previously reported 65% increase in plasma exposure in these patients can be associated with a 63% increase in bECF. The PBPK simulations showed that there is a significant correlation between sATX and bECF in human.

CONCLUSION

Saliva levels can be used to predict atomoxetine pharmacological response.

Determination of In Vivo Biotransformation Kinetics Using Early-Time Biota Concentrations

Technical challenges have hampered the characterization of biotransformation kinetics-a critical link in understanding and predicting the toxicokinetics and ecotoxicology of organic compounds. A shortcut approach to characterize the in vivo biotransformation rate constant (k_M ) with incomplete pathway or metabolite details was proposed. The value of k_M can be derived as 2tln1fPC(t)) , with f_PC (t) being the molar equivalent fraction of the parent compound (PC) at an early time t in both constant exposure and decay source chemical uptake scenarios. The approximation-based k_M values agreed well with k_M values derived from rigorous fitting or toxicokinetic modeling (n = 42, root mean square error = 0.30) with accuracy exceeding those of typical toxicokinetic or partitioning models. The method is accurate when sampling time is adequately resolved (i.e., t < ln(2)/k_M ) but will likely produce biased k_M values with improper time-averaging. The approximate equation yields consistent theoretical expectations for fast and slow biotransformation reactions and is fully compatible with standard bioaccumulation and toxicity testing protocols. The simplification strategy circumvents statistical complications and numerical issues inherent in regressing or modeling the toxicokinetics of multimetabolite systems and may be adapted to similar problems at other physiological scales or ecotoxicological contexts. The method can help advance interspecies comparison of chemical metabolism and support the development of in vitro-in vivo extrapolations and in silico models needed for building next-generation ecological and health risk-assessment practices. Environ Toxicol Chem 2022;41:148-158. © 2021 SETAC.

Pseudoelementary Steps: A Key Concept and Tool for Studying the Kinetics and Mechanisms of Complex Chemical Systems

The concept of a pseudoelementary step (PEStep) is reviewed, a key concept for approaching the analysis of kinetics data and associated, underlying mechanisms of complex chemical systems. Following a brief Introduction, a definition of a PEStep is given: a PEStep is an initial building block for more complex reactions, that is a starting point for the initial analysis of the observed kinetics and then constructing initial, deliberately minimalistic mechanistic models for complex reactions. PESteps are, therefore and typically, composites of underlying elementary step reactions and can be very useful if not required for the inverse problem of discovering mechanisms from experimental observables for complex reactions. It is the use of PESteps in the inverse problem of mechanism determination that is a primary focus of this review. After a section detailing the results of a literature search of "pseudoelementary step" and related terms such as "pseudoelementary process", pedagogically illustrative examples are given of the use of the PEStep concept in approaching and elucidating the mechanisms of complex reactions. This review shows how the underlying elementary steps of a catalytic cycle were successfully uncovered via a PEStep approach, addresses the classic case of the use of PESteps in determining the mechanisms of oscillating reactions, and examines a well-studied case of an Ir(0)_n nanoparticle formation reaction. This latter example is illustrative in that the Ir(0)_n nanoparticle formation reaction consisting of thousands of underlying elementary steps that, however, can be treated initially kinetically as just two PESteps, a reduction in complexity of 3 orders of magnitude. Known weaknesses and caveats of the PEStep approach are also summarized and discussed. A short summary of Horituti's "Stoichiometric Number" concept is provided, a concept that would appear to merit further investigation and use in the study of complex reactions. Finally, a section is provided that lists a few, selected areas where the PEStep concept and methodology are expected to prove especially important in the future, and a Conclusions section is provided that lists 11 bullet points. The latter serves as a summary of this first review of the PEStep concept and its importance in dealing with the kinetics and in elucidating the mechanisms of more complex, multistep reactions.

CLUE: Exact maximal reduction of kinetic models by constrained lumping of differential equations

MOTIVATION

Detailed mechanistic models of biological processes can pose significant challenges for analysis and parameter estimations due to the large number of equations used to track the dynamics of all distinct configurations in which each involved biochemical species can be found. Model reduction can help tame such complexity by providing a lower-dimensional model in which each macro-variable can be directly related to the original variables.

RESULTS

We present CLUE, an algorithm for exact model reduction of systems of polynomial differential equations by constrained linear lumping. It computes the smallest dimensional reduction as a linear mapping of the state space such that the reduced model preserves the dynamics of user-specified linear combinations of the original variables. Even though CLUE works with nonlinear differential equations, it is based on linear algebra tools, which makes it applicable to high-dimensional models. Using case studies from the literature, we show how CLUE can substantially lower model dimensionality and help extract biologically intelligible insights from the reduction.

AVAILABILITY

An implementation of the algorithm and relevant resources to replicate the experiments herein reported are freely available for download at https://github.com/pogudingleb/CLUE.

SUPPLEMENTARY INFORMATION

Supplementary data are available at Bioinformatics online.

New insights on the ESIPT process based on solid-state data and state-of-the-art computational methods

Benzothiazole derivatives were used as models to study the excited-state intramolecular proton transfer (ESIPT) from an experimental and theoretical point of view. The experimental electronic and vibrational results were compared with a comprehensive selection of state-of-the-art computational methods in a workflow approach. The latter were performed based on modern techniques, such as DLPNO-CCSD(T), which gives the reference energies and current methodologies for ESIPT analysis, such as molecular dynamics and charge density difference testing. The theoretical vibrational results were focused on the stretch vibrational-mode of the hydroxyl group, which indicated a large increase in the intramolecular hydrogen bond strength, which facilitates the ESIPT process. Theoretically, the optimization of a large number of molecules shows that π-stacking plays a fundamental role in benzothiazole stabilization, with a remarkably strong intramolecular hydrogen bond. The potential energy surface of the ESIPT reactive benzothiazole (4HBS) has a clear transition state where ESIPT is easily observed with a large difference in energy between the enol and keto tautomer. Additionally, molecular dynamics showed that the ESIPT process occurs very fast. The tautomer appears around 8.7 fs and the enolic form is regenerated in just 24 fs, closing the Förster cycle. The calculated Stokes shift could be related to the ESIPT process and the experimental solid-state emission spectrum matched almost perfectly with the theoretical one. In contrast, for the non-ESIPT benzothiazole (4HBSN), the agreement between theory and experiment was limited, probably due to intermolecular interaction effects that are not considered in these calculations.

Sensitivity analysis of the reaction occurrence and recurrence times in steady-state biochemical networks

Continuous-time stationary Markov jump processes among discrete sites are encountered in disparate biochemical ambits. Sites and connecting dynamical events form a 'network' in which the sites are the available system's states, and the events are site-to-site transitions, or even neutral processes in which the system does not change site but the event is however detectable. Examples include conformational transitions in single biomolecules, stochastic chemical kinetics in the space of the molecules copy numbers, and even macroscopic steady-state reactive mixtures if one adopts the viewpoint of a tagged molecule (or even of a molecular moiety) whose state may change when it is involved in a chemical reaction. Among the variety of dynamical descriptors, here we focus on the first occurrence times (starting from a given site) and on the recurrence times of an event of interest. We develop the sensitivity analysis for the lowest moments of the statistical distribution of such times with respect to the rate constants of the network. In particular, simple expressions and inequalities allow us to establish a direct relationship between selective variation of rate constants and effect on average times and variances. As illustrative cases we treat the substrate inhibition in enzymatic catalysis in which a tagged enzyme molecule jumps between three states, and the basic two-site model of stochastic gene expression in which the single gene switches between active and inactive forms.

Sensitivity Based Order Reduction of a Chemical Membrane Degradation Model for Low-Temperature Proton Exchange Membrane Fuel Cells

The chemical degradation of the perfluorinated sulfonic acid (PFSA) ion-exchange membrane as a result of an attack from a radical species, originating as a by-product of the oxygen reduction reaction, represents a significant limiting factor in a wider adoption of low-temperature proton exchange membrane fuel cells (LT-PEMFCs). The efficient mathematical modeling of these processes is therefore a crucial step in the further development of proton exchange membrane fuel cells. Starting with an extensive kinetic modeling framework, describing the whole range of chemical processes leading to the membrane degradation, we use the mathematical method of sensitivity analysis to systematically reduce the number of both chemical species and reactions needed to efficiently and accurately describe the chemical degradation of the membrane. The analysis suggests the elimination of chemical reactions among the radical species, which is supported by the physicochemical consideration of the modeled reactions, while the degradation of Nafion backbone can be significantly simplified by lumping several individual species concentrations. The resulting reduced model features only 12 species coupled by 8 chemical reactions, compared to 19 species coupled by 23 reactions in the original model. The time complexity of the model, analyzed on the basis of its stiffness, however, is not significantly improved in the process. Nevertheless, the significant reduction in the model system size and number of parameters represents an important step in the development of a computationally efficient coupled model of various fuel cell degradation processes. Additionally, the demonstrated application of sensitivity analysis method shows a great potential for further use in the optimization of models of operation and degradation of fuel cell components.

Representation of Kinetics Models in Batch Flotation as Distributed First-Order Reactions

Four kinetic models are studied as first-order reactions with flotation rate distribution f(k): (i) deterministic nth-order reaction, (ii) second-order with Rectangular f(k), (iii) Rosin–Rammler, and (iv) Fractional kinetics. These models are studied because they are considered as alternatives to the first-order reactions. The first-order representation leads to the same recovery R(t) as in the original domain. The first-order R∞-f(k) are obtained by inspection of the R(t) formulae or by inverse Laplace Transforms. The reaction orders of model (i) are related to the shape parameters of first-order Gamma f(k)s. Higher reaction orders imply rate concentrations at k ≈ 0 in the first-order domain. Model (ii) shows reverse J-shaped first-order f(k)s. Model (iii) under stretched exponentials presents mounded first-order f(k)s, whereas model (iv) with derivative orders lower than 1 shows from reverse J-shaped to mounded first-order f(k)s. Kinetic descriptions that lead to the same R(t) cannot be differentiated between each other. However, the first-order f(k)s can be studied in a comparable domain.

Model Reduction for Kinetic Models of Biological Systems

High dimensionality continues to be a challenge in computational systems biology. The kinetic models of many phenomena of interest are high-dimensional and complex, resulting in large computational effort in the simulation. Model order reduction (MOR) is a mathematical technique that is used to reduce the computational complexity of high-dimensional systems by approximation with lower dimensional systems, while retaining the important information and properties of the full order system. Proper orthogonal decomposition (POD) is a method based on Galerkin projection that can be used for reducing the model order. POD is considered an optimal linear approach since it obtains the minimum squared distance between the original model and its reduced representation. However, POD may represent a restriction for nonlinear systems. By applying the POD method for nonlinear systems, the complexity to solve the nonlinear term still remains that of the full order model. To overcome the complexity for nonlinear terms in the dynamical system, an approach called the discrete empirical interpolation method (DEIM) can be used. In this paper, we discuss model reduction by POD and DEIM to reduce the order of kinetic models of biological systems and illustrate the approaches on some examples. Additional computational costs for setting up the reduced order system pay off for large-scale systems. In general, a reduced model should not be expected to yield good approximations if different initial conditions are used from that used to produce the reduced order model. We used the POD method of a kinetic model with different initial conditions to compute the reduced model. This reduced order model is able to predict the full order model for a variety of different initial conditions.

Reducing a model of sugar metabolism in peach to catch different patterns among genotypes

Several studies have been conducted to understand the dynamic of primary metabolisms in fruit by translating them into mathematics models. An ODE kinetic model of sugar metabolism has been developed by Desnoues et al. (2018) to simulate the accumulation of different sugars during peach fruit development. Two major drawbacks of this model are (a) the number of parameters to calibrate and (b) its integration time that can be long due to non-linearity and time-dependent input functions. Together, these issues hamper the use of the model for a large panel of genotypes, for which few data are available. In this paper, we present a model reduction scheme that explicitly addresses the specificity of genetic studies in that: (i) it yields a reduced model that is adapted to the whole expected genetic diversity (ii) it maintains network structure and variable identity, in order to facilitate biological interpretation. The proposed approach is based on the combination and the systematic evaluation of different reduction methods. Thus, we combined multivariate sensitivity analysis, structural simplification and timescale-based approaches to simplify the number and the structure of ordinary differential equations of the model. The original and reduced models were compared based on three criteria, namely the corrected Aikake Information Criterion (AIC_C), the calibration time and the expected error of the reduced model over a progeny of virtual genotypes. The resulting reduced model not only reproduces the predictions of the original one but presents many advantages including a reduced number of parameters to be estimated and shorter calibration time, opening new promising perspectives for genetic studies and virtual breeding. The validity of the reduced model was further evaluated by calibration on 30 additional genotypes of an inter-specific peach progeny for which few data were available.

Automated proper lumping for simplification of linear physiologically based pharmacokinetic systems

Physiologically based pharmacokinetic (PBPK) models are an important type of systems model used commonly in drug development before commencement of first-in-human studies. Due to structural complexity, these models are not easily utilised for future data-driven population pharmacokinetic (PK) analyses that require simpler models. In the current study we aimed to explore and automate methods of simplifying PBPK models using a proper lumping technique. A linear 17-state PBPK model for fentanyl was identified from the literature. Four methods were developed to search the optimal lumped model, including full enumeration (the reference method), non-adaptive random search (NARS), scree plot plus NARS, and simulated annealing (SA). For exploratory purposes, it was required that the total area under the fentanyl arterial concentration–time curve (AUC) between the lumped and original models differ by 0.002% at maximum. In full enumeration, a 4-state lumped model satisfying the exploratory criterion was found. In NARS, a lumped model with the same number of lumped states was found, requiring a large number of random samples. The scree plot provided a starting lumped model to NARS and the search completed within a short time. In SA, a 4-state lumped model was consistently delivered. In simplify an existing linear fentanyl PBPK model, SA was found to be robust and the most efficient and may be suitable for general application to other larger-scale linear systems. Ultimately, simplified PBPK systems with fundamental mechanisms may be readily used for data-driven PK analyses.

Usefulness of continuous probability distributions of rates for modelling radionuclide biokinetics in humans and animals

Modelling the biokinetics of radionuclide excretion or retention is important in nuclear medicine and following accidental/malicious radioactivity releases. Sums of discrete exponential decay rates are often used, but we hypothesized that continuous probability distributions (CPD) of decay rates can describe the data more parsimoniously and robustly. We tested this hypothesis on diverse human and animal data sets involving various radionuclides (including plutonium, strontium, caesium) measured in the laboratory and in regions contaminated by the Fukushima and Chernobyl nuclear accidents. We used four models on each data set: mono-exponential (ME) with one discrete decay rate, bi-exponential (BE) with two rates, gamma-exponential (GE) with a Gamma distribution of stretched-exponential rates, and power-decay (PD) with a Gamma distribution of power-decay rates. Information-theoretic model selection suggested that radionuclide biokinetics, e.g. for plutonium in humans, are often better described by CPD models like GE and PD, than by discrete rates (ME and BE). Extrapolation of models fitted to data at short times to longer times was frequently more robust for CPD formalisms. We suggest that using a set of several CPD and discrete-rate models, and comparing them by information-theoretic methods, is a promising strategy to enhance the analysis of radionuclide excretion and retention kinetics.

Insights into Dynamic Network States Using Metabolomic Data

Metabolomic data is the youngest of the high-throughput data types; however, it is potentially one of the most informative, as it provides a direct, quantitative biochemical phenotype. There are a number of ways in which metabolomic data can be analyzed in systems biology; however, the thermodynamic and kinetic relevance of these data cannot be overstated. Genome-scale metabolic network reconstructions provide a natural context to incorporate metabolomic data in order to provide insight into the condition-specific kinetic characteristics of metabolic networks. Herein we discuss how metabolomic data can be incorporated into constraint-based models in a flexible framework that enables scaling from small pathways to cell-scale models, while being able to accommodate coarse-grained to more detailed, allosteric interactions, all using the well-known principle of mass action.

Evaluating model reduction under parameter uncertainty

BACKGROUND

The dynamics of biochemical networks can be modelled by systems of ordinary differential equations. However, these networks are typically large and contain many parameters. Therefore model reduction procedures, such as lumping, sensitivity analysis and time-scale separation, are used to simplify models. Although there are many different model reduction procedures, the evaluation of reduced models is difficult and depends on the parameter values of the full model. There is a lack of a criteria for evaluating reduced models when the model parameters are uncertain.

RESULTS

We developed a method to compare reduced models and select the model that results in similar dynamics and uncertainty as the original model. We simulated different parameter sets from the assumed parameter distributions. Then, we compared all reduced models for all parameter sets using cluster analysis. The clusters revealed which of the reduced models that were similar to the original model in dynamics and variability. This allowed us to select the smallest reduced model that best approximated the full model. Through examples we showed that when parameter uncertainty was large, the model should be reduced further and when parameter uncertainty was small, models should not be reduced much.

CONCLUSIONS

A method to compare different models under parameter uncertainty is developed. It can be applied to any model reduction method. We also showed that the amount of parameter uncertainty influences the choice of reduced models.

Model reduction in mathematical pharmacology : Integration, reduction and linking of PBPK and systems biology models

In this paper we present a framework for the reduction and linking of physiologically based pharmacokinetic (PBPK) models with models of systems biology to describe the effects of drug administration across multiple scales. To address the issue of model complexity, we propose the reduction of each type of model separately prior to being linked. We highlight the use of balanced truncation in reducing the linear components of PBPK models, whilst proper lumping is shown to be efficient in reducing typically nonlinear systems biology type models. The overall methodology is demonstrated via two example systems; a model of bacterial chemotactic signalling in Escherichia coli and a model of extracellular regulatory kinase activation mediated via the extracellular growth factor and nerve growth factor receptor pathways. Each system is tested under the simulated administration of three hypothetical compounds; a strong base, a weak base, and an acid, mirroring the parameterisation of pindolol, midazolam, and thiopental, respectively. Our method can produce up to an 80% decrease in simulation time, allowing substantial speed-up for computationally intensive applications including parameter fitting or agent based modelling. The approach provides a straightforward means to construct simplified Quantitative Systems Pharmacology models that still provide significant insight into the mechanisms of drug action. Such a framework can potentially bridge pre-clinical and clinical modelling - providing an intermediate level of model granularity between classical, empirical approaches and mechanistic systems describing the molecular scale.

Memory functions reveal structural properties of gene regulatory networks

Gene regulatory networks (GRNs) control cellular function and decision making during tissue development and homeostasis. Mathematical tools based on dynamical systems theory are often used to model these networks, but the size and complexity of these models mean that their behaviour is not always intuitive and the underlying mechanisms can be difficult to decipher. For this reason, methods that simplify and aid exploration of complex networks are necessary. To this end we develop a broadly applicable form of the Zwanzig-Mori projection. By first converting a thermodynamic state ensemble model of gene regulation into mass action reactions we derive a general method that produces a set of time evolution equations for a subset of components of a network. The influence of the rest of the network, the bulk, is captured by memory functions that describe how the subnetwork reacts to its own past state via components in the bulk. These memory functions provide probes of near-steady state dynamics, revealing information not easily accessible otherwise. We illustrate the method on a simple cross-repressive transcriptional motif to show that memory functions not only simplify the analysis of the subnetwork but also have a natural interpretation. We then apply the approach to a GRN from the vertebrate neural tube, a well characterised developmental transcriptional network composed of four interacting transcription factors. The memory functions reveal the function of specific links within the neural tube network and identify features of the regulatory structure that specifically increase the robustness of the network to initial conditions. Taken together, the study provides evidence that Zwanzig-Mori projections offer powerful and effective tools for simplifying and exploring the behaviour of GRNs.

Gene regulatory networks are essential for cell fate specification and function. But the recursive links that comprise these networks often make determining their properties and behaviour complicated. Computational models of these networks can also be difficult to decipher. To reduce the complexity of such models we employ a Zwanzig-Mori projection approach. This allows a system of ordinary differential equations, representing a network, to be reduced to an arbitrary subnetwork consisting of part of the initial network, with the rest of the network (bulk) captured by memory functions. These memory functions account for the bulk by describing signals that return to the subnetwork after some time, having passed through the bulk. We show how this approach can be used to simplify analysis and to probe the behaviour of a gene regulatory network. Applying the method to a transcriptional network in the vertebrate neural tube reveals previously unappreciated properties of the network. By taking advantage of the structure of the memory functions we identify interactions within the network that are unnecessary for sustaining correct patterning. Upon further investigation we find that these interactions are important for conferring robustness to variation in initial conditions. Taken together we demonstrate the validity and applicability of the Zwanzig-Mori projection approach to gene regulatory networks.

Understanding and reducing complex systems pharmacology models based on a novel input-response index

A growing understanding of complex processes in biology has led to large-scale mechanistic models of pharmacologically relevant processes. These models are increasingly used to study the response of the system to a given input or stimulus, e.g., after drug administration. Understanding the input-response relationship, however, is often a challenging task due to the complexity of the interactions between its constituents as well as the size of the models. An approach that quantifies the importance of the different constituents for a given input-output relationship and allows to reduce the dynamics to its essential features is therefore highly desirable. In this article, we present a novel state- and time-dependent quantity called the input-response index that quantifies the importance of state variables for a given input-response relationship at a particular time. It is based on the concept of time-bounded controllability and observability, and defined with respect to a reference dynamics. In application to the brown snake venom-fibrinogen (Fg) network, the input-response indices give insight into the coordinated action of specific coagulation factors and about those factors that contribute only little to the response. We demonstrate how the indices can be used to reduce large-scale models in a two-step procedure: (i) elimination of states whose dynamics have only minor impact on the input-response relationship, and (ii) proper lumping of the remaining (lower order) model. In application to the brown snake venom-fibrinogen network, this resulted in a reduction from 62 to 8 state variables in the first step, and a further reduction to 5 state variables in the second step. We further illustrate that the sequence, in which a recursive algorithm eliminates and/or lumps state variables, has an impact on the final reduced model. The input-response indices are particularly suited to determine an informed sequence, since they are based on the dynamics of the original system. In summary, the novel measure of importance provides a powerful tool for analysing the complex dynamics of large-scale systems and a means for very efficient model order reduction of nonlinear systems.

Maximal aggregation of polynomial dynamical systems

Ordinary differential equations (ODEs) with polynomial derivatives are a fundamental tool for understanding the dynamics of systems across many branches of science, but our ability to gain mechanistic insight and effectively conduct numerical evaluations is critically hindered when dealing with large models. Here we propose an aggregation technique that rests on two notions of equivalence relating ODE variables whenever they have the same solution (backward criterion) or if a self-consistent system can be written for describing the evolution of sums of variables in the same equivalence class (forward criterion). A key feature of our proposal is to encode a polynomial ODE system into a finitary structure akin to a formal chemical reaction network. This enables the development of a discrete algorithm to efficiently compute the largest equivalence, building on approaches rooted in computer science to minimize basic models of computation through iterative partition refinements. The physical interpretability of the aggregation is shown on polynomial ODE systems for biochemical reaction networks, gene regulatory networks, and evolutionary game theory.

Statistical physics approaches to subnetwork dynamics in biochemical systems

We apply a Gaussian variational approximation to model reduction in large biochemical networks of unary and binary reactions. We focus on a small subset of variables (subnetwork) of interest, e.g. because they are accessible experimentally, embedded in a larger network (bulk). The key goal is to write dynamical equations reduced to the subnetwork but still retaining the effects of the bulk. As a result, the subnetwork-reduced dynamics contains a memory term and an extrinsic noise term with non-trivial temporal correlations. We first derive expressions for this memory and noise in the linearized (Gaussian) dynamics and then use a perturbative power expansion to obtain first order nonlinear corrections. For the case of vanishing intrinsic noise, our description is explicitly shown to be equivalent to projection methods up to quadratic terms, but it is applicable also in the presence of stochastic fluctuations in the original dynamics. An example from the epidermal growth factor receptor signalling pathway is provided to probe the increased prediction accuracy and computational efficiency of our method.

Fluctuation relations between hierarchical kinetically equivalent networks with Arrhenius-type transitions and their roles in systems and structural biology

The legality of using simple kinetic schemes to determine the stochastic properties of a complex system depends on whether the fluctuations generated from hierarchical equivalent schemes are consistent with one another. To analyze this consistency, we perform lumping processes on the stochastic differential equations and the generalized fluctuation-dissipation theorem and apply them to networks with the frequently encountered Arrhenius-type transition rates. The explicit Langevin force derived from those networks enables us to calculate the state fluctuations caused by the intrinsic and extrinsic noises on the free energy surface and deduce their relations between kinetically equivalent networks. In addition to its applicability to wide classes of network related systems, such as those in structural and systems biology, the result sheds light on the fluctuation relations for general physical variables in Keizer's canonical theory.

Methods of Model Reduction for Large-Scale Biological Systems: A Survey of Current Methods and Trends

Complex models of biochemical reaction systems have become increasingly common in the systems biology literature. The complexity of such models can present a number of obstacles for their practical use, often making problems difficult to intuit or computationally intractable. Methods of model reduction can be employed to alleviate the issue of complexity by seeking to eliminate those portions of a reaction network that have little or no effect upon the outcomes of interest, hence yielding simplified systems that retain an accurate predictive capacity. This review paper seeks to provide a brief overview of a range of such methods and their application in the context of biochemical reaction network models. To achieve this, we provide a brief mathematical account of the main methods including timescale exploitation approaches, reduction via sensitivity analysis, optimisation methods, lumping, and singular value decomposition-based approaches. Methods are reviewed in the context of large-scale systems biology type models, and future areas of research are briefly discussed.

QSP Toolbox: Computational Implementation of Integrated Workflow Components for Deploying Multi-Scale Mechanistic Models

Quantitative systems pharmacology (QSP) modeling has become increasingly important in pharmaceutical research and development, and is a powerful tool to gain mechanistic insights into the complex dynamics of biological systems in response to drug treatment. However, even once a suitable mathematical framework to describe the pathophysiology and mechanisms of interest is established, final model calibration and the exploration of variability can be challenging and time consuming. QSP models are often formulated as multi-scale, multi-compartment nonlinear systems of ordinary differential equations. Commonly accepted modeling strategies, workflows, and tools have promise to greatly improve the efficiency of QSP methods and improve productivity. In this paper, we present the QSP Toolbox, a set of functions, structure array conventions, and class definitions that computationally implement critical elements of QSP workflows including data integration, model calibration, and variability exploration. We present the application of the toolbox to an ordinary differential equations-based model for antibody drug conjugates. As opposed to a single stepwise reference model calibration, the toolbox also facilitates simultaneous parameter optimization and variation across multiple in vitro, in vivo, and clinical assays to more comprehensively generate alternate mechanistic hypotheses that are in quantitative agreement with available data. The toolbox also includes scripts for developing and applying virtual populations to mechanistic exploration of biomarkers and efficacy. We anticipate that the QSP Toolbox will be a useful resource that will facilitate implementation, evaluation, and sharing of new methodologies in a common framework that will greatly benefit the community.

A combined model reduction algorithm for controlled biochemical systems

BACKGROUND

Systems Biology continues to produce increasingly large models of complex biochemical reaction networks. In applications requiring, for example, parameter estimation, the use of agent-based modelling approaches, or real-time simulation, this growing model complexity can present a significant hurdle. Often, however, not all portions of a model are of equal interest in a given setting. In such situations methods of model reduction offer one possible approach for addressing the issue of complexity by seeking to eliminate those portions of a pathway that can be shown to have the least effect upon the properties of interest.

METHODS

In this paper a model reduction algorithm bringing together the complementary aspects of proper lumping and empirical balanced truncation is presented. Additional contributions include the development of a criterion for the selection of state-variable elimination via conservation analysis and use of an 'averaged' lumping inverse. This combined algorithm is highly automatable and of particular applicability in the context of 'controlled' biochemical networks.

RESULTS

The algorithm is demonstrated here via application to two examples; an 11 dimensional model of bacterial chemotaxis in Escherichia coli and a 99 dimensional model of extracellular regulatory kinase activation (ERK) mediated via the epidermal growth factor (EGF) and nerve growth factor (NGF) receptor pathways. In the case of the chemotaxis model the algorithm was able to reduce the model to 2 state-variables producing a maximal relative error between the dynamics of the original and reduced models of only 2.8% whilst yielding a 26 fold speed up in simulation time. For the ERK activation model the algorithm was able to reduce the system to 7 state-variables, incurring a maximal relative error of 4.8%, and producing an approximately 10 fold speed up in the rate of simulation. Indices of controllability and observability are additionally developed and demonstrated throughout the paper. These provide insight into the relative importance of individual reactants in mediating a biochemical system's input-output response even for highly complex networks.

CONCLUSIONS

Through application, this paper demonstrates that combined model reduction methods can produce a significant simplification of complex Systems Biology models whilst retaining a high degree of predictive accuracy. In particular, it is shown that by combining the methods of proper lumping and empirical balanced truncation it is often possible to produce more accurate reductions than can be obtained by the use of either method in isolation.

Lumping evolutionary game dynamics on networks

We study evolutionary game dynamics on networks (EGN), where players reside in the vertices of a graph, and games are played between neighboring vertices. The model is described by a system of ordinary differential equations which depends on players payoff functions, as well as on the adjacency matrix of the underlying graph. Since the number of differential equations increases with the number of vertices in the graph, the analysis of EGN becomes hard for large graphs. Building on the notion of lumpability for Markov chains, we identify conditions on the network structure allowing to reduce the original graph. In particular, we identify a partition of the vertex set of the graph and show that players in the same block of a lumpable partition have equivalent dynamical behaviors, whenever their payoff functions and initial conditions are equivalent. Therefore, vertices belonging to the same partition block can be merged into a single vertex, giving rise to a reduced graph and consequently to a simplified system of equations. We also introduce a tighter condition, called strong lumpability, which can be used to identify dynamical symmetries in EGN which are related to the interchangeability of players in the system.

Michaelis-Menten dynamics in protein subnetworks

To understand the behaviour of complex systems, it is often necessary to use models that describe the dynamics of subnetworks. It has previously been established using projection methods that such subnetwork dynamics generically involves memory of the past and that the memory functions can be calculated explicitly for biochemical reaction networks made up of unary and binary reactions. However, many established network models involve also Michaelis-Menten kinetics, to describe, e.g., enzymatic reactions. We show that the projection approach to subnetwork dynamics can be extended to such networks, thus significantly broadening its range of applicability. To derive the extension, we construct a larger network that represents enzymes and enzyme complexes explicitly, obtain the projected equations, and finally take the limit of fast enzyme reactions that gives back Michaelis-Menten kinetics. The crucial point is that this limit can be taken in closed form. The outcome is a simple procedure that allows one to obtain a description of subnetwork dynamics, including memory functions, starting directly from any given network of unary, binary, and Michaelis-Menten reactions. Numerical tests show that this closed form enzyme elimination gives a much more accurate description of the subnetwork dynamics than the simpler method that represents enzymes explicitly and is also more efficient computationally.

A Six-Stage Workflow for Robust Application of Systems Pharmacology

Quantitative and systems pharmacology (QSP) is increasingly being applied in pharmaceutical research and development. One factor critical to the ultimate success of QSP is the establishment of commonly accepted language, technical criteria, and workflows. We propose an integrated workflow that bridges conceptual objectives with underlying technical detail to support the execution, communication, and evaluation of QSP projects.