A demonstration of the inhomogeneity of the local dielectric response of proteins by molecular dynamics simulations

Interactive patches over amyloid-β oligomers mediate fractal self-assembly

The monomeric units of intrinsically disordered proteins self-assemble into oligomers, protofilaments, and eventually fibrils which may turn into amyloid. The aggregation of these proteins is primarily studied in bulk with no restriction on their degrees of freedom. Herein we experimentally demonstrate that amyloid-β (Aβ) aggregation under diffusion-limited conditions leads to its fractal self-assembly. Confocal microscopy and scanning electron microscopy with energy dispersion x-ray analysis were used to confirm that the fractal self-assemblies were formed from Aβ rather than the salt present in the two supporting media: deionized water and phosphate buffered saline. The results from the molecular docking experiments implicated that electrostatic and hydrophobic patches on the solvent-accessible surface area of the Aβ oligomers mediate the fractal self-assembly. These implications were tested with laser light scattering experiments on the oligomers formed by breaking mature fibrils of Aβ through sonication, which were observed to self-assemble into fractals when sonicated solutions were drop casted. The electrostatic interactions modulate the fractal morphologies with pH of the solution, which leads to a morphological phase transition observed through the variation in their fractal dimension. These transitions provide experimental evidence for the existing theoretical framework in terms of different kinetic models. The higher surface-to-volume ratio of these fractal self-assemblies may have applications in drug delivery, biosensing, and other biomedical applications.

Multiscale modeling shows that dielectric differences make NaV channels faster than KV channels

The generation of action potentials in excitable cells requires different activation kinetics of voltage-gated Na (Na_V) and K (K_V) channels. Na_V channels activate much faster and allow the initial Na⁺ influx that generates the depolarizing phase and propagates the signal. Recent experimental results suggest that the molecular basis for this kinetic difference is an amino acid side chain located in the gating pore of the voltage sensor domain, which is a highly conserved isoleucine in K_V channels but an equally highly conserved threonine in Na_V channels. Mutagenesis suggests that the hydrophobicity of this side chain in Shaker K_V channels regulates the energetic barrier that gating charges cross as they move through the gating pore and control the rate of channel opening. We use a multiscale modeling approach to test this hypothesis. We use high-resolution molecular dynamics to study the effect of the mutation on polarization charge within the gating pore. We then incorporate these results in a lower-resolution model of voltage gating to predict the effect of the mutation on the movement of gating charges. The predictions of our hierarchical model are fully consistent with the tested hypothesis, thus suggesting that the faster activation kinetics of Na_V channels comes from a stronger dielectric polarization by threonine (Na_V channel) produced as the first gating charge enters the gating pore compared with isoleucine (K_V channel).

Influence of Dielectric Environment upon Isotope Effects on Glycoside Heterolysis: Computational Evaluation and Atomic Hessian Analysis

Isotope effects depend upon the polarity of the bulk medium in which a chemical process occurs. Implicit solvent calculations with molecule-shaped cavities show that the equilibrium isotope effect (EIE) for heterolysis of the glycosidic bonds in 5'-methylthioadenosine and in 2-(p-nitrophenoxy)tetrahydropyran, both in water, are very sensitive in the range 2 ≤ ε ≤ 10 to the relative permittivity of the continuum surrounding the oxacarbenium ion. However, different implementations of nominally the same PCM method can lead to opposite trends being predicted for the same molecule. Computational modeling of the influence of the inhomogeneous effective dielectric surrounding a substrate within the protein environment of an enzymic reaction requires an explicit treatment. The EIE (K_H/K_D) for transfer of cyclopentyl, cyclohexyl, tetrahydrofuranyl and tetrahydropyranyl cations from water to cyclohexane is predicted by B3LYP/6-31+G(d) calculations with implicit solvation and confirmed by B3LYP/6-31+G(d)/OPLS-AA calculations with averaging over many explicit solvation configurations. Atomic Hessian analysis, whereby the full Hessian is reduced to the elements belonging to a single atom at the site of isotopic substitution, reveals a remarkable result for both implicit and explicit solvation: the influence of the solvent environment on these EIEs is essentially captured completely by only a 3 × 3 block of the Hessian, although these values must correctly reflect the influence of the whole environment. QM/MM simulation with ensemble averaging has an important role to play in assisting the meaningful interpretation of observed isotope effects for chemical reactions both in solution and catalyzed by enzymes.

