Lattice microbes: High‐performance stochastic simulation method for the reaction‐diffusion master equation

Spatial stochastic simulation is a valuable technique for studying reactions in biological systems. With the availability of high‐performance computing (HPC), the method is poised to allow integration of data from structural, single‐molecule and biochemical studies into coherent computational models of cells. Here, we introduce the Lattice Microbes software package for simulating such cell models on HPC systems. The software performs either well‐stirred or spatially resolved stochastic simulations with approximated cytoplasmic crowding in a fast and efficient manner. Our new algorithm efficiently samples the reaction‐diffusion master equation using NVIDIA graphics processing units and is shown to be two orders of magnitude faster than exact sampling for large systems while maintaining an accuracy of ～ 0.1%. Display of cell models and animation of reaction trajectories involving millions of molecules is facilitated using a plug‐in to the popular VMD visualization platform. The Lattice Microbes software is open source and available for download at http://www.scs.illinois.edu/schulten/lm © 2012 Wiley Periodicals, Inc.     

Various methods for enhancement in heterogeneous catalysed biodiesel production reaction rate

The biodiesel can be produced by transesterification and esterification reaction with edible and non-edible oil respectively. These reactions are catalysed by both homogeneous and heterogenous catalyst. The transesterification reactions for heterogeneous catalysts proceed at a relatively slow rate. The heterogeneous reaction mixture constitutes a three-phase system, oil/alcohol/catalyst therefore the mass transfer limitation controls the reaction rate. In the present study Tetra-hydrofuran, Hexane and Heptane as co-solvent have been tested for the transesterification and esterification reaction. Tetra-hydrofuran, assists in decreasing the mass transfer between the oil and methanol phases. Tetra-hydrofuran was found to be the best co-solvent. It also represents that the presence of Tetra-hydrofuran only enhance the solubility of phases not final equilibrium yield. Results are compared with homogeneous and heterogeneous catalysed reaction. Ultrasonic irradiation influenced the mixing of the reaction mixture and resulted in higher yields, for each co-solvent case.     

Criteria for diffusion limitation in coordination polymerization

The following criteria are proposed to judge whether a coordination polymerization may be diffusion controlled or not: (1) If the number-average molecular weight and polydispersity of the polymer calculated from kinetic rate constants as a function of time agree with the experimental values, the polymerization is not diffusion controlled. (2) The polymerization may be diffusion controlled if the Thiele modulus, the ratio of the characteristic diffusion time to the characteristic reaction time, is much greater than unity; if it is much smaller than unity, the polymerization is reaction controlled. (3) If an initial linear dependence of rate of polymerization on catalyst concentration changes over to a square-root dependence, the polymerization may be diffusion limited. (4) The polymerization is likely to be diffusion limited if the instantaneous rate of polymerization is proportional to the rate of particle growth when the proportionality coefficient is the surface area of the particle. Criterion (1) is a necessary and sufficient condition as stated, as its converse is not true. All the other criteria are merely necessary but not sufficient conditions. The established Ziegler–Natta catalysts have activities too low to cause diffusion limitation; the Phillips catalyst system is likely to be diffusion limited. The polydispersity of polyolefins produced with Ziegler–Natta catalysts are not the consequence of diffusion control but are the characteristics of the catalysts in their kinetics of initiation, propagation, chain transfer, and termination.     

A comparative study of the master equation approach and the discrete galerkin method for the simulation of free radical polymerization

The present article deals with the mathematical treatment of free radical polymerization reactions. As a typical example the synthesis of poly(methyl methacrylate) under realistic experimental conditions is investigated. Since the mathematical treatment of the kinetic rate equations raises severe numerical problems, alternative approaches are required. In this paper two of these methods, i.e. the discrete Galerkin method and the master equation approach, are compared. The discrete Galerkin method circumvents difficulties encountered by the direct integration of the kinetic rate equations but requires much a priori knowledge of the chemical reaction system. Within the framework of the master equation approach the polymerization reaction is regarded as a stochastic process. For the simulation of this stochastic process a modified algorithm is presented. The example of the polymerization of methyl methacrylate shows that the master equation approach is an efficient tool in the simulation of free radical polymerization reactions.     

Hybrid modeling and simulation of stochastic effects on progression through the eukaryotic cell cycle

