A generalized Langevin dynamics approach to model solvent dynamics effects on proteins via a solvent-accessible surface. The carboxypeptidase A inhibitor protein as a model

A generalized Langevin dynamics (GLD) scheme is derived for (bio)macromolecules having internal structure, arbitrary shapes and a size larger than solvent molecules (i.e. proteins). The concept of solvent-accessible surface area (SASA) is used to incorporate solvent effects via external forces thereby avoiding its explicit molecular representation. A simulation algorithm is implemented in the GROMOS molecular dynamics (MD) program including random forces and memory effects, while solvation effects enter via derivatives of the surface area. The potato carboxypeptidase inhibitor (PCI), a small protein, is used to numerically test the approach. This molecule has N- and C-terminal tails whose structure and fluctuations are solvent dependent. A 1-ns MD trajectory was analyzed in depth. X-ray and NMR structures are used in conjunction with MD simulations with and without explicit solvent to gauge the quality of the results. All the analyses showed that the GLD simulation approached the results obtained for the MD simulation with explicit simple-point-charge-model water molecules. The SASAs of the polar atoms show a natural exposure towards the solvent direction. A FLS solvent simulation was completed in order to sense memory effects. The approach and results presented here could be of great value for developing alternatives to the use of explicit solvent molecules in the MD simulation of proteins, expanding its use and the time-scale explored. Received: 2 February 2000 / Revised: 12 March 2000 / Accepted: 26 May 2000 / Published online: 2 November 2000     

Ensembles of microelectrodes: A digital- simulation

The rate of diffusion to ensembles of microelectrodes is calculated by digital simulation, based on the explicit finite-difference technique. Two limiting cases are shown in the current-time relationship at short times, when the diffusion layer is small in comparison to the radius of each microelectrode, the current behaves as in the case of semi-infinite planar diffusion to the active area on the surface. At long times, when the diffusion layers of adjacent electrodes overlap significantly, the current approaches the case of planar diffusion to the total area of the surface, including the non-active regions. The results are compared with previous calculations based on analytical solutions, which were derived on the basis of different approximations. Good agreement between reported experimental data and results based on simulation is observed.Edge effects arising as a result of contributions of radial components to the total diffusion rate are shown to be significant, even on relatively large electrodes. The errors which arise due to such effects are estimated as a function of the size of the electrode and the duration of measurement.     

Polymerization molecular dynamics simulations. I. Cross-linked atomistic models for poly(methacrylate) networks

A methodology to build cross-linked atomistic structures for poly(methacrylates) is presented. The methodology is based on a polymerization molecular dynamics (MD) scheme in which monomers are allowed to react with each other during an MD simulation. The criteria for forming chemical bonds is based on distances between prespecified reactive centres on the monomers and growing polymer chains. The simulated rate of monomer diffusion is increased through the use of scaled charges and approximations in non-bond energy calculations. Local energetic stresses caused by bond formation are relaxed through the use of an explicit velocity rescaling molecular dynamics algorithm. An example network consisting of a typical dental resin is constructed and presented.     

A modified next reaction method for simulating chemical systems with time dependent propensities and delays

Chemical reaction systems with a low to moderate number of molecules are typically modeled as discrete jump Markov processes. These systems are oftentimes simulated with methods that produce statistically exact sample paths such as the Gillespie algorithm or the next reaction method. In this paper we make explicit use of the fact that the initiation times of the reactions can be represented as the firing times of independent, unit rate Poisson processes with internal times given by integrated propensity functions. Using this representation we derive a modified next reaction method and, in a way that achieves efficiency over existing approaches for exact simulation, extend it to systems with time dependent propensities as well as to systems with delays.     

Adaptive finite element simulation of chronoamperometry at microdisc electrodes

In this paper we shall introduce a transient finite element algorithm by considering the simplest problem of numerical simulation of the chronoamperometric current for an E reaction mechanism at a microdisc electrode. Such a numerical simulation is made difficult by the presence of a boundary singularity where the electrode meets the insulator (often known as the ‘edge-effect’) and by a time-singularity caused by the impulsive start of the experiment. We attempt to overcome these problems by using an adaptive finite element algorithm in which we derive an a posteriori bound on the error in the computed current. This is used to guide mesh refinement and adaptive time-stepping, resulting in a fully automated algorithm which is both accurate and efficient.     

Introduction

