Mathematical analysis of multienzyme systems. I. Modelling of the glycolysis of human erythrocytes.

A mathematical model for the glycolysis of human erythrocytes is presented which takes into account ATP-synthesis and -consumption. A set of three differential equations describes the steady states and the time-dependent changes of the metabolite concentrations under blood storage conditions. For a given parameter combination there are in general three stationary points of the system, one of which is unstable. At a low ATP-need the ATP-level is relatively constant for variations in the rate constant of the ATP-consuming processes. Above a critical level of the energy consumption the system breaks down. An important role of the 2.3P2G-bypass of the erythrocytes is its action as an "energy buffer", wasting ATP in case of ATP-overproduction and producing ATP in case of underproduction. A parameter combination consistent with the data on the isolated enzymes was found which gives a good agreement of theoretical predictions with the measured metabolite concentrations. Under blood preservation conditions the difference of the rates of ATP-production and -consumption is the most important factor for a high ATP-level over long periods.     

An integral equation approach to measuring turnover in nonsteady compartmental and distributed systems

Physiological systems are often modelled by a set of compartments. Alternatively they can be described by the diffusion-convection-reaction equations governing distributed systems. The problem considered here is that of identifying a continuously changing input of some metabolite )tracee), endogenous to the system and hence inaccessible, when a nonlinear or time-varying component is also introduced into the loss parameter, as for example through feedback mechanisms. A tracer is used to determine the steady-state impulse response under time-invariant, linear conditions. A known input of tracer is also administered when the system is driven out of steady state. The integral equations developed utilize the predetermined impulse response, the measured concentrations of both tracer and tracee (output) in some region of the system to estimate the changing loss parameter and the unknown input in a continuous fashion.     

Non-Lethal Control of the Cariogenic Potential of an Agent-Based Model for Dental Plaque

Dental caries or tooth decay is a prevalent global disease whose causative agent is the oral biofilm known as plaque. According to the ecological plaque hypothesis, this biofilm becomes pathogenic when external challenges drive it towards a state with a high proportion of acid-producing bacteria. Determining which factors control biofilm composition is therefore desirable when developing novel clinical treatments to combat caries, but is also challenging due to the system complexity and the existence of multiple bacterial species performing similar functions. Here we employ agent-based mathematical modelling to simulate a biofilm consisting of two competing, distinct types of bacterial populations, each parameterised by their nutrient uptake and aciduricity, periodically subjected to an acid challenge resulting from the metabolism of dietary carbohydrates. It was found that one population was progressively eliminated from the system to give either a benign or a pathogenic biofilm, with a tipping point between these two fates depending on a multiplicity of factors relating to microbial physiology and biofilm geometry. Parameter sensitivity was quantified by individually varying the model parameters against putative experimental measures, suggesting non-lethal interventions that can favourably modulate biofilm composition. We discuss how the same parameter sensitivity data can be used to guide the design of validation experiments, and argue for the benefits of in silico modelling in providing an additional predictive capability upstream from in vitro experiments.     

Isotope-assisted metabolic flux analysis as an equality-constrained nonlinear program for improved scalability and robustness

