Molecular alignment in a liquid induced by a nonresonant laser field: Molecular dynamics simulation

We carried out molecular dynamics (MD) simulations for a dilute aqueous solution of pyrimidine in order to investigate the mechanisms of field-induced molecular alignment in a liquid phase. An anisotopically polarizable molecule can be aligned in a liquid phase by the interaction with a nonresonant intense laser field. We derived the effective forces induced by a nonresonant field on the basis of the concept of the average of the total potential over one optical cycle. The results of MD simulations show that a pyrimidine molecule is aligned in an aqueous solution by a linearly polarized field of light intensity I approximately 10(13) W/cm2 and wavelength lambda = 800 nm. The temporal behavior of field-induced alignment is adequately reproduced by the solution of the Fokker-Planck equation for a model system in which environmental fluctuations are represented by Gaussian white noise. From this analysis, we have revealed that the time required for alignment in a liquid phase is in the order of the reciprocals of rotational diffusion coefficients of a solute molecule. The degree of alignment is determined by the anisotropy of the polarizability of a molecule, light intensity, and temperature. We also discuss differences between the mechanisms of optical alignment in a gas phase and a liquid phase.     

Gene regulatory networks: a coarse-grained, equation-free approach to multiscale computation

We present computer-assisted methods for analyzing stochastic models of gene regulatory networks. The main idea that underlies this equation-free analysis is the design and execution of appropriately initialized short bursts of stochastic simulations; the results of these are processed to estimate coarse-grained quantities of interest, such as mesoscopic transport coefficients. In particular, using a simple model of a genetic toggle switch, we illustrate the computation of an effective free energy Phi and of a state-dependent effective diffusion coefficient D that characterize an unavailable effective Fokker-Planck equation. Additionally we illustrate the linking of equation-free techniques with continuation methods for performing a form of stochastic "bifurcation analysis"; estimation of mean switching times in the case of a bistable switch is also implemented in this equation-free context. The accuracy of our methods is tested by direct comparison with long-time stochastic simulations. This type of equation-free analysis appears to be a promising approach to computing features of the long-time, coarse-grained behavior of certain classes of complex stochastic models of gene regulatory networks, circumventing the need for long Monte Carlo simulations.     

A multiscale model for kinetics of formation and disintegration of spherical micelles

Dynamics of self-assembly and structural transitions in surfactant systems often involve a large span of length and time scales. A comprehensive understanding of these processes requires development of models connecting phenomena taking place on different scales. In this paper, we develop a multiscale model for formation and disintegration of spherical nonionic micelles. The study is performed under the assumption that the dominant mechanism of micelle formation (disintegration) is a stepwise addition (removal) of single monomers to (from) a surfactant aggregate. Different scales of these processes are investigated using a combination of coarse-grained molecular dynamics simulations, analytical and numerical solution of stochastic differential equations, and a numerical solution of kinetic equations. The removal of a surfactant from an aggregate is modeled by a Langevin equation for a single reaction coordinate, the distance between the centers of mass of the surfactant and the aggregate, with parameters obtained from a series of constrained molecular dynamics simulations. We demonstrate that the reverse process of addition of a surfactant molecule to an aggregate involves at least two additional degrees of freedom, orientation of the surfactant molecule and micellar microstructure. These additional degrees of freedom play an active role in the monomer addition process and neglecting their contribution leads to qualitative discrepancies in predicted surfactant addition rates. We propose a stochastic model for the monomer addition which takes the two additional degrees of freedom into account and extracts the model parameters from molecular dynamics simulations. The surfactant addition rates are determined from Brownian dynamics simulations of this model. The obtained addition and removal rates are then incorporated into the kinetic model of micellar formation and disintegration.     

Dynamic reaction coordinate in thermally fluctuating environment in the framework of the multidimensional generalized Langevin equations