The formalism for a configuration interaction approach is presented, in which explicit introduction of interelectronic coordinates into individual configurations is accomplished through the use of spherical Gaussian correlation factors. Formulas for all required matrix elements over these configurations are derived, and a computationally convenient algorithm to evaluate the matrix elements is presented. In addition, analysis of the form of the correlation factor shows that both angular and radial correlation can be obtained using spherical Gaussian correlation factors.     

Fast second order electron transfer reactions coupled to redox protein electrochemistry: Experiment and digital simulation

Chemical modification of metal surfaces by chemisorption provides a versatile method for the production of electrode interfaces which can be selective for the direct electrochemistry of one redox protein over another. The electrochemistry of a mixture of horse heart cytochrome c and spinach plastocyanin has been investigated at gold surfaces made selective for first one and then the other protein. The resulting cyclic voltammetry is quite unusual, containing pre-shoulders to both reduction current and reoxidation current peaks. The results have been interpreted in terms of fast second order electron transfer reactions taking place between the two proteins in homogeneous solution. This rationalisation has been corroborated by an explicit digital simulation of the proposed reaction scheme, using second order RKI. There are three independently variable parameters to the simulation: forward kinetic parameter, reverse kinetic parameter, and concentration ratio of non-electrode-active species to electrode-active species. The simulation has been used to explore a number of interesting trends in these parameters. Five such sequences of simulated cyclic voltammograms are reported, together with peak current and potential data in most cases. Attention is drawn to the possibility for further interesting experimental mixed redox protein electrochemistry at selective surfaces.     

Application of Parallel Hybrid Algorithm in Massively Parallel GPGPU—The Improved Effective and Efficient Method for Calculating Coulombic Interactions in Simulations of Many Ions with SIMION

In our previous study, we introduced a new hybrid approach to effectively approximate the total force on each ion during a trajectory calculation in mass spectrometry device simulations, and the algorithm worked successfully with SIMION. We took one step further and applied the method in massively parallel general-purpose computing with GPU (GPGPU) to test its performance in simulations with thousands to over a million ions. We took extra care to minimize the barrier synchronization and data transfer between the host (CPU) and the device (GPU) memory, and took full advantage of the latency hiding. Parallel codes were written in CUDA C++ and implemented to SIMION via the user-defined Lua program. In this study, we tested the parallel hybrid algorithm with a couple of basic models and analyzed the performance by comparing it to that of the original, fully-explicit method written in serial code. The Coulomb explosion simulation with 128,000 ions was completed in 309?s, over 700 times faster than the 63?h taken by the original explicit method in which we evaluated two-body Coulomb interactions explicitly on one ion with each of all the other ions. The simulation of 1,024,000 ions was completed in 2650?s. In another example, we applied the hybrid method on a simulation of ions in a simple quadrupole ion storage model with 100,000 ions, and it only took less than 10 d. Based on our estimate, the same simulation is expected to take 5-7 y by the explicit method in serial code.     

Avoiding negative populations in explicit Poisson tau-leaping

The explicit tau-leaping procedure attempts to speed up the stochastic simulation of a chemically reacting system by approximating the number of firings of each reaction channel during a chosen time increment tau as a Poisson random variable. Since the Poisson random variable can have arbitrarily large sample values, there is always the possibility that this procedure will cause one or more reaction channels to fire so many times during tau that the population of some reactant species will be driven negative. Two recent papers have shown how that unacceptable occurrence can be avoided by replacing the Poisson random variables with binomial random variables, whose values are naturally bounded. This paper describes a modified Poisson tau-leaping procedure that also avoids negative populations, but is easier to implement than the binomial procedure. The new Poisson procedure also introduces a second control parameter, whose value essentially dials the procedure from the original Poisson tau-leaping at one extreme to the exact stochastic simulation algorithm at the other; therefore, the modified Poisson procedure will generally be more accurate than the original Poisson procedure.     

Extending the horizon: towards the efficient modeling of large biomolecular complexes in atomic detail

The application of evaluation of implicit solvent methods for the simulation of biomolecules is described. Detailed comparisons with explicit solvent are described for the modeling of peptide and proteins in continuum aqueous solvent. In addition, we are presenting new data on the simulation of DNA with implicit solvent and describe the development of a heterogeneous dielectric model for the simulation of integral membranes. The performance of implicit solvent simulations based on the GBMV generalized Born method is compared with explicit solvent simulations, and implications for the simulation of very large biomolecular complexes is discussed. We are anticipating that the work described herein will lead to new, efficient modeling tools that will allow the simulation of longer timescales and larger system sizes in order to meet current and future challenges by the experimental community.     