Computational simulation of mechanism and isotope effects on acetal heterolysis as a model for glycoside hydrolysis

DFT calculations for the equilibrium isotope effect for deuterium substitution at the anomeric centre Cα in 2-(p-nitrophenoxy)tetrahydropyran with continuum solvation show significant variation in the range of relative permittivity 2 ≤ ε ≤ 10. One-dimensional scans of potential energy (with implicit solvation by water) or of free energy (from QM/MM potentials of mean force with explicit aqueous solvation with a hybrid AM1/OPLS method) for heterolysis of the bond between Cα and the nucleofuge do not show a transition state. A two-dimensional free-energy surface that considers also the distance between Cα and a nucleophilic water indicates a pre-association DN*ANint ‡ mechanism with a transition state involving nucleophilic attack upon an ion-pair intermediate, and this is supported by good agreement between the mean values of the calculated and experimental α-D KIEs. However, the magnitudes of the standard deviations about the mean values for the making and breaking C–O bonds suggest that the transition state is rather plastic, with Cα–Onu≈2 ± 0.4 Å and Cα–Olg≈3 ± 0.5 Å. Not only is nucleophilic solvent assistance necessary, but there is also evidence for electrophilic assistance through specific hydrogen bonding to the nucleofuge.

Protein-Assisted Formation and Stabilization of Catalytically Active Polyoxometalate Species

The effect of the protein environment on the formation and stabilization of an elusive catalytically active polyoxometalate (POM) species, K₆ [Hf(α₂ -P₂ W₁₇ O₆₁ )] (1), is reported. In the co-crystal of hen egg-white lysozyme (HEWL) with 1, the catalytically active monomeric species is observed, originating from the dimeric 1:2 POM form, while it is intrinsically unstable under physiological pH conditions. The protein-assisted dissociation of the dimeric POM was rationalized by means of DFT calculations. The dissociation process is unfavorable in bulk water, but becomes favorable in the protein-POM complex due to the low dielectric response at the protein surface. The crystal structure shows that the monomeric form is stabilized by electrostatic and water-mediated hydrogen bonding interactions with the protein. It interacts at three distinct sites, close to the aspartate-containing hydrolysis sites, demonstrating high selectivity towards peptide bonds containing this residue.

Model for screened, charge-regulated electrostatics of an eye lens protein: Bovine gammaB-crystallin

We model screened, site-specific charge regulation of the eye lens protein bovine gammaB-crystallin (γB) and study the probability distributions of its proton occupancy patterns. Using a simplified dielectric model, we solve the linearized Poisson-Boltzmann equation to calculate a 54×54 work-of-charging matrix, each entry being the modeled voltage at a given titratable site, due to an elementary charge at another site. The matrix quantifies interactions within patches of sites, including γB charge pairs. We model intrinsic pK values that would occur hypothetically in the absence of other charges, with use of experimental data on the dependence of pK values on aqueous solution conditions, the dielectric model, and literature values. We use Monte Carlo simulations to calculate a model grand-canonical partition function that incorporates both the work-of-charging and the intrinsic pK values for isolated γB molecules and we calculate the probabilities of leading proton occupancy configurations, for 4<pH<8 and Debye screening lengths from 6 to 20 Å. We select the interior dielectric value to model γB titration data. At pH 7.1 and Debye length 6.0 Å, on a given γB molecule the predicted top occupancy pattern is present nearly 20% of the time, and 90% of the time one or another of the first 100 patterns will be present. Many of these occupancy patterns differ in net charge sign as well as in surface voltage profile. We illustrate how charge pattern probabilities deviate from the multinomial distribution that would result from use of effective pK values alone and estimate the extents to which γB charge pattern distributions broaden at lower pH and narrow as ionic strength is lowered. These results suggest that for accurate modeling of orientation-dependent γB-γB interactions, consideration of numerous pairs of proton occupancy patterns will be needed.