The eukaryotic cell cycle is regulated by a complicated chemical reaction network. Although many deterministic models have been proposed, stochastic models are desired to capture noise in the cell resulting from low numbers of critical species. However, converting a deterministic model into one that accurately captures stochastic effects can result in a complex model that is hard to build and expensive to simulate. In this paper, we first apply a hybrid (mixed deterministic and stochastic) simulation method to such a stochastic model. With proper partitioning of reactions between deterministic and stochastic simulation methods, the hybrid method generates the same primary characteristics and the same level of noise as Gillespie's stochastic simulation algorithm, but with better efficiency. By studying the results generated by various partitionings of reactions, we developed a new strategy for hybrid stochastic modeling of the cell cycle. The new approach is not limited to using mass-action rate laws. Numerical experiments demonstrate that our approach is consistent with characteristics of noisy cell cycle progression, and yields cell cycle statistics in accord with experimental observations.     

Efficient stochastic simulation of simultaneous reaction and diffusion in a gas-liquid interface

A stochastic simulation of simultaneous reaction and diffusion is proposed for the gas-liquid interface formed in the surface of a gas bubble within a liquid. The interface between a carbon dioxide bubble and an aqueous solution of calcium hydroxide was simulated as an application example, taken from the integrated production of calcium carbonate. First Gillespie’s stochastic simulation algorithm was applied in separate reaction and diffusion simulations. The results from these simulations were consistent with deterministic solutions based on differential equations. However it was observed that stochastic diffusion simulations are extremely slow. The sampling of diffusion events was accelerated applying a group molecule transfer scheme based on the binomial distribution function. Simulations of the reaction-diffusion in the gas-liquid interface based on the standard Gillespie’s stochastic algorithm were also slow. However the application of the binomial distribution function scheme allowed to compute the concentration profiles in the gas-liquid interface in a fraction of the time required with the standard Gillespie’s stochastic algorithm.     

Binomial tau-leap spatial stochastic simulation algorithm for applications in chemical kinetics

In cell biology, cell signaling pathway problems are often tackled with deterministic temporal models, well mixed stochastic simulators, and/or hybrid methods. But, in fact, three dimensional stochastic spatial modeling of reactions happening inside the cell is needed in order to fully understand these cell signaling pathways. This is because noise effects, low molecular concentrations, and spatial heterogeneity can all affect the cellular dynamics. However, there are ways in which important effects can be accounted without going to the extent of using highly resolved spatial simulators (such as single-particle software), hence reducing the overall computation time significantly. We present a new coarse grained modified version of the next subvolume method that allows the user to consider both diffusion and reaction events in relatively long simulation time spans as compared with the original method and other commonly used fully stochastic computational methods. Benchmarking of the simulation algorithm was performed through comparison with the next subvolume method and well mixed models (MATLAB), as well as stochastic particle reaction and transport simulations (CHEMCELL, Sandia National Laboratories). Additionally, we construct a model based on a set of chemical reactions in the epidermal growth factor receptor pathway. For this particular application and a bistable chemical system example, we analyze and outline the advantages of our presented binomial tau-leap spatial stochastic simulation algorithm, in terms of efficiency and accuracy, in scenarios of both molecular homogeneity and heterogeneity.     

Accurate hybrid stochastic simulation of a system of coupled chemical or biochemical reactions

The dynamical solution of a well-mixed, nonlinear stochastic chemical kinetic system, described by the Master equation, may be exactly computed using the stochastic simulation algorithm. However, because the computational cost scales with the number of reaction occurrences, systems with one or more "fast" reactions become costly to simulate. This paper describes a hybrid stochastic method that partitions the system into subsets of fast and slow reactions, approximates the fast reactions as a continuous Markov process, using a chemical Langevin equation, and accurately describes the slow dynamics using the integral form of the "Next Reaction" variant of the stochastic simulation algorithm. The key innovation of this method is its mechanism of efficiently monitoring the occurrences of slow, discrete events while simultaneously simulating the dynamics of a continuous, stochastic or deterministic process. In addition, by introducing an approximation in which multiple slow reactions may occur within a time step of the numerical integration of the chemical Langevin equation, the hybrid stochastic method performs much faster with only a marginal decrease in accuracy. Multiple examples, including a biological pulse generator and a large-scale system benchmark, are simulated using the exact and proposed hybrid methods as well as, for comparison, a previous hybrid stochastic method. Probability distributions of the solutions are compared and the weak errors of the first two moments are computed. In general, these hybrid methods may be applied to the simulation of the dynamics of a system described by stochastic differential, ordinary differential, and Master equations.     

Reactions on cell membranes: comparison of continuum theory and Brownian dynamics simulations

Biochemical transduction of signals received by living cells typically involves molecular interactions and enzyme-mediated reactions at the cell membrane, a problem that is analogous to reacting species on a catalyst surface or interface. We have developed an efficient Brownian dynamics algorithm that is especially suited for such systems and have compared the simulation results with various continuum theories through prediction of effective enzymatic rate constant values. We specifically consider reaction versus diffusion limitation, the effect of increasing enzyme density, and the spontaneous membrane association/dissociation of enzyme molecules. In all cases, we find the theory and simulations to be in quantitative agreement. This algorithm may be readily adapted for the stochastic simulation of more complex cell signaling systems.     

