Efficient parametric analysis of the chemical master equation through model order reduction

ABSTRACT: BACKGROUND: Stochastic biochemical reaction networks are commonly modelled by the chemical master equation, and can be simulated as first order linear differential equations through a finite state projection. Due to the very high state space dimension of these equations, numerical simulations are computationally expensive. This is a particular problem for analysis tasks requiring repeated simulations for different parameter values. Such tasks are computationally expensive to the point of infeasibility with the chemical master equation. RESULTS: In this article, we apply parametric model order reduction techniques in order to construct accurate low-dimensional parametric models of the chemical master equation. These surrogate models can be used in various parametric analysis task such as identifiability analysis, parameter estimation, or sensitivity analysis. As biological examples, we consider two models for gene regulation networks, a bistable switch and a network displaying stochastic oscillations. CONCLUSIONS: The results show that the parametric model reduction yields efficient models of stochastic biochemical reaction networks, and that these models can be useful for systems biology applications involving parametric analysis problems such as parameter exploration, optimization, estimation or sensitivity analysis.     

Salting-out in the aqueous single-protein solution: the effect of shape factor

A molecular-thermodynamic model is developed to describe salt-induced protein precipitation. The protein-protein interaction goes through the potential of mean force. An equation of state is derived based on the generalized van der Waals partition function. The attractive term including the potential of mean force is perturbed by the statistical mechanical perturbation theory. The precipitation behaviors are studied by calculating the partition coefficient with various conditions such as the ionic strength and the shape of protein. Our results show that the protein shape plays a significant role in the protein precipitation behavior.     

Sensitivity function-based model reduction: A bacterial gene expression case study

Mathematical models used to predict the behavior of genetically modified organisms require 1). a (rather) large number of state variables, and 2). complicated kinetic expressions containing a large number of parameters. Since these models are hardly identifiable and of limited use in model-based optimization and control strategies, a generic methodology based on sensitivity function analysis is presented to reduce the model complexity at the level of the kinetics, while maintaining high prediction power. As a case study to illustrate the method and results obtained, the influence of the dissolved oxygen concentration on the cytN gene expression in the bacterium Azospirillum brasilense Sp7 is modeled. As a first modeling approach, available mechanistic knowledge is incorporated into a mass balance equation model with 3 states and 14 parameters. The large differences in order of magnitude of the model parameters identified on the available experimental data indicate 1). possible structural problems in the kinetic model and, associated with this, 2). a possibly too high number of model parameters. A careful sensitivity function analysis reveals that a reduced model with only seven parameters is almost as accurate as the original model.     

Theoretical nuclear magnetic resonance-lineshape study of a two-site chemical exchange model of sodium ions (I = 3/2) in an intracellular environment.

We present a theoretical calculation of the lineshape function based on the solution of the semiclassical Liouville equation, of a two-site chemical exchange model of biological relevance. The bound site is allowed to be in the slow region regime that is the inverse quadrupole interaction of one bound site is in the same range as the reorientational correlation time. We compare different chemical exchange models, and several different physical situations are investigated. The variation of the width at half height (WHH) and the relative intensity (l/lo) is shown to be important, experimentally accessible quantities that are useful in order to discriminate between different model systems.     

A note on a possible mathematical approach to the theory of individual freedom

As suggested in previous publications, freedom may be defined quantitatively as a restriction upon the choice of a number of activities. If the choice is determined by maximizing the satisfaction function, it is suggested that freedom may be defined in terms of the satisfaction function. If an individual is isolated and no physical restrictions limit his choice of activities, he is free to choose any activity in an amount which maximizes his satisfaction. This isolated state may be considered therefore as that of maximum freedom. If the individual interacts with another, he will choose different amounts of his object of satisfaction depending on whether he behaves egoistically or altruistically. But in either case the value chosen will not maximize his satisfaction function considered alone. A simple analytical expression is suggested as a measure of freedom in this case, and some problems which arise from this suggestion are mentioned.     

