REACH coarse-grained biomolecular simulation: transferability between different protein structural classes

Coarse graining of protein interactions provides a means of simulating large biological systems. The REACH (Realistic Extension Algorithm via Covariance Hessian) coarse-graining method, in which the force constants of a residue-scale elastic network model are calculated from the variance-covariance matrix obtained from atomistic molecular dynamics (MD) simulation, involves direct mapping between scales without the need for iterative optimization. Here, the transferability of the REACH force field is examined between protein molecules of different structural classes. As test cases, myoglobin (all α), plastocyanin (all β), and dihydrofolate reductase (α/β) are taken. The force constants derived are found to be closely similar in all three proteins. An MD version of REACH is presented, and low-temperature coarse-grained (CG) REACH MD simulations of the three proteins are compared with atomistic MD results. The mean-square fluctuations of the atomistic MD are well reproduced by the CGMD. Model functions for the CG interactions, derived by averaging over the three proteins, are also shown to produce fluctuations in good agreement with the atomistic MD. The results indicate that, similarly to the use of atomistic force fields, it is now possible to use a single, generic REACH force field for all protein studies, without having first to derive parameters from atomistic MD simulation for each individual system studied. The REACH method is thus likely to be a reliable way of determining spatiotemporal motion of a variety of proteins without the need for expensive computation of long atomistic MD simulations.     

REACH Coarse-Grained Simulation of a Cellulose Fiber

A molecular level understanding of the structure, dynamics and mechanics of cellulose fibers can aid in understanding the recalcitrance of biomass to hydrolysis in cellulosic biofuel production. Here, a residue-scale REACH (Realistic Extension Algorithm via Covariance Hessian) coarse-grained force field was derived from all-atom molecular dynamics (MD) simulations of the crystalline Iβ cellulose fibril. REACH maps the atomistic covariance matrix onto coarse-grained elastic force constants. The REACH force field was found to reproduce the positional fluctuations and low-frequency vibrational spectra from the all-atom model, allowing elastic properties of the cellulose fibril to be characterized using the coarse-grained force field with a speedup of >20 relative to atomistic MD on systems of the same size. The calculated longitudinal/transversal Young's modulus and the velocity of sound are in agreement with experiment. The persistence length of a 36-chain cellulose microcrystal was estimated to be ～380 μm. Finally, the normal-mode analysis with the REACH force field suggests that intrinsic dynamics might facilitate the deconstruction of the cellulose fibril from the hydrophobic surface.     

REACH Coarse-Grained Normal Mode Analysis of Protein Dimer Interaction Dynamics

The REACH (realistic extension algorithm via covariance Hessian) coarse-grained biomolecular simulation method is a self-consistent multiscale approach directly mapping atomistic molecular dynamics simulation results onto a residue-scale model. Here, REACH is applied to calculate the dynamics of protein-protein interactions. The intra- and intermolecular fluctuations and the intermolecular vibrational densities of states derived from atomistic molecular dynamics are well reproduced by the REACH normal modes. The phonon dispersion relations derived from the REACH lattice dynamics model of crystalline ribonuclease A are also in satisfactory agreement with the corresponding all-atom results. The REACH model demonstrates that increasing dimer interaction strength decreases the translational and rotational intermolecular vibrational amplitudes, while their vibrational frequencies are relatively unaffected. A comparative study of functionally interacting biological dimers with crystal dimers, which are formed artificially via crystallization, reveals a relation between their static structures and the interprotein dynamics: i.e., the consequence of the extensive interfaces of biological dimers is reduction of the intermonomer translational and rotational amplitudes, but not the frequencies.     

Spectral signatures of heterogeneous protein ensembles revealed by MD Simulations of 2DIR spectra

A model for the calculation of amide I FTIR and 2DIR spectra taking into account fluctuations in hydrogen bonding and structure from molecular dynamics (MD) simulations is tested on three systems. It is found that although the homogeneous lineshape approximation yields satisfactory FTIR spectra, 2DIR spectra are sensitive to the inhomogeneity naturally present in any solvated protein and the common approximations of a static structure and averaged-effect solvent are invalid. By building on the local amide Hamiltonian and incorporating site energy variation with electrostatic-based models and disorder from MD trajectories, good agreement is obtained between calculated and measured 2DIR spectra. The largest contribution to the observed inhomogeneity is found to be the fluctuating site energies, which in turn are most sensitive to the water solvent. With the ability to accurately predict 2DIR spectra from atomistic simulations, new opportunities to test force fields and mechanistic predictions from MD are revealed.     

