A system architecture for parallel analysis of flux-balanced metabolic pathways

Elementary flux mode (EFM) analysis is a well-studied method in constraint-based modeling of metabolic networks. In EFM analysis, a network is decomposed into minimal functional pathways based on the assumption of balanced metabolic fluxes. In this paper, a system architecture is proposed that approximately models the functionality of metabolic networks. The AND/OR graph model is used to represent the metabolic network and each processing element in the system emulates the functionality of a metabolite. The system is implemented on a graphics processing unit (GPU) as the hardware platform using CUDA environment. The proposed architecture takes advantage of the inherent parallelism in the network structure in terms of both pathway and metabolite traversal. The function of each element is defined such that it can find flux-balanced pathways. Pathways in both small and large metabolic networks are applied to the proposed architecture and the results are discussed.     

Polypeptide folding on a conformational‐space network: Dependence of network topology on the structural discretization procedure

Mapping the conformational space of a polypeptide onto a network of conformational states involves a number of subjective choices, mostly in relation to the definition of conformation and its discrete nature in a network framework. Here, we evaluate the robustness of the topology of conformational‐space networks derived from Molecular Dynamics (MD) simulations with respect to the use of different discretization (clustering) methods, variation of their parameters, simulation length and analysis time‐step, and removing high‐frequency motions from the coordinate trajectories. In addition, we investigate the extent to which polypeptide dynamics can be reproduced on the resulting networks when assuming Markovian behavior. The analysis is based on eight 500 ns and eight 400 ns MD simulations in explicit water of two 10‐residue peptides. Three clustering algorithms were used, two of them based on the pair‐wise root‐mean‐square difference between structures and one on dihedral‐angle patterns. A short characteristic path length and a power‐law behavior of the probability distribution of the node degree are obtained irrespective of the clustering method or the value of any of the tested parameters. The average cliquishness is consistently one or two orders of magnitude larger than that of a random realization of a network of corresponding size and connectivity. The cliquishness as function of node degree and the kinetic properties of the networks are found to be most dependent on clustering method and/or parameters. Although Markovian simulations on the networks reproduce cluster populations accurately, their kinetic properties most often differ from those observed in the MD simulations. © 2010 Wiley Periodicals, Inc. J Comput Chem, 2010     

A Systematic Petri Net Approach for Multiple-Scale Modeling and Simulation of Biochemical Processes

A method to exploit hybrid Petri nets for modeling and simulating biochemical processes in a systematic way was introduced. Both molecular biology and biochemical engineering aspects are manipulated. With discrete and continuous elements, the hybrid Petri nets can easily handle biochemical factors such as metabolites concentration and kinetic behaviors. It is possible to translate both molecular biological behavior and biochemical processes workflow into hybrid Petri nets in a natural manner. As an example, penicillin production bioprocess is modeled to illustrate the concepts of the methodology. Results of the dynamic of production parameters in the bioprocess were simulated and observed diagrammatically. Current problems and post-genomic perspectives were also discussed.     

High Accuracy Protein Structures from Minimal Sparse Paramagnetic Solid‐State NMR Restraints

There is a pressing need for new computational tools to integrate data from diverse experimental approaches in structural biology. We present a strategy that combines sparse paramagnetic solid‐state NMR restraints with physics‐based atomistic simulations. Our approach explicitly accounts for uncertainty in the interpretation of experimental data through the use of a semi‐quantitative mapping between the data and the restraint energy that is calibrated by extensive simulations. We apply our approach to solid‐state NMR data for the model protein GB1 labeled with Cu²⁺‐EDTA at six different sites. We are able to determine the structure to 0.9 Å accuracy within a single day of computation on a GPU cluster. We further show that in some cases, the data from only a single paramagnetic tag are sufficient for accurate folding.     

A constant-time kinetic Monte Carlo algorithm for simulation of large biochemical reaction networks