Modelling of reaction–diffusion processes: the theory of catalytic electrode processes at hemispheroidal ultramicroelectrodes

The transient chronoamperometric current for a catalytic electrode reaction (EC′ reaction) at a hemispheroidal Ultramicroelectrode is reported. This simple new approximate multidimensional analytical expression is valid for all values of time and all values of the rate constant. It is more accurate for all times and large K than for short times and small K. The analytical results are compared with recently presented digital simulation results (Galceran and co-workers, J. Electroanal. Chem 466 (1999) 15) for a disc electrode and are found to be in good agreement.     

Kinetic models for mechanoenzymes: structural aspects under large loads

A broad class of chemical kinetic model for mechanoenzymes is analyzed theoretically in order to uncover structural aspects of the underlying free-energy landscape that determine the behavior under large resisting and assisting loads, specifically the turnover rate or, for a translocatory motor protein, the mean velocity, say, V. A systematic graphical reduction algorithm is presented that provides explicit analytical expressions for mean occupation times in individual biomechanochemical states, for the splitting or backward/forward fractions, for the overall mean dwell time, and for the turnover rate. Application to the previously studied N-state sequential and (N alpha,N beta)-parallel-chain models provides explicit structural criteria (independent of the zero-load transition rates) that determine whether mid /V/ diverges to large values or, conversely, exhibits extrema and converges to a vanishing value as the externally imposed load grows. Closed-form analytical extensions accommodate side-chain and looped side-chain reaction sequences in the enzymatic cycle. A general divided-pathway model is analyzed in detail.     

Comparison of various implicit solvent models in molecular dynamics simulations of immunoglobulin G light chain dimer

The present study tests performance of different solvation models applied to molecular dynamics simulation of a large, dimeric protein molecule. Analytical Continuum Electrostatics (ACE) with two different parameter sets, older V98 and new V01, and Effective Energy Function (EEF) are employed in molecular dynamics simulation of immunoglobulin G (IgG) light chain dimer and variable domain of IgG light chain. Results are compared with explicit solvent and distance dependent dielectric constant (DDE) calculations. The overall analysis shows that the EEF method yields results comparable to explicit solvent simulations; however, the stability of simulations is lower. On the other hand, the ACE_V98 model does not seem to achieve the accuracy or stability expected in nanosecond timescale MD simulation for the studied systems. The ACE_V01 model greatly improves stability of the calculation; nonetheless, changes in radius of gyration and solvent accessible surface of the studied systems may indicate that the parameter set still needs to be improved if the method is supposed to be used for simulations of large, polymeric proteins. Additionally, electrostatic contribution to the solvation free energy calculated in the ACE model is compared with a numerical treatment of the dielectric continuum model. Wall clock time of all simulations is compared. It shows that EEF calculation is six times faster than corresponding ACE and 50 times faster than explicit solvent simulations.     

Surface hopping simulation of the vibrational relaxation of I2 in liquid xenon using the collective probabilities algorithm

A surface hopping simulation of the vibrational relaxation of highly excited I(2) in liquid xenon is presented. The simulation is performed by using the collective probabilities algorithm which assures the coincidence of the classical and quantum populations. The agreement between the surface hopping simulation results and the experimental measurements for the vibrational energy decay curves at different solvent densities and temperatures is shown to be good. The overlap of the decay curves when the time axis is linearly scaled is explained in terms of the perturbative theory for the rate constants. The contribution of each solvent atom to the change of the quantum populations of the solute molecule is used to analyze the mechanism of the relaxation process     

Explicit symplectic integrators of molecular dynamics algorithms for rigid-body molecules in the canonical, isobaric-isothermal, and related ensembles

The authors propose explicit symplectic integrators of molecular dynamics (MD) algorithms for rigid-body molecules in the canonical and isobaric-isothermal ensembles. They also present a symplectic algorithm in the constant normal pressure and lateral surface area ensemble and that combined with the Parrinello-Rahman algorithm. Employing the symplectic integrators for MD algorithms, there is a conserved quantity which is close to Hamiltonian. Therefore, they can perform a MD simulation more stably than by conventional nonsymplectic algorithms. They applied this algorithm to a TIP3P pure water system at 300 K and compared the time evolution of the Hamiltonian with those by the nonsymplectic algorithms. They found that the Hamiltonian was conserved well by the symplectic algorithm even for a time step of 4 fs. This time step is longer than typical values of 0.5-2 fs which are used by the conventional nonsymplectic algorithms.     