Modeling isotopomer distributions in biochemical networks using isotopomer mapping matrices

Within the last decades NMR spectroscopy has undergone tremendous development and has become a powerful analytical tool for the investigation of intracellular flux distributions in biochemical networks using (13)C-labeled substrates. Not only are the experiments much easier to conduct than experiments employing radioactive tracer elements, but NMR spectroscopy also provides additional information on the labeling pattern of the metabolites. Whereas the maximum amount of information obtainable with (14)C-labeled substrates is the fractional enrichment in the individual carbon atom positions, NMR spectroscopy can also provide information on the degree of labeling at neighboring carbon atom positions by analyzing multiplet patterns in NMR spectra or using 2-dimensional NMR spectra. It is possible to quantify the mole fractions of molecules that show a specific labeling pattern, i.e., information of the isotopomer distribution in metabolite pools can be obtained. The isotopomer distribution is the maximum amount of information that in theory can be obtained from (13)C-tracer studies. The wealth of information contained in NMR spectra frequently leads to overdetermined algebraic systems. Consequently, fluxes must be estimated by nonlinear least squares analysis, in which experimental labeling data is compared with simulated steady state isotopomer distributions. Hence, mathematical models are required to compute the steady state isotopomer distribution as a function of a given set of steady state fluxes. Because 2(n) possible labeling patterns exist in a molecule of n carbon atoms, and each pattern corresponds to a separate state in the isotopomer model, these models are inherently complex. Model complexity, so far, has restricted usage of isotopomer information to relatively small metabolic networks. A general methodology for the formulation of isotopomer models is described. The model complexity of isotopomer models is reduced to that of classical metabolic models by expressing the 2(n) isotopomer mass balances of a metabolite pool in a single matrix equation. Using this approach an isotopomer model has been implemented that describes label distribution in primary carbon metabolism, i.e., in a metabolic network including the Embden-Meyerhof-Parnas and pentose phosphate pathway, the tricarboxylic acid cycle, and selected anaplerotic reaction sequences. The model calculates the steady state label distribution in all metabolite pools as a function of the steady state fluxes and is applied to demonstrate the effect of selected anaplerotic fluxes on the labeling pattern of the pathway intermediates. (c) 1997 John Wiley & Sons, Inc. Biotechnol Bioeng 55:831-840, 1997.     

Chemosensory responses of microorganisms: Asymmetric diffusion as the limit of a biased random walk

A macroscopic asymmetric diffusion equation to model the responses of microbial populations to chemical attractants and repellents is derived from a biased random walk model for the motion of individual cells. The equation is different from the well-known Keller-Segel equation which contains a Fickian diffusion term. The implications of this difference for selected problems of biological interest are considered. In particular the aggregation of a population of microorganisms in a region of high concentration of attractant is discussed. Some similarities and limitations of both models are noted.     

A Partition Function Approximation Using Elementary Symmetric Functions

In statistical mechanics, the canonical partition function can be used to compute equilibrium properties of a physical system. Calculating however, is in general computationally intractable, since the computation scales exponentially with the number of particles in the system. A commonly used method for approximating equilibrium properties, is the Monte Carlo (MC) method. For some problems the MC method converges slowly, requiring a very large number of MC steps. For such problems the computational cost of the Monte Carlo method can be prohibitive. Presented here is a deterministic algorithm – the direct interaction algorithm (DIA) – for approximating the canonical partition function in operations. The DIA approximates the partition function as a combinatorial sum of products known as elementary symmetric functions (ESFs), which can be computed in operations. The DIA was used to compute equilibrium properties for the isotropic 2D Ising model, and the accuracy of the DIA was compared to that of the basic Metropolis Monte Carlo method. Our results show that the DIA may be a practical alternative for some problems where the Monte Carlo method converge slowly, and computational speed is a critical constraint, such as for very large systems or web-based applications.     