The time evolution of species concentrations in biochemical reaction networks is often modeled using the stochastic simulation algorithm (SSA) [Gillespie, J. Phys. Chem. 81, 2340 (1977)]. The computational cost of the original SSA scaled linearly with the number of reactions in the network. Gibson and Bruck developed a logarithmic scaling version of the SSA which uses a priority queue or binary tree for more efficient reaction selection [Gibson and Bruck, J. Phys. Chem. A 104, 1876 (2000)]. More generally, this problem is one of dynamic discrete random variate generation which finds many uses in kinetic Monte Carlo and discrete event simulation. We present here a constant-time algorithm, whose cost is independent of the number of reactions, enabled by a slightly more complex underlying data structure. While applicable to kinetic Monte Carlo simulations in general, we describe the algorithm in the context of biochemical simulations and demonstrate its competitive performance on small- and medium-size networks, as well as its superior constant-time performance on very large networks, which are becoming necessary to represent the increasing complexity of biochemical data for pathways that mediate cell function.     

Kinetic Modeling of the Decomposition of Carbon Tetrachloride in Thermal Plasma

Decomposition of carbon tetrachloride in a RF thermal plasma reactor was investigated in argon atmosphere. The net conversion of CCl₄ and the main products of its decomposition were determined from the mass spectrometric analysis of outlet gases. Flow and temperature profiles in the reactor were calculated and concentration profiles of the species along the axis of the reactor were estimated using a newly developed chemical kinetic mechanism, containing 12 species and 34 reaction steps. The simulations indicated that all carbon tetrachloride decomposed within a few microseconds. However, CCl₄ was partly recombined from its decomposition products. The calculations predicted 70\% net conversion of CCl₄, which was close to the experimentally determined value of 60\%. A thermodynamic equilibrium model also simulated the decomposition. Results of the kinetic and thermodynamic simulations agreed well above 2000 K. However, below 2000 K the thermodynamic equilibrium model gave wrong predictions. Therefore, application of detailed kinetic mechanisms is recommended for modeling CCl₄ decomposition under thermal plasma conditions.     

Comparing normal modes across different models and scales: Hessian reduction versus coarse‐graining

Dimension reduction is often necessary when attempting to reach longer length and time scales in molecular simulations. It is realized by constraining degrees of freedom or by coarse‐graining the system. When evaluating the accuracy of a dimensional reduction, there is a practical challenge: the models yield vectors with different lengths, making a comparison by calculating their dot product impossible. This article investigates mapping procedures for normal mode analysis. We first review a horizontal mapping procedure for the reduced Hessian techniques, which projects out degrees of freedom. We then design a vertical mapping procedure for the “implosion” of the all‐atom (AA) Hessian to a coarse‐grained scale that is based upon vibrational subsystem analysis. This latter method derives both effective force constants and an effective kinetic tensor. Next, a series of metrics is presented for comparison across different scales, where special attention is given to proper mass‐weighting. The dimension‐dependent metrics, which require prior mapping for proper evaluation, are frequencies, overlap of normal mode vectors, probability similarity, Hessian similarity, collectivity of modes, and thermal fluctuations. The dimension‐independent metrics are shape derivatives, elastic modulus, vibrational free energy differences, heat capacity, and projection on a predefined basis set. The power of these metrics to distinguish between reasonable and unreasonable models is tested on a toy alpha helix system and a globular protein; both are represented at several scales: the AA scale, a Gō‐like model, a canonical elastic network model, and a network model with intentionally unphysical force constants. Published 2012 Wiley Periodicals, Inc.     

Efficient gaussian density formulation of volume and surface areas of macromolecules on graphical processing units

We present an algorithm to efficiently compute accurate volumes and surface areas of macromolecules on graphical processing unit (GPU) devices using an analytic model which represents atomic volumes by continuous Gaussian densities. The volume of the molecule is expressed by means of the inclusion–exclusion formula, which is based on the summation of overlap integrals among multiple atomic densities. The surface area of the molecule is obtained by differentiation of the molecular volume with respect to atomic radii. The many‐body nature of the model makes a port to GPU devices challenging. To our knowledge, this is the first reported full implementation of this model on GPU hardware. To accomplish this, we have used recursive strategies to construct the tree of overlaps and to accumulate volumes and their gradients on the tree data structures so as to minimize memory contention. The algorithm is used in the formulation of a surface area‐based non‐polar implicit solvent model implemented as an open source plug‐in (named GaussVol) for the popular OpenMM library for molecular mechanics modeling. GaussVol is 50 to 100 times faster than our best optimized implementation for the CPUs, achieving speeds in excess of 100 ns/day with 1 fs time‐step for protein‐sized systems on commodity GPUs. © 2017 Wiley Periodicals, Inc.     