Lid opening and conformational stability of T1 Lipase is mediated by increasing chain length polar solvents

The dynamics and conformational landscape of proteins in organic solvents are events of potential interest in nonaqueous process catalysis. Conformational changes, folding transitions, and stability often correspond to structural rearrangements that alter contacts between solvent molecules and amino acid residues. However, in nonaqueous enzymology, organic solvents limit stability and further application of proteins. In the present study, molecular dynamics (MD) of a thermostable Geobacillus zalihae T1 lipase was performed in different chain length polar organic solvents (methanol, ethanol, propanol, butanol, and pentanol) and water mixture systems to a concentration of 50%. On the basis of the MD results, the structural deviations of the backbone atoms elucidated the dynamic effects of water/organic solvent mixtures on the equilibrium state of the protein simulations in decreasing solvent polarity. The results show that the solvent mixture gives rise to deviations in enzyme structure from the native one simulated in water. The drop in the flexibility in H₂O, MtOH, EtOH and PrOH simulation mixtures shows that greater motions of residues were influenced in BtOH and PtOH simulation mixtures. Comparing the root mean square fluctuations value with the accessible solvent area (SASA) for every residue showed an almost correspondingly high SASA value of residues to high flexibility and low SASA value to low flexibility. The study further revealed that the organic solvents influenced the formation of more hydrogen bonds in MtOH, EtOH and PrOH and thus, it is assumed that increased intraprotein hydrogen bonding is ultimately correlated to the stability of the protein. However, the solvent accessibility analysis showed that in all solvent systems, hydrophobic residues were exposed and polar residues tended to be buried away from the solvent. Distance variation of the tetrahedral intermediate packing of the active pocket was not conserved in organic solvent systems, which could lead to weaknesses in the catalytic H-bond network and most likely a drop in catalytic activity. The conformational variation of the lid domain caused by the solvent molecules influenced its gradual opening. Formation of additional hydrogen bonds and hydrophobic interactions indicates that the contribution of the cooperative network of interactions could retain the stability of the protein in some solvent systems. Time-correlated atomic motions were used to characterize the correlations between the motions of the atoms from atomic coordinates. The resulting cross-correlation map revealed that the organic solvent mixtures performed functional, concerted, correlated motions in regions of residues of the lid domain to other residues. These observations suggest that varying lengths of polar organic solvents play a significant role in introducing dynamic conformational diversity in proteins in a decreasing order of polarity.

An improved quantum biochemistry description of the glutamate–GluA2 receptor binding within an inhomogeneous dielectric function framework

A new methodology to define the inhomogeneous dielectric constant of protein residues, to apply to the calculation of protein–ligand properties such as the electrostatic interaction.

Explaining RANKL inhibition by OPG through quantum biochemistry computations and insights into peptide-design for the treatment of osteoporosis

Quantum biochemistry computations are applied to precisely describe important protein–protein interactions, providing a basis for the design of inhibitory peptides against osteoporosis.

Characterization of the Dielectric Constant in the Trichoderma reesei Cel7B Active Site