Thermodynamics of information transfer between subunits in oligomeric enzymes and kinetic cooperativity. 1. Thermodynamics of subunit interactions, partition functions and enzyme reaction rate

The strength of quaternary constraints between two subunits of a polymeric enzyme depends upon the number of neighboring subunits and upon whether these subunits are liganded or not. These quaternary constraints between two subunits of a complex polymeric enzyme may be expressed, however, in terms of quaternary constraints that exist within ideal dimers. The influence of quaternary constraints on the reaction rate of a complex polymeric enzyme may thus be expressed in terms of the intersubunit strain that exists within dimers. This conclusion, that was far from evident, appears to be the consequence of the postulates of structural kinetics, and derive as well from usual thermodynamic principles. The structural steady-state equations may be expressed in terms of partition and sub-partition functions. As applied to structural kinetic models, a partition function expresses how, during the steady state, the energy of a population of enzyme molecules is distributed over n states. Similarly a sub-partition function describes how, during the steady state, the energy of these enzyme molecules is partitioned among only n-k of these states. Although the concept of partition function was initially formulated for equilibrium processes, it may be extended without any loss of generality to non-equilibrium processes. Moreover it is reminiscent of the concept of binding polynomial presented some years ago by Wyman for the equilibrium binding of a ligand to a protein. With this formalism derived from statistical mechanics, a structural rate equation may be derived from the ratio of a sub-partition function of degree n-1 and of a partition function of degree n. Again these properties are the consequence of the postulates of structural kinetics associated with simple ideas derived from statistical thermodynamics.     

Studies of the unfolding of an unstable mutant of staphylococcal nuclease: Evidence for low temperature unfolding and compactness of the high temperature unfolded state

Fluorescence and circular dichroism data as a function of temperature were obtained to characterize the unfolding of nuclease A and two of its less stable mutants. These spectroscopic data were obtained with a modified instrument that enables the nearly simultaneous detection of both fluorescence and CD data on the same sample. A global analysis of these multiple datasets yielded an excellent fit of a model that includes a change in the heat capacity change, ΔC_p, between the unfolded and native states. This analysis gives a ΔC_p of 2.2 kcal/mol/·K for thermal unfolding of the WT protein and 1.3 and 1.8 kcal/mol/K for the two mutants. These ΔC_p values are consistent with significant population of the cold unfolded state at ∼0°C. Independent evidence for the existence of a cold unfolded state is the observation of a separately migrating peak in size exclusion chromatography. The new chromatographic peak is seen near 0°C, has a partition coefficient corresponding to a larger hydrodynamic radius, and shows a red-shifted fluorescence spectrum, as compared to the native protein. Data also indicate that the high-temperature unfolded form of mutant nuclease is relatively compact. Size exclusion chromatography shows the high temperature unfolded form to have a hydrodynamic radius that is larger than that for the native form, but smaller than that for the urea or pH-induced unfolded forms. Addition of chemical denaturants to the high-temperature unfolded form causes a further unfolding of the protein, as indicated by an increase in the apparent hydrodynamic radius and a decrease in the rotational correlation time for Trp140 (as determined by fluorescence anisotropy decay measurements). Proteins 28:227–240, 1997 © 1997 Wiley-Liss Inc.     

Modeling quantitative trait Loci and interpretation of models

A quantitative genetic model relates the genotypic value of an individual to the alleles at the loci that contribute to the variation in a population in terms of additive, dominance, and epistatic effects. This partition of genetic effects is related to the partition of genetic variance. A number of models have been proposed to describe this relationship: some are based on the orthogonal partition of genetic variance in an equilibrium population. We compare a few representative models and discuss their utility and potential problems for analyzing quantitative trait loci (QTL) in a segregating population. An orthogonal model implies that estimates of the genetic effects are consistent in a full or reduced model in an equilibrium population and are directly related to the partition of the genetic variance in the population. Linkage disequilibrium does not affect the estimation of genetic effects in a full model, but would in a reduced model. Certainly linkage disequilibrium would complicate the detection of QTL and epistasis. Using different models does not influence the detection of QTL and epistasis. However, it does influence the estimation and interpretation of genetic effects.     

