RPYFMM: Parallel Adaptive Fast Multipole Method for Rotne-Prager-Yamakawa Tensor in Biomolecular Hydrodynamics Simulations.

RPYFMM is a software package for the efficient evaluation of the potential field governed by the Rotne-Prager-Yamakawa (RPY) tensor interactions in biomolecular hydrodynamics simulations. In our algorithm, the RPY tensor is decomposed as a linear combination of four Laplace interactions, each of which is evaluated using the adaptive fast multipole method (FMM) [1] where the exponential expansions are applied to diagonalize the multipole-to-local translation operators. RPYFMM offers a unified execution on both shared and distributed memory computers by leveraging the DASHMM library [2, 3]. Preliminary numerical results show that the interactions for a molecular system of 15 million particles (beads) can be computed within one second on a Cray XC30 cluster using 12, 288 cores, while achieving approximately 54% strong-scaling efficiency.

Enzyme localization, crowding, and buffers collectively modulate diffusion-influenced signal transduction: Insights from continuum diffusion modeling.

Biochemical reaction networks consisting of coupled enzymes connect substrate signaling events with biological function. Substrates involved in these reactions can be strongly influenced by diffusion "barriers" arising from impenetrable cellular structures and macromolecules, as well as interactions with biomolecules, especially within crowded environments. For diffusion-influenced reactions, the spatial organization of diffusion barriers arising from intracellular structures, non-specific crowders, and specific-binders (buffers) strongly controls the temporal and spatial reaction kinetics. In this study, we use two prototypical biochemical reactions, a Goodwin oscillator, and a reaction with a periodic source/sink term to examine how a diffusion barrier that partitions substrates controls reaction behavior. Namely, we examine how conditions representative of a densely packed cytosol, including reduced accessible volume fraction, non-specific interactions, and buffers, impede diffusion over nanometer length-scales. We find that diffusion barriers can modulate the frequencies and amplitudes of coupled diffusion-influenced reaction networks, as well as give rise to "compartments" of decoupled reactant populations. These effects appear to be intensified in the presence of buffers localized to the diffusion barrier. These findings have strong implications for the role of the cellular environment in tuning the dynamics of signaling pathways.

Predicting the influence of long-range molecular interactions on macroscopic-scale diffusion by homogenization of the Smoluchowski equation.

The macroscopic diffusion constant for a charged diffuser is in part dependent on (1) the volume excluded by solute "obstacles" and (2) long-range interactions between those obstacles and the diffuser. Increasing excluded volume reduces transport of the diffuser, while long-range interactions can either increase or decrease diffusivity, depending on the nature of the potential. We previously demonstrated [P. M. Kekenes-Huskey et al., Biophys. J. 105, 2130 (2013)] using homogenization theory that the configuration of molecular-scale obstacles can both hinder diffusion and induce diffusional anisotropy for small ions. As the density of molecular obstacles increases, van der Waals (vdW) and electrostatic interactions between obstacle and a diffuser become significant and can strongly influence the latter's diffusivity, which was neglected in our original model. Here, we extend this methodology to include a fixed (time-independent) potential of mean force, through homogenization of the Smoluchowski equation. We consider the diffusion of ions in crowded, hydrophilic environments at physiological ionic strengths and find that electrostatic and vdW interactions can enhance or depress effective diffusion rates for attractive or repulsive forces, respectively. Additionally, we show that the observed diffusion rate may be reduced independent of non-specific electrostatic and vdW interactions by treating obstacles that exhibit specific binding interactions as "buffers" that absorb free diffusers. Finally, we demonstrate that effective diffusion rates are sensitive to distribution of surface charge on a globular protein, Troponin C, suggesting that the use of molecular structures with atomistic-scale resolution can account for electrostatic influences on substrate transport. This approach offers new insight into the influence of molecular-scale, long-range interactions on transport of charged species, particularly for diffusion-influenced signaling events occurring in crowded cellular environments.

Dewetting-controlled binding of ligands to hydrophobic pockets.

We report on a combined atomistic molecular dynamics simulation and implicit solvent analysis of a generic hydrophobic pocket-ligand (host-guest) system. The approaching ligand induces complex wetting-dewetting transitions in the weakly solvated pocket. The transitions lead to bimodal solvent fluctuations which govern magnitude and range of the pocket-ligand attraction. A recently developed implicit water model, based on the minimization of a geometric functional, captures the sensitive aqueous interface response to the concave-convex pocket-ligand configuration semiquantitatively.

Computer-aided molecular design.