An attempt is made to evaluate the dielectric constant of the Trichoderma reesei Cel7B active site. Through kinetic measurements, the pKa value of the catalytic acid E201 is determined. Mutations (away from E201) with net charge changes are introduced to perturb the E201 pKa. It is shown that the mutation with a +1 charge change (including G225R, G230R, and A335R) decreases the pKa of E201, whereas the mutation with a -1 charge change (including Q149E, A222D, G225D, and G230D) increases the pKa. This effect is consistent with the electrostatic interaction between the changed charge and the E201 side chain. The fitting of the experimental data yields an apparent dielectric constant of 25-80. Molecular dynamics simulations with explicit water molecules indicate that the high solvent accessibility of the active site contributes largely to the high dielectric constant. ONIOM calculations show that high dielectric constant benefits the catalysis through decreasing the energy of the transition state relative to that of the enzyme substrate complex.

Unstable, metastable, or stable halogen bonding interaction involving negatively charged donors? A statistical and computational chemistry study

The noncovalent halogen bonding could be attributed to the attraction between the positively charged σ-hole and a nucleophile. Quantum mechanics (QM) calculation indicated that the negatively charged organohalogens have no positively charged σ-hole on their molecular surface, leading to a postulation of repulsion between negatively charged organohalogens and nucleophiles in vacuum. However, PDB survey revealed that 24% of the ligands with halogen bonding geometry could be negatively charged. Moreover, 36% of ionizable drugs in CMC (Comprehensive Medicinal Chemistry) are possibly negatively charged at pH 7.0. QM energy scan showed that the negatively charged halogen bonding is probably metastable in vacuum. However, the QM calculated bonding energy turned negative in various solvents, suggesting that halogen bonding with negatively charged donors should be stable in reality. Indeed, QM/MM calculation on three crystal structures with negatively charged ligands revealed that the negatively charged halogen bonding was stable. Hence, we concluded that halogen bonding with negatively charged donors is unstable or metastable in vacuum but stable in protein environment, and possesses similar geometric and energetic characteristics as conventional halogen bonding. Therefore, negatively charged organohalogens are still effective halogen bonding donors for medicinal chemistry and other applications.

Solvent effects on isotope effects: methyl cation as a model system

The isotopic sensitivity (CH3(+) vs CD3(+)) of the equilibrium between the methyl cation in vacuum and in solution has been investigated. Two alternative options for describing the shape of the solute cavity within the widely used polarized continuum model for implicit solvation were compared; the UFF and UA0 methods give equilibrium isotope effects (EIEs) that vary as a function of the dielectric constant in opposite directions. The same isotope effect was also obtained as the average over 40 structures from a hybrid quantum mechanical/molecular mechanical molecular dynamics simulation for the methyl cation explicitly solvated by many water molecules; the inverse value of the EIE agrees with UFF but not UA0. The opposing trends may be satisfactorily explained in terms of the different degrees of exposure of the atomic charges to the dielectric continuum in cavities of different shapes.

Study of pK values and effective dielectric constants of ionizable residues in pentapeptides and in staphylococcal nuclease (SNase) using a mean-field approach

The determination of pK values of amino acid residues as a function of temperature and ionic concentration is crucial to understanding the dynamics of various biological processes such as adsorption of peptides and their interactions with active sites of enzymes. In this study we developed a mean-field model to calculate the position-dependent dielectric constants of ionizable groups and the mean electrostatic potential on the surface. Such potential, which takes into account the contributions exerted by neighboring groups and ions in solution, is responsible for the fine-tuning of the pK value of each residue. The proposed model was applied to the amino acids Asp, Glu, Lys, His, Tyr, and Cys, and since the results were consistent with experimentally obtained values, the model was extended and applied to computation of pK values of Gly and Ala pentapeptides and of ionizable residues of the enzyme staphylococcal nuclease (SNase). In this latter case, we used an approach similar to a first-neighbors approximation, and the results turned out to be in good agreement with previously reported data when considering only the interactions of charged groups located at distances of maximally 20 Å. These considerations and the little computational cost involved turn the suggested approach into a promising tool for the modeling of force fields in computational simulations.

L-Asparagine crystals with wide gap semiconductor features: optical absorption measurements and density functional theory computations