A kinetic and equilibrium study of ligand binding to the monomeric and dimeric haem-containing globins of two chitons.

The radular muscles of the amphineuran molluscs Amaurochiton glaucus and Sipharochiton pelliserpentis contain both a dimeric and a monomeric form of myoglobin. The dimeric form of the protein is composed of two polypeptide chains covalently linked to each other via one or more disulphide bonds. The dimeric protein shows co-operative O2-binding curves. Kinetic investigations indicate that CO binding is co-operative in the dimeric protein, subsequent to full photolysis, but mono-exponential following 10% photolysis. O2 recombination following part photolysis is mono-exponential in the dimeric form, whereas O2 dissociation kinetics indicates the presence of chain heterogeneity. The monomeric form of the protein exhibits mono-exponential time courses in all the experimental situations explored. Although the rate constants associated with the reactions of individual dimer and monomer molecular species are very different, the two species of chiton investigated show remarkably similar properties when compared with each other. All the reactions studied are pH-independent in the range pH 6-8. Amino acid analysis indicates that the monomeric units that combine to form the dimeric species are not identical with the naturally occurring monomeric form. A comparison is made between the chiton myoglobins and other similar O2-binding proteins.     

Role of N‐glycosylation in EGFR ectodomain ligand binding

The epidermal growth factor receptor (EGFR) is a tyrosine kinase protein, overexpressed in several cancers. The extracellular domain of EGFR is known to be heavily glycosylated. Growth factor (mostly epidermal growth factor or EGF) binding activates EGFR. This occurs by inducing the transition from the autoinhibited tethered conformation to an extended conformation of the monomeric form of EGFR and by stabilizing the flexible preformed dimer. Activated EGFR adopts a back‐to‐back dimeric conformation after binding of another homologous receptor to its extracellular domain as the dimeric partner. Several antibodies inhibit EGFR by targeting the growth factor binding site or the dimeric interfaces. Glycosylation has been shown to be important for modulating the stability and function of EGFR. Here, atomistic MD simulations show that N‐glycosylation of the EGFR extracellular domain plays critical roles in the binding of growth factors, monoclonal antibodies, and the dimeric partners to the monomeric EGFR extracellular domain. N‐glycosylation results in the formation of several noncovalent interactions between the glycans and EGFR extracellular domain near the EGF binding site. This stabilizes the growth factor binding site, resulting in stronger interactions (electrostatic) between the growth factor and EGFR. N‐glycosylation also helps maintain the dimeric interface and plays distinct roles in binding of antibodies to spatially separated epitopes of the EGFR extracellular domain. Analysis of SNP data suggests the possibility of altered glycosylation with functional consequences. Proteins 2017; 85:1529–1549. © 2017 Wiley Periodicals, Inc.     

Systematic multiscale parameterization of heterogeneous elastic network models of proteins

We present a method to parameterize heterogeneous elastic network models (heteroENMs) of proteins to reproduce the fluctuations observed in atomistic simulations. Because it is based on atomistic simulation, our method allows the development of elastic coarse-grained models of proteins under different conditions or in different environments. The method is simple and applicable to models at any level of coarse-graining. We validated the method in three systems. First, we computed the persistence length of ADP-bound F-actin, using a heteroENM model. The value of 6.1 ± 1.6 μm is consistent with the experimentally measured value of 9.0 ± 0.5 μm. We then compared our method to a uniform elastic network model and a realistic extension algorithm via covariance Hessian (REACH) model of carboxy myoglobin, and found that the heteroENM method more accurately predicted mean-square fluctuations of α-carbon atoms. Finally, we showed that the method captures critical differences in effective harmonic interactions for coarse-grained models of the N-terminal Bin/amphiphysin/Rvs (N-BAR) domain of amphiphysin, by building models of N-BAR both bound to a membrane and free in solution.     

Coarse-grained modeling of the actin filament derived from atomistic-scale simulations