A framework recently developed for the extraction of a dynamic reaction coordinate to mediate reactions buried in a multidimensional Langevin equation is extended to the generalized Langevin equations without a priori assumption of the forms of the potential (in general, nonlinearly coupled systems) and the friction kernel. The equation of motion with memory effect can be transformed into an equation without memory at the cost of an increase in the dimensionality of the system, and hence the theoretical framework developed for the (nonlinear) Langevin formulation can be generalized to the non-Markovian process with colored noise. It is found that the increased dimension can be physically interpreted as effective modes of the fluctuating environment. As an illustrative example, we apply this theory to a multidimensional generalized Langevin equation for motion on the Müller-Brown potential surface with an exponential friction kernel. Numerical simulations find a boundary between the highly reactive region and the less reactive region in the space of initial conditions. The location of the boundary is found to depend significantly on both the memory kernel and the nonlinear couplings. The theory extracts a reaction coordinate whose sign determines the fate of the reaction taking into account thermally fluctuating environments, memory effect, and nonlinearities. It is found that the location of the boundary of reactivity is satisfactorily reproduced as the zero of the statistical average of the new reaction coordinate, which is an analytical functional of both the original position coordinates and velocities of the system, and of the properties of the environment.     

Coarse-grained model for phospholipid/cholesterol bilayer

We construct a coarse-grained (CG) model for dipalmitoylphosphatidylcholine (DPPC)/cholesterol bilayers and apply it to large-scale simulation studies of lipid membranes. Our CG model is a two-dimensional representation of the membrane, where the individual lipid and sterol molecules are described by pointlike particles. The effective intermolecular interactions used in the model are systematically derived from detailed atomic-scale molecular dynamics simulations using the Inverse Monte Carlo technique, which guarantees that the radial distribution properties of the CG model are consistent with those given by the corresponding atomistic system. We find that the coarse-grained model for the DPPC/cholesterol bilayer is substantially more efficient than atomistic models, providing a speedup of approximately eight orders of magnitude. The results are in favor of formation of cholesterol-rich and cholesterol-poor domains at intermediate cholesterol concentrations, in agreement with the experimental phase diagram of the system. We also explore the limits of the coarse-grained model, and discuss the general validity and applicability of the present approach.     

Comparison between theoretical values and simulation results of viscosity for the dissipative particle dynamics method

In order to investigate the validity of the dissipative particle dynamics method, which is a mesoscopic simulation technique, we have derived an expression for viscosity from the equation of motion of dissipative particles. In the concrete, we have shown the Fokker-Planck equation in phase space, and macroscopic conservation equations such as the equation of continuity and the equation of momentum conservation. The basic equations of the single-particle and pair distribution functions have been derived using the Fokker-Planck equation. The solutions of these distribution functions have approximately been solved by the perturbation method under the assumption of molecular chaos. The expressions of the viscosity due to momentum and dissipative forces have been obtained using the approximate solutions of the distribution functions. Also, we have conducted nonequilibrium dynamics simulations to investigate the influence of the parameters, which have appeared in defining the equation of motion in the dissipative particle dynamics method. The theoretical values of the viscosity due to dissipative forces in the Hoogerbrugge-Koelman theory are in good agreement with the simulation results obtained by the nonequilibrium dynamics method, except in the range of small number densities. There are restriction conditions for taking appropriate values of the number density, number of particles, time interval, shear rate, etc., to obtain physically reasonable results by means of dissipative particle dynamics simulations.     

Stochastic theory of large-scale enzyme-reaction networks: finite copy number corrections to rate equation models

Chemical reactions inside cells occur in compartment volumes in the range of atto- to femtoliters. Physiological concentrations realized in such small volumes imply low copy numbers of interacting molecules with the consequence of considerable fluctuations in the concentrations. In contrast, rate equation models are based on the implicit assumption of infinitely large numbers of interacting molecules, or equivalently, that reactions occur in infinite volumes at constant macroscopic concentrations. In this article we compute the finite-volume corrections (or equivalently the finite copy number corrections) to the solutions of the rate equations for chemical reaction networks composed of arbitrarily large numbers of enzyme-catalyzed reactions which are confined inside a small subcellular compartment. This is achieved by applying a mesoscopic version of the quasisteady-state assumption to the exact Fokker-Planck equation associated with the Poisson representation of the chemical master equation. The procedure yields impressively simple and compact expressions for the finite-volume corrections. We prove that the predictions of the rate equations will always underestimate the actual steady-state substrate concentrations for an enzyme-reaction network confined in a small volume. In particular we show that the finite-volume corrections increase with decreasing subcellular volume, decreasing Michaelis-Menten constants, and increasing enzyme saturation. The magnitude of the corrections depends sensitively on the topology of the network. The predictions of the theory are shown to be in excellent agreement with stochastic simulations for two types of networks typically associated with protein methylation and metabolism.     

