Molecular dynamics simulation of highly charged proteins: comparison of the particle-particle particle-mesh and reaction field methods for the calculation of electrostatic interactions

Computational studies on the molecular insights of aptamer induced poly(N-isopropylacrylamide)-graft-graphene oxide for on/off- switchable whole-cell cancer diagnostics

This work deals with first-principles and in silico studies of graphene oxide-based whole-cell selective aptamers for cancer diagnostics utilising a tunable-surface strategy. Herein, graphene oxide (GO) was constructed as a surface-based model with poly(N-isopropylacrylamide) (PNIPAM) covalently grafted as an "on/off"-switch in triggering interactions with the cancer-cell protein around its lower critical solution temperature. The atomic building blocks of the aptamer and the PNIPAM adsorbed onto the GO was investigated at the density functional theory (DFT) level. The presence of the monomer of PNIPAM stabilised the system's π-π interaction between GO and its nucleobases as confirmed by higher bandgap energy, satisfying the eigenvalues of the single-point energy observed rather than the nucleobase and the GO complex independently. The unaltered geometrical structures of the surface emphasise the physisorption type interaction between the nucleobase and the GO/NIPAM surface. The docking result for the aptamer and the protein, highlighted the behavior of the PNIPAM-graft-GO  is exhibiting globular and extended conformations, further supported by molecular dynamics (MD) simulations. These studies enabled a better understanding of the thermal responsive behavior of the polymer-enhanced GO complex for whole-cell protein interactions through computational methods.

All-Atom Continuous Constant pH Molecular Dynamics With Particle Mesh Ewald and Titratable Water

Development of a pH stat to properly control solution pH in biomolecular simulations has been a long-standing goal in the community. Toward this goal recent years have witnessed the emergence of the so-called constant pH molecular dynamics methods. However, the accuracy and generality of these methods have been hampered by the use of implicit-solvent models or truncation-based electrostatic schemes. Here we report the implementation of the particle mesh Ewald (PME) scheme into the all-atom continuous constant pH molecular dynamics (CpHMD) method, enabling CpHMD to be performed with a standard MD engine at a fractional added computational cost. We demonstrate the performance using pH replica-exchange CpHMD simulations with titratable water for a stringent test set of proteins, HP36, BBL, HEWL, and SNase. With the sampling time of 10 ns per replica, most pK_a's are converged, yielding the average absolute and root-mean-square deviations of 0.61 and 0.77, respectively, from experiment. Linear regression of the calculated vs experimental pK_a shifts gives a correlation coefficient of 0.79, a slope of 1, and an intercept near 0. Analysis reveals inadequate sampling of structure relaxation accompanying a protonation-state switch as a major source of the remaining errors, which are reduced as simulation prolongs. These data suggest PME-based CpHMD can be used as a general tool for pH-controlled simulations of macromolecular systems in various environments, enabling atomic insights into pH-dependent phenomena involving not only soluble proteins but also transmembrane proteins, nucleic acids, surfactants, and polysaccharides.

Scaling of Multimillion-Atom Biological Molecular Dynamics Simulation on a Petascale Supercomputer

A strategy is described for a fast all-atom molecular dynamics simulation of multimillion-atom biological systems on massively parallel supercomputers. The strategy is developed using benchmark systems of particular interest to bioenergy research, comprising models of cellulose and lignocellulosic biomass in an aqueous solution. The approach involves using the reaction field (RF) method for the computation of long-range electrostatic interactions, which permits efficient scaling on many thousands of cores. Although the range of applicability of the RF method for biomolecular systems remains to be demonstrated, for the benchmark systems the use of the RF produces molecular dipole moments, Kirkwood G factors, other structural properties, and mean-square fluctuations in excellent agreement with those obtained with the commonly used Particle Mesh Ewald method. With RF, three million- and five million-atom biological systems scale well up to ∼30k cores, producing ∼30 ns/day. Atomistic simulations of very large systems for time scales approaching the microsecond would, therefore, appear now to be within reach.