Classification and stability of global inhomogeneous solutions of a macroscopic model of cell motion

Many micro-organisms use chemotaxis for aggregation, resulting in stable patterns. In this paper, the amoeba Dictyostelium discoideum serves as a model organism for understanding the conditions for aggregation and classification of resulting patterns. To accomplish this, a 1D nonlinear diffusion equation with chemotaxis that models amoeba behavior is analyzed. A classification of the steady state solutions is presented, and a Lyapunov functional is used to determine conditions for stability of inhomogenous solutions. Changing the chemical sensitivity, production rate of the chemical attractant, or domain length can cause the system to transition from having an asymptotic steady state, to having asymptotically stable single-step solution and multi-stepped stable plateau solutions.     

Electronegativity—a perspective

Electronegativity is a very useful concept but it is not a physical observable; it cannot be determined experimentally. Most practicing chemists view it as the electron-attracting power of an atom in a molecule. Various formulations of electronegativity have been proposed on this basis, and predictions made using different formulations generally agree reasonably well with each other and with chemical experience. A quite different approach, loosely linked to density functional theory, is based on a ground-state free atom or molecule, and equates electronegativity to the negative of an electronic chemical potential. A problem that is encountered with this approach is the differentiation of a noncontinuous function. We show that this approach leads to some results that are not chemically valid. A formulation of atomic electronegativity that does prove to be effective is to express it as the average local ionization energy on an outer contour of the atom’s electronic density.     

Formation of unique structure in polypeptide chains. Theoretical investigation with the aid of a replica approach

A replica approach analogous to that used in spin glass systems is implemented to study the configurational space of a heteropolymeric model of protein with a quenched, disordered sequence of links in the limit of a large number of link types. It is shown that there exists a threshold value of chain heterogeneity which separates two qualitatively different types of behavior. For a low degree of heterogeneity the protein globule is like a homopolymer in a collapsed state without definite chain folds: an exponentially large number of folds make a significant contribution to the partition function in this regime. After the threshold heterogeneity has been overcome, the chain freezes drastically but without latent heat; few (approx. 1) frozen states with definite chain folds are thermodynamically dominant in this state. The relation of these results to thermodynamic aspects of protein folding is discussed.     

Matrix notation for efficient development of first-principles models within PAT applications: integrated modeling of antibiotic production with Streptomyces coelicolor

A matrix notation coupled to macroscopic principles is introduced as a means to develop first- principles models in an efficient and structured way within PAT applications. The notation was evaluated for developing an integrated biological, chemical (pH modeling) and physical (gas-liquid exchange) model for describing antibiotic production with Streptomyces coelicolor in batch fermentations. The model provided statistically adequate fits to all the monitored macroscopic biological, chemical and physical data of the process, except the phosphate uptake dynamics. This phosphate discrepancy is hypothesized to result from the internal storage of phosphate as polyphosphate prior to the exponential growth phase. The antibiotic production was associated with the stationary phase and its kinetics was adequately described using a modified Luedeking-Piret equation. Further, the maintenance was best described by employing a combination of Pirt and Herbert models, a result that was supported by a model-based hypothesis testing. Overall the process knowledge currently incorporated in the model is believed to be useful both for process optimization purposes and for further testing of hypotheses aiming at improving the mechanistic understanding of antibiotic production with S. coelicolor. Last but not least, the matrix notation is believed to be a promising supporting tool for efficient development and communication of complex dynamic models within a PAT framework.     

Physical invariant strain energy function for passive myocardium

