Search

Atomic-level information is essential to explain the specific interactions governing protein-protein recognition in terms of structure and dynamics. Of particular interest is a characterization of the time-dependent kinetic aspects of protein-protein association and dissociation. A powerful framework to characterize the dynamics of complex molecular systems is provided by Markov State Models (MSMs). The central idea is to construct a reduced stochastic model of the full system by defining a set of conformational featured microstates and determining the matrix of transition probabilities between them. While a MSM framework can sometimes be very effective, different combinations of input featurization and simulation methods can significantly affect the robustness and the quality of the information generated from MSMs in the context of protein association. Here, a systematic examination of a variety of MSMs methodologies is undertaken to clarify these issues. To circumvent the uncertainties caused by sampling issues, we use a simplified coarse-grained model of the barnase-barstar protein complex. A sensitivity analysis is proposed to identify the microstates of an MSM that contribute most to the error in conjunction with the transition-based reweighting analysis method for a more efficient and accurate MSM construction.

Computational determination of the equilibrium state of heterogeneous phospholipid mem-branes is a significant challenge. We wish to explore the rich phase diagram of these multi-component systems. However, the diffusion and mixing times in membranes are long com-pared to typical times of computer simulations. To speed up the relaxation times, advanced simulation methods are used. We evaluate the combination of enhanced sampling tech-niques such as MDAS (Molecular Dynamics with Alchemical Steps) and MC-MD (Monte Carlo with Molecular Dynamics) with a coarse-grained model of membranes (Martini) to re-duce the number of steps and force evaluations that are needed to reach equilibrium. We illustrate a significant gain compared to straightforward Molecular Dynamics of the Martini model by factors between three to ten. The combination is a useful tool to enhance the study of phase separation and the formation of domains in biological membranes.

Determining the drug-target residence time (RT) is of major interest in drug discovery given that this kinetic parameter often represents a better indicator of in vivo drug efficacy than binding affinity. However, obtaining drug-target unbinding rates poses significant challenges, both computationally and experimentally. This is particularly palpable for complex systems like G Protein-Coupled Receptors (GPCRs) whose ligand unbinding typically requires very long timescales oftentimes inaccessible by standard molecular dynamics simulations. Enhanced sampling methods offer a useful alternative, and their efficiency can be further improved by using machine learning tools to identify optimal reaction coordinates. Here, we test the combination of two machine learning techniques, automatic mutual information noise omission and reweighted autoencoded variational Bayes for enhanced sampling, with infrequent metadynamics to efficiently study the unbinding kinetics of two classical drugs with different RTs in a prototypic GPCR, the μ-opioid receptor. Dissociation rates derived from these computations are within one order of magnitude from experimental values. We also use the simulation data to uncover the dissociation mechanisms of these drugs, shedding light on the structures of rate-limiting transition states, which, alongside metastable poses, are difficult to obtain experimentally but important to visualize when designing drugs with a desired kinetic profile.

State space truncation methods are widely used to approximate solutions of the chemical master equation. While most methods of this kind focus on truncating the state space directly, in this work we propose modifying the underlying chemical reaction network by introducing \emph{slack reactants} that indirectly truncate the state space. More specifically, slack reactants introduce an expanded chemical reaction network and impose a truncation scheme based on user defined properties, such as mass-conservation. This network structure allows us to prove inheritance of special properties of the original model, such as irreducibility and complex balancing. We use the network structure imposed by slack reactants to prove the convergence of the stationary distribution and first arrival times. We then provide examples comparing our method with the stationary finite state projection and finite buffer methods. Our slack reactant system appears to be more robust than some competing methods with respect to calculating first arrival times., 29 pages, 5 Figures

We have developed an algorithm coupling mesoscopic simulations on different levels in a hierarchy of Cartesian meshes. Based on the multiscale nature of the chemical reactions, some molecules in the system will live on a fine-grained mesh, while others live on a coarse-grained mesh. By allowing molecules to transfer from the fine levels to the coarse levels when appropriate, we show that we can save up to three orders of magnitude of computational time compared to microscopic simulations or highly resolved mesoscopic simulations, without losing significant accuracy. We demonstrate this in several numerical examples with systems that cannot be accurately simulated with a coarse-grained mesoscopic model.

TURBOMOLE is a collaborative, multi-national software development project aiming to provide highly efficient and stable computational tools for quantum chemical simulations of molecules, clusters, periodic systems, and solutions. The TURBOMOLE software suite is optimized for widely available, inexpensive, and resource-efficient hardware such as multi-core workstations and small computer clusters. TURBOMOLE specializes in electronic structure methods with outstanding accuracy-cost ratio, such as density functional theory including local hybrids and the random phase approximation (RPA), GW-Bethe-Salpeter methods, second-order Møller-Plesset theory, and explicitly correlated coupled-cluster methods. TURBOMOLE is based on Gaussian basis sets and has been pivotal for the development of many fast and low-scaling algorithms in the past three decades, such as integral-direct methods, fast multipole methods, the resolution-of-the-identity approximation, imaginary frequency integration, Laplace transform, and pair natural orbital methods. This review focuses on recent additions to TURBOMOLE's functionality, including excited-state methods, RPA and Green's function methods, relativistic approaches, high-order molecular properties, solvation effects, and periodic systems. A variety of illustrative applications along with accuracy and timing data are discussed. Moreover, available interfaces to users as well as other software are summarized. TURBOMOLE's current licensing, distribution, and support model are discussed, and an overview of TURBOMOLE's development workflow is provided. Challenges such as communication and outreach, software infrastructure, and funding are highlighted.