Toward the correction of effective electrostatic forces in explicit-solvent molecular dynamics simulations: restraints on solvent-generated electrostatic potential and solvent polarization

Despite considerable advances in computing power, atomistic simulations under nonperiodic boundary conditions, with Coulombic electrostatic interactions and in systems large enough to reduce finite-size associated errors in thermodynamic quantities to within the thermal energy, are still not affordable. As a result, periodic boundary conditions, systems of microscopic size and effective electrostatic interaction functions are frequently resorted to. Ensuing artifacts in thermodynamic quantities are nowadays routinely corrected a posteriori, but the underlying configurational sampling still descends from spurious forces. The present study addresses this problem through the introduction of on-the-fly corrections to the physical forces during an atomistic molecular dynamics simulation. Two different approaches are suggested, where the force corrections are derived from special potential energy terms. In the first approach, the solvent-generated electrostatic potential sampled at a given atom site is restrained to a target value involving corrections for electrostatic artifacts. In the second approach, the long-range regime of the solvent polarization around a given atom site is restrained to the Born polarization, i.e., the solvent polarization corresponding to the ideal situation of a macroscopic system under nonperiodic boundary conditions and governed by Coulombic electrostatic interactions. The restraints are applied to the explicit-water simulation of a hydrated sodium ion, and the effect of the restraints on the structural and energetic properties of the solvent is illustrated. Furthermore, by means of the calculation of the charging free energy of a hydrated sodium ion, it is shown how the electrostatic potential restraint translates into the on-the-fly consideration of the corresponding free-energy correction terms. It is discussed how the restraints can be generalized to situations involving several solute particles. Although the present study considers a very simple system only, it is an important step toward the on-the-fly elimination of finite-size and approximate-electrostatic artifacts during atomistic molecular dynamics simulations.

The zero-multipole summation method for estimating electrostatic interactions in molecular dynamics: analysis of the accuracy and application to liquid systems

In the preceding paper [I. Fukuda, J. Chem. Phys. 139, 174107 (2013)], the zero-multipole (ZM) summation method was proposed for efficiently evaluating the electrostatic Coulombic interactions of a classical point charge system. The summation takes a simple pairwise form, but prevents the electrically non-neutral multipole states that may artificially be generated by a simple cutoff truncation, which often causes large energetic noises and significant artifacts. The purpose of this paper is to judge the ability of the ZM method by investigating the accuracy, parameter dependencies, and stability in applications to liquid systems. To conduct this, first, the energy-functional error was divided into three terms and each term was analyzed by a theoretical error-bound estimation. This estimation gave us a clear basis of the discussions on the numerical investigations. It also gave a new viewpoint between the excess energy error and the damping effect by the damping parameter. Second, with the aid of these analyses, the ZM method was evaluated based on molecular dynamics (MD) simulations of two fundamental liquid systems, a molten sodium-chlorine ion system and a pure water molecule system. In the ion system, the energy accuracy, compared with the Ewald summation, was better for a larger value of multipole moment l currently induced until l ≲ 3 on average. This accuracy improvement with increasing l is due to the enhancement of the excess-energy accuracy. However, this improvement is wholly effective in the total accuracy if the theoretical moment l is smaller than or equal to a system intrinsic moment L. The simulation results thus indicate L ∼ 3 in this system, and we observed less accuracy in l = 4. We demonstrated the origins of parameter dependencies appearing in the crossing behavior and the oscillations of the energy error curves. With raising the moment l we observed, smaller values of the damping parameter provided more accurate results and smoother behaviors with respect to cutoff length were obtained. These features can be explained, on the basis of the theoretical error analyses, such that the excess energy accuracy is improved with increasing l and that the total accuracy improvement within l ⩽ L is facilitated by a small damping parameter. Although the accuracy was fundamentally similar to the ion system, the bulk water system exhibited distinguishable quantitative behaviors. A smaller damping parameter was effective in all the practical cutoff distance, and this fact can be interpreted by the reduction of the excess subset. A lower moment was advantageous in the energy accuracy, where l = 1 was slightly superior to l = 2 in this system. However, the method with l = 2 (viz., the zero-quadrupole sum) gave accurate results for the radial distribution function. We confirmed the stability in the numerical integration for MD simulations employing the ZM scheme. This result is supported by the sufficient smoothness of the energy function. Along with the smoothness, the pairwise feature and the allowance of the atom-based cutoff mode on the energy formula lead to the exact zero total-force, ensuring the total-momentum conservations for typical MD equations of motion.