Analysis of mass transport limitation during emulsion polymerization using reaction calorimetry and conductivity measurement

The seeded emulsion polymerization of styrene was investigated by reaction calorimetry in combination with conductivity measurements using the RC1e calorimeter (Mettler Toledo). Varying the stirrer speed showed the existence of mass transport limitation in the lower speed range. For the case of using the autoclave HP60 with an anchor stirrer a minimum stirrer speed of 400 rpm is required to be sure of adequate mixing in the reactor.     

A distributed parameter diffusion-reaction model for the alcoholic fermentation process

A distributed parameter model is developed for the yeast floc in the alcoholic fermentation process. The model takes into consideration the external mass transfer resistances, the mass transfer resistance through the cellular membrane, and the diffusion resistances inside the floc. The two-point boundary value differential equations for the membrane are manipulated analytically, whereas the nonlinear twopoint boundary value differential equations of diffusion and reaction inside the floc have been approximated using the orthogonal collocation technique. The evaluation of the necessary diffusion coefficients have involved a relatively large number of assumptions because of the present limited knowledge regarding the complex process of diffusion and biochemical reactions in these systems.     

Modelling the free-radical polymerization of methyl methacrylate over the complete range of conversion

A model is proposed for the radical-induced polymerization of methyl methacrylate which combines concepts given in the literature. Consecutive steps for diffusion and reaction are used for both the propagation and the termination process, and the termination at higher degrees of conversion is dominated by reaction diffusion as described by Buback. The diffusional contributions to the processes are semi-empirically based on the free-volume theory, following a suggestion made by Marten and Hamielec, and the diffusion of the macroradicals is assumed to depend on the instantaneous molar mass (as representative for their chain length) and on the cumulative molar mass (representing the matrix of the dead polymer). With one fixed and four adjustable parameters, the model successfully not only describes the bulk polymerization but also takes account of the presence of solvent, chain transfer agent and pre-polymer.     

Calculating the True and Observed Rates of Complex Heterogeneous Catalytic Reactions

Equations of the theory of steady-state complex reactions are considered in matrix form. A set of stage stationarity equations is given, and an algorithm is described for deriving the canonic set of stationarity equations with appropriate corrections for the existence of fast stages in a mechanism. A formula for calculating the number of key compounds is presented. The applicability of the Gibbs rule to estimating the number of independent compounds in a complex reaction is analyzed. Some matrix equations relating the rates of dependent and key substances are derived. They are used as a basis to determine the general diffusion stoichiometry relationships between temperature, the concentrations of dependent reaction participants, and the concentrations of key reaction participants in a catalyst grain. An algorithm is described for calculating heat and mass transfer in a catalyst grain with respect to arbitrary complex heterogeneous catalytic reactions.     

Kinetics of Ge(IV) Extraction Using a Microstructured Membrane Contactor

For hydrometallurgical metal extraction, the mass transfer by diffusion is additionally coupled with a chemical reaction of the metal ion with the extractant molecule, usually at the interface. Consequently, the kinetics is either limited by diffusion or chemical reaction (or both, respectively, in a mixed regime). Conventional methods for determining the extraction kinetics often lead to a misinterpretation, especially for fast reactions, and an isolated view of the interfacial reaction is restricted. With a new microcontactor setup, it is possible to perform a comprehensive kinetic analysis with very low sample volumes compared to established methods. Additionally, it is possible to quantify all individual mass transfer resistances and identify the extraction regime with the developed mass transfer model. The chemical reaction part is investigated isolated, to derive rate laws and kinetic constants. The methodology is discussed for the extraction of Ge(IV) from aqueous solutions with the two extractants 5,8‐diethyl‐7‐hydroxydodecan‐6‐oxime (LIX 63) and 7‐(4‐ethyl‐1‐methyloctyl)‐8‐hydroxyquinoline (Kelex 100). It was found that the extraction with LIX 63 is reaction limited and with Kelex 100 the reaction resistance is in the same order of magnitude as the diffusion resistances. The obtained results provide fundamental kinetic data for the design of solvent extraction equipment.     

Numerical simulation of polymerization in interdigital multilamination micromixers