Conventional molecular dynamics simulations are incapable of sampling many important interactions in biomolecular systems due to their high dimensionality and rough energy landscapes. To observe rare events and calculate transition rates in these systems, enhanced sampling is a necessity. In particular, the study of ligand-protein interactions necessitates a diverse ensemble of protein conformations and transition states, and for many systems, this occurs on prohibitively long time scales. Previous strategies such as WExplore that can be used to determine these types of ensembles are hindered by problems related to the regioning of conformational space. Here, we propose a novel, regionless, enhanced sampling method that is based on the weighted ensemble framework. In this method, a value referred to as "trajectory variation" is optimized after each cycle through cloning and merging operations. This method allows for a more consistent measurement of observables and broader sampling resulting in the efficient exploration of previously unexplored conformations. We demonstrate the performance of this algorithm with the N-dimensional random walk and the unbinding of the trypsin-benzamidine system. The system is analyzed using conformation space networks, the residence time of benzamidine is confirmed, and a new unbinding pathway for the trypsin-benzamidine system is found. We expect that resampling of ensembles by variation optimization will be a useful general tool to broadly explore free energy landscapes.

Spatial stochastic models of single cell kinetics are capable of capturing both fluctuations in molecular numbers and the spatial dependencies of the key steps of intracellular regulatory networks. The spatial stochastic model can be simulated both on a detailed microscopic level using particle tracking and on a mesoscopic level using the reaction–diffusion master equation. However, despite substantial progress on simulation efficiency for spatial models in the last years, the computational cost quickly becomes prohibitively expensive for tasks that require repeated simulation of thousands or millions of realizations of the model. This limits the use of spatial models in applications such as multicellular simulations, likelihood-free parameter inference, and robustness analysis. Further approximation of the spatial dynamics is needed to accelerate such computational engineering tasks. We here propose a multiscale model where a compartment-based model approximates a detailed spatial stochastic model. The compartment model is constructed via a first-exit time analysis on the spatial model, thus capturing critical spatial aspects of the fine-grained simulations, at a cost close to the simple well-mixed model. We apply the multiscale model to a canonical model of negative-feedback gene regulation, assess its accuracy over a range of parameters, and demonstrate that the approximation can yield substantial speedups for likelihood-free parameter inference. eSSENCE

Computational strategies aimed at unveiling the thermodynamic and kinetic properties of G Protein-Coupled Receptor (GPCR) activation require extensive molecular dynamics simulations of the receptor embedded in an explicit lipid-water environment. A possible method for efficiently sampling the conformational space of such a complex system is metadynamics (MetaD) with path collective variables (CV). Here, we applied well-tempered MetaD with path CVs to one of the few GPCRs for which both inactive and fully active experimental structures are available, the μ-opioid receptor (MOR), and assessed the ability of this enhanced sampling method to estimate thermodynamic properties of receptor activation in line with those obtained by more computationally expensive adaptive sampling protocols. While n-body information theory (nBIT) analysis of these simulations confirmed that MetaD can efficiently characterize ligand-induced allosteric communication across the receptor, standard MetaD cannot be used directly to derive kinetic rates because transitions are accelerated by a bias potential. Applying the principle of Maximum Caliber (MaxCal) to the free-energy landscape of morphine-bound MOR reconstructed from MetaD, we obtained Markov State Models (MSMs) that yield kinetic rates of MOR activation in agreement with those obtained by adaptive sampling. Taken together, these results suggest that the MetaD-MaxCal combination creates an efficient strategy for estimating thermodynamic and kinetic properties of GPCR activation at an affordable computational cost.

We present a comprehensive benchmark study of the adsorption energy of a single water molecule on the (001) LiH surface using periodic coupled cluster and quantum Monte Carlo theories. We benchmark and compare different implementations of quantum chemical wave function based theories in order to verify the reliability of the predicted adsorption energies and the employed approximations. Furthermore we compare the predicted adsorption energies to those obtained employing widely used van der Waals density-functionals. Our findings show that quantum chemical approaches are becoming a robust and reliable tool for condensed phase electronic structure calculations, providing an additional tool that can also help in potentially improving currently available van der Waals density-functionals.

Evaluating solvation entropies directly and combining with direct energy calculations is one way of calculating free energies of solvation and is used by Inhomogeneous Fluid Solvation Theory (IFST). The configurational entropy of a fluid is a function of the interatomic correlations and can thus be expressed in terms of correlation functions. The entropies in this work are directly calculated from a truncated series of integrals over these correlation functions. Many studies truncate all terms higher than the solvent-solute correlations. This study includes an additional solvent-solvent correlation term and assesses the associated free energy when IFST is applied to a fixed Lennard-Jones particle solvated in neon. The strength of the central potential is varied to imitate larger solutes. Average free energy estimates with both levels of IFST are able to reproduce the estimate made using the Free energy Perturbation (FEP) to within 0.16 kcal/mol. We find that the signal from the solvent-solvent correlations is very weak. Our conclusion is that for monatomic fluids simulated by pairwise classical potentials the correction term is relatively small in magnitude. This study shows it is possible to reproduce the free energy from a path based method like FEP, by only considering the endpoints of the path. This method can be directly applied to more complex solutes which break the spherical symmetry of this study.

Single cells exhibit a significant amount of variability in transcript levels, which arises from slow, stochastic transitions between gene expression states. Elucidating the nature of these states and understanding how transition rates are affected by different regulatory mechanisms require state-of-the-art methods to infer underlying models of gene expression from single cell data. A Bayesian approach to statistical inference is the most suitable method for model selection and uncertainty quantification of kinetic parameters using small data sets. However, this approach is impractical because current algorithms are too slow to handle typical models of gene expression. To solve this problem, we first show that time-dependent mRNA distributions of discrete-state models of gene expression are dynamic Poisson mixtures, whose mixing kernels are characterized by a piecewise deterministic Markov process. We combined this analytical result with a kinetic Monte Carlo algorithm to create a hybrid numerical method that accelerates the calculation of time-dependent mRNA distributions by 1000-fold compared to current methods. We then integrated the hybrid algorithm into an existing Monte Carlo sampler to estimate the Bayesian posterior distribution of many different, competing models in a reasonable amount of time. We demonstrate that kinetic parameters can be reasonably constrained for modestly sampled data sets if the model is known a priori. If there are many competing models, Bayesian evidence can rigorously quantify the likelihood of a model relative to other models from the data. We demonstrate that Bayesian evidence selects the true model and outperforms approximate metrics typically used for model selection.