Zero-multipole summation method for efficiently estimating electrostatic interactions in molecular system

The zero-multipole summation method has been developed to efficiently evaluate the electrostatic Coulombic interactions of a point charge system. This summation prevents the electrically non-neutral multipole states that may artificially be generated by a simple cutoff truncation, which often causes large amounts of energetic noise and significant artifacts. The resulting energy function is represented by a constant term plus a simple pairwise summation, using a damped or undamped Coulombic pair potential function along with a polynomial of the distance between each particle pair. Thus, the implementation is straightforward and enables facile applications to high-performance computations. Any higher-order multipole moment can be taken into account in the neutrality principle, and it only affects the degree and coefficients of the polynomial and the constant term. The lowest and second moments correspond respectively to the Wolf zero-charge scheme and the zero-dipole summation scheme, which was previously proposed. Relationships with other non-Ewald methods are discussed, to validate the current method in their contexts. Good numerical efficiencies were easily obtained in the evaluation of Madelung constants of sodium chloride and cesium chloride crystals.

Dependence of protein crystal stability on residue charge states and ion content of crystal solvent

Protein crystallization is frequently induced by the addition of various precipitants, which directly affect protein solubility. In addition to organic cosolvents and long-chain polymers, salts belong to the most widely used precipitants in protein crystallography. However, despite such widespread usage, their mode of action at the atomistic level is still largely unknown. Here, we perform extensive molecular dynamics simulations of the villin headpiece crystal unit cell to examine its stability at different concentrations of sodium sulfate. We show that the inclusion of ions in crystal solvent at high concentration can prevent large rearrangements of the asymmetric units and a loss of symmetry of the unit cell without significantly affecting protein dynamics. Of importance, a similar result can be achieved by neutralizing several specific charged residues suggesting that they may play an active role in crystal destabilization due to unfavorable electrostatic interactions. Our results provide a microscopic picture behind salt-induced stabilization of a protein crystal and further suggest that adequate modeling of realistic crystallization conditions may be necessary for successful molecular dynamics simulations of protein crystals.

Charge-leveling and proper treatment of long-range electrostatics in all-atom molecular dynamics at constant pH

Recent development of constant pH molecular dynamics (CpHMD) methods has offered promise for adding pH-stat in molecular dynamics simulations. However, until now the working pH molecular dynamics (pHMD) implementations are dependent in part or whole on implicit-solvent models. Here we show that proper treatment of long-range electrostatics and maintaining charge neutrality of the system are critical for extending the continuous pHMD framework to the all-atom representation. The former is achieved here by adding forces to titration coordinates due to long-range electrostatics based on the generalized reaction field method, while the latter is made possible by a charge-leveling technique that couples proton titration with simultaneous ionization or neutralization of a co-ion in solution. We test the new method using the pH-replica-exchange CpHMD simulations of a series of aliphatic dicarboxylic acids with varying carbon chain length. The average absolute deviation from the experimental pK(a) values is merely 0.18 units. The results show that accounting for the forces due to extended electrostatics removes the large random noise in propagating titration coordinates, while maintaining charge neutrality of the system improves the accuracy in the calculated electrostatic interaction between ionizable sites. Thus, we believe that the way is paved for realizing pH-controlled all-atom molecular dynamics in the near future.