Computation of binding free energy with molecular dynamics and grand canonical Monte Carlo simulations

The binding of a ligand to a receptor is often associated with the displacement of a number of bound water molecules. When the binding site is exposed to the bulk region, this process may be sampled adequately by standard unbiased molecular dynamics trajectories. However, when the binding site is deeply buried and the exchange of water molecules with the bulk region may be difficult to sample, the convergence and accuracy in free energy perturbation (FEP) calculations can be severely compromised. These problems are further compounded when a reduced system including only the region surrounding the binding site is simulated. To address these issues, we couple molecular dynamics (MD) with grand canonical Monte Carlo (GCMC) simulations to allow the number of water to fluctuate during an alchemical FEP calculation. The atoms in a spherical inner region around the binding pocket are treated explicitly while the influence of the outer region is approximated using the generalized solvent boundary potential (GSBP). At each step during thermodynamic integration, the number of water in the inner region is equilibrated with GCMC and energy data generated with MD is collected. Free energy calculations on camphor binding to a deeply buried pocket in cytochrome P450cam, which causes about seven water molecules to be expelled, are used to test the method. It concluded that solvation free energy calculations with the GCMC/MD method can greatly improve the accuracy of the computed binding free energy compared to simulations with fixed number of water.     

Molecular simulations of solute transport in xylose isomerase crystals

Cross-linked enzyme crystals (CLECs) enclose an extensive regular matrix of chiral solvent-filled nanopores, via which ions and solutes travel in and out. Several cross-linked enzyme crystals have recently been used for chiral separation and as biocatalysts. We studied the dynamics of solute transport in orthorhombic d-xylose isomerase (XI) crystals by means of Brownian dynamics (BD) and molecular dynamics (MD) simulations, which show how the protein residues influence the dynamics of solute molecules in confined regions inside the lattice. In the BD simulations, coarse-grained beads represent solutes of different sizes. The diffusion of S-phenylglycine molecules inside XI crystals is investigated by long-time MD simulations. The computed diffusion coefficients within a crystal are found to be orders of magnitude lower than in bulk water. The simulation results are compared to the recent experimental studies of diffusion and reaction inside XI crystals. The insights obtained from simulations allow us to understand the nature of solute-protein interactions and transport phenomena in CLECs, which is useful for the design of novel nanoporous biocatalysts and bioseparations based on CLECs.     

Multidimensional vibrational spectroscopy for tunneling processes in a dissipative environment

Simulating tunneling processes as well as their observation are challenging problems for many areas. In this study, we consider a double-well potential system coupled to a heat bath with a linear-linear (LL) and square-linear (SL) system-bath interactions. The LL interaction leads to longitudinal (T1) and transversal (T2) homogeneous relaxations, whereas the SL interaction leads to the inhomogeneous dephasing (T2*) relaxation in the white noise limit with a rotating wave approximation. We discuss the dynamics of the double-well system under infrared (IR) laser excitations from a Gaussian-Markovian quantum Fokker-Planck equation approach, which was developed by generalizing Kubo's stochastic Liouville equation. Analytical expression of the Green function is obtained for a case of two-state-jump modulation by performing the Fourier-Laplace transformation. We then calculate a two-dimensional infrared signal, which is defined by the four-body correlation function of optical dipole, for various noise correlation time, system-bath coupling parameters, and temperatures. It is shown that the bath-induced vibrational excitation and relaxation dynamics between the tunneling splitting levels can be detected as the isolated off-diagonal peaks in the third-order two-dimensional infrared (2D-IR) spectroscopy for a specific phase matching condition. Furthermore, this spectroscopy also allows us to directly evaluate the rate constants for tunneling reactions, which relates to the coherence between the splitting levels; it can be regarded as a novel technique for measuring chemical reaction rates. We depict the change of reaction rates as a function of system-bath coupling strength and a temperature through the 2D-IR signal.     

Kramers problem for nonequilibrium current-induced chemical reactions