The probability density functions (PDFs) of the local measure of pressure as a function of the sampling volume are computed for a model Lennard-Jones (LJ) fluid using the Method of Planes (MOP) and Volume Averaging (VA) techniques. This builds on the study of Heyes, Dini, and Smith [J. Chem. Phys. 145, 104504 (2016)] which only considered the VA method for larger subvolumes. The focus here is typically on much smaller subvolumes than considered previously, which tend to the Irving-Kirkwood limit where the pressure tensor is defined at a point. The PDFs from the MOP and VA routes are compared for cubic subvolumes, V=ℓ3. Using very high grid-resolution and box-counting analysis, we also show that any measurement of pressure in a molecular system will fail to exactly capture the molecular configuration. This suggests that it is impossible to obtain the pressure in the Irving-Kirkwood limit using the commonly employed grid based averaging techniques. More importantly, below ℓ≈3 in LJ reduced units, the PDFs depart from Gaussian statistics, and for ℓ=1.0, a double peaked PDF is observed in the MOP but not VA pressure distributions. This departure from a Gaussian shape means that the average pressure is not the most representative or common value to arise. In addition to contributing to our understanding of local pressure formulas, this work shows a clear lower limit on the validity of simply taking the average value when coarse graining pressure from molecular (and colloidal) systems. EPSRC

The problem of predicting polymorphism in atomic and molecular crystals constitutes a significant challenge both experimentally and theoretically. From the theoretical viewpoint, polymorphism prediction falls into the general class of problems characterized by an underlying rough energy landscape, and consequently, free energy based enhanced sampling approaches can be brought to bear on the problem. In this paper, we build on a scheme previously introduced by two of the authors in which the lengths and angles of the supercell are targeted for enhanced sampling via temperature accelerated adiabatic free energy dynamics [T. Q. Yu and M. E. Tuckerman, Phys. Rev. Lett. 107, 015701 (2011)]. Here, that framework is expanded to include general order parameters that distinguish different crystalline arrangements as target collective variables for enhanced sampling. The resulting free energy surface, being of quite high dimension, is nontrivial to reconstruct, and we discuss one particular strategy for performing the free energy analysis. The method is applied to the study of polymorphism in xenon crystals at high pressure and temperature using the Steinhardt order parameters without and with the supercell included in the set of collective variables. The expected fcc and bcc structures are obtained, and when the supercell parameters are included as collective variables, we also find several new structures, including fcc states with hcp stacking faults. We also apply the new method to the solid-liquid phase transition in copper at 1300 K using the same Steinhardt order parameters. Our method is able to melt and refreeze the system repeatedly, and the free energy profile can be obtained with high efficiency.

Delay is an important and ubiquitous aspect of many biochemical processes. For example, delay plays a central role in the dynamics of genetic regulatory networks as it stems from the sequential assembly of first mRNA and then protein. Genetic regulatory networks are therefore frequently modeled as stochastic birth-death processes with delay. Here we examine the relationship between delay birth-death processes and their appropriate approximating delay chemical Langevin equations. We prove that the distance between these two descriptions, as measured by expectations of functionals of the processes, converges to zero with increasing system size. Further, we prove that the delay birth-death process converges to the thermodynamic limit as system size tends to infinity. Our results hold for both fixed delay and distributed delay. Simulations demonstrate that the delay chemical Langevin approximation is accurate even at moderate system sizes. It captures dynamical features such as the spatial and temporal distributions of transition pathways in metastable systems, oscillatory behavior in negative feedback circuits, and cross-correlations between nodes in a network. Overall, these results provide a foundation for using delay stochastic differential equations to approximate the dynamics of birth-death processes with delay.

Analysis of complex networks has been widely used as a powerful tool for investigating various physical, chemical, and biological processes. To understand the emergent properties of these complex systems, one of the most basic issues is to determine the structure and topology of the underlying networks. Recently, a new theoretical approach based on first-passage analysis has been developed for investigating the relationship between structure and dynamic properties for network systems with exponential dwell time distributions. However, many real phenomena involve transitions with non-exponential waiting times. We extend the first-passage method to uncover the structure of distinct pathways in complex networks with non-exponential dwell time distributions. It is found that the analysis of early time dynamics provides explicit information on the length of the pathways associated to their dynamic properties. It reveals a universal relationship that we have condensed in one general equation, which relates the number of intermediate states on the shortest path to the early time behavior of the first-passage distributions. Our theoretical predictions are confirmed by extensive Monte Carlo simulations.

We investigate the reversible diffusion-influenced reaction of an isolated pair in the context of the area reactivity model that describes the reversible binding of a single molecule in the presence of a binding site in terms of a generalized version of the Feynman-Kac equation in two dimensions. We compute the corresponding exact Green's function in the Laplace domain for both the initially unbound and bound molecule. We discuss convolution relations that facilitate the calculation of the binding and survival probabilities. Furthermore, we calculate an exact analytical expression for the Green's function in the time domain by inverting the Laplace transform via the Bromwich contour integral and derive expressions for the binding and survival probability in the time domain as well. We numerically confirm the accuracy of the obtained expressions by propagating the generalized Feynman-Kac equation in the time domain. Our results should be useful for comparing the area reactivity model with the contact reactivity model.

The nucleophilic attack of a chloride ion on methyl chloride is an important prototype SN2 reaction in organic chemistry that is known to be sensitive to the effects of the surrounding solvent. Herein, we develop a highly accurate Specific Reaction Parameter (SRP) model based on the Austin Model 1 Hamiltonian for chlorine to study the effects of solvation into an aqueous environment on the reaction mechanism. To accomplish this task, we apply high-level quantum mechanical calculations to study the reaction in the gas phase and combined quantum mechanical/molecular mechanical simulations with TIP3P and TIP4P-ew water models and the resulting free energy profiles are compared with those determined from simulations using other fast semi-empirical quantum models. Both gas phase and solution results with the SRP model agree very well with experiment and provide insight into the specific role of solvent on the reaction coordinate. Overall, the newly parameterized SRP Hamiltonian is able to reproduce both the gas phase and solution phase barriers, suggesting it is an accurate and robust model for simulations in the aqueous phase at greatly reduced computational cost relative to comparably accurate ab initio and density functional models.