Mathematical treatment of concentration profiles and anodic current of amperometric multilayer enzyme electrodes

Formulae are presented for explicit solutions of partial differential equations describing the transient behaviour of concentration profiles and the anodic current in amperometric multilayer enzyme electrodes. The mathematical treatment is based on reaction/diffusion models with irreversible first-order catalytic reactions. Numerical results are presented to demonstrate the feasibility of the given formulae. The expression for the anodic current is used to fit kinetic parameters to measured current/time data.     

Molecular dynamics extended for fluctuating networks: Application to water

Molecular simulation models are increasingly important tools in efforts to understand the role that water plays in biochemical processes. However, existing models of water have limited capacity to deal with the characteristics of hydrogen bond networks. This article proposes a new fluctuating network (FN) algorithm as an extension of the standard molecular dynamics algorithm. The new algorithm allows for the simulation of a molecular system based on an underlying network, such as the hydrogen bond network in water. This algorithm distinguishes strong from weak network connections, applying a potential that best describes the specific connection behavior. We model liquid water with this new technique using a single‐site, isotropic, short‐range potential. We successfully reproduce liquid water's signature molecular spacing (as represented by the radial distribution function) and characterize its dynamic properties including the exponential hydrogen bond lifetime distribution, diffusion rate, and average hydrogen bonds per molecule. The FN algorithm allows exploration of the behavior of networked systems where explicit coordination limits are required. As such it could also be used to model covalent interactions, reaction dynamics, and applied to simulation of cellular networks. © 2012 Wiley Periodicals, Inc.     

Helix formation in a pentapeptide: experiment and force-field dependent dynamics

We used a combined approach of experiment and simulation to determine the helical population and folding pathway of a small helix forming blocked pentapeptide, Ac-(Ala)(5)-NH(2). Experimental structural characterization of this blocked peptide was carried out with far UV circular dichroism spectroscopy, FTIR, and NMR measurements. These measurements confirm the presence of the α-helical state in a buffer solution. Direct molecular dynamics and replica-exchange simulations of the pentapeptide were performed using several popular force fields with explicit solvent. The simulations yielded statistically reliable estimates of helix populations, melting curves, folding, and nucleation times. The distributions of conformer populations are used to measure folding cooperativity. Finally, a statistical analysis of the sample of helix-coil transition paths was performed. The details of the calculated helix populations, folding kinetics and pathways vary with the employed force field. Interestingly, the helix populations, folding, and unfolding times obtained from most of the studied force fields are in qualitative agreement with each other and with available experimental data, with the deviations corresponding to several kcal/mol in energy at 300 K. Most of the force fields also predict qualitatively similar transition paths, with unfolding initiated at the C-terminus. Accuracy of potential energy parameters, rather than conformational sampling may be the limiting factor in current molecular simulations.     

Adaptative finite element simulation of currents at microelectrodes to a guaranteed accuracy. Application to a simple model problem

In this series of papers we consider the general problem of numerical simulation of the currents at microelectrodes using an adaptive finite element approach. Microelectrodes typically consist of an electrode embedded (or recessed) in an insulating material. For all such electrodes, numerical simulation is made difficult by the presence of a boundary singularity at the electrode edge (where the electrode meets the insulator), manifested by the large increase in the current density at this point, often referred to as the ‘edge-effect’. Our approach to overcoming this problem involves the derivation of an a posteriori bound on the error in the numerical approximation for the current that can be used to drive an adaptive mesh-generation algorithm. This allows us to calculate the current to within a prescribed tolerance. We begin by demonstrating the power of the method for a simple model problem — an E reaction mechanism at a microdisc electrode — for which the analytical solution is known. In this paper we give the background to the problem, and show how an a posteriori error bound can be used to drive an adaptive mesh-generation algorithm. We then use the algorithm to solve our model problem and obtain very accurate results on comparatively coarse meshes in minimal computing time. We give the technical details of the background theory and the derivation of the error bound in the accompanying paper.