A coarse-grained (CG) procedure that incorporates the information obtained from all-atom molecular dynamics (MD) simulations is presented and applied to actin filaments (F-actin). This procedure matches the averaged values and fluctuations of the effective internal coordinates that are used to define a CG model to the values extracted from atomistic MD simulations. The fluctuations of effective internal coordinates in a CG model are computed via normal-mode analysis (NMA), and the computed fluctuations are matched with the atomistic MD results in a self-consistent manner. Each actin monomer (G-actin) is coarse-grained into four sites, and each site corresponds to one of the subdomains of G-actin. The potential energy of a CG G-actin contains three bonds, two angles, and one dihedral angle; effective harmonic bonds are used to describe the intermonomer interactions in a CG F-actin. The persistence length of a CG F-actin was found to be sensitive to the cut-off distance of assigning intermonomer bonds. Effective harmonic bonds for a monomer with its third nearest neighboring monomers are found to be necessary to reproduce the values of persistence length obtained from all-atom MD simulations. Compared to the elastic network model, incorporating the information of internal coordinate fluctuations enhances the accuracy and robustness for a CG model to describe the shapes of low-frequency vibrational modes. Combining the fluctuation-matching CG procedure and NMA, the achievable time- and length scales of modeling actin filaments can be greatly enhanced. In particular, a method is described to compute the force-extension curve using the CG model developed in this work and NMA. It was found that F-actin is easily buckled under compressive deformation, and a writhing mode is developed as a result. In addition to the bending and twisting modes, this novel writhing mode of F-actin could also play important roles in the interactions of F-actin with actin-binding proteins and in the force-generation process via polymerization.     

Modeling the dynamics of the solvated SL1 domain of HIV-1 genomic RNA

A mode-coupling solution of the Smoluchowski diffusion equation (MCD theory), designed to describe the dynamics of wobbling macromolecules in water, is applied to a macromolecular bead model including water beads in the nearest layers. The necessary statistical averages are evaluated by time averaging along a molecular dynamics (MD) trajectory where both solute and water are introduced as atomistic models. The cross peaks in (1)H nuclear Overhauser effect spectroscopy (NOESY) NMR spectra that are routinely measured to determine biological structures are here calculated for the mutated 23 nucleotides stem-loop fragment of the SL1 domain in the HIV-1(Lai) genomic RNA. The calculations are in acceptable agreement with experiments without requiring any screening of the hydrodynamic interactions. The screening of hydrodynamics was necessary in previous MCD calculations obtained by using the same full atomistic MD trajectory, but a nonsolvated frictional model.     

Role of pH on dimeric interactions for DENV envelope protein: an insight from molecular dynamics study

The entry of dengue viruses is mediated by pH triggering in the host cells. Here we have studied the DENV E protein stability and binding of its units at low and normal pH using MD and MM-PB/SA method for the first time. To investigate the role of pH in dissociation of dimeric protein, we have performed a concise study of hydrogen bonding and other interactions between units of dimer at low and normal pH. The Generalized Born calculation connotes that dimeric unit was relatively less stable and less proned for dimerisation at low pH. Our results provide a theoretical verification for previous assumptions of pH triggering mechanism of dengue envelope protein. During the pH alteration, we found a large decrement in salt bridges which were observed at normal pH. We also compared the flexibility of each unit and found that they exhibit different fluctuations during molecular dynamics simulations.     

Monomeric state of S100P protein: Experimental and molecular dynamics study