The chemical master equation (CME) and the mathematically equivalent stochastic simulation algorithm (SSA) assume that the reactant molecules in a chemically reacting system are "dilute"and "well- mixed"throughout the containing volume. Here we clarify what those two conditions mean, and we show why their satisfaction is necessary in order for bimolecular reactions to physically occur in the manner assumed by the CME and the SSA. We prove that these conditions are closely connected, in that a system will stay well-mixed if and only if it is dilute. We explore the implications of these validity conditions for the reaction-diffusion (or spatially inhomogeneous) extensions of the CME and the SSA to systems whose containing volumes are not necessarily well-mixed, but can be partitioned into cubical subvolumes (voxels) that are. We show that the validity conditions, together with an additional condition that is needed to ensure the physical validity of the diffusion-induced jump probability rates of molecules between voxels, require the voxel edge length to have a strictly positive lower bound. We prove that if the voxel edge length is steadily decreased in a way that respects that lower bound, the average rate at which bimolecular reactions occur in the reaction-diffusion CME and SSA will remain constant, while the average rate of diffusive transfer reactions will increase as the inverse square of the voxel edge length. We conclude that even though the reaction-diffusion CME and SSA are inherently approximate, and cannot be made exact by shrinking the voxel size to zero, they should nevertheless be useful in many practical situations. © 2014 AIP Publishing LLC.

The technique of constraining high frequency modes of molecular motion is an effective way to increase simulation time scale and improve conformational sampling in molecular dynamics simulations. However, it has been shown that constraints on higher frequency modes such as bond lengths and bond angles stiffen the molecular model, thereby introducing systematic biases in the statistical behavior of the simulations. Fixman proposed a compensating potential to remove such biases in the thermodynamic and kinetic properties calculated from dynamics simulations. Previous implementations of the Fixman potential have been limited to only short serial chain systems. In this paper, we present a spatial operator algebra based algorithm to calculate the Fixman potential and its gradient within constrained dynamics simulations for branched topology molecules of any size. Our numerical studies on molecules of increasing complexity validate our algorithm by demonstrating recovery of the dihedral angle probability distribution function for systems that range in complexity from serial chains to protein molecules. We observe that the Fixman compensating potential recovers the free energy surface of a serial chain polymer, thus annulling the biases caused by constraining the bond lengths and bond angles. The inclusion of Fixman potential entails only a modest increase in the computational cost in these simulations. We believe that this work represents the first instance where the Fixman potential has been used for general branched systems, and establishes the viability for its use in constrained dynamics simulations of proteins and other macromolecules.

Markov state models (MSMs) have been demonstrated to be a powerful method for computationally studying intramolecular processes such as protein folding and macromolecular conformational changes. In this article, we present a new approach to construct MSMs that is applicable to modeling a broad class of multi-molecular assembly reactions. Distinct structures formed during assembly are distinguished by their undirected graphs, which are defined by strong subunit interactions. Spatial inhomogeneities of free subunits are accounted for using a recently developed Gaussian-based signature. Simplifications to this state identification are also investigated. The feasibility of this approach is demonstrated on two different coarse-grained models for virus self-assembly. We find good agreement between the dynamics predicted by the MSMs and long, unbiased simulations, and that the MSMs can reduce overall simulation time by orders of magnitude., Comment: 12 pages, 11 figures

A theory of nonlinear response of chemical kinetics, in which multiple perturbations are used to probe the time evolution of nonlinear chemical systems, is developed. Expressions for nonlinear chemical response functions and susceptibilities, which can serve as multidimensional measures of the kinetic pathways and rates, are derived. A new class of multidimensional measures that combine multiple perturbations and measurements is also introduced. Nonlinear fluctuation-dissipation relations for steady-state chemical systems, which replace operations of concentration measurement and perturbations, are proposed. Several applications to the analysis of complex reaction mechanisms are provided. © 2014 AIP Publishing LLC.

In this work we examine and extend the theory of proximal radial distribution functions for molecules in solution. We point out two formal extensions, the first of which generalizes the proximal distribution function hierarchy approach to the complete, angularly dependent molecular pair distribution function. Second, we generalize from the traditional right-handed solute-solvent proximal distribution functions to the left-handed distributions. The resulting neighbor hierarchy convergence is shown to provide a measure of the coarse-graining of the internal solute sites with respect to the solvent. Simulation of the test case of a deca-alanine peptide shows that this coarse-graining measure converges at a length scale of approximately 5 amino acids for the system considered.