Theoretical chemistry, as implemented on fast computers, is beginning to yield accurate predictions of the thermodynamic and kinetic properties of large molecular assemblies. In addition to providing detailed insights into the origins of molecular activity, theoretical calculations can be used to design new molecules with specific properties. This article describes two types of calculations that show special promise as design tools, the thermodynamic cycle-perturbation method and the Brownian reactive dynamics method. These methods can be applied to calculate equilibrium and rate constants that describe many aspects of molecular recognition, stability, and reactivity.

Phenylalanine transfer RNA: molecular dynamics simulation.

Yeast phenylalanine transfer RNA was subjected to a 12-picosecond molecular dynamics simulation. The principal features of the x-ray crystallographic analysis are reproduced, and the amplitudes of atomic displacements appear to be determined by the degree of exposure of the atoms. An analysis of the hydrogen bonds shows a correlation between the average length of a bond and the fluctuation in that length and reveals a rocking motion of bases in Watson-Crick guanine X cytosine base pairs. The in-plane motions of the bases are generally of larger amplitude than the out-of-plane motions, and there are correlations in the motions of adjacent bases.

Open "back door" in a molecular dynamics simulation of acetylcholinesterase.

The enzyme acetylcholinesterase generates a strong electrostatic field that can attract the cationic substrate acetylcholine to the active site. However, the long and narrow active site gorge seems inconsistent with the enzyme's high catalytic rate. A molecular dynamics simulation of acetylcholinesterase in water reveals the transient opening of a short channel, large enough to pass a water molecule, through a thin wall of the active site near tryptophan-84. This simulation suggests that substrate, products, or solvent could move through this "back door," in addition to the entrance revealed by the crystallographic structure. Electrostatic calculations show a strong field at the back door, oriented to attract the substrate and the reaction product choline and to repel the other reaction product, acetate. Analysis of the open back door conformation suggests a mutation that could seal the back door and thus test the hypothesis that thermal motion of this enzyme may open multiple routes of access to its active site.

Acetylcholinesterase: mechanisms of covalent inhibition of H447I mutant determined by computational analyses.

The reaction mechanisms of two inhibitor TFK(+) and TFK(0) binding to H447I mutant mouse acetylcholinesterase (mAChE) have been investigated by using a combined ab initio quantum mechanical/molecular mechanical (QM/MM) approach and classical molecular dynamics (MD) simulations. TFK(+) binding to the H447I mutant may proceed with a different reaction mechanism from the wild-type. A water molecule takes over the role of His447 and participates in the bond breaking and forming as a "charge relayer". Unlike in the wild-type mAChE case, Glu334, a conserved residue from the catalytic triad, acts as a catalytic base in the reaction. The calculated energy barrier for this reaction is about 8kcal/mol. These predictions await experimental verification. In the case of the neutral ligand TFK(0), however, multiple MD simulations on the TFK(0)/H447I complex reveal that none of the water molecules can be retained in the active site as a "catalytic" water. Taken together our computational studies confirm that TFK(0) is almost inactive in the H447I mutant, and also provide detailed mechanistic insights into the experimental observations.

Theory of biomolecular recognition.

Specific, noncovalent binding of biomolecules can only be understood by considering structural, thermodynamic, and kinetic issues. The theoretical foundations for such analyses have been clarified in the past year. Computational techniques for both particle-based and continuum models continue to improve and to yield useful insights into an ever wider range of biomolecular systems.

Molecular modeling methods. Basic techniques and challenging problems.

An overview is presented of computer modeling and simulation methods that play an increasing role in drug design: quantum chemical methods, molecular mechanics, molecular dynamics and Brownian dynamics. The application of molecular dynamics for the prediction of thermodynamic properties like free energy differences and binding constants is discussed. The Brownian dynamics method is presented in connection with the calculation of effective electrostatic forces using the Poisson-Boltzmann equation, which allows one to sample ligand-binding geometries and to predict the kinetics of diffusion-limited enzyme reactions. New techniques that have recently been extensively developed, such as the global energy minimization and quantum-classical dynamics methods, are also introduced. The molecular modeling methods are illustrated with selected examples.

Binding of an antiviral agent to a sensitive and a resistant human rhinovirus. Computer simulation studies with sampling of amino acid side-chain conformation. I. Mapping the rotamers of residue 188 of viral protein 1.