Free radical polymerization in microfluidic devices modeled with the help of numerical simulations is discussed. The simulation method used allows the simultaneous solvation of partial differential equations resulting from the hydrodynamics, thermal and mass transfer (convection, diffusion and chemical reaction). Three microfluidic devices are modeled, two interdigital multilamination micromixers respectively with a large and short focusing section, and a simple T-junction followed by a microtube reactor together considered as a bilamination micromixer with a large focusing section. The simulations show that in spite of the heat released by the polymerization reaction, the thermal transfer in such microfluidic devices is high enough to ensure isothermal conditions. Moreover, for low radial Peclet number, microfluidic devices with a large focusing section can achieve better control over the polymerization than a laboratory scale reactor as the polydispersity index obtained is very close to the theoretical limiting value. As the characteristic dimension of the microfluidic device increases, i.e. for high radial Peclet number, the reactive medium cannot be fully homogenized by the diffusion transport before leaving the system resulting in a high polydispersity index and a loss in the control of the polymerization.     

Preparation of Polymer Microparticles through Nonaqueous Suspension Polycondensations: Part V—Modeling and Parameter Estimation for Poly(butylene succinate) Polycondensations

Polycondensation polymers are normally produced through bulk and solution polymerization processes, which are characterized by significant mass and heat transfer constraints and difficult polymer purification (when prepared in solution). Therefore, it is desirable to develop industrial processes that can circumvent some of these limitations. Recently, a suspension polycondensation process has been developed, rendering the industrial process simpler and enabling the manufacture of polycondensation polymer microparticles. For this reason, the present work builds a phenomenological model to describe the analyzed suspension polycondensation reactions and estimate the model parameters required to simulate poly(butylene succinate) suspension polycondensations. It is shown that both the suspending medium and the reaction conditions can affect the mass transfer resistance for removal of water and that mass transfer rate coefficients are controlled mainly by reaction temperature and solubility of water in the suspending medium, leading to higher mass transfer rates when polymerizations are carried out in soybean oil (when compared to paraffin) at higher temperatures.     

HZSM-5分子筛上甲醇制烯烃典型产物的传质行为研究

分子筛催化甲醇制烯烃反应(MTO)是典型的扩散主导反应过程,运用频率响应技术系统研究了几种典型产物分子(乙烯/乙烷、丙烯/丙烷、苯)在HZSM-5分子筛上的扩散行为。结果表明,频率响应法成功辨析了不同产物分子的传质规律,证实C2和C3烃分子在HZSM-5微孔孔道内具有相近的扩散速率,但由于受晶体表面阻碍效应影响不同,乙烷分子可自由进出HZSM-5分子筛孔道,而丙烷分子则受到较显著的微孔孔道扩散限制。另外,苯分子的扩散速率显著小于C2和C3分子,且苯分子受晶体表面阻抗效应的影响较小。本研究结果可用于解释HZSM-5分子筛在MTO反应中产物选择性的特点及表面结焦原因,进而从传质角度为高活性、选择性以及稳定性的高效甲醇转化制烃催化剂的定向开发提供理论指导。     

Understanding supported reactions in spherical compartments: a general algorithm to model and determine rate constants, diffusion coefficients, and spatial product distributions

A general algorithm allowing the numerical modeling of the time and space dependence of product formation in spherical reaction volumes is described. The algorithm is described by the complete set of mass balance equations. On the basis of these equations, the effects of the diffusion coefficient, reaction rate, bead size, reagent excess, and packing density of the resin beads on the overall reaction rates are determined for second-order reactions. Experimental data of reaction progress are employed to calculate reaction rates and diffusion coefficients in polymer-supported reactions. In addition, the conditions for shell-like product formation are determined, and various strategies for the radial patterning of resin beads are compared. The effect of diffusion on polymer-supported enzyme-catalyzed reactions of the Michaelis-Menten type is treated, as well. Finally, the effects of typical nonideal solid-phase phenomena, namely, the inhomogeneity of rate constants and the concentration dependence of diffusion coefficients, on overall rates are discussed.     

Hybrid stochastic simulations of intracellular reaction–diffusion systems

With the observation that stochasticity is important in biological systems, chemical kinetics have begun to receive wider interest. While the use of Monte Carlo discrete event simulations most accurately capture the variability of molecular species, they become computationally costly for complex reaction–diffusion systems with large populations of molecules. On the other hand, continuous time models are computationally efficient but they fail to capture any variability in the molecular species. In this study a hybrid stochastic approach is introduced for simulating reaction–diffusion systems. We developed an adaptive partitioning strategy in which processes with high frequency are simulated with deterministic rate-based equations, and those with low frequency using the exact stochastic algorithm of Gillespie. Therefore the stochastic behavior of cellular pathways is preserved while being able to apply it to large populations of molecules. We describe our method and demonstrate its accuracy and efficiency compared with the Gillespie algorithm for two different systems. First, a model of intracellular viral kinetics with two steady states and second, a compartmental model of the postsynaptic spine head for studying the dynamics of Ca+2 and NMDA receptors.