The calculation of a protein-ligand binding free energy based on molecular dynamics (MD) simulations generally relies on a thermodynamic cycle in which the ligand is alchemically inserted into the system, both in the solvated protein and free in solution. The corresponding ligand-insertion free energies are typically calculated in nanoscale computational boxes simulated under periodic boundary conditions and considering electrostatic interactions defined by a periodic lattice-sum. This is distinct from the ideal bulk situation of a system of macroscopic size simulated under non-periodic boundary conditions with Coulombic electrostatic interactions. This discrepancy results in finite-size effects, which affect primarily the charging component of the insertion free energy, are dependent on the box size, and can be large when the ligand bears a net charge, especially if the protein is charged as well. This article investigates finite-size effects on calculated charging free energies using as a test case the binding of the ligand 2-amino-5-methylthiazole (net charge +1 e) to a mutant form of yeast cytochrome c peroxidase in water. Considering different charge isoforms of the protein (net charges -5, 0, +3, or +9 e), either in the absence or the presence of neutralizing counter-ions, and sizes of the cubic computational box (edges ranging from 7.42 to 11.02 nm), the potentially large magnitude of finite-size effects on the raw charging free energies (up to 17.1 kJ mol(-1)) is demonstrated. Two correction schemes are then proposed to eliminate these effects, a numerical and an analytical one. Both schemes are based on a continuum-electrostatics analysis and require performing Poisson-Boltzmann (PB) calculations on the protein-ligand system. While the numerical scheme requires PB calculations under both non-periodic and periodic boundary conditions, the latter at the box size considered in the MD simulations, the analytical scheme only requires three non-periodic PB calculations for a given system, its dependence on the box size being analytical. The latter scheme also provides insight into the physical origin of the finite-size effects. These two schemes also encompass a correction for discrete solvent effects that persists even in the limit of infinite box sizes. Application of either scheme essentially eliminates the size dependence of the corrected charging free energies (maximal deviation of 1.5 kJ mol(-1)). Because it is simple to apply, the analytical correction scheme offers a general solution to the problem of finite-size effects in free-energy calculations involving charged solutes, as encountered in calculations concerning, e.g., protein-ligand binding, biomolecular association, residue mutation, pKa and redox potential estimation, substrate transformation, solvation, and solvent-solvent partitioning.

The dependence of various types of linear signals on the phase profile of broadband optical pulses is examined using fundamental time translation invariance symmetry of multipoint correlation functions. The frequency-domain wave-mixing analysis presented here unifies several arguments made earlier with respect to the conditions whereby coherent control schemes may be used. © 2013 AIP Publishing LLC.

Excited-state vibrational dynamics in molecules can be studied by an electronically off-resonant Raman process induced by a probe pulse with variable delay with respect to an actinic pulse. We establish the connection between several variants of the technique that involve either spontaneous or stimulated Raman detection and different pulse configurations. By using loop diagrams in the frequency domain, we show that all signals can be described as six wave mixing which depend on the same four point molecular correlation functions involving two transition dipoles and two polarizabilities and accompanied by a different gating. Simulations for the stochastic two-state-jump model illustrate the origin of the absorptive and dispersive features observed experimentally. © 2013 AIP Publishing LLC.

The thermodynamic perturbation theory was tested against newly obtained Monte Carlo computer simulations to describe the major features of the hydrophobic effect in a simple 3D-Mercedes-Benz water model: the temperature and hydrophobe size dependence on entropy, enthalpy, and free energy of transfer of a simple hydrophobic solute into water. An excellent agreement was obtained between the theoretical and simulation results. Further, the thermodynamic perturbation theory qualitatively correctly (with respect to the experimental data) describes the solvation thermodynamics under conditions where the simulation results are difficult to obtain with good enough accuracy, e.g., at high pressures.

The diffusion tensor of complex macromolecules in Stokes flow is often approximated by the bead models. The bead models are known to reproduce the experimental diffusion coefficients of a single macromolecule, but the accuracy of their calculation of the whole multi-body diffusion tensor, which is important for Brownian dynamics simulations, has not been closely investigated. As a first step, we assess the accuracy of the bead model calculated diffusion tensor of two spheres. Our results show that the bead models produce very accurate diffusion tensors for two spheres where a reasonable number of beads are used and there is no bead overlap.

Effective coarse-grained representations of protein-protein interaction potentials are vital in the modeling of large scale systems. We develop a method to fit an arbitrary number of effective charges to approximate the electrostatic potential of a protein at a given pH in an ionic solution. We find that the effective charges can reproduce an input potential calculated from a high resolution Poisson-Boltzmann calculation. Since the effective charges used in this model are not constrained to the locations of the original charged groups, the extra degrees of freedom allows us to reproduce the field anisotropy with fewer charges. The fitting procedure uses a number of approximations in the charge magnitudes, initial conditions, and multipoles to speed convergence. The most significant gains are found by fitting the multipole moments of the effective charge potential to the moments of the original field. We show that the Yukawa potential is not only sufficient as a pairwise summation in reproducing the potential, but comes naturally from the linearized expansion of the Poisson-Boltzmann equation. We compute interaction energies and find excellent agreement to the original potential. From the effective charge model we compute the electrostatic contribution to the second virial coefficient.

The observables of multidimensional infrared spectroscopy may be calculated from nonlinear vibrational response functions. Fully quantum dynamical calculations of vibrational response functions are generally impractical, while completely classical calculations are qualitatively incorrect at long times. These challenges motivate the development of semiclassical approximations to quantum mechanics, which use classical mechanical information to reconstruct quantum effects. The mean-trajectory (MT) approximation is a semiclassical approach to quantum vibrational response functions employing classical trajectories linked by deterministic transitions representing the effects of the radiation-matter interaction. Previous application of the MT approximation to the third-order response function R(3)(t3, t2, t1) demonstrated that the method quantitatively describes the coherence dynamics of the t3 and t1 evolution times, but is qualitatively incorrect for the waiting-time t2 period. Here we develop an optimized version of the MT approximation by elucidating the connection between this semiclassical approach and the double-sided Feynman diagrams (2FD) that represent the quantum response. Establishing the direct connection between 2FD and semiclassical paths motivates a systematic derivation of an optimized MT approximation (OMT). The OMT uses classical mechanical inputs to accurately reproduce quantum dynamics associated with all three propagation times of the third-order vibrational response function.

The interpretation of protein amide I infrared spectra has been greatly assisted by the observation that the vibrational frequency of a peptide unit reports on its local electrostatic environment. However, the interpretation of spectra remains largely qualitative due to a lack of direct quantitative connections between computational models and experimental data. Here, we present an empirical parameterization of an electrostatic amide I frequency map derived from the infrared absorption spectra of 28 dipeptides. The observed frequency shifts are analyzed in terms of the local electrostatic potential, field, and field gradient, evaluated at sites near the amide bond in molecular dynamics simulations. We find that the frequency shifts observed in experiment correlate very well with the electric field in the direction of the C=O bond evaluated at the position of the amide oxygen atom. A linear best-fit mapping between observed frequencies and electric field yield sample standard deviations of 2.8 and 3.7 cm(-1) for the CHARMM27 and OPLS-AA force fields, respectively, and maximum deviations (within our data set) of 9 cm(-1). These results are discussed in the broader context of amide I vibrational models and the effort to produce quantitative agreement between simulated and experimental absorption spectra.