Stable isotope-assisted metabolic flux analysis (MFA) is a powerful method to estimate carbon flow and partitioning in metabolic networks. At its core, MFA is a parameter estimation problem wherein the fluxes and metabolite pool sizes are model parameters that are estimated, via optimization, to account for measurements of steady-state or isotopically-nonstationary isotope labeling patterns. As MFA problems advance in scale, they require efficient computational methods for fast and robust convergence. The structure of the MFA problem enables it to be cast as an equality-constrained nonlinear program (NLP), where the equality constraints are constructed from the MFA model equations, and the objective function is defined as the sum of squared residuals (SSR) between the model predictions and a set of labeling measurements. This NLP can be solved by using an algebraic modeling language (AML) that offers state-of-the-art optimization solvers for robust parameter estimation and superior scalability to large networks. When implemented in this manner, the optimization is performed with no distinction between state variables and model parameters. During each iteration of such an optimization, the system state is updated instead of being calculated explicitly from scratch, and this occurs concurrently with improvement in the model parameter estimates. This optimization approach starkly contrasts with traditional “shooting” methods where the state variables and model parameters are kept distinct and the system state is computed afresh during each iteration of a stepwise optimization. Our NLP formulation uses the MFA modeling framework of Wiechert et al. [1], which is amenable to incorporation of the model equations into an NLP. The NLP constraints consist of balances on either elementary metabolite units (EMUs) or cumomers. In this formulation, both the steady-state and isotopically-nonstationary MFA (inst-MFA) problems may be solved as an NLP. For the inst-MFA case, the ordinary differential equation (ODE) system describing the labeling dynamics is transcribed into a system of algebraic constraints for the NLP using collocation. This large-scale NLP may be solved efficiently using an NLP solver implemented on an AML. In our implementation, we used the reduced gradient solver CONOPT, implemented in the General Algebraic Modeling System (GAMS). The NLP framework is particularly advantageous for inst-MFA, scaling well to large networks with many free parameters, and having more robust convergence properties compared to the shooting methods that compute the system state and sensitivities at each iteration. Additionally, this NLP approach supports the use of tandem-MS data for both steady-state and inst-MFA when the cumomer framework is used. We assembled a software, eiFlux, written in Python and GAMS that uses the NLP approach and supports both steady-state and inst-MFA. We demonstrate the effectiveness of the NLP formulation on several examples, including a genome-scale inst-MFA model, to highlight the scalability and robustness of this approach. In addition to typical inst-MFA applications, we expect that this framework and our associated software, eiFlux, will be particularly useful for applying inst-MFA to complex MFA models, such as those developed for eukaryotes (e.g. algae) and co-cultures with multiple cell types.     

High-pressure freezing of soybean nodules leads to an improved preservation of ultrastructure

High-pressure freezing of chemically untreated nodules of soybean (Glycine max (L.) Merr.), in sharp contrast to chemical fixation and prefixation, appears to preserve the ultrastructure close to the native state. This is supported by the observation that the peribacteroid membrane of high-pressure-frozen samples is tightly wrapped around the bacteroids, a finding that is fully consistent with the current views on the physiology of oxygen and metabolite transport between plant cytosol and bacteroids. In soybean root nodules, the plant tissue and the enclosed bacteria are so dissimilar that conventional aldehyde-fixation procedures are unable to preserve the overall native ultrastructure. This was demonstrated by high-pressure freezing of nodules that had been pre-fixed in glutaraldehyde at various buffer molalities: no buffer strength tested preserved all ultrastructural aspects that could be seen after high-pressure freezing of chemically untreated nodules.     

Numerical methods and parameter estimation of a structured population model with discrete events in the life history

We consider two numerical methods for the solution of a physiologically structured population (PSP) model with multiple life stages and discrete event reproduction. The model describes the dynamic behaviour of a predator-prey system consisting of rotifers predating on algae. The nitrate limited algal prey population is modelled unstructured and described by an ordinary differential equation (ODE). The formulation of the rotifer dynamics is based on a simple physiological model for their two life stages, the egg and the adult stage. An egg is produced when an energy buffer reaches a threshold value. The governing equations are coupled partial differential equations (PDE) with initial and boundary conditions. The population models together with the equation for the dynamics of the nutrient result in a chemostat model. Experimental data are used to estimate the model parameters. The results obtained with the explicit finite difference (FD) technique compare well with those of the Escalator Boxcar Train (EBT) method. This justifies the use of the fast FD method for the parameter estimation, a procedure which involves repeated solution of the model equations.     

Modular synthesis of Petri nets for modeling flexible manufacturing systems

This paper proposes a modular Petri net synthesis method for modeling flexible manufacturing systems based on synchronization among control processes of the manufacturing resources (such as robots and machines). In the method, the target system is modeled in a bottom-up and uniform manner by first describing the system's control processes using strongly connected state machines (SCSMs) as the basic modules. Each SCSM may contain multiple tokens to represent resources from the same type such as spaces in a buffer. Next, the common transitions and common transition subnets of the modules are merged to represent their synchronization. The system model constructed is proven to be conservative and thus bounded. Moreover, a restricted class of merged nets is proven to be live and reversible. For general classes of merged nets, this paper shows theorems that easily calculateP-invariants of the final net without solving the linear system equations. TheseP-invariants can be used to help in verifying the model's qualitative properties such as liveness.     

Entropy consumption in primary photosynthesis