S100 proteins constitute a large subfamily of the EF-hand superfamily of calcium binding proteins. They possess one classical EF-hand Ca²⁺-binding domain and an atypical EF-hand domain. Most of the S100 proteins form stable symmetric homodimers. An analysis of literature data on S100 proteins showed that their physiological concentrations could be much lower than dissociation constants of their dimeric forms. It means that just monomeric forms of these proteins are important for their functioning. In the present work, thermal denaturation of apo-S100P protein monitored by intrinsic tyrosine fluorescence has been studied at various protein concentrations within the region from 0.04–10 μM. A transition from the dimeric to monomeric form results in a decrease in protein thermal stability shifting the mid-transition temperature from 85 to 75 °C. Monomeric S100P immobilized on the surface of a sensor chip of a surface plasmon resonance instrument forms calcium dependent 1 to 1 complexes with human interleukin-11 (equilibrium dissociation constant 1.2 nM). In contrast, immobilized interleukin-11 binds two molecules of dimeric S100P with dissociation constants of 32 nM and 288 nM. Since effective dissociation constant of dimeric S100P protein is very low (0.5 μM as evaluated from our data) the sensitivity of the existing physical methods does not allow carrying out a detailed study of S100P monomer properties. For this reason, we have used molecular dynamics methods to evaluate structural changes in S100P upon its transition from the dimeric to monomeric state. 80-ns molecular dynamics simulations of kinetics of formation of S100P, S100B and S100A11 monomers from the corresponding dimers have been carried out. It was found that during the transition from the homo-dimer to monomer form, the three S100 monomer structures undergo the following changes: (1) the helices in the four-helix bundles within each monomer rotate in order to shield the exposed non-polar residues; (2) almost all lost contacts at the dimer interface are substituted with equivalent and newly formed interactions inside each monomer, and new stabilizing interactions are formed; and (3) all monomers recreate functional hydrophobic cores. The results of the present study show that both dimeric and monomeric forms of S100 proteins can be functional.     

Protein molecular dynamics with electrostatic force entirely determined by a single Poisson-Boltzmann calculation

The electrostatic force including the intramolecular Coulombic interactions and the electrostatic contribution of solvation effect were entirely calculated by using the finite difference Poisson-Boltzmann method (FDPB), which was incorporated into the GROMOS96 force field to complete a new finite difference stochastic dynamics procedure (FDSD). Simulations were performed on an insulin dimer. Different relative dielectric constants were successively assigned to the protein interior; a value of 17 was selected as optimal for our system. The simulation data were analyzed and compared with those obtained from 500-ps molecular dynamics (MD) simulation with explicit water and a 500-ps conventional stochastic dynamics (SD) simulation without the mean solvent force. The results indicate that the FDSD method with GROMOS96 force field is suitable to study the dynamics and structure of proteins in solution if used with the optimal protein dielectric constant.     

Amino acid sequence of myoglobin from the chitonLiolophura japonica and a phylogenetic tree for molluscan globins

Myoglobin was isolated from the radular muscle of the chitonLiolophura japonica, a primitive archigastropodic mollusc.Liolophura contains three monomeric myoglobins (I, II, and III), and the complete amino acid sequence of myoglobin I has been determined. It is composed of 145 amino acid residues, and the molecular mass was calculated to be 16,070 D. The E7 distal histidine, which is replaced by valine or glutamine in several molluscan globins, is conserved inLiolophura myoglobin. The autoxidation rate at physiological conditions indicated thatLiolophura oxymyoglobin is fairly stable when compared with other molluscan myoglobins. The amino acid sequence ofLiolophura myoglobin shows low homology (11–21%) with molluscan dimeric myoglobins and hemoglobins, but shows higher homology (26–29%) with monomeric myoglobins from the gastropodic molluscsAplysia, Dolabella, andBursatella. A phylogenetic tree was constructed from 19 molluscan globin sequences. The tree separated them into two distinct clusters, a cluster for muscle myoglobins and a cluster for erythrocyte or gill hemoglobins. The myoglobin cluster is divided further into two subclusters, corresponding to monomeric and dimeric myoglobins, respectively.Liolophura myoglobin was placed on the branch of monomeric myoglobin lineage, showing that it diverged earlier from other monomeric myoglobins. The hemoglobin cluster is also divided into two subclusters. One cluster contains homodimeric, heterodimeric, tetrameric, and didomain chains of erythrocyte hemoglobins of the blood clamsAnadara, Scapharca, andBarbatia. Of special interest is the other subcluster. It consists of three hemoglobin chains derived from the bacterial symbiont-harboring clamsCalyptogena andLucina, in which hemoglobins are supposed to play an important role in maintaining the symbiosis with sulfide bacteria.     

Molecular dynamics guided study of salt bridge length dependence in both fluorinated and non-fluorinated parallel dimeric coiled-coils