Results of optical absorption measurements are presented together with calculated structural, electronic, and optical properties for the anhydrous monoclinic L-asparagine crystal. Density functional theory (DFT) within the generalized gradient approximation (GGA) including dispersion effects (TS, Grimme) was employed to perform the calculations. The optical absorption measurements revealed that the anhydrous monoclinic L-asparagine crystal is a wide band gap material with 4.95 eV main gap energy. DFT-GGA+TS simulations, on the other hand, produced structural parameters in very good agreement with X-ray data. The lattice parameter differences Δa, Δb, Δc between theory and experiment were as small as 0.020, 0.051, and 0.022 Å, respectively. The calculated band gap energy is smaller than the experimental data by about 15%, with a 4.23 eV indirect band gap corresponding to Z → Γ and Z → β transitions. Three other indirect band gaps of 4.30 eV, 4.32 eV, and 4.36 eV are assigned to α3 → Γ, α1 → Γ, and α2 → Γ transitions, respectively. Δ-sol computations, on the other hand, predict a main band gap of 5.00 eV, just 50 meV above the experimental value. Electronic wavefunctions mainly originating from O 2p-carboxyl, C 2p-side chain, and C 2p-carboxyl orbitals contribute most significantly to the highest valence and lowest conduction energy bands, respectively. By varying the lattice parameters from their converged equilibrium values, we show that the unit cell is less stiff along the b direction than for the a and c directions. Effective mass calculations suggest that hole transport behavior is more anisotropic than electron transport, but the mass values allow for some charge mobility except along a direction perpendicular to the molecular layers of L-asparagine which form the crystal, so anhydrous monoclinic L-asparagine crystals could behave as wide gap semiconductors. Finally, the calculations point to a high degree of optical anisotropy for the absorption and complex dielectric function, with more structured curves for incident light polarized along the 100 and 101 directions.

Electrostatics controls the formation of amyloid superstructures in protein aggregation

The possibility for proteins to aggregate in different superstructures, i.e. large-scale polymorphism, has been widely observed, but an understanding of the physicochemical mechanisms behind it is still out of reach. Here we present a theoretical model for the description of a generic aggregate formed from an ensemble of charged proteins. The model predicts the formation of multifractal structures with the geometry of the growth determined by the electrostatic interactions between single proteins. The model predictions are successfully verified in comparison with experimental curves for aggregate growth allowing us to reveal the mechanism of formation of such complex structures. The model is general and is able to predict aggregate morphologies occurring both in vivo and in vitro. Our findings provide a framework where the physical interactions between single proteins, the aggregate morphology, and the growth kinetics are connected into a single model in agreement with the experimental data.

Calculating the sensitivity and robustness of binding free energy calculations to force field parameters

Binding free energy calculations offer a thermodynamically rigorous method to compute protein-ligand binding, and they depend on empirical force fields with hundreds of parameters. We examined the sensitivity of computed binding free energies to the ligand's electrostatic and van der Waals parameters. Dielectric screening and cancellation of effects between ligand-protein and ligand-solvent interactions reduce the parameter sensitivity of binding affinity by 65%, compared with interaction strengths computed in the gas-phase. However, multiple changes to parameters combine additively on average, which can lead to large changes in overall affinity from many small changes to parameters. Using these results, we estimate that random, uncorrelated errors in force field nonbonded parameters must be smaller than 0.02 e per charge, 0.06 Å per radius, and 0.01 kcal/mol per well depth in order to obtain 68% (one standard deviation) confidence that a computed affinity for a moderately-sized lead compound will fall within 1 kcal/mol of the true affinity, if these are the only sources of error considered.

The stabilization effect of dielectric constant and acidic amino acids on arginine-arginine (Arg-Arg) pairings: database survey and computational studies