Molecular dynamics simulations of aqueous ions at the liquid–vapor interface accelerated using graphics processors

Molecular dynamics (MD) simulations are a vital tool in chemical research, as they are able to provide an atomistic view of chemical systems and processes that is not obtainable through experiment. However, large‐scale MD simulations require access to multicore clusters or supercomputers that are not always available to all researchers. Recently, scientists have returned to exploring the power of graphics processing units (GPUs) for various applications, such as MD, enabled by the recent advances in hardware and integrated programming interfaces such as NVIDIA's CUDA platform. One area of particular interest within the context of chemical applications is that of aqueous interfaces, the salt solutions of which have found application as model systems for studying atmospheric process as well as physical behaviors such as the Hoffmeister effect. Here, we present results of GPU‐accelerated simulations of the liquid–vapor interface of aqueous sodium iodide solutions. Analysis of various properties, such as density and surface tension, demonstrates that our model is consistent with previous studies of similar systems. In particular, we find that the current combination of water and ion force fields coupled with the ability to simulate surfaces of differing area enabled by GPU hardware is able to reproduce the experimental trend of increasing salt solution surface tension relative to pure water. In terms of performance, our GPU implementation performs equivalent to CHARMM running on 21 CPUs. Finally, we address possible issues with the accuracy of MD simulaions caused by nonstandard single‐precision arithmetic implemented on current GPUs. © 2010 Wiley Periodicals, Inc. J Comput Chem, 2011     

Development of Petri Net-Based Dynamic Model for Improved Production of Farnesyl Pyrophosphate by Integrating Mevalonate and Methylerythritol Phosphate Pathways in Yeast

In this case study, we designed a farnesyl pyrophosphate (FPP) biosynthetic network using hybrid functional Petri net with extension (HFPNe) which is derived from traditional Petri net theory and allows easy modeling with graphical approach of various types of entities in the networks together. Our main objective is to improve the production of FPP in yeast, which is further converted to amorphadiene (AD), a precursor of artemisinin (antimalarial drug). Natively, mevalonate (MEV) pathway is present in yeast. Methyl erythritol phosphate pathways (MEP) are present only in higher plant plastids and eubacteria, but not present in yeast. IPP and DAMP are common isomeric intermediate in these two pathways, which immediately yields FPP. By integrating these two pathways in yeast, we augmented the FPP synthesis approximately two folds higher (431.16 U/pt) than in MEV pathway alone (259.91 U/pt) by using HFPNe technique. Further enhanced FPP levels converted to AD by amorphadiene synthase gene yielding 436.5 U/pt of AD which approximately two folds higher compared to the AD (258.5 U/pt) synthesized by MEV pathway exclusively. Simulation and validation processes performed using these models are reliable with identified biological information and data.     

Effects of Tao‐Hong‐Si‐Wu decoction on acute blood stasis in rats based on a LC‐Q/TOF‐MS metabolomics and network approach

A novel approach using metabolomics coupled with a metabolic network was used to investigate the effects of Tao‐Hong‐Si‐Wu decoction (THSWD) on the rat model of acute blood stasis syndrome. Acute blood stasis syndrome was induced by placing the rats in ice‐cold water following two injections with epinephrine. The hemorheological indicators [whole blood viscosity (WBV) and plasma viscosity (PV)] and the blood coagulation indicators [thrombin time (TT), prothrombin time (PT), activated partial thromboplastin time (APTT) and fibrinogen (FIB)] were detected. The nonparametric univariate method and multivariate statistical analysis were performed for determining the potential biomarkers. A correlation map was structured between biochemical indicators and hub metabolites to explain the effects mechanism of THSWD. After the administration of THSWD, the levels of WBV, PV, TT, APTT and FIB returned to levels observed in the control group. According to metabolomics coupled with metabolic network analysis, the intervention of THSWD in rats with acute blood stasis syndrome induced substantial and characteristic changes in their metabolic profiles. Fifteen metabolites were screened, which mainly involved 10 pathways and five hub metabolites, namely, l ‐glutamate, l ‐phenylalanine, N‐acylsphingosine, arachidonic acid and phosphatidate. The biochemical indicators and hub metabolites could be adjusted to close to normal levels by THSWD. Therefore, combining metabolomics and metabolic network helped to evaluate the effects of THSWD on acute blood stasis.     