The alpha-helical coiled-coil is one of the most common oligomerization motifs found in both native and engineered proteins. To better understand the stability and dynamics of the coiled-coil motifs, including those modified by fluorination, several fluorinated and nonfluorinated parallel dimeric coiled-coil protein structures were designed and modeled. We also attempt to investigate how changing the length and geometry of the important stabilizing salt bridges influences the coiled-coil protein structure. Molecular dynamics (MD) and free energy simulations with AMBER used a particle mesh Ewald treatment of the electrostatics in explicit TIP3P solvent with balanced force field treatments. Preliminary studies with legacy force fields (ff94, ff96, and ff99) show a profound instability of the coiled-coil structures in short MD simulation. Significantly, better behavior is evident with the more balanced ff99SB and ff03 protein force fields. Overall, the results suggest that the coiled-coil structures can readily accommodate the larger acidic arginine or S-2,7-diaminoheptanedoic acid mutants in the salt bridge, whereas substitution of the smaller L-ornithine residue leads to rapid disruption of the coiled-coil structure on the MD simulation time scale. This structural distortion of the secondary structure allows both the formation of large hydration pockets proximal to the charged groups and within the hydrophobic core. Moreover, the increased structural fluctuations and movement lead to a decrease in the water occupancy lifetimes in the hydration pockets. In contrast, analysis of the hydration in the stable dimeric coiled-coils shows high occupancy water sites along the backbone residues with no water occupancy in the hydrophobic core, although transitory water interactions with the salt bridge residues are evident. The simulations of the fluorinated coiled-coils suggest that in some cases fluorination electrostatically stabilizes the intermolecular coiled-coil salt bridges. Structural analyses also reveal different side chain rotamer preferences for leucine when compared with 5,5,5,5',5',5'-hexafluoroleucine mutants. These observed differences in the side chain rotamer populations suggest differential changes in the side chain conformational entropy upon coiled-coil formation when the protein is fluorinated. The free energy of hydration of the isolated 5,5,5,5',5',5'-hexafluoroleucine amino acid is calculated to be 1.1 kcal/mol less stable than leucine; this hydrophobic penalty in the monomer may provide a driving force for coiled-coil dimer formation. Estimation of the ellipticity at 222 nm from a series of snapshots from the MD simulations with DicroCalc shows distinct increases in the ellipticity when the coiled-coil is fluorinated, which suggests that the helicity in the folded coiled-coils is greater when fluorinated.     

Amino acid sequence of myoglobin from the chitonLiolophura japonica and a phylogenetic tree for molluscan globins

Myoglobin was isolated from the radular muscle of the chitonLiolophura japonica, a primitive archigastropodic mollusc.Liolophura contains three monomeric myoglobins (I, II, and III), and the complete amino acid sequence of myoglobin I has been determined. It is composed of 145 amino acid residues, and the molecular mass was calculated to be 16,070 D. The E7 distal histidine, which is replaced by valine or glutamine in several molluscan globins, is conserved inLiolophura myoglobin. The autoxidation rate at physiological conditions indicated thatLiolophura oxymyoglobin is fairly stable when compared with other molluscan myoglobins. The amino acid sequence ofLiolophura myoglobin shows low homology (11–21%) with molluscan dimeric myoglobins and hemoglobins, but shows higher homology (26–29%) with monomeric myoglobins from the gastropodic molluscsAplysia, Dolabella, andBursatella. A phylogenetic tree was constructed from 19 molluscan globin sequences. The tree separated them into two distinct clusters, a cluster for muscle myoglobins and a cluster for erythrocyte or gill hemoglobins. The myoglobin cluster is divided further into two subclusters, corresponding to monomeric and dimeric myoglobins, respectively.Liolophura myoglobin was placed on the branch of monomeric myoglobin lineage, showing that it diverged earlier from other monomeric myoglobins. The hemoglobin cluster is also divided into two subclusters. One cluster contains homodimeric, heterodimeric, tetrameric, and didomain chains of erythrocyte hemoglobins of the blood clamsAnadara, Scapharca, andBarbatia. Of special interest is the other subcluster. It consists of three hemoglobin chains derived from the bacterial symbiont-harboring clamsCalyptogena andLucina, in which hemoglobins are supposed to play an important role in maintaining the symbiosis with sulfide bacteria.     