Knox and Parson have objected to our previous conclusion on possible negative entropy production during primary photochemistry, i.e., from photon absorption to primary charge separation, by considering a pigment system in which primary photochemistry is not specifically considered. This approach does not address our proposal. They suggest that when a pigment absorbs light and passes to an excited state, its entropy increases by hν/T. This point is discussed in two ways: (i) from considerations based on the energy gap law for excited state relaxation; (ii) using classical thermodynamics, in which free energy is introduced into the pigment (antenna) system by photon absorption. Both approaches lead us to conclude that the excited state and the ground state are isoentropic, in disagreement with Knox and Parson. A discussion on total entropy changes specifically during the charge separation process itself indicates that this process may be almost isoentropic and thus our conclusions on possible negentropy production associated with the sequence of reactions which go from light absorption to the first primary charge separation event, due to its very high thermodynamic efficiency, remain unchanged.     

The use of different test systems with yeasts for the evaluation of chemically induced gene conversions and gene mutations

The wide variety of genetic alterations that can be induced in human populations when exposed to chemical genotoxic substances present in our environment may be predictable using laboratory organisms such as yeasts.In the present paper methodologies are described for analysing the genetic effects induced by a well known chemical mutagen, ethyl methanesulfonate (EMS). Haploid or diploid yeast cells have been treated in vitro, in buffer or in the presence of mouse liver microsomes, and in vivo, in the peritoneum of the mouse (host-mediated assay).With these different methods of assaying the genetic activity of a compound, its metabolic activation occuring in the mammalian body is taken into account: this might lead to a more reliable extrapolation of data from laboratory experiments to man.The relationships between doses and frequencies of the induced genetic effects are described by equations obtained after regression analysis of the data, thus allowing a quantitative comparison among different methodologies and different genetic systems.One genetic system analyzed is represented by forward-mutations scored phenotypically on a non selective medium. Mutations induced in five loci with different sensitivity and average data of mutation-induction per locus have been derived. The second genetic system was provided by scoring on a selective medium mitotic gene-conversions induced in two loci with different kinetics.Haploid cells of Schizosaccharomyces pombe and diploid cells of Saccharomyces cerevisiae were submitted to analysis for the evaluation of gene-mutations and gene-conversions respectively.     

Advantage of storage in a fluctuating environment