We discuss the use of tunneling electron current to control and catalyze chemical reactions. Assuming the separation of time scales for electronic and nuclear dynamics we employ Langevin equation for a reaction coordinate. The Langevin equation contains nonconservative current-induced forces and gives nonequilibrium, effective potential energy surface for current-carrying molecular systems. The current-induced forces are computed via Keldysh nonequilibrium Green's functions. Once a nonequilibrium, current-depended potential energy surface is defined, the chemical reaction is modeled as an escape of a Brownian particle from the potential well. We demonstrate that the barrier between the reactant and the product states can be controlled by the bias voltage. When the molecule is asymmetrically coupled to the electrodes, the reaction can be catalyzed or stopped depending on the polarity of the tunneling current.     

Coarse-graining kinetic networks in photosynthetic energy transport

Photosynthetic systems utilize hundreds of chlorophylls to collect sunlight and transport the energy to the reaction center with remarkably high quantum efficiency, however, the large size of the system together with the complex interactions among the components make it extremely challenging to understand the dynamics of light harvesting in large photosynthetic systems. To shed light on this problem, we present a structure-based theoretical framework that can be used to calculate transition rate matrix describing energy transport in photosynthetic systems and network clustering methods that provide simplified coarse-grained model revealing key structures guiding the light harvesting process. We constructed an effective model for energy transport in a Photosystem II supercomplex and applied several network clustering methods to generate coarse-grained kinetic cluster models for the system. Furthermore, we evaluated the performances of the network clustering methods, and show that a spectral clustering method and a minimum cut approach produce accurate coarse-grained models for the PSII-sc system. The results indicate that finding bottlenecks of energy transport is a crucial factor for reduced representations of photosynthetic light harvesting, and the overall work presented in this paper should provide a comprehensive theoretical framework to elucidate the dynamics of light harvesting in photosynthetic systems.     

How accurate are the nonlinear chemical Fokker-Planck and chemical Langevin equations?

The chemical Fokker-Planck equation and the corresponding chemical Langevin equation are commonly used approximations of the chemical master equation. These equations are derived from an uncontrolled, second-order truncation of the Kramers-Moyal expansion of the chemical master equation and hence their accuracy remains to be clarified. We use the system-size expansion to show that chemical Fokker-Planck estimates of the mean concentrations and of the variance of the concentration fluctuations about the mean are accurate to order Ω(-3∕2) for reaction systems which do not obey detailed balance and at least accurate to order Ω(-2) for systems obeying detailed balance, where Ω is the characteristic size of the system. Hence, the chemical Fokker-Planck equation turns out to be more accurate than the linear-noise approximation of the chemical master equation (the linear Fokker-Planck equation) which leads to mean concentration estimates accurate to order Ω(-1∕2) and variance estimates accurate to order Ω(-3∕2). This higher accuracy is particularly conspicuous for chemical systems realized in small volumes such as biochemical reactions inside cells. A formula is also obtained for the approximate size of the relative errors in the concentration and variance predictions of the chemical Fokker-Planck equation, where the relative error is defined as the difference between the predictions of the chemical Fokker-Planck equation and the master equation divided by the prediction of the master equation. For dimerization and enzyme-catalyzed reactions, the errors are typically less than few percent even when the steady-state is characterized by merely few tens of molecules.     

Eliminating fast reactions in stochastic simulations of biochemical networks: a bistable genetic switch

In many stochastic simulations of biochemical reaction networks, it is desirable to "coarse grain" the reaction set, removing fast reactions while retaining the correct system dynamics. Various coarse-graining methods have been proposed, but it remains unclear which methods are reliable and which reactions can safely be eliminated. We address these issues for a model gene regulatory network that is particularly sensitive to dynamical fluctuations: a bistable genetic switch. We remove protein-DNA and/or protein-protein association-dissociation reactions from the reaction set using various coarse-graining strategies. We determine the effects on the steady-state probability distribution function and on the rate of fluctuation-driven switch flipping transitions. We find that protein-protein interactions may be safely eliminated from the reaction set, but protein-DNA interactions may not. We also find that it is important to use the chemical master equation rather than macroscopic rate equations to compute effective propensity functions for the coarse-grained reactions.