The mutation of valine 188 to leucine in the viral protein 1 of human rhinovirus 14 renders the virus resistant to certain antiviral compounds. Thermodynamic-cycle perturbation theory provides a means of calculating the difference in the binding free energies of an antiviral compound to the wild-type virus and to the mutant virus. In calculating the relevant free-energy differences in molecular dynamics simulations, it is important to sample the multiple rotational isomers of residue 188 correctly. In general, these rotamers will not be fully sampled during a single molecular dynamics simulation. However, the contributions of all the rotamers to the free-energy differences associated with mutation of residue 188 may be considered explicitly once they have been identified and their relative free energies determined. Therefore, we describe here the mapping of the rotamers of residue 188 by steric-bump search and energy minimization techniques, and by the computation of potentials of mean force (p.m.f.s.) using umbrella sampling. The usefulness, validity and efficiency of these methods of examining rotameric states is discussed. Adiabatic mapping by energy minimization was found to be unreliable for this residue due to the small magnitude of its interactions with the surrounding protein atoms. Ambiguities in the adiabatic maps were resolved by computing p.m.f.s. The p.m.f. for valine 188 in the unliganded wild-type virus shows a minimum corresponding to the crystallographically observed conformation of valine 188. The p.m.f.s. for valine 188 in the liganded virus and for leucine 188 in the unliganded mutant virus suggest that the experimentally observed conformations may be interpreted as averages of a number of conformations corresponding to those at the minima in the p.m.f.s. The calculations suggest also that the conformation of leucine 188 may change when the ligand binds. The use of the calculated p.m.f.s. to compute the difference in the free energy of binding of an antiviral compound to the wild-type and mutant rhinoviruses is described in the accompanying article.

Binding of an antiviral agent to a sensitive and a resistant human rhinovirus. Computer simulation studies with sampling of amino acid side-chain conformations. II. Calculation of free-energy differences by thermodynamic integration.

Thermodynamic-cycle perturbation theory and molecular dynamics simulations were used to calculate the difference in the free energy of binding of the antiviral compound WIN53338 to the wild-type human rhinovirus 14 and to a drug-resistant mutant of the virus in which valine 188 of the viral protein 1 is mutated to leucine. Because of the difficulty of achieving adequate sampling of all of the rotational isomers of amino acid side-chains in molecular dynamics simulations, an explicit treatment of the effects of the existence of multiple rotational isomers of residue 188 on the calculated free energies was used. The rotamers of residue 188 were first mapped by steric and energetic techniques as described in the accompanying article. Thermodynamic integration was then carried out during simulations of the virus, both with and without the antiviral compound bound, by mutating residue 188 while restraining its side-chain to one conformation. The contributions of the other rotamers of residue 188 to the free-energy changes for this mutation were then added to those calculated by thermodynamic integration as correction factors. Binding of WIN53338 to the wild-type virus was calculated to be favored over binding to the mutant virus by 1.7(+/- 3.0) kcal/mol. This is consistent with experimental data which, if differences in activity are assumed to be due to differences in binding, indicate that the binding affinity of WIN53338 for the wild-type virus is at least 0.15 to 1.7 kcal/mol greater than for the mutant virus. Thermodynamic integration was also performed in the conventional manner without restraints and was found to give less accurate results.

Electrostatic contributions to the stability of halophilic proteins.

Examination of the first crystal structures of proteins from a halophilic organism suggests that an abundance of acidic residues distributed over the protein surface is a key determinant of adaptation to high-salt conditions. Although one extant theory suggests that acidic residues are favored because of their superior water-binding capacity, it is clear that extensive repulsive electrostatic interactions will also be present in such proteins at physiological pH. To investigate the magnitude and importance of such electrostatic interactions, we conducted a theoretical analysis of their contributions to the salt and pH-dependence of stability of two halophilic proteins. Our approach centers on use of the Poisson-Boltzmann equation of classical electrostatics, applied at an atomic level of detail to crystal structures of the proteins. We first show that in using the method, it is important to account for the fact that the dielectric constant of water decreases at high salt concentrations, in order to reproduce experimental changes in pKa values of small acids and bases. We then conduct a comparison of salt and pH effects on the stability of 2Fe-2S ferredoxins from the halophile Haloarcula marismortui and the non-halophile anabaena. In both proteins, substantial upward shifts in pKa accompany protein folding, though shifts are considerably larger, on average, in the halophile. Upward shifts for basic residues occur because of favorable salt-bridge interactions, whilst upward shifts for acidic residues result from unfavorable electrostatic interactions with other acidic groups. Our calculations suggest that at pH 7 the stability of the halophilic protein is decreased by 18.2 kcal/mol on lowering the salt concentration from 5 M to 100 mM, a result that is in line with the fact that halophilic proteins generally unfold at low salt concentrations. For comparison, the non-halophilic ferredoxin is calculated to be destabilized by only 5.1 kcal/mol over the same range. Analysis of the pH stability curve suggests that lowering the pH should increase the intrinsic stability of the halophilic protein at low salt concentrations, although in practice this is not observed because of aggregation effects. We report the results of a similar analysis carried out on the tetrameric malate dehydrogenase from H. marismortui. In this case, we investigated the salt and pH dependence of the various monomer-monomer interactions present in the tetramer. All monomer-monomer interactions are found to make substantial contributions to the salt-dependence of stability of the tetramer. Excellent agreement is obtained between our calculated results for the stability of the tetramer and experimental results. In particular, the finding that at 4 M NaCl, the tetramer is stable only between pH 4.8 and 10 is accurately reproduced. Taken together, our results suggest that repulsive electrostatic interactions between acidic residues are a major factor in the destabilization of halophilic proteins in low-salt conditions, and that these interactions remain destabilizing even at high salt concentrations. As a consequence, the role of acidic residues in halophilic proteins may be more to prevent aggregation than to make a positive contribution to intrinsic protein stability.