Rate‐Reactivity Model: A New Theoretical Basis for Systematic Kinetic Characterization of Heterogeneous Catalysts

Non‐steady‐state kinetic measurements contain a wealth of information about catalytic reactions and other gas–solid chemical interactions, which is extracted from experimental data via kinetic models. The standard mathematical framework of microkinetic models, which are typically used in computational catalysis and for advanced modeling of steady‐state data, encounters multiple challenges when applied to non‐steady‐state data. Robust phenomenological models, such as the steady‐state Langmuir–Hinshelwood–Hougen–Watson equations, are presently unavailable for non‐steady‐state data. Herein, a novel modeling framework is proposed to fulfill this need. The rate‐reactivity model (RRM) is formulated in terms of experimentally observable quantities including the gaseous transformation rates, concentrations, and surface uptakes. The model is linear with respect to these quantities and their pairwise products, and it is also linear in terms of its parameters (reactivities). The RRM parameters have a clear physicochemical meaning and fully characterize the kinetic behavior of a specific catalyst state, but unlike microkinetic models that rely on hypothetical surface intermediates and specific reaction networks, the RRM does not require any assumptions regarding the underlying mechanism. The systematic RRM‐based procedure outlined in this paper enables an effective comparison of various catalysts and the construction of more detailed microkinetic models in a rational manner. The model was applied to temporal analysis of products pulse‐response data as an example, but it is more generally applicable to other non‐steady‐state techniques that provide time‐resolved rates and concentrations. Several numerical examples are given to illustrate the application of the model to simple model reactions.     

Using gene expression programming to infer gene regulatory networks from time-series data

Gene regulatory networks inference is currently a topic under heavy research in the systems biology field. In this paper, gene regulatory networks are inferred via evolutionary model based on time-series microarray data. A non-linear differential equation model is adopted. Gene expression programming (GEP) is applied to identify the structure of the model and least mean square (LMS) is used to optimize the parameters in ordinary differential equations (ODEs). The proposed work has been first verified by synthetic data with noise-free and noisy time-series data, respectively, and then its effectiveness is confirmed by three real time-series expression datasets. Finally, a gene regulatory network was constructed with 12 Yeast genes. Experimental results demonstrate that our model can improve the prediction accuracy of microarray time-series data effectively.     

Probing the molecular interactions in the diffusion of water through epoxy and epoxy–bismaleimide networks

Fourier transform infrared spectroscopy in the near‐infrared (NIR) frequency range was used to investigate the molecular interactions occurring between absorbed water molecules and networks based on a tetrafunctional epoxy resin. One of these networks was a typical formulation containing 4,4′‐diamino diphenylsulfone as a hardener, and the other was a modified resin containing 4,4′‐bismaleimide‐diphenylmethane (BMI) as a coreactive monomer. Molecular spectroscopy analysis confirmed the existence of mobile water localized into network defects (microvoids) that did not interact with the networks and water molecules bound to the networks through hydrogen‐bonding interactions. In the BMI‐containing system, the fraction of bound water decreased significantly with respect to the unmodified epoxy resin. This was a relevant result because the bound water was primarily responsible for the plasticization of the network and for the consequent worsening of mechanical performance. Water diffusion was investigated with gravimetric sorption measurements and time‐resolved Fourier transform NIR spectroscopy measurements. These showed that the presence of BMI decreased the water uptake at equilibrium, enhanced the diffusivity, and reduced the activation energy for diffusion. A dual‐mode model for diffusion was found to be suitable for accurately describing the mass‐transport process in both investigated systems. The results of the model simulations allowed us to estimate the ratio of free and bound water, which was in good agreement with that obtained from the spectroscopic analysis. © 2002 Wiley Periodicals, Inc. J Polym Sci Part B: Polym Phys 40: 922–938, 2002     

The Viedma Deracemization of Racemic Conglomerate Mixtures as a Paradigm of Spontaneous Mirror Symmetry Breaking in Aggregation and Polymerization