Effect of the Reaction Field on Molecular Forces and Torques Revealed by an Image-Charge Solvation Model

We recently developed the Image-Charge Solvation Model (ICSM), which is an explicit/implicit hybrid model to accurately account for long-range electrostatic forces in molecular dynamics simulations [Lin et al., J. Chem. Phys., 131, 154103, 2009]. The ICSM has a productive spherical volume within the simulation cell for which key physical properties of bulk water are reproduced, such as density, radial distribution function, diffusion constants and dielectric properties. Although the reaction field (RF) is essential, it typically accounts for less than 2% of the total electrostatic force on a water molecule. This observation motivates investigating further the role of the RF within the ICSM. In this report we focus on distributions of forces and torques on water molecules as a function of distance from the origin and make extensive tests over a range of model parameters where Coulomb forces are decomposed into direct interactions from waters modeled explicitly and the RF. Molecular torques due to the RF typically account for 20% of the total torque, revealing why the RF plays an important role in the dielectric properties of simulated water. Moreover, it becomes clear that the buffer layer in the ICSM is essential to mitigate artifacts caused by the discontinuous change in dielectric constants at the explicit/implicit interface.

Non-Ewald methods: theory and applications to molecular systems

Several non-Ewald methods for calculating electrostatic interactions have recently been developed, such as the Wolf method, the reaction field method, the pre-averaging method, and the zero-dipole summation method, for molecular dynamics simulations of various physical systems, including biomolecular systems. We review the theories of these approaches and their potential applications to molecular simulations, and discuss their relationships.

Assessment of Common Simulation Protocols for Simulations of Nanopores, Membrane Proteins, and Channels

Molecular dynamics (MD) simulation has become a common technique to study biological systems. Transport of small molecules through carbon nanotubes and membrane proteins has been an intensely studied topic, and MD simulations have been able to provide valuable predictions, many of which have later been experimentally proven. Simulations of such systems pose challenges, and unexpected problems in commonly used protocols and methods have been found in the past few years. The two main reasons why some were not found before are that most of these newly discovered errors do not lead to unstable simulations. Furthermore, some of them manifest themselves only after relatively long simulation times. We assessed the reliability of the most common simulations protocols by MD and stochastic dynamics (SD) or Langevin dynamics, simulations of an alpha hemolysin nanochannel embedded in a palmitoyloleoylphosphatidylcholine (POPC) lipid bilayer. Our findings are that (a) reaction field electrostatics should not be used in simulations of such systems, (b) local thermostats should be preferred over global ones since the latter may lead to an unphysical temperature distribution,

A simple polarizable continuum solvation model for electrolyte solutions

We propose a Debye-Hückel-like screening model (DESMO) that generalizes the familiar conductor-like screening model (COSMO) to solvents with non-zero ionic strength and furthermore provides a numerical generalization of the Debye-Hückel model that is applicable to non-spherical solute cavities. The numerical implementation of DESMO is based upon the switching/Gaussian (SWIG) method for smooth cavity discretization, which we have recently introduced in the context of polarizable continuum models (PCMs). This approach guarantees that the potential energy is a smooth function of the solute geometry and analytic gradients for DESMO are reported here. The SWIG formalism also facilitates analytic implementation of two other PCMs that are based on a screened Coulomb potential: the "integral equation formalism" (IEF-PCM) and the "surface and simulation of volume polarization for electrostatics" [SS(V)PE] method. Fully analytic implementations of these screened PCMs are reported here for the first time. Numerical results, for model systems where an exact solution of the linearized Poisson-Boltzmann equation is available, demonstrate that these screened PCMs are highly accurate. In realistic test cases, they are as accurate as the best available three-dimensional finite-difference methods. In polar solvents, DESMO is nearly as accurate as more sophisticated screened PCMs, but is significantly simpler and more efficient.