Prediction of pH-dependent properties of proteins.

We describe what may be the most accurate approach currently available for the calculation of the pKas of ionizable groups in proteins. The accuracy is assessed by comparison of computed pKas with 60 measured pKas in a total of seven proteins. The overall root-mean-square error is 0.89 pKa units. Linear regression analysis of computed versus measured pKas yields a slope of 0.95, y-intercept of -0.02 and a correlation coefficient of 0.96. The proposed approach also picks out many of the shifted pKas of groups in enzyme active sites and special salt bridges. However, it does yield several over-shifted pKas and tends to underestimate pKa shifts which result from desolvation effects. We examine the ability of the new approach to reproduce the dependence of protein stability upon pH, using the ionization polynomial formalism. Overall features of the stability curves are reproduced, but the quantitative agreement is not particularly good. The reasons for the disagreement may have to do both with insufficient accuracy in the theory and with uncertainty in the nature of the unfolded state of proteins. The methodology described here is based upon finite difference solutions of the Poisson-Boltzmann equation. Its success depends upon the use of the rather high protein dielectric constant of 20. However, theoretical considerations and the fact that pKa shifts which result from desolvation are underestimated here imply that the dielectric constant of the protein interior actually is lower than 20. We suggest that the high protein dielectric constant improves the overall agreement with experiment because it accounts approximately for phenomena which tend to mitigate pKa shifts and which are not specifically included in the model. These include conformational relaxation and specific ion-binding. Future models based upon a low protein dielectric constant and treating such phenomena explicitly might yield improved agreement with experiment.

Computer simulations of actin polymerization can explain the barbed-pointed end asymmetry.

Computer simulations of actin polymerization were performed to investigate the role of electrostatic interactions in determining polymerization rates. Atomically detailed models of actin monomers and filaments were used in conjunction with a Brownian dynamics method. The simulations were able to reproduce the measured barbed end association rates over a range of ionic strengths and predicted a slower growing pointed end, in agreement with experiment. Similar simulations neglecting electrostatic interactions indicate that configurational and entropic factors may actually favor polymerization at the pointed end, but electrostatic interactions remove this trend. This result would indicate that polymerization at the pointed end is not only limited by diffusion, but faces electrostatic forces that oppose binding. The binding of the actin depolymerizing factor (ADF) and G-actin complex to the end of a filament was also simulated. In this case, electrostatic steering effects lead to an increase in the simulated association rate. Together, the results indicate that simulations provide a realistic description of both polymerization and the binding of more complex structures to actin filaments.

Combined conformational search and finite-difference Poisson-Boltzmann approach for flexible docking. Application to an operator mutation in the lambda repressor-operator complex.

The N-terminal domain of the phage lambda repressor binds as a dimer to its palindromic DNA operator sequence. In addition to a helix-turn-helix DNA recognition motif, the first six amino acids of the phage lambda repressor form a flexible peptide segment which wraps around DNA. Site-directed mutagenesis studies have shown that amino acid replacements or partial removal of the arm structure, or changes in the DNA sequence contacting the N-terminal arm, can lower the repressor-operator binding affinity by several orders of magnitude. The finite-difference Poisson-Boltzmann approach in combination with a conformational search procedure was used to study energetic contributions of the lambda arm to repressor-operator recognition based on the high resolution X-ray structure. It allows for the local relaxation of the structure upon changing the DNA sequence in the lambda arm binding region. A simplified potential energy function including torsional, truncated Lennard-Jones and approximate electrostatic terms is used in the initial step to screen out energetically unfavorable structures. The electrostatic energy of selected conformations is subsequently calculated more accurately using the finite-difference Poisson-Boltzmann approach. The method was applied to study the effect of a C-->T mutation at position 6 of the consensus half-site of the operator. This base-pair contacts Lys4 which is part of the arm segment. Keeping only the Lys4 side-chain mobile and with the wild-type DNA operator sequence, several conformations close to the X-ray structure were identified as those with lowest energy. In the case of the DNA mutation, lowest energy conformations differed significantly from those selected for the wild-type sequence. These initial calculations indicate that the approach might be a useful tool to estimate conformational and energetic effects upon mutagenesis of protein-DNA complexes.