Molecular dynamics simulation of dimeric and monomeric forms of human prion protein: insight into dynamics and properties

A central theme in prion protein research is the detection of the process that underlies the conformational transition from the normal cellular prion form (PrP(C)) to its pathogenic isoform (PrP(Sc)). Although the three-dimensional structures of monomeric and dimeric human prion protein (HuPrP) have been revealed by NMR spectroscopy and x-ray crystallography, the process underlying the conformational change from PrP(C) to PrP(Sc) and the dynamics and functions of PrP(C) remain unknown. The dimeric form is thought to play an important role in the conformational transition. In this study, we performed molecular dynamics (MD) simulations on monomeric and dimeric HuPrP at 300 K and 500 K for 10 ns to investigate the differences in the properties of the monomer and the dimer from the perspective of dynamic and structural behaviors. Simulations were also undertaken with Asp178Asn and acidic pH, which is known as a disease-associated factor. Our results indicate that the dynamics of the dimer and monomer were similar (e.g., denaturation of helices and elongation of the beta-sheet). However, additional secondary structure elements formed in the dimer might result in showing the differences in dynamics and properties between the monomer and dimer (e.g., the greater retention of dimeric than monomeric tertiary structure).     

Electrostatic Interaction of Globins with Phospholipid Membranes

Biophysics - For the structures of sperm whale myoglobin (swMb), horse heart myoglobin (hhMb), hemoglobin I (HbI) from botfly Gasterophilus intestinalis (giHbI), and monomeric and dimeric...     

Kinetic analysis of the binding of monomeric and dimeric ephrins to Eph receptors: correlation to function in a growth cone collapse assay

Eph receptors and ephrins play important roles in regulating cell migration and positioning during both normal and oncogenic tissue development. Using a surface plasma resonance (SPR) biosensor, we examined the binding kinetics of representative monomeric and dimeric ephrins to their corresponding Eph receptors and correlated the apparent binding affinity with their functional activity in a neuronal growth cone collapse assay. Our results indicate that the Eph receptor binding of dimeric ephrins, formed through fusion with disulfide-linked Fc fragments, is best described using a bivalent analyte model as a two-step process involving an initial monovalent 2:1 binding followed by a second bivalent 2:2 binding. The bivalent binding dramatically decreases the apparent dissociation rate constants with little effect on the initial association rate constants, resulting in a 30- to 6000-fold decrease in apparent equilibrium dissociation constants for the binding of dimeric ephrins to Eph receptors relative to their monomeric counterparts. Interestingly, the change was more prominent in the A-class ephrin/Eph interactions than in the B-class of ephrins to Eph receptors. The increase in apparent binding affinities correlated well with increased activation of Eph receptors and the resulting growth cone collapse. Our kinetic analysis and correlation of binding affinity with function helped us better understand the interactions between ephrins and Eph receptors and should be useful in the design of inhibitors that interfere with the interactions.     

Application of the theory of enzyme subunit interactions to ATP-hydrolyzing enzymes. The case of Na,K-ATPase.

The theory developed by T. L. Hill (1977, Proc. Natl. Acad. Sci. USA, 74:3632-3636) for enzyme interactions is applied to a dimeric enzyme, the subunits of which may each exist in three distinct states (as in a uni-bi kinetic mechanism). It is shown that when simultaneous binding of substrate to both subunits is excluded, the complex kinetic mechanism of the dimer reduces to a simpler scheme with two distinct, but analogous, cycles that are in principle separately observable in kinetic experiments. Because of the intersubunit interactions, which are explicitly taken into account, the two cycles have different Michaelis constants and maximal velocities. The model exhibits negative cooperativity and enhanced reactivity, relative to a monomeric enzyme. The theory is applied to Na,K-ATPase for which a complete, bicyclic, kinetic mechanism and rate constants are available. When taken together with other evidence, structural as well as functional, the striking similarity of the observed kinetics with that developed for a dimeric enzyme strongly suggests that the functional unit of Na,K-ATPase is a dimer. The free energy differences (calculated from the known rate constants) between intermediates are 6-16 kJ/mol, comparable, for example, to the free energy associated with the formation of a base pair in a nucleic acid double helix. The possible relevance of these results for other ATPases is briefly discussed.