We investigate the reversible diffusion-influenced reaction of an isolated pair in the presence of a non-Markovian generalization of the backreaction boundary condition in two space dimensions. Following earlier work by Agmon and Weiss, we consider residence time probability densities that decay slower than an exponential and that are characterized by a single parameter 0σ ≤ 1. We calculate an exact expression for a Green's function of the two-dimensional diffusion equation subject to a non-Markovian backreaction boundary condition that is valid for arbitrary σ and for all times. We use the obtained expression to derive the survival probability for the initially unbound pair and we calculate an exact expression for the probability S(t[line]*) that the initially bound particle is unbound. Finally, we obtain an approximate solution for long times. In particular, we show that the ultimate fate of the bound state is complete dissociation, as in the Markovian case. However, the limiting value is approached quite differently: Instead of a ~t(-1) decay, we obtain 1 - S(t[line]*) ~ t(-σ)ln t. The derived expressions should be relevant for a better understanding of reversible membrane-bound reactions in cell biology.

Increasingly, theoretical studies of proteins focus on large systems. This trend demands the development of computational models that are fast, to overcome the growing complexity, and accurate, to capture the physically relevant features. To address this demand, we introduce a protein model that uses all-atom architecture to ensure the highest level of chemical detail while employing effective pair potentials to represent the effect of solvent to achieve the maximum speed. The effective potentials are derived for amino acid residues based on the condition that the solvent-free model matches the relevant pair-distribution functions observed in explicit solvent simulations. As a test, the model is applied to alanine polypeptides. For the chain with 10 amino acid residues, the model is found to reproduce properly the native state and its population. Small discrepancies are observed for other folding properties and can be attributed to the approximations inherent in the model. The transferability of the generated effective potentials is investigated in simulations of a longer peptide with 25 residues. A minimal set of potentials is identified that leads to qualitatively correct results in comparison with the explicit solvent simulations. Further tests, conducted for multiple peptide chains, show that the transferable model correctly reproduces the experimentally observed tendency of polyalanines to aggregate into β-sheets more strongly with the growing length of the peptide chain. Taken together, the reported results suggest that the proposed model could be used to succesfully simulate folding and aggregation of small peptides in atomic detail. Further tests are needed to assess the strengths and limitations of the model more thoroughly.

We investigate reversible diffusion-influenced reactions of an isolated pair in two dimensions. To this end, we employ convolution relations that permit deriving the survival probability of the reversible reaction directly in terms of the survival probability of the irreversible reaction. Furthermore, we make use of the mean reaction time approximation to write the irreversible survival probability in restricted spaces as a single exponential. In this way, we obtain exact and approximative expressions in the time domain for the reversible survival probability for three different two dimensional diffusion spaces: The infinite and restricted plane as well as the surface of a sphere. Our obtained results should prove useful in the context of membrane-bound reversible diffusion-influenced reactions in cell biology., Comment: 3 figures

We have developed a method for modeling spatial stochastic biochemical reactions in complex, three-dimensional, and time-dependent domains using the reaction-diffusion master equation formalism. In particular, we look to address the fully coupled problems that arise in systems biology where the shape and mechanical properties of a cell are determined by the state of the biochemistry and vice versa. To validate our method and characterize the error involved, we compare our results for a carefully constructed test problem to those of a microscale implementation. We demonstrate the effectiveness of our method by simulating a model of polarization and shmoo formation during the mating of yeast. The method is generally applicable to problems in systems biology where biochemistry and mechanics are coupled, and spatial stochastic effects are critical.

In the present work, we address a major challenge facing the modeling of biochemical reaction networks: when using stochastic simulations, the computational load and number of unknown parameters may dramatically increase with system size and complexity. A proposed solution to this challenge is the reduction of models by utilizing nonlinear reaction rate laws in place of a complex multi-reaction mechanism. This type of model reduction in stochastic systems often fails when applied outside of the context in which it was initially conceived. We hypothesize that the use of nonlinear rate laws fails because a single reaction is inherently Poisson distributed and cannot match higher order statistics. In this study we explore the use of Hill-type rate laws as an approximation for gene regulation, specifically transcription repression. We matched output data for several simple gene networks to determine Hill-type parameters. We show that the models exhibit inaccuracies when placed into a simple feedback repression model. By adding an additional abstract reaction to the models we account for second-order statistics. This split Hill rate law matches higher order statistics and demonstrates that the new model is able to more accurately describe the mean protein output. Finally, the modified Hill model is shown to be modular and models retain accuracy when placed into a larger multi-gene network. The work as presented may be used in gene regulatory or cell-signaling networks, where multiple binding events can be captured by Hill kinetics. The added benefit of the proposed split-Hill kinetics is the improved accuracy in modeling stochastic effects. We demonstrate these benefits with a few specific reaction network examples.

A new algorithm, "HiER-leap", is derived which improves on the computational properties of the ER-leap algorithm for exact accelerated simulation of stochastic chemical kinetics. Unlike ER-leap, HiER-leap utilizes a hierarchical or divide-and-conquer organization of reaction channels into tightly coupled "blocks" and is thereby able to speed up systems with many reaction channels. Like ER-leap, HiER-leap is based on the use of upper and lower bounds on the reaction propensities to define a rejection sampling algorithm with inexpensive early rejection and acceptance steps. But in HiER-leap, large portions of intra-block sampling may be done in parallel. An accept/reject step is used to synchronize across blocks. This method scales well when many reaction channels are present and has desirable asymptotic properties. The algorithm is exact, parallelizable and achieves a significant speedup over SSA and ER-leap on certain problems. This algorithm offers a potentially important step towards efficient in silico modeling of entire organisms., Comment: 22 pages, 3 figures