Molecular dynamics simulations of the hyperthermophilic protein sac7d from Sulfolobus acidocaldarius: contribution of salt bridges to thermostability.

Hyperthermophilic proteins often possess an increased number of surface salt bridges compared with their mesophilic homologues. However, salt bridges are generally thought to be of minor importance in protein stability at room temperature. In an effort to understand why this may no longer be true at elevated temperatures, we performed molecular dynamics simulations of the hyperthermophilic protein Sac7d at 300 K, 360 K, and 550 K. The three trajectories are stable on the nanosecond timescale, as evidenced by the analysis of several time-resolved properties. The simulations at 300 K and (to a lesser extent) 360 K are also compatible with nuclear Overhauser effect-derived distances. Raising the temperature from 300 K to 360 K results in a less favourable protein-solvent interaction energy, and a more favourable intraprotein interaction energy. Both effects are almost exclusively electrostatic in nature and dominated by contributions due to charged side-chains. The reduced solvation is due to a loss of spatial and orientational structure of water around charged side-chains, which is a consequence of the increased thermal motion in the solvent. The favourable change in the intraprotein Coulombic interaction energy is essentially due to the tightening of salt bridges. Assuming that charged side-chains are on average more distant from one another in the unfolded state than in the folded state, it follows that salt bridges may contribute to protein stability at elevated temperatures because (i) the solvation free energy of charged side-chains is more adversely affected in the unfolded state than in the folded state by an increase in temperature, and (ii) due to the tightening of salt bridges, unfolding implies a larger unfavourable increase in the intraprotein Coulombic energy at higher temperature. Possible causes for the unexpected stability of the protein at 550 K are also discussed.

Self-organizing neural networks bridge the biomolecular resolution gap.

Topology-representing neural networks are employed to generate pseudo-atomic structures of large-scale protein assemblies by combining high-resolution data with volumetric data at lower resolution. As an application example, actin monomers and structural subdomains are located in a three-dimensional (3D) image reconstruction from electron micrographs. To test the reliability of the method, the resolution of the atomic model of an actin polymer is lowered to a level typically encountered in electron microscopic reconstructions. The atomic model is restored with a precision nine times the nominal resolution of the corresponding low-resolution density. The presented self-organizing computing method may be used as an information-processing tool for the synthesis of structural data from a variety of biophysical sources.

Computer simulation of protein-protein association kinetics: acetylcholinesterase-fasciculin.

Computer simulations were performed to investigate the role of electrostatic interactions in promoting fast association of acetylcholinesterase with its peptidic inhibitor, the neurotoxin fasciculin. The encounter of the two macromolecules was simulated with the technique of Brownian dynamics (BD), using atomically detailed structures, and association rate constants were calculated for the wild-type and a number of mutant proteins. In a first set of simulations, the ordering of the experimental rate constants for the mutant proteins was correctly reproduced, although the absolute values of the rate constants were overestimated by a factor of around 30. Rigorous calculations of the full electrostatic interaction energy between the two proteins indicate that this overestimation of association rates results at least in part from approximations made in the description of interaction energetics in the BD simulations. In particular, the initial BD simulations neglect the unfavourable electrostatic desolvation effects that result from the exclusion of high dielectric solvent that accompanies the approach of the two low dielectric proteins. This electrostatic desolvation component is so large that the overall contribution of electrostatics to the binding energy of the complex is unlikely to be strongly favourable. Nevertheless, electrostatic interactions are still responsible for increased association rates, because even if they are unfavourable in the fully formed complex, they are still favourable at intermediate protein-protein separation distances. It therefore appears possible for electrostatic interactions to promote the kinetics of binding even if they do not make a strongly favourable contribution to the thermodynamics of binding. When an approximate description of these electrostatic desolvation effects is included in a second set of BD simulations, the relative ordering of the mutant proteins is again correctly reproduced, but now association rate constants that are much closer in magnitude to the experimental values are obtained. Inclusion of electrostatic desolvation effects also improves reproduction of the experimental ionic strength dependence of the wild-type association rate.