Database survey in this study revealed that about one-third of the protein structures deposited in the Protein Data Bank (PDB) contain arginine-arginine (Arg-Arg) pairing with a carbon···carbon (CZ···CZ) interaction distance less than 5 Å. All the Arg-Arg pairings were found to bury in a polar environment composed of acidic residues, water molecules, and strong polarizable or negatively charged moieties from binding site or bound ligand. Most of the Arg-Arg pairings are solvent exposed and 68.3% Arg-Arg pairings are stabilized by acidic residues, forming Arg-Arg-Asp/Glu clusters. Density functional theory (DFT) was then employed to study the effect of environment on the pairing structures. It was revealed that Arg-Arg pairings become thermodynamically stable (about -1 kcal/mol) as the dielectric constant increases to 46.8 (DMSO), in good agreement with the results of the PDB survey. DFT calculations also demonstrated that perpendicular Arg-Arg pairing structures are favorable in low dielectric constant environment, while in high dielectric constant environment parallel structures are favorable. Additionally, the acidic residues can stabilize the Arg-Arg pairing structures to a large degree. Energy decomposition analysis of Arg-Arg pairings and Arg-Arg-Asp/Glu clusters showed that both solvation and electrostatic energies contribute significantly to their stability. The results reported herein should be very helpful for understanding Arg-Arg pairing and its application in drug design.

Anisotropy of the Coulomb interaction between folded proteins: consequences for mesoscopic aggregation of lysozyme

Toward quantitative description of protein aggregation, we develop a computationally efficient method to evaluate the potential of mean force between two folded protein molecules that allows for complete sampling of their mutual orientation. Our model is valid at moderate ionic strengths and accounts for the actual charge distribution on the surface of the molecules, the dielectric discontinuity at the protein-solvent interface, and the possibility of protonation or deprotonation of surface residues induced by the electric field due to the other protein molecule. We apply the model to the protein lysozyme, whose solutions exhibit both mesoscopic clusters of protein-rich liquid and liquid-liquid separation; the former requires that protein form complexes with typical lifetimes of approximately milliseconds. We find the electrostatic repulsion is typically lower than the prediction of the Derjaguin-Landau-Verwey-Overbeek theory. The Coulomb interaction in the lowest-energy docking configuration is nonrepulsive, despite the high positive charge on the molecules. Typical docking configurations barely involve protonation or deprotonation of surface residues. The obtained potential of mean force between folded lysozyme molecules is consistent with the location of the liquid-liquid coexistence, but produces dimers that are too short-lived for clusters to exist, suggesting lysozyme undergoes conformational changes during cluster formation.

The dynamics of Ca2+ ions within the solvation shell of calbindin D9k

The encounter of a Ca²⁺ ion with a protein and its subsequent binding to specific binding sites is an intricate process that cannot be fully elucidated from experimental observations. We have applied Molecular Dynamics to study this process with atomistic details, using Calbindin D9k (CaB) as a model protein. The simulations show that in most of the time the Ca²⁺ ion spends within the Debye radius of CaB, it is being detained at the 1st and 2nd solvation shells. While being detained near the protein, the diffusion coefficient of the ion is significantly reduced. However, due to the relatively long period of detainment, the ion can scan an appreciable surface of the protein. The enhanced propagation of the ion on the surface has a functional role: significantly increasing the ability of the ion to scan the protein's surface before being dispersed to the bulk. The contribution of this mechanism to Ca²⁺ binding becomes significant at low ion concentrations, where the intervals between successive encounters with the protein are getting longer. The efficiency of the surface diffusion is affected by the distribution of charges on the protein's surface. Comparison of the Ca²⁺ binding dynamics in CaB and its E60D mutant reveals that in the wild type (WT) protein the carboxylate of E60 function as a preferred landing-site for the Ca²⁺ arriving from the bulk, followed by delivering it to the final binding site. Replacement of the glutamate by aspartate significantly reduced the ability to transfer Ca²⁺ ions from D60 to the final binding site, explaining the observed decrement in the affinity of the mutated protein to Ca²⁺.