Simulations of a chemical kinetics model, based on the free‐energy relationships of classical primary nucleation theory, show that the deracemization phenomenon in systems of achiral or fast racemizing compounds yielding enantiopure crystals as the more stable solid phase is a true spontaneous mirror symmetry breaking process (SMSB). That is, the achievement of a stationary chiral state is more stable than the racemic one. The model translates the free‐energy relationships determined by the existence of a critical size cluster to a chemical kinetics model, in which the consideration of forward and backward reaction rate constants avoids the misuse of network parameters that violate thermodynamic constraints (microreversibility principle), which would lead to apparent in silico SMSB. The model provides qualitative agreement for deracemizations by mechanical attrition of visible crystals, as well as for those obtained under temperature gradients. The analysis of the effect of the system parameters to obtain a SMSB scenario shows that the network possesses the principal characteristics of SMSB networks: 1) an enantioselective autocatalytic stage, corresponding to the non‐linear kinetics of enantioselective (homochiral) cluster‐to‐cluster growth, and 2) the mutual inhibition step originating in the backward flow of chiral clusters towards smaller achiral clusters, or even to a racemizing monomer. The application of such a SMSB kinetic model to enantioselective polymerizations and to chiral biopolymers is discussed.     

Understanding the elastomeric properties of polymer networks

It is shown that Monte‐Carlo (MC) simulations of the elastic behaviour of chains in networks using realistic rotational‐isomeric‐state (RIS) chain models are able to reproduce experimentally observed deviations from Gaussian network behaviour in uniaxial extension and also to interpret, quantitatively, stress‐optical properties. In stress‐strain behaviour, an increase in the proportion of fully extended chains with increasing macroscopic strain gives rise to a steady decrease in the rate of change of the Helmholtz energy of a network, causing a reduction in network modulus at moderate macroscopic strains. There is no need to invoke a transition from affine to phantom chain behaviour as deformation increases. To evaluate stress‐optical properties, the average orientation of segments with respect to the deformation axis is calculated, taking into account the interdependence of segment orientation and chain orientation as chains become more extended and aligned under uniaxial stress. The MC method gives, in agreement with experiment, values of stress‐optical coefficient that are dependent upon both deformation ratio and network‐chain length. The method highlights serious shortcomings in the classical Gaussian model of stress‐optical behaviour. Applications of the simulation methods to the quantitative modelling of the stress‐strain behaviour of poly(dimethyl siloxane) networks and the stress‐optical behaviour of polyethylene networks are described.     

The ORAC Assay: Mathematical Analysis of the Rate Equations and Some Practical Considerations

ORAC (oxygen radical absorbance capacity), a method widely used for measuring the total antioxidant capacity of biological samples, can also be used for the determination of the relative reactivity of an antioxidant compound (XH) by examining the dependence of the rate of consumption of the probe (PH) on the concentration of XH; initial conditions are chosen in such a way that the rate of consumption of the starting reactants may be assumed to follow a drastically simplified kinetic scheme, and the steady‐state approximation for the concentration of the azo compound peroxyl (ROO^•) radical is invoked to simplify the analysis. Here we first attempted to find an analytical solution to the coupled first‐order ordinary differential equations (ODEs) of the minimal ORAC kinetic system, applying Lie symmetry group theory without any precondition. However, the Lie symmetry transformations applied to the Chini equation, which appeared after mathematical transformations, showed that the form of the coefficients of the Chini equation precluded the analytical solution of the minimal ORAC kinetic system through symmetry reduction. Consequently, an approximate analytical solution was sought, valid for the case when the bimolecular rate constant of XH with ROO^• (i.e., k_x ) was much larger than that of PH with ROO^• (i.e., k_p ). Using numerical solutions of the original set of ODEs of the ORAC kinetic system, the quality of the approximate solution was inspected under conditions that correspond to those employed in several ORAC methods together with a low initial concentration of the azo compound radical initiator. The simulations allowed us to conclude that the approximate analytical solution of the ODEs of the minimal ORAC kinetic system was not entirely devoid of academic interest, but its applicability was restricted to conditions where both k_x ≫ k_p and the initial concentration of XH was higher than that of PH.