Prediction of a new class of RNA recognition motif

The observation that activation domains (AD) of procarboxypeptidases are rather long compared to the pro-regions of other zymogens raises the possibility that they could play additional roles apart from precluding enzymatic activity within the proenzyme and helping in its folding process. In the present work, we compared the overall pro-domain tertiary structure with several proteins belonging to the same fold in the structural classification of proteins (SCOP) database by using structure and sequence comparisons. The best score obtained was between the activation domain of human procarboxypeptidase A4 (ADA4h) and the human U1A protein from the U1 snRNP. Structural alignment revealed the existence of RNP1- and RNP2-related sequences in ADA4h. After modeling ADA4h on U1A, the new structure was used to extract a new sequence pattern characteristic for important residues at key positions. The new sequence pattern allowed scanning protein sequences to predict the RNA-binding function for 32 sequences undetected by PFAM. Unspecific RNA electrophoretic mobility shift assays experimentally supported the prediction that ADA4h binds an RNA motif similar to the U1A binding-motif of stem-loop II of U1 small nuclear RNA. The experiments carried out with ADA4h in the present work suggest the sharing of a common ancestor with other RNA recognition motifs. However, the fact that key residues preventing activity within the proenzyme are also key residues for RNA binding might have induced the activation domains of procarboxypeptidases to evolve from the canonical RNP1 and RNP2 sequences.

Molecular dynamics simulations of a reversibly folding beta-heptapeptide in methanol: influence of the treatment of long-range electrostatic interactions

Eight 100-ns molecular dynamics simulations of a beta-heptapeptide in methanol at 340 K (within cubic periodic computational boxes of about 6-nm edge) are reported and compared. These simulations were performed with three different charge-state combinations at the peptide termini, one of them with or without a neutralizing chloride counterion, and using either the lattice-sum (LS) or reaction-field (RF) scheme to handle electrostatic interactions. The choice of the electrostatic scheme has essentially no influence on the folding-unfolding equilibrium when the peptide termini are uncharged and only a small influence when the peptide is positively charged at its N-terminus (with or without inclusion of a neutralizing chloride counterion). However, when the peptide is zwitterionic, the LS scheme leads to preferential sampling of the high-dipole folded helical state, whereas the RF scheme leads to preferential sampling of a low-dipole unfolded salt-bridged state. A continuum electrostatics analysis based on the sampled configurations (zwitterionic case) suggests that the LS scheme stabilizes the helical state through artificial periodicity, but that the magnitude of this perturbation is essentially negligible (compared to the thermal energy) for the large box size and relatively polar solvent considered. The results thus provide clear evidence (continuum electrostatics analysis) for the absence of LS artifacts and some indications (still not definitive because of the limited sampling of the folding-unfolding transition) for the presence of RF artifacts in this specific system.

Influence of NH-Sgamma bonding interactions on the structure and dynamics of metallothioneins

Mammalian metallothioneins ([Formula: see text]) show a clustered arrangement of the metal ions and a nonregular protein structure. The solution structures of Cd(3)-thiolate cluster containing beta-domain of mouse beta-MT-1 and rat beta-MT-2 show high structural similarities, but widely differing structure dynamics. Molecular dynamics simulations revealed a substantially increased number of NH-Sgamma hydrogen bonds in beta-MT-2, features likely responsible for the increased stability of the Cd(3)-thiolate cluster and the enfolding protein domain. Alterations in the NH-Sgamma hydrogen-bonding network may provide a rationale for the differences in dynamic properties encountered in the beta-domains of MT-1, -2, and -3 isoforms, believed to be essential for their different biological function.

Coarse-graining in interaction space: a systematic approach for replacing long-range electrostatics with short-range potentials