How to describe nondynamic electron correlation is still a major challenge to density functional theory (DFT). Recent models designed particularly for this problem, such as Becke'05 (B05) and Perdew-Staroverov-Tao-Scuseria (PSTS) functionals employ the exact-exchange density, the efficient calculation of which is technically quite challenging. We have recently implemented self-consistently the B05 functional based on an efficient resolution-identity (RI) technique. In this study, we report a self-consistent RI implementation of the PSTS functional. In contrast to its original implementation, our version brings no limitation on the choice of the basis set. We have also implemented the Mori-Sanchez-Cohen-Yang-2 (MCY2) functional, another recent DFT method that includes full exact exchange. The performance of PSTS, B05, and MCY2 is validated on thermochemistry, reaction barriers, and dissociation energy curves, with an emphasis on nondynamic correlation effects in the discussion. All three methods perform rather well in general, B05 and MCY2 being on average somewhat better than PSTS. We include also results with other functionals that represent various aspects of the development in this field in recent years, including B3LYP, M06-HF, M06-2X, ωB97X, and TPSSh. The performance of the heavy-parameterized functionals M06-2X and ωB97X is on average better than that of B05, MCY2, and PSTS for standard thermodynamic properties and reactions, while the latter functionals do better in hydrogen abstraction reactions and dissociation processes. In particular, B05 is found to be the only functional that yields qualitatively correct dissociation curves for two-center symmetric radicals like \documentclass[12pt]{minimal}\begin{document}${\rm He}_{2}^{+}$\end{document} He 2+. Finally, we compare the performance of all these functionals on a strongly correlated exemplary case system, the NO dimer. Only PSTS, B05, and MCY2 describe the system qualitatively correctly. Overall, this new type of functionals show good promise of overcoming some of the difficulties DFT encounters for systems with strong nondynamic correlation.

A self-consistent method is presented for the calculation of the local dielectric permittivity and electrostatic potential generated by a solute of arbitrary shape and charge distribution in a polar and polarizable liquid. The structure and dynamics behavior of the liquid at the solute/liquid interface determine the spatial variations of the density and the dielectric response. Emphasis here is on the treatment of the interface. The method is an extension of conventional methods used in continuum protein electrostatics, and can be used to estimate changes in the static dielectric response of the liquid as it adapts to charge redistribution within the solute. This is most relevant in the context of polarizable force fields, during electron structure optimization in quantum chemical calculations, or upon charge transfer. The method is computationally efficient and well suited for code parallelization, and can be used for on-the-fly calculations of the local permittivity in dynamics simulations of systems with large and heterogeneous charge distributions, such as proteins, nucleic acids, and polyelectrolytes. Numerical calculation of the system free energy is discussed for the general case of a liquid with field-dependent dielectric response.

We derive an exact Green's function of the diffusion equation for a pair of disk-shaped interacting particles in two dimensions subject to a backreaction boundary condition. Furthermore, we use the obtained function to calculate exact expressions for the survival probability and the time-dependent rate coefficient for the initially unbound pair and the survival probability of the bound state. The derived expressions will be of particular utility for the description of reversible membrane-bound reactions in cell biology.

Many biochemical networks have complex multidimensional dynamics and there is a long history of methods that have been used for dimensionality reduction for such reaction networks. Usually a deterministic mass action approach is used; however, in small volumes, there are significant fluctuations from the mean which the mass action approach cannot capture. In such cases stochastic simulation methods should be used. In this paper, we evaluate the applicability of one such dimensionality reduction method, the quasi-steady state approximation (QSSA) [L. Menten and M. Michaelis, “Die kinetik der invertinwirkung,” Biochem. Z 49, 333369 (1913)] for dimensionality reduction in case of stochastic dynamics. First, the applicability of QSSA approach is evaluated for a canonical system of enzyme reactions. Application of QSSA to such a reaction system in a deterministic setting leads to Michaelis-Menten reduced kinetics which can be used to derive the equilibrium concentrations of the reaction species. In the case of stochastic simulations, however, the steady state is characterized by fluctuations around the mean equilibrium concentration. Our analysis shows that a QSSA based approach for dimensionality reduction captures well the mean of the distribution as obtained from a full dimensional simulation but fails to accurately capture the distribution around that mean. Moreover, the QSSA approximation is not unique. We have then extended the analysis to a simple bistable biochemical network model proposed to account for the stability of synaptic efficacies; the substrate of learning and memory [J. E. Lisman, “A mechanism of memory storage insensitive to molecular turnover: A bistable autophosphorylating kinase,” Proc. Natl. Acad. Sci. U.S.A. 82, 3055–3057 (1985)]10.1073/pnas.82.9.3055. Our analysis shows that a QSSA based dimensionality reduction method results in errors as big as two orders of magnitude in predicting the residence times in the two stable states.