Principal axis formulations are regularly used in isotropic elasticity, but they are not often used in dealing with anisotropic problems. In this paper, based on a principal axis technique, we develop a physical invariant orthotropic constitutive equation for incompressible solids, where it contains only a one variable (general) function. The corresponding strain energy function depends on six invariants that have immediate physical interpretation. These invariants are useful in facilitating an experiment to obtain a specific constitutive equation for a particular type of materials. The explicit appearance of the classical ground-state constants in the constitutive equation simplifies the calculation for their admissible values. A specific constitutive model is proposed for passive myocardium, and the model fits reasonably well with existing simple shear and biaxial experimental data. It is also able to predict a set of data from a simple shear experiment.     

Protein folding failure sets high-temperature limit on growth of phage P22 in Salmonella enterica serovar Typhimurium

The high-temperature limit for growth of microorganisms differs greatly depending on their species and habitat. The importance of an organism's ability to manage thermal stress is reflected in the ubiquitous distribution of the heat shock chaperones. Although many chaperones function to reduce protein folding defects, it has been difficult to identify the specific protein folding pathways that set the high-temperature limit of growth for a given microorganism. We have investigated this for a simple system, phage P22 infection of Salmonella enterica serovar Typhimurium. Production of infectious particles exhibited a broad maximum of 150 phage per cell when host cells were grown at between 30 and 39 degrees C in minimal medium. Production of infectious phage declined sharply in the range of 40 to 41 degrees C, and at 42 degrees C, production had fallen to less than 1% of the maximum rate. The host cells maintained optimal division rates at these temperatures. The decrease in phage infectivity was steeper than the loss of physical particles, suggesting that noninfectious particles were formed at higher temperatures. Sodium dodecyl sulfate-polyacrylamide gel electrophoresis revealed a decrease in the tailspike adhesins assembled on phage particles purified from cultures incubated at higher temperatures. The infectivity of these particles was restored by in vitro incubation with soluble tailspike trimers. Examination of tailspike folding and assembly in lysates of phage-infected cells confirmed that the fraction of polypeptide chains able to reach the native state in vivo decreased with increasing temperature, indicating a thermal folding defect rather than a particle assembly defect. Thus, we believe that the folding pathway of the tailspike adhesin sets the high-temperature limit for P22 formation in Salmonella serovar Typhimurium.     

A Quasistationary Analysis of a Stochastic Chemical Reaction: Keizer’s Paradox

For a system of biochemical reactions, it is known from the work of T.G. Kurtz [J. Appl. Prob. 8, 344 (1971)] that the chemical master equation model based on a stochastic formulation approaches the deterministic model based on the Law of Mass Action in the infinite system-size limit in finite time. The two models, however, often show distinctly different steady-state behavior. To further investigate this “paradox,” a comparative study of the deterministic and stochastic models of a simple autocatalytic biochemical reaction, taken from a text by the late J. Keizer, is carried out. We compute the expected time to extinction, the true stochastic steady state, and a quasistationary probability distribution in the stochastic model. We show that the stochastic model predicts the deterministic behavior on a reasonable time scale, which can be consistently obtained from both models. The transition time to the extinction, however, grows exponentially with the system size. Mathematically, we identify that exchanging the limits of infinite system size and infinite time is problematic. The appropriate system size that can be considered sufficiently large, an important parameter in numerical computation, is also discussed.     

The Gibbs-Duhem equation in membrane transport

It is argued that the Gibbs-Duhem equation alone cannot be used for deriving conclusions about the pressure gradient in membrane permeation. Statements regarding spatial variation of pressure in conjunction with chemical potential gradients of the components can legitimately be drawn from an equation that results from a combination of the G-D equation and the mechanical equilibrium equation. The derived equation has been applied here for explaining the mechanics of osmosis. In a further application, the frictional model has been improved here because the driving force also includes the membrane-solute potential interaction, thus allowing the solute partition coefficient to appear in the calculations naturally. By recognizing that because of the membrane-solution interaction, external forces of both potential and frictional character are present, the dissipation function is shown to depend explicitly on the centre-of-mass velocity. Thus the reference velocity for diffusive fluxes cannot be chosen arbitrarily, making Prigogine's theorem invalid in this approach to describing membrane permeation.