We will elaborate the evolutionary course of an ecosystem consisting of a population in a chemostat environment with periodically fluctuating nutrient supply. The organisms that make up the population consist of structural biomass and energy storage compartments. In a constant chemostat environment a species without energy storage always out-competes a species with energy reserves. This hinders evolution of species with storage from those without storage. Using the adaptive dynamics approach for non-equilibrium ecological systems we will show that in a fluctuating environment there are multiple stable evolutionary singular strategies (ss's): one for a species without, and one for a species with energy storage. The evolutionary end-point depends on the initial evolutionary state. We will formulate the invasion fitness in terms of Floquet multipliers for the oscillating non-autonomous system. Bifurcation theory is used to study points where due to evolutionary development by mutational steps, the long-term dynamics of the ecological system changes qualitatively. To that end, at the ecological time scale, the trait value at which invasion of a mutant into a resident population becomes possible can be calculated using numerical bifurcation analysis where the trait is used as the free parameter, because it is just a bifurcation point. In a constant environment there is a unique stable equilibrium for one species following the “competitive exclusion” principle. In contrast, due to the oscillatory dynamics on the ecological time scale two species may coexist. That is, non-equilibrium dynamics enhances biodiversity. However, we will show that this coexistence is not stable on the evolutionary time scale and always one single species survives.     

Transport of proteins into mitochondria: a potassium diffusion potential is able to drive the import of ADP/ATP carrier.

The transfer of cytoplasmically synthesized precursor proteins into or across the inner mitochondrial membrane is dependent on energization of the membrane. To investigate the role of this energy requirement, a buffer system was developed in which efficient import of ADP/ATP carrier into mitochondria from the receptor-bound state occurred. This import was rapid and was dependent on divalent cations, whereas the binding of precursor proteins to the mitochondrial surface was slow and was independent of added divalent cations. Using this buffer system, the import of ADP/ATP carrier could be driven by a valinomycin-induced potassium diffusion potential. The protonophore carbonylcyanide m-chlorophenyl-hydrazone was not able to abolish this import. Imposition of a delta pH did not stimulate the import. We conclude that the membrane potential delta psi itself and not the total protonmotive force delta p is the required energy source.     

GEPASI: a software package for modelling the dynamics, steady states and control of biochemical and other systems

GEPASI is a software system for modelling chemical and biochemicalreaction networks on computers running Microsoft Windows. Forany system of up to 45 metabolites and 45 reactions, each withany user-defined or one of 35 predefined rate equations, onecan produce trajectories of the metabolite concentrations andobtain a steady state (if it does exist). When steady-statesolutions are produced, elasticity and control coefficients,as defined in metabolic control analysis, are calculated. GEPASIalso allows the automatic generation of a sequence of simulationswith different combinations of parameter values, effectivelyscanning a hyper-solid in parameter space. Together with theability to produce user-defined columnar data files, these featuresallow for both very quick and systematic study of biochemicalpathway models. The source code (in C) is available on requestfrom the author, and while the user interface is dependent onhaving MS-Windows as the operating system, the numerical partis portable to other operating systems. GEPASI is suitable bothfor research and educational purposes. Although GEPASI was writtenwith biochemical pathways in mind, it can equally be used tosimulate other dynamical systems.     

Interrogating the effect of enzyme kinetics on metabolism using differentiable constraint-based models

Metabolic models are typically characterized by a large number of parameters. Traditionally, metabolic control analysis is applied to differential equation-based models to investigate the sensitivity of predictions to parameters. A corresponding theory for constraint-based models is lacking, due to their formulation as optimization problems. Here, we show that optimal solutions of optimization problems can be efficiently differentiated using constrained optimization duality and implicit differentiation. We use this to calculate the sensitivities of predicted reaction fluxes and enzyme concentrations to turnover numbers in an enzyme-constrained metabolic model of Escherichia coli. The sensitivities quantitatively identify rate limiting enzymes and are mathematically precise, unlike current finite difference based approaches used for sensitivity analysis. Further, efficient differentiation of constraint-based models unlocks the ability to use gradient information for parameter estimation. We demonstrate this by improving, genome-wide, the state-of-the-art turnover number estimates for E. coli. Finally, we show that this technique can be generalized to arbitrarily complex models. By differentiating the optimal solution of a model incorporating both thermodynamic and kinetic rate equations, the effect of metabolite concentrations on biomass growth can be elucidated. We benchmark these metabolite sensitivities against a large experimental gene knockdown study, and find good alignment between the predicted sensitivities and in vivo metabolome changes. In sum, we demonstrate several applications of differentiating optimal solutions of constraint-based metabolic models, and show how it connects to classic metabolic control analysis.     

A possible mechanism of the generation of singlet molecular oxygen in NADPH-dependent microsomal lipid peroxidation

A simplified system, consisting of NADPH, Fe³⁺-ADP, EDTA, liposomes, NADPH-cytochrome c reductase and Tris · HCl buffer (pH 6.8), has been employed in studies of the generation of singlet oxygen in NADPH-dependent microsomal lipid peroxidation.The light emitted by the system involves ¹Δg type molecular oxygen identifiable by its characteristic emission spectrum and its behavior with β-carotene. The generation of another excited species (a compound in the triplet state) could be demonstrated in this system by changes of light intensity and emission spectra which arise from photosensitizer (9, 10-dibromoanthracene sulfonate, eosin, Rose-Bengal)-mediated energy transfers.Chemiluminescence in the visible region was markedly quenched by various radical trappers and by an inhibitor of NADPH-cytochrome c reductase, but not by superoxide dismutase. During the early stage of lipid peroxidation, the intensity of chemiluminescence was proportional to the square of the concentration of lipid peroxide.These characteristics suggest that singlet oxygen and a compound in the triplet state (probably a carbonyl compound) are generated by a self-reaction of lipid peroxy radicals.     

Irreversible transitions in the 6-phosphofructokinase/fructose 1,6-bisphosphatase cycle.

The dynamics of the fructose 6-phosphate fructose-1,6-bisphosphate cycle operating in an open and homogeneous system reconstituted from purified enzymes was extensively studied. In addition to 6-phosphofructokinase and fructose-1,6-bisphosphatase, pyruvate kinase, adenylate kinae and glucose-6-phosphate isomerase were involved. In that multi-enzyme system, the main source of non-linearity is the reciprocal effect of AMP on the activities of 6-phosphofructokinase and fructose-1,6-bisphosphatase. Depending upon the experimental parameter values, stable attractors, various types of multiple states and sustained oscillations were shown to occur. In the present report we show that irreversible transitions are also likely to occur for realistic operating conditions. Two parameters of the system, that is the adenylate energy charge of the influx and the fructose-1,6-bisphosphatase maximal activity, are potential candidates to provoke such irreversible transitions from one steady state to the other: (a) when varying the maximal activity of fructose-1,6-bisphosphatase, the system can jump irreversibly from a low to a high stable steady state, and (b) when the adenylate energy charge of the influx is the changing parameter, irreversible transitions occur from a high stable steady state to a stable oscillatory state (limit cycle motion). This behavior can be predicted by constructing the loci of limit points and Hopf bifurcation points.     

Modeling Granulomas in Response to Infection in the Lung

Alveolar macrophages play a large role in the innate immune response of the lung. However, when these highly immune-regulatory cells are unable to eradicate pathogens, the adaptive immune system, which includes activated macrophages and lymphocytes, particularly T cells, is called upon to control the pathogens. This collection of immune cells surrounds, isolates and quarantines the pathogen, forming a small tissue structure called a granuloma for intracellular pathogens like Mycobacterium tuberculosis (Mtb). In the present work we develop a mathematical model of the dynamics of a granuloma by a system of partial differential equations. The ‘strength’ of the adaptive immune response to infection in the lung is represented by a parameter α, the flux rate by which T cells and M1 macrophages that immigrated from the lymph nodes enter into the granuloma through its boundary. The parameter α is negatively correlated with the ‘switching time’, namely, the time it takes for the number of M1 type macrophages to surpass the number of infected, M2 type alveolar macrophages. Simulations of the model show that as α increases the radius of the granuloma and bacterial load in the granuloma both decrease. The model is used to determine the efficacy of potential host-directed therapies in terms of the parameter α, suggesting that, with fixed dosing level, an infected individual with a stronger immune response will receive greater benefits in terms of reducing the bacterial load.     

Advantage of storage in a fluctuating environment

We will elaborate the evolutionary course of an ecosystem consisting of a population in a chemostat environment with periodically fluctuating nutrient supply. The organisms that make up the population consist of structural biomass and energy storage compartments. In a constant chemostat environment a species without energy storage always out-competes a species with energy reserves. This hinders evolution of species with storage from those without storage. Using the adaptive dynamics approach for non-equilibrium ecological systems we will show that in a fluctuating environment there are multiple stable evolutionary singular strategies (ss's): one for a species without, and one for a species with energy storage. The evolutionary end-point depends on the initial evolutionary state. We will formulate the invasion fitness in terms of Floquet multipliers for the oscillating non-autonomous system. Bifurcation theory is used to study points where due to evolutionary development by mutational steps, the long-term dynamics of the ecological system changes qualitatively. To that end, at the ecological time scale, the trait value at which invasion of a mutant into a resident population becomes possible can be calculated using numerical bifurcation analysis where the trait is used as the free parameter, because it is just a bifurcation point. In a constant environment there is a unique stable equilibrium for one species following the "competitive exclusion" principle. In contrast, due to the oscillatory dynamics on the ecological time scale two species may coexist. That is, non-equilibrium dynamics enhances biodiversity. However, we will show that this coexistence is not stable on the evolutionary time scale and always one single species survives.     

Bifurcation analysis of a neural network model

This paper describes the analysis of the well known neural network model by Wilson and Cowan. The neural network is modeled by a system of two ordinary differential equations that describe the evolution of average activities of excitatory and inhibitory populations of neurons. We analyze the dependence of the model's behavior on two parameters. The parameter plane is partitioned into regions of equivalent behavior bounded by bifurcation curves, and the representative phase diagram is constructed for each region. This allows us to describe qualitatively the behavior of the model in each region and to predict changes in the model dynamics as parameters are varied. In particular, we show that for some parameter values the system can exhibit long-period oscillations. A new type of dynamical behavior is also found when the system settles down either to a stationary state or to a limit cycle depending on the initial point.