We calculate the pair diffusion coefficient D(r) as a function of the distance r between two hard-sphere particles in a dense monodisperse suspension. The distance-dependent pair diffusion coefficient describes the hydrodynamic interactions between particles in a fluid that are central to theories of polymer and colloid dynamics. We determine D(r) from the propagators (Green's functions) of particle pairs obtained from discontinuous molecular dynamics simulations. At distances exceeding 3 molecular diameters, the calculated pair diffusion coefficients are in excellent agreement with predictions from exact macroscopic hydrodynamic theory for large Brownian particles suspended in a solvent bath, as well as the Oseen approximation. However, the asymptotic 1/r distance dependence of D(r) associated with hydrodynamic effects emerges only after the pair distance dynamics has been followed for relatively long times, indicating non-negligible memory effects in the pair diffusion at short times. Deviations of the calculated D(r) from the hydrodynamic models at short distances r reflect the underlying many-body fluid structure, and are found to be correlated to differences in the local available volume. The procedure used here to determine the pair diffusion coefficients can also be used for single-particle diffusion in confinement with spherical symmetry., 6 pages, 5 figures

The generalized Fick-Jacobs equation is widely used to study diffusion of Brownian particles in three-dimensional tubes and quasi-two-dimensional channels of varying constraint geometry. We show how this equation can be applied to study the slowdown of unconstrained diffusion in the presence of obstacles. Specifically, we study diffusion of a point Brownian particle in the presence of identical cylindrical obstacles arranged in a square lattice. The focus is on the effective diffusion coefficient of the particle in the plane perpendicular to the cylinder axes, as a function of the cylinder radii. As radii vary from zero to one half of the lattice period, the effective diffusion coefficient decreases from its value in the obstacle free space to zero. Using different versions of the generalized Fick-Jacobs equation, we derive simple approximate formulas, which give the effective diffusion coefficient as a function of the cylinder radii, and compare their predictions with the values of the effective diffusion coefficient obtained from Brownian dynamics simulations. We find that both Reguera-Rubi and Kalinay-Percus versions of the generalized Fick-Jacobs equation lead to quite accurate predictions of the effective diffusion coefficient (with maximum relative errors below 4% and 7%, respectively) over the entire range of the cylinder radii from zero to one half of the lattice period.

In many physical, chemical, and biological systems energy and charge transfer processes are of utmost importance. To determine the influence of the environment on these transport processes, equilibrium molecular dynamics simulations become more and more popular. From these simulations, one usually determines the thermal fluctuations of certain energy gaps, which are then either used to perform ensemble-averaged wave packet simulations, also called Ehrenfest dynamics, or to employ a density matrix approach via spectral densities. These two approaches are analyzed through energy gap fluctuations that are generated to correspond to a predetermined spectral density. Subsequently, density matrix and wave packet simulations are compared through population dynamics and absorption spectra for different parameter regimes. Furthermore, a previously proposed approach to enforce the correct long-time behavior in the wave packet simulations is probed and an improvement is proposed.

Diffusion is often the rate determining step in many biological processes. Currently, the two main computational methods for studying diffusion are stochastic methods, such as Brownian dynamics, and continuum methods, such as the finite element method. This paper proposes a new hybrid diffusion method that couples the strengths of each of these two methods. The method is derived for a general multidimensional system, and is presented using a basic test case for 1D linear and radially symmetric diffusion systems.

The weighted histogram analysis method (WHAM) is routinely used for computing free energies and expectations from multiple ensembles. Existing derivations of WHAM require observations to be discretized into a finite number of bins. Yet, WHAM formulas seem to hold even if the bin sizes are made arbitrarily small. The purpose of this article is to demonstrate both the validity and value of the multi-state Bennet acceptance ratio (MBAR) method seen as a binless extension of WHAM. We discuss two statistical arguments to derive the MBAR equations, in parallel to the self-consistency and maximum likelihood derivations already known for WHAM. We show that the binless method, like WHAM, can be used not only to estimate free energies and equilibrium expectations, but also to estimate equilibrium distributions. We also provide a number of useful results from the statistical literature, including the determination of MBAR estimators by minimization of a convex function. This leads to an approach to the computation of MBAR free energies by optimization algorithms, which can be more effective than existing algorithms. The advantages of MBAR are illustrated numerically for the calculation of absolute protein-ligand binding free energies by alchemical transformations with and without soft-core potentials. We show that binless statistical analysis can accurately treat sparsely distributed interaction energy samples as obtained from unmodified interaction potentials that cannot be properly analyzed using standard binning methods. This suggests that binless multi-state analysis of binding free energy simulations with unmodified potentials offers a straightforward alternative to the use of soft-core potentials for these alchemical transformations.

We present here the implementation of a self-consistent approach to the calculation of excitation energies within regular Kohn-Sham density functional theory. The method is based on the n-order constricted variational density functional theory (CV(n)-DFT) [T. Ziegler, M. Seth, M. Krykunov, J. Autschbach, and F. Wang, J. Chem. Phys. 130, 154102 (2009)] doi: 10.1063/1.3114988. and its self-consistent formulation (SCF-CV(∞)-DFT) [J. Cullen, M. Krykunov, and T. Ziegler, Chem. Phys. 391, 11 (2011)] doi: 10.1016/j.chemphys.2011.05.021.. A full account is given of the way in which SCF-CV(∞)-DFT is implemented. The SCF-CV(∞)-DFT scheme is further applied to transitions from occupied π orbitals to virtual π* orbitals. The same series of transitions has been studied previously by high-level ab initio methods. We compare here the performance of SCF-CV(∞)-DFT to that of time dependent density functional theory (TD-DFT), CV(n)-DFT and ΔSCF-DFT, with the ab initio results as a benchmark standard. It is finally demonstrated how adiabatic TD-DFT and ΔSCF-DFT are related through different approximations to SCF-CV(∞)-DFT.

When a diffusing particle escapes from a spherical cavity through a narrow, not too long tunnel, the escape kinetics is essentially single-exponential. The presence of reversible binding sites on the cavity wall leads to retention of the particle in the system and converts the single-exponential kinetics into bi-exponential. We develop a theory that describes these effects. The theory shows how the delay time and the average number of binding events depend on the geometric and kinetic parameters of the system. To study the effect of the cavity shape, we also analyze the kinetics when the particle escapes from a cylindrical cavity with reversible binding sites.

The Fourier transform Coulomb (FTC) method has been shown to be effective for the fast and accurate calculation of long-range Coulomb interactions between diffuse (low-energy cutoff) densities in quantum mechanical (QM) systems. In this work, we split the potential of a compact (high-energy cutoff) density into short-range and long-range components, similarly to how point charges are handled in the Ewald mesh methods in molecular mechanics simulations. With this linear scaling QM Ewald mesh method, the long-range potential of compact densities can be represented on the same grid as the diffuse densities that are treated by the FTC method. The new method is accurate and significantly reduces the amount of computational time on short-range interactions, especially when it is compared to the continuous fast multipole method.