A short-range effective potential for long-range electrostatic interactions in homogeneously disordered condensed phase systems has been determined with a novel approach to coarse-graining in interaction space. As opposed to coarse-graining the system resolution, this approach "coarsens" the system's interactions by mapping multiple configurations of an accurate long-range atomistic potential onto a more efficient, short-range effective potential with a force-matching (FM) method. Developing an empirical potential in this manner is fundamentally different from existing strategies because it utilizes condensed-phase (as opposed to gas-phase) atomistic interactions to determine general pair potentials defined on distance meshes (as opposed to fitting predetermined functional forms). The resulting short-range ( approximately 10 A) effective potential reproduces structural, dynamical, and many thermodynamic properties of liquid water, ions in water, and hydrophobes in water, with unprecedented accuracy. The effective potential is also shown to be transferable to a nonaqueous molten salt system. With continued development, such effective potentials may provide an accurate and highly efficient alternative to Ewald-based long-range electrostatics methods.

Combining the lattice-sum and reaction-field approaches for evaluating long-range electrostatic interactions in molecular simulations

A new scheme, the lattice-sum-emulated reaction-field (LSERF) method, is presented that combines the lattice-sum (LS) and reaction-field (RF) approaches for evaluating electrostatic interactions in molecular simulations. More precisely, the LSERF scheme emulates a RF calculation (based on an atomic cutoff) via the LS machinery. This is achieved by changing the form of the electrostatic interactions in a standard LS calculation (Coulombic) to the form corresponding to RF electrostatics (Coulombic plus quadratic reaction-field correction term, truncated at the cutoff distance). It is shown (both analytically and numerically) that in the limit of infinite reciprocal-space accuracy, (i) the LSERF scheme with a finite reaction-field cutoff and a given reaction-field permittivity is identical to the RF scheme with the same parameters (and an atomic cutoff), and (ii) the LSERF scheme is identical to the LS scheme in the limit of an infinite reaction-field cutoff, irrespective of the reaction-field permittivity. This new scheme offers two key advantages: (i) from a conceptual point of view, it shows that there is a continuity between the RF and LS schemes and unifies them into a common framework; (ii) from a practical point of view, it allows us to perform RF calculations with arbitrarily large reaction-field cutoff distances for the same computational costs as a corresponding LS calculation. The optimal choice for the cutoff will be the one that achieves the best compromise between artifacts arising from the dielectric heterogeneity of the system (short cutoff) and its artificial periodicity (long cutoff). The implementation of the LSERF method is extremely easy, requiring only very limited modifications of any standard LS code. For practical applications to biomolecular systems, the use of the LSERF scheme with large reaction-field cutoff distances is expected to represent a significant improvement over the current RF simulations involving comparatively much shorter cutoffs.

Positive selection in MAOA gene is human exclusive: determination of the putative amino acid change selected in the human lineage

Monoamine oxidase A (MAOA) is the X-linked gene responsible for deamination and subsequent degradation of several neurotransmitters and other amines. Among other activities, the gene has been shown to play a role in locomotion, circadian rhythm, and pain sensitivity and to have a critical influence on behavior and cognition. Previous studies have reported a non-neutral evolution of the gene attributable to positive selection in the human lineage. To determine whether this selection was human-exclusive or shared with other species, we performed a population genetic analysis of the pattern of nucleotide variation in non-human species, including bonobo, chimpanzee, gorilla, and orangutan. Footprints of positive selection were absent in all analyzed species, suggesting that positive selection has been recent and unique to humans. To determine which human-unique genetic changes could have been responsible for this differential evolution, the coding region of the gene was compared between human, chimpanzee, and gorilla. Only one human exclusive non-conservative change is present in the gene: Glu151Lys. This human substitution affects protein dimerization according to a three-dimensional structural model that predicts a non-negligible functional shift. This is the only candidate position at present to have been selected to fixation in humans during an episode of positive selection. Divergence analysis among species has shown that, even under positive selection in the human lineage, the MAOA gene did not experience accelerated evolution in any of the analyzed lineages, and that tools such as K(a)/ K(s) would not have detected the selective history of the gene.