On the Use of Different Dielectric Constants for Computing Individual and Pairwise Terms in Poisson−Boltzmann Studies of Protein Ionization Equilibrium

Delineating redox cooperativity in water-soluble and membrane multiheme cytochromes through protein design

Nature has evolved diverse electron transport proteins and multiprotein assemblies essential to the generation and transduction of biological energy. However, substantially modifying or adapting these proteins for user-defined applications or to gain fundamental mechanistic insight can be hindered by their inherent complexity. De novo protein design offers an attractive route to stripping away this confounding complexity, enabling us to probe the fundamental workings of these bioenergetic proteins and systems, while providing robust, modular platforms for constructing completely artificial electron-conducting circuitry. Here, we use a set of de novo designed mono-heme and di-heme soluble and membrane proteins to delineate the contributions of electrostatic micro-environments and dielectric properties of the surrounding protein medium on the inter-heme redox cooperativity that we have previously reported. Experimentally, we find that the two heme sites in both the water-soluble and membrane constructs have broadly equivalent redox potentials in isolation, in agreement with Poisson-Boltzmann Continuum Electrostatics calculations. BioDC, a Python program for the estimation of electron transfer energetics and kinetics within multiheme cytochromes, also predicts equivalent heme sites, and reports that burial within the low dielectric environment of the membrane strengthens heme-heme electrostatic coupling. We conclude that redox cooperativity in our diheme cytochromes is largely driven by heme electrostatic coupling and confirm that this effect is greatly strengthened by burial in the membrane. These results demonstrate that while our de novo proteins present minimalist, new-to-nature constructs, they enable the dissection and microscopic examination of processes fundamental to the function of vital, yet complex, bioenergetic assemblies.

High Charge Density in Peptide Dendrimers is Required to Destabilize Membranes: Insights into Endosome Evasion

Peptide dendrimers are a type of branched, symmetric, and topologically well-defined molecule that have already been used as delivery systems for nucleic acid transfection. Several of the most promising sequences showed high efficiency in many key steps of transfection, namely, binding siRNA, entering cells, and evading the endosome. However, small changes to the peptide dendrimers, such as in the hydrophobic core, the amino acid chirality, or the total available charges, led to significantly different experimental results with unclear mechanistic insights. In this work, we built a computational model of several of those peptide dendrimers (MH18, MH13, and MH47) and some of their variants to study the molecular details of the structure and function of these molecules. We performed CpHMD simulations in the aqueous phase and in interaction with a lipid bilayer to assess how conformation and protonation are affected by pH in different environments. We found that while the different peptide dendrimer sequences lead to no substantial structural differences in the aqueous phase, the total charge and, more importantly, the total charge density are key for the capacity of the dendrimer to interact and destabilize the membrane. These dendrimers become highly charged when the pH changes from 7.5 to 4.5, and the presence of a high charge density, which is decreased for MH47 that has four fewer titratable lysines, is essential to trigger membrane destabilization. These findings are in excellent agreement with the experimental data and help us to understand the high efficiency of some dendrimers and why the dendrimer MH47 is unable to complete the transfection process. This evidence provides further understanding of the mode of action of these peptide dendrimers and will be pivotal for the future design of new sequences with improved transfection capabilities.

Fluctuation Relations to Calculate Protein Redox Potentials from Molecular Dynamics Simulations

The tunable design of protein redox potentials promises to open a range of applications in biotechnology and catalysis. Here, we introduce a method to calculate redox potential changes by combining fluctuation relations with molecular dynamics simulations. It involves the simulation of reduced and oxidized states, followed by the instantaneous conversion between them. Energy differences introduced by the perturbations are obtained using the Kubo-Onsager approach. Using a detailed fluctuation relation coupled with Bayesian inference, these are postprocessed into estimates for the redox potentials in an efficient manner. This new method, denoted MD + CB, is tested on a de novo four-helix bundle heme protein (the m4D2 "maquette") and five designed mutants, including some mutants characterized experimentally in this work. The MD + CB approach is found to perform reliably, giving redox potential shifts with reasonably good correlation (0.85) to the experimental values for the mutants. The MD + CB approach also compares well with redox potential shift predictions using a continuum electrostatic method. The estimation method employed within the MD + CB approach is straightforwardly transferable to standard equilibrium MD simulations and holds promise for redox protein engineering and design applications.

An expandable, modular de novo protein platform for precision redox engineering

The electron-conducting circuitry of life represents an as-yet untapped resource of exquisite, nanoscale biomolecular engineering. Here, we report the characterization and structure of a de novo diheme "maquette" protein, 4D2, which we subsequently use to create an expanded, modular platform for heme protein design. A well-folded monoheme variant was created by computational redesign, which was then utilized for the experimental validation of continuum electrostatic redox potential calculations. This demonstrates how fundamental biophysical properties can be predicted and fine-tuned. 4D2 was then extended into a tetraheme helical bundle, representing a 7 nm molecular wire. Despite a molecular weight of only 24 kDa, electron cryomicroscopy illustrated a remarkable level of detail, indicating the positioning of the secondary structure and the heme cofactors. This robust, expressible, highly thermostable and readily designable modular platform presents a valuable resource for redox protein design and the future construction of artificial electron-conducting circuitry.

Computational Study of the pH-Dependent Ionic Environment around β-Lactoglobulin

Ions are involved in multiple biological processes and may exist bound to biomolecules or may be associated with their surface. Although the presence of ions in nucleic acids has traditionally gained more interest, ion-protein interactions, often with a marked dependency on pH, are beginning to gather attention. Here we present a detailed analysis on the binding and distribution of ions around β-lactoglobulin using a constant-pH MD (CpHMD) method, at a pH range 3-8, and compare it with the more traditional Poisson-Boltzmann (PB) model and the existing experimental data. Most analyses used ion concentration maps built around the protein, obtained from either the CpHMD simulations or PB calculations. The requirements of approximate charge neutrality and ionic strength equal to bulk, imposed on the MD box, imply that the absolute value of the ion excess should be half the protein charge, which is in agreement with experimental observation on other proteins ( Proc. Natl. Acad. Sci. U.S.A. 2021, 118, e2015879118) and lends support to this protocol. In addition, the protein total charge (including territorially bound ions) estimated with MD is in excellent agreement with electrophoretic measurements. Overall, the CpHMD simulations show good agreement with the nonlinear form of the PB (NLPB) model but not with its linear form, which involves a theoretical inconsistency in the calculation of the concentration maps. In several analyses, the observed pH-dependent trends for the counterions and co-ions are those generally expected, and the ion concentration maps correctly converge to the bulk ionic strength as one moves away from the protein. Despite the overall similarity, the CpHMD and NLPB approaches show some discrepancies when analyzed in more detail, which may be related to an apparent overestimation of counterion excess and underestimation of co-ion exclusion by the NLPB model, particularly at short distances from the protein.

Increasing the Realism of in Silico pHLIP Peptide Models with a Novel pH Gradient CpHMD Method

The pH-low insertion peptides (pHLIP) are pH-dependent membrane inserting peptides, whose function depends on the cell microenvironment acidity. Several peptide variants have been designed to improve upon the wt-sequence, particularly the state transition kinetics and the selectivity for tumor pH. The variant 3 (Var3) peptide is a 27 residue long peptide, with a key titrating residue (Asp-13) that, despite showing a modest performance in liposomes (pK^ins ∼ 5.0), excelled in tumor cell experiments. To help rationalize these results, we focused on the pH gradient in the cell membrane, which is one of the crucial properties that are not present in liposomes. We extended our CpHMD-L method and its pH replica-exchange (pHRE) implementation to include a pH gradient and mimic the pHLIP-membrane microenvironment in a cell where the internal pH is fixed (pH 7.2) and the external pH is allowed to change. We showed that, by properly modeling the pH-gradient, we can correctly predict the experimentally observed loss and gain of performance in tumor cells experiments by the wt and Var3 sequences, respectively. In sum, the pH gradient implementation allowed for more accurate and realistic pK_a estimations and was a pivotal step in bridging the in silico data and the in vivo cell experiments.

Extending the Stochastic Titration CpHMD to CHARMM36m

The impact of pH on proteins is significant but often neglected in molecular dynamics simulations. Constant-pH Molecular Dynamics (CpHMD) is the state-of-the-art methodology to deal with these effects. However, it still lacks widespread adoption by the scientific community. The stochastic titration CpHMD is one of such methods that, until now, only supported the GROMOS force field family. Here, we extend this method's implementation to include the CHARMM36m force field available in the GROMACS software package. We test this new implementation with a diverse group of proteins, namely, lysozyme, Staphylococcal nuclease, and human and E. coli thioredoxins. All proteins were conformationally stable in the simulations, even at extreme pH values. The RMSE values (pK_a prediction vs experimental) obtained were very encouraging, in particular for lysozyme and human thioredoxin. We have also identified a few residues that challenged the CpHMD simulations, highlighting scenarios where the method still needs improvement independently of the force field. The CHARMM36m all-atom implementation was more computationally efficient when compared with the GROMOS 54A7, taking advantage of a shorter nonbonded interaction cutoff and a less frequent neighboring list update. The new extension will allow the study of pH effects in many systems for which this force field is particularly suited, i.e., proteins, membrane proteins, lipid bilayers, and nucleic acids.

Unveiling molecular details behind improved activity at neutral to alkaline pH of an engineered DyP-type peroxidase

DyP-type peroxidases (DyPs) are microbial enzymes that catalyze the oxidation of a wide range of substrates, including synthetic dyes, lignin-derived compounds, and metals, such as Mn²⁺ and Fe²⁺, and have enormous biotechnological potential in biorefineries. However, many questions on the molecular basis of enzyme function and stability remain unanswered. In this work, high-resolution structures of PpDyP wild-type and two engineered variants (6E10 and 29E4) generated by directed evolution were obtained. The X-ray crystal structures revealed the typical ferredoxin-like folds, with three heme access pathways, two tunnels, and one cavity, limited by three long loops including catalytic residues. Variant 6E10 displays significantly increased loops' flexibility that favors function over stability: despite the considerably higher catalytic efficiency, this variant shows poorer protein stability compared to wild-type and 29E4 variants. Constant-pH MD simulations revealed a more positively charged microenvironment near the heme pocket of variant 6E10, particularly in the neutral to alkaline pH range. This microenvironment affects enzyme activity by modulating the pK _a of essential residues in the heme vicinity and should account for variant 6E10 improved activity at pH 7-8 compared to the wild-type and 29E4 that show optimal enzymatic activity close to pH 4. Our findings shed light on the structure-function relationships of DyPs at the molecular level, including their pH-dependent conformational plasticity. These are essential for understanding and engineering the catalytic properties of DyPs for future biotechnological applications.

Novel US-CpHMD Protocol to Study the Protonation-Dependent Mechanism of the ATP/ADP Carrier

We have designed a protocol combining constant-pH molecular dynamics (CpHMD) simulations with an umbrella sampling (US) scheme (US-CpHMD) to study the mechanism of ADP/ATP transport (import and export) by their inner mitochondrial membrane carrier protein [ADP/ATP carrier (AAC)]. The US scheme helped overcome the limitations of sampling the slow kinetics involved in these substrates' transport, while CpHMD simulations provided an unprecedented realism by correctly capturing the associated protonation changes. The import of anionic substrates along the mitochondrial membrane has a strong energetic disadvantage due to a smaller substrate concentration and an unfavorable membrane potential. These limitations may have created an evolutionary pressure on AAC to develop specific features benefiting the import of ADP. In our work, the potential of mean force profiles showed a clear selectivity in the import of ADP compared to ATP, while in the export, no selectivity was observed. We also observed that AAC sequestered both substrates at longer distances in the import compared to the export process. Furthermore, only in the import process do we observe transient protonation of both substrates when going through the AAC cavity, which is an important advantage to counteract the unfavorable mitochondrial membrane potential. Finally, we observed a substrate-induced disruption of the matrix salt-bridge network, which can promote the conformational transition (from the C- to M-state) required to complete the import process. This work unraveled several important structural features where the complex electrostatic interactions were pivotal to interpreting the protein function and illustrated the potential of applying the US-CpHMD protocol to other transport processes involving membrane proteins.

Approach to Study pH-Dependent Protein Association Using Constant-pH Molecular Dynamics: Application to the Dimerization of β-Lactoglobulin

Protein-protein association is often mediated by electrostatic interactions and modulated by pH. However, experimental and computational studies have often overlooked the effect of association on the protonation state of the protein. In this work, we present a methodological approach based on constant-pH molecular dynamics (MD), which aims to provide a detailed description of a pH-dependent protein-protein association, and apply it to the dimerization of β-lactoglobulin (BLG). A selection of analyses is performed using the data generated by constant-pH MD simulations of monomeric and dimeric forms of bovine BLG, in the pH range 3-8. First, we estimate free energies of dimerization using a computationally inexpensive approach based on the Wyman-Tanford linkage theory, calculated in a new way through the use of thermodynamically based splines. The individual free energy contribution of each titratable site is also calculated, allowing for identification of relevant residues. Second, the correlations between the proton occupancies of pairs of sites are calculated (using the Pearson coefficient), and extensive networks of correlated sites are observed at acidic pH values, sometimes involving distant pairs. In general, strongly correlated sites are also slow proton exchangers and contribute significantly to the pH-dependency of the dimerization free energy. Third, we use ionic density as a fingerprint of protein charge distribution and observe electrostatic complementarity between the monomer faces that form the dimer interface, more markedly at the isoionic point (where maximum dimerization occurs) than at other pH values, which might contribute to guide the association. Finally, the pH-dependent dimerization modes are inspected using PCA, among other analyses, and two states are identified: a relaxed state at pH 4-8 (with the typical alignment of the crystallographic structure) and a compact state at pH 3-4 (with a tighter association and rotated alignment). This work shows that an approach based on constant-pH MD simulations can produce rich detailed pictures of pH-dependent protein associations, as illustrated for BLG dimerization.

Prediction of protein pK a with representation learning

The behavior of proteins is closely related to the protonation states of the residues. Therefore, prediction and measurement of pK a are essential to understand the basic functions of proteins. In this work, we develop a new empirical scheme for protein pK a prediction that is based on deep representation learning. It combines machine learning with atomic environment vector (AEV) and learned quantum mechanical representation from ANI-2x neural network potential (J. Chem. Theory Comput. 2020, 16, 4192). The scheme requires only the coordinate information of a protein as the input and separately estimates the pK a for all five titratable amino acid types. The accuracy of the approach was analyzed with both cross-validation and an external test set of proteins. Obtained results were compared with the widely used empirical approach PROPKA. The new empirical model provides accuracy with MAEs below 0.5 for all amino acid types. It surpasses the accuracy of PROPKA and performs significantly better than the null model. Our model is also sensitive to the local conformational changes and molecular interactions.

Predicting stable binding modes from simulated dimers of the D76N mutant of β 2-microglobulin

•

β2m D76N mutant populates an aggregation-prone monomer (I2) with unstructured termini.

•

MD and MM-PBSA indicate that I2 dimers are stabilized by hydrophobic interactions.

•

The termini regions and BC- and DE-loops are prevalent in the most stable interfaces.

•

The most stable dimer has a limited growth potential without structural rearrangement.

The D76N mutant of the β2m protein is a biologically motivated model system to study protein aggregation. There is strong experimental evidence, supported by molecular simulations, that D76N populates a highly dynamic conformation (which we originally named I2) that exposes aggregation-prone patches as a result of the detachment of the two terminal regions. Here, we use Molecular Dynamics simulations to study the stability of an ensemble of dimers of I2 generated via protein–protein docking. MM-PBSA calculations indicate that within the ensemble of investigated dimers the major contribution to interface stabilization at physiological pH comes from hydrophobic interactions between apolar residues. Our structural analysis also reveals that the interfacial region associated with the most stable binding modes are particularly rich in residues pertaining to both the N- and C-terminus, as well residues from the BC- and DE-loops. On the other hand, the less stable interfaces are stabilized by intermolecular interactions involving residues from the CD- and EF-loops. By focusing on the most stable binding modes, we used a simple geometric rule to propagate the corresponding dimer interfaces. We found that, in the absence of any kind of structural rearrangement occurring at an early stage of the oligomerization pathway, some interfaces drive a self-limited growth process, while others can be propagated indefinitely allowing the formation of long, polymerized chains. In particular, the interfacial region of the most stable binding mode reported here falls in the class of self-limited growth.

pK a Calculations in Membrane Proteins from Molecular Dynamics Simulations

The conformational changes of membrane proteins are crucial to their function and usually lead to fluctuations in the electrostatic environment of the protein surface. A very effective way to quantify these changes is by calculating the pK _a values of the protein's titratable residues, which can be regarded as electrostatic probes. To achieve this, we need to take advantage of the fast and reliable pK _a calculators developed for globular proteins and adapt them to include the explicit effects of membranes. Here, we provide a detailed linear response approximation protocol that uses our own software (PypKa) to calculate reliable pK _a values from short MD simulations of membrane proteins.

Improved Protocol to Tackle the pH Effects on Membrane-Inserting Peptides

Many important biological pathways rely on membrane-interacting peptides or proteins, which can alter the biophysical properties of the cell membrane by simply adsorbing to its surface to undergo a full insertion process. To study these phenomena with atomistic detail, model peptides have been used to refine the current computational methodologies. Improvements have been made with force-field parameters, enhanced sampling techniques to obtain faster sampling, and the addition of chemical-physical properties, such as pH, whose influence dramatically increases at the water/membrane interface. The pH (low) insertion peptide (pHLIP) is a peptide that inserts across a membrane bilayer depending on the pH due to the presence of a key residue (Asp14) whose acidity-induced protonation triggers the whole process. The complex nature of these peptide/membrane interactions resulted in sampling limitations of the protonation and configurational space albeit using state-of-the-art methods such as the constant-pH molecular dynamics. To address this issue and circumvent those limitations, new simulations were performed with our newly developed pH-replica exchange method using wild-type (wt)-pHLIP in different 2-oleoyl-1-palmitoyl-sn-glycero-3-phosphocholine membrane sizes. This technique provided enhanced sampling and allowed for the calculation of more complete Asp14 pK_a profiles. The conformational heterogeneity derived from strong electrostatic interactions between Asp14 and the lipid phosphate groups was identified as the source of most pK_a variability. In spite of these persistent and harder-to-equilibrate phosphate interactions, the pK_a values at deeper regions (6.0-6.2) still predicted the experimental pK of insertion (6.0) since the electrostatic perturbation decays as the residue inserts further into the membrane. We also observed that reducing the system size leads to membrane deformations where it increasingly loses the ability to accommodate the pHLIP-induced perturbations. This indicates that large membrane patches, such as 256 or even 352 lipids, are needed to obtain stable and more realistic pHLIP/membrane systems. These results strengthen our method pK_a predictive and analytical capabilities to study the intricate play of electrostatic effects of the peptide/membrane interface, granting confidence for future applications in similar systems.

Improved GROMOS 54A7 Charge Sets for Phosphorylated Tyr, Ser, and Thr to Deal with pH-Dependent Binding Phenomena

Phosphorylation is a ubiquitous post-translational modification in proteins, and the phosphate group is present constitutively or transiently in most biological building blocks. These phosphorylated biomolecules are involved in many high-affinity binding/unbinding events that rely predominantly on electrostatic interactions. To build accurate models of these molecules, we need an improved description of the atomic partial charges for all relevant protonation states. In this work, we showed that the commonly used protocols to derive atomic partial charges using well-solvated molecules are inadequate to model the protonation equilibria in binding events. We introduced a protocol based on PB/MC calculations with a single representative conformation (of both protonation states) and used the resulting pK_a estimations to help manually curate the atomic partial charges. The final charge set, which is fully compatible with the GROMOS 54A7 force field, proved to be very effective in modeling the protonation equilibrium in different phosphorylated peptides in the free (tetrapeptides, pentapeptides, and pY1021) and protein-complexed forms (pY1021/PLC-γ1 complex). This was particularly important in the case of the pY1021 bound to the SH2 domain of PLC-γ1, where only our curated charge set captured the correct protonation equilibrium at the neutral to slightly acidic pH range. The binding/unbinding phenomena in that pH range are biologically relevant, and to improve our models, we need to go beyond the commonly used protocols and obtain revised force field parameters for these molecules.

The Lysosomotropic Activity of Hydrophobic Weak Base Drugs is Mediated via Their Intercalation into the Lysosomal Membrane

Lipophilic weak base therapeutic agents, termed lysosomotropic drugs (LDs), undergo marked sequestration and concentration within lysosomes, hence altering lysosomal functions. This lysosomal drug entrapment has been described as luminal drug compartmentalization. Consistent with our recent finding that LDs inflict a pH-dependent membrane fluidization, we herein demonstrate that LDs undergo intercalation and concentration within lysosomal membranes. The latter was revealed experimentally and computationally by (a) confocal microscopy of fluorescent compounds and drugs within lysosomal membranes, and (b) molecular dynamics modeling of the pH-dependent membrane insertion and accumulation of an assortment of LDs, including anticancer drugs. Based on the multiple functions of the lysosome as a central nutrient sensory hub and a degradation center, we discuss the molecular mechanisms underlying the alteration of morphology and impairment of lysosomal functions as consequences of LDs’ intercalation into lysosomes. Our findings bear important implications for drug design, drug induced lysosomal damage, diseases and pertaining therapeutics.

The Early Phase of β2m Aggregation: An Integrative Computational Study Framed on the D76N Mutant and the ΔN6 Variant

Human β2-microglobulin (b2m) protein is classically associated with dialysis-related amyloidosis (DRA). Recently, the single point mutant D76N was identified as the causative agent of a hereditary systemic amyloidosis affecting visceral organs. To get insight into the early stage of the β2m aggregation mechanism, we used molecular simulations to perform an in depth comparative analysis of the dimerization phase of the D76N mutant and the ΔN6 variant, a cleaved form lacking the first six N-terminal residues, which is a major component of ex vivo amyloid plaques from DRA patients. We also provide first glimpses into the tetramerization phase of D76N at physiological pH. Results from extensive protein–protein docking simulations predict an essential role of the C- and N-terminal regions (both variants), as well as of the BC-loop (ΔN6 variant), DE-loop (both variants) and EF-loop (D76N mutant) in dimerization. The terminal regions are more relevant under acidic conditions while the BC-, DE- and EF-loops gain importance at physiological pH. Our results recapitulate experimental evidence according to which Tyr10 (A-strand), Phe30 and His31 (BC-loop), Trp60 and Phe62 (DE-loop) and Arg97 (C-terminus) act as dimerization hot-spots, and further predict the occurrence of novel residues with the ability to nucleate dimerization, namely Lys-75 (EF-loop) and Trp-95 (C-terminus). We propose that D76N tetramerization is mainly driven by the self-association of dimers via the N-terminus and DE-loop, and identify Arg3 (N-terminus), Tyr10, Phe56 (D-strand) and Trp60 as potential tetramerization hot-spots.

Identification of the Initial Steps in Signal Transduction in the α4β2 Nicotinic Receptor: Insights from Equilibrium and Nonequilibrium Simulations

Nicotinic acetylcholine receptors (nAChRs) modulate synaptic transmission in the nervous system. These receptors have emerged as therapeutic targets in drug discovery for treating several conditions, including Alzheimer's disease, pain, and nicotine addiction. In this in silico study, we use a combination of equilibrium and nonequilibrium molecular dynamics simulations to map dynamic and structural changes induced by nicotine in the human α4β2 nAChR. They reveal a striking pattern of communication between the extracellular binding pockets and the transmembrane domains (TMDs) and show the sequence of conformational changes associated with the initial steps in this process. We propose a general mechanism for signal transduction for Cys-loop receptors: the mechanistic steps for communication proceed firstly through loop C in the principal subunit, and are subsequently transmitted, gradually and cumulatively, to loop F of the complementary subunit, and then to the TMDs through the M2-M3 linker.

PKA17-A Coarse-Grain Grid-Based Methodology and Web-Based Software for Predicting Protein pK a Shifts

We have developed and tested PKA17, a coarse-grain grid-based model for predicting protein pK _a shifts. Our pK _a predictor is currently deployed via a website interface. We have carried out parameter fitting using 442 Asp, Glu, His, Lys, and Arg residues for which experimental results are available in the literature. PROPKA software has been used for benchmarking. The average unsigned error and root-mean-square deviation (RMSD) have been found to be 0.628 and 0.831 pH units, respectively, for PKA17. The corresponding results with PROPKA are 0.761 and 1.063 units. We have assessed the robustness of the developed PKA17 methodology with a number of tests and have also explored the possibility of using a combination of PROPKA and PKA17 calculations in order to improve the accuracy of predicted pK _a values for protein residues. We have also once again confirmed that protein acidity constants are influenced almost entirely by residues in the immediate spatial proximity of the ionizable amino acids. The resulting PKA17 software has been deployed online with a web-based interface at http://users.wpi.edu/~jpcvitkovic/pka_calc.html. © 2019 Wiley Periodicals, Inc.

Role of Counterions in Constant-pH Molecular Dynamics Simulations of PAMAM Dendrimers

Electrostatic interactions play a pivotal role in the structure and mechanism of action of most biomolecules. There are several conceptually different methods to deal with electrostatics in molecular dynamics simulations. Ionic strength effects are usually introduced using such methodologies and can have a significant impact on the quality of the final conformation space obtained. We have previously shown that full system neutralization can lead to wrong lipidic phases in the 25% PA/PC bilayer (J. Chem. Theory Comput. 2014,10, 5483-5492). In this work, we investigate how two limit approaches to the ionic strength treatment (implicitly with GRF or using full system neutralization with either GRF or PME) can influence the conformational space of the second-generation PAMAM dendrimer. Constant-pH MD simulations were used to map PAMAM's conformational space at its full pH range (from 2.5 to 12.5). Our simulations clearly captured the coupling between protonation and conformation in PAMAM. Interestingly, the dendrimer conformational distribution was almost independent of the ionic strength treatment methods, which is in contrast to what we have observed in charged lipid bilayers. Overall, our results confirm that both GRF with implicit ionic strength and a fully neutralized system with PME are valid approaches to model charged globular systems, using the GROMOS 54A7 force field.

Coupled molecular dynamics and continuum electrostatic method to compute the ionization pKa's of proteins as a function of pH. Test on a large set of proteins

A computational method, to predict the pKa values of the ionizable residues Asp, Glu, His, Tyr, and Lys of proteins, is presented here. Calculation of the electrostatic free-energy of the proteins is based on an efficient version of a continuum dielectric electrostatic model. The conformational flexibility of the protein is taken into account by carrying out molecular dynamics simulations of 10 ns in implicit water. The accuracy of the proposed method of calculation of pKa values is estimated from a test set of experimental pKa data for 297 ionizable residues from 34 proteins. The pKa-prediction test shows that, on average, 57, 86, and 95% of all predictions have an error lower than 0.5, 1.0, and 1.5 pKa units, respectively. This work contributes to our general understanding of the importance of protein flexibility for an accurate computation of pKa, providing critical insight about the significance of the multiple neutral states of acid and histidine residues for pKa-prediction, and may spur significant progress in our effort to develop a fast and accurate electrostatic-based method for pKa-predictions of proteins as a function of pH.

The key role of glutamate 172 in the mechanism of type II NADH:quinone oxidoreductase of Staphylococcus aureus

Type II NADH:quinone oxidoreductases (NDH-2s) are membrane bound enzymes that deliver electrons to the respiratory chain by oxidation of NADH and reduction of quinones. In this way, these enzymes also contribute to the regeneration of NAD⁺, allowing several metabolic pathways to proceed. As for the other members of the two-Dinucleotide Binding Domains Flavoprotein (tDBDF) superfamily, the enzymatic mechanism of NDH-2s is still little explored and elusive. In this work we addressed the role of the conserved glutamate 172 (E172) residue in the enzymatic mechanism of NDH-2 from Staphylococcus aureus. We aimed to test our earlier hypothesis that E172 plays a key role in proton transfer to allow the protonation of the quinone. For this we performed a complete biochemical characterization of the enzyme's variants E172A, E172Q and E172S. Our steady state kinetic measurements show a clear decrease in the overall reaction rate, and our substrate interaction studies indicate the binding of the two substrates is also affected by these mutations. Interestingly our fast kinetic results show quinone reduction is more affected than NADH oxidation. We have also determined the X-ray crystal structure of the E172S mutant (2.55Ǻ) and compared it with the structure of the wild type (2.32Ǻ). Together these results support our hypothesis for E172 being of central importance in the catalytic mechanism of NDH-2, which may be extended to other members of the tDBDF superfamily.

Structural and Functional insights into the catalytic mechanism of the Type II NADH:quinone oxidoreductase family

Type II NADH:quinone oxidoreductases (NDH-2s) are membrane proteins involved in respiratory chains. These proteins contribute indirectly to the establishment of the transmembrane difference of electrochemical potential by catalyzing the reduction of quinone by oxidation of NAD(P)H. NDH-2s are widespread enzymes being present in the three domains of life. In this work, we explored the catalytic mechanism of NDH-2 by investigating the common elements of all NDH-2s, based on the rationale that conservation of such elements reflects their structural/functional importance. We observed conserved sequence motifs and structural elements among 1762 NDH-2s. We identified two proton pathways possibly involved in the protonation of the quinone. Our results led us to propose the first catalytic mechanism for NDH-2 family, in which a conserved glutamate residue, E₁₇₂ (in NDH-2 from Staphylococcus aureus) plays a key role in proton transfer to the quinone pocket. This catalytic mechanism may also be extended to the other members of the two-Dinucleotide Binding Domains Flavoprotein (tDBDF) superfamily, such as sulfide:quinone oxidoreductases.

Effect of a pH Gradient on the Protonation States of Cytochrome c Oxidase: A Continuum Electrostatics Study

Cytochrome c oxidase (CcO) couples the reduction of dioxygen to water with transmembrane proton pumping, which leads to the generation of an electrochemical gradient. In this study we analyze how one of the components of the electrochemical gradient, the difference in pH across the membrane, or ΔpH, influences the protonation states of residues in CcO. We modified our continuum electrostatics/Monte Carlo (CE/MC) method in order to include the ΔpH and applied it to the study of CcO, in what is, to our best knowledge, the first CE/MC study of CcO in the presence of a pH gradient. The inclusion of a transmembrane pH gradient allows for the identification of residues whose titration behavior depends on the pH on both sides of the membrane. Among the several residues with unusual titration profiles, three are well-known key residues in the proton transfer process of CcO: E286_I, Y288_I, and K362_I. All three residues have been previously identified as being critical for the catalytic or proton pumping functions of CcO. Our results suggest that when the pH gradient increases, these residues may be part of a regulatory mechanism to stem the proton flow.

The role of electrostatics in TrxR electron transfer mechanism: A computational approach

Thioredoxin reductase (TrxR) is an important enzyme in the control of the intracellular reduced redox environment. It transfers electrons from NADPH to several molecules, including its natural partner, thioredoxin. Although there is a generally accepted model describing how the electrons are transferred along TrxR, which involves a flexible arm working as a "shuttle," the molecular details of such mechanism are not completely understood. In this work, we use molecular dynamics simulations with Poisson-Boltzmann/Monte Carlo pKa calculations to investigate the role of electrostatics in the electron transfer mechanism. We observed that the combination of redox/protonation states of the N-terminal (FAD and Cys59/64) and C-terminal (Cys497/Selenocysteine498) redox centers defines the preferred relative positions and allows for the flexible arm to work as the desired "shuttle." Changing the redox/ionization states of those key players, leads to electrostatic triggers pushing the arm into the pocket when oxidized, and pulling it out, once it has been reduced. The calculated pKa values for Cys497 and Selenocysteine498 are 9.7 and 5.8, respectively, confirming that the selenocysteine is indeed deprotonated at physiological pH. This can be an important advantage in terms of reactivity (thiolate/selenolate are more nucleophilic than thiol/selenol) and ability to work as an electrostatic trigger (the "shuttle" mechanism) and may be the reason why TrxR uses selenium instead of sulfur. Proteins 2016; 84:1836-1843. © 2016 Wiley Periodicals, Inc.

Structure of Escherichia coli Flavodiiron Nitric Oxide Reductase

Flavodiiron proteins (FDPs) are present in organisms from all domains of life and have been described so far to be involved in the detoxification of oxygen or nitric oxide (NO), acting as O₂ and/or NO reductases. The Escherichia coli FDP, named flavorubredoxin (FlRd), is the most extensively studied FDP. Biochemical and in vivo studies revealed that FlRd is involved in NO detoxification as part of the bacterial defense mechanisms against reactive nitrogen species. E. coli FlRd has a clear preference for NO as a substrate in vitro, exhibiting a very low reactivity toward O₂. To contribute to the understanding of the structural features defining this substrate selectivity, we determined the crystallographic structure of E. coli FlRd, both in the isolated and reduced states. The overall tetrameric structure revealed a highly conserved flavodiiron core domain, with a metallo-β-lactamase-like domain containing a diiron center, and a flavodoxin domain with a flavin mononucleotide cofactor. The metal center in the oxidized state has a μ-hydroxo bridge coordinating the two irons, while in the reduced state, this moiety is not detected. Since only the flavodiiron domain was observed in these crystal structures, the structure of the rubredoxin domain was determined by NMR. Tunnels for the substrates were identified, and through molecular dynamics simulations, no differences for O₂ or NO permeation were found. The present data represent the first structure for a NO-selective FDP.

Coupling between protonation and conformation in cytochrome c oxidase: Insights from constant-pH MD simulations

Cytochrome c oxidases (CcOs) are the terminal enzymes of the respiratory chain in mitochondria and most bacteria. These enzymes reduce dioxygen (O(2)) to water and, simultaneously, generate a transmembrane electrochemical proton gradient. Despite their importance in the aerobic metabolism and the large amount of structural and biochemical data available for the A1-type CcO family, there is still no consensually accepted description of the molecular mechanisms operating in this protein. A substantial number of questions about the CcO's working mechanism remain to be answered, including how the protonation behavior of some key residues is modulated during a reduction cycle and how is the conformation of the protein affected by protonation. The main objective of this work was to study the protonation-conformation coupling in CcOs and identify the molecular factors that control the protonation state of some key residues. In order to directly capture the interplay between protonation and conformational effects, we have performed constant-pH MD simulations of an A1-type CcO inserted into a lipid bilayer in two redox states (oxidized and reduced) at physiological pH. From the simulations, we were able to identify several groups with unusual titration behavior that are highly dependent on the protein redox state, including the A-propionate from heme a and the D-propionate from heme a3, two key groups possibly involved in proton pumping. The protonation state of these two groups is heavily influenced by subtle conformational changes in the protein (notably of R481(I) and R482(I)) and by small changes in the hydrogen bond network.

Constant-pH Molecular Dynamics Study of Kyotorphin in an Explicit Bilayer

To our knowledge, we present the first constant-pH molecular dynamics study of the neuropeptide kyotorphin in the presence of an explicit lipid bilayer. The overall conformation freedom of the peptide was found to be affected by the interaction with the membrane, in accordance with previous results using different methodologies. Analysis of the interactions between the N-terminus amine group of the peptide and several lipid atoms shows that the membrane is able to stabilize both ionized and neutral forms of kyotorphin, resulting in a pKa value that is similar to the one obtained in water. This illustrates how a detailed molecular model of the membrane leads to rather different results than would be expected from simply regarding it as a low-dielectric slab.

Constant-pH MD Simulations Portray the Protonation and Structural Behavior of Four Decapeptides Designed to Coordinate Cu(2+)

The cyclic decapeptide C-Asp, containing one Asp residue and three His residues, was designed by Fragoso et al. (Chem. Eur. J. 2013, 19, 2076) to bind Cu(2+) exclusively through the side chain groups and mimic copper coordination in metalloproteins. A variant of the cyclodecapeptide where Asp is substituted by Asn (C-Asn) has also been synthesized in addition to the linear ("open") counterparts of both forms (O-Asp and O-Asn), testing the importance of cyclization and the presence of Asp in Cu(2+) coordination (Chem. Eur. J. 2013, 19, 2076; Dalton Trans. 2013, 42, 6182). All peptides formed a major species at neutral pH that was able to coordinate Cu(2+) exclusively through the neutral imidazole groups and the Asp side chain, when present, with C-Asp being the most effective. A detailed description of the protonation behavior of each histidine could help understanding the coordination species being formed in the pH range and eventually further optimizing the peptide's design. However, the standard current methods (NMR titrations) are not very suited for proximal groups titrating in the same pH range. In this work, we used the stochastic titration constant-pH molecular dynamics method to calculate the protonation curves and pKa of each titrable residue in the four decapeptides, in the absence of Cu(2+) ions. The global protonation curves obtained in our simulations are in very good agreement with the existing potentiometric titration curves. The histidines are titrating very closely, and the Asp forms abundant salt bridges with the basic residues, displaying an unusually low pKa value. In addition, we could observe that the four peptides are very unstructured in the absence of copper, and not even the cyclic forms exhibit a significant β-sheet, unlike what could be expected from the presence of β-turn inducer units in this type of scaffold.

The Effect of Membrane Environment on Surfactant Protein C Stability Studied by Constant-pH Molecular Dynamics

Pulmonary surfactant protein C (SP-C) is a small peptide with two covalently linked fatty acyl chains that plays a crucial role in the formation and stabilization of the pulmonary surfactant reservoirs during the compression and expansion steps of the respiratory cycle. Although its function is known to be tightly related to its highly hydrophobic character and key interactions maintained with specific lipid components, much is left to understand about its molecular mechanism of action. Also, although it adopts a mainly helical structure while associated with the membrane, factors as pH variation and deacylation have been shown to affect its stability and function. In this work, the conformational behavior of both the acylated and deacylated SP-C isoforms was studied in a DPPC bilayer under different pH conditions using constant-pH molecular dynamics simulations. Our findings show that both protein isoforms are remarkably stable over the studied pH range, even though the acylated isoform exhibits a labile helix-turn-helix motif rarely observed in the other isoform. We estimate similar tilt angles for the two isoforms over the studied pH range, with a generally higher degree of internalization of the basic N-terminal residues in the deacylated case, and observe and discuss some protonation-conformation coupling effects. Both isoforms establish contacts with the surrounding lipid molecules (preferentially with the sn-2 ester bonds) and have a local effect on the conformational behavior of the surrounding lipid molecules, the latter being more pronounced for acylated SP-C.

Understanding and Manipulating Electrostatic Fields at the Protein-Protein Interface Using Vibrational Spectroscopy and Continuum Electrostatics Calculations

Biological function emerges in large part from the interactions of biomacromolecules in the complex and dynamic environment of the living cell. For this reason, macromolecular interactions in biological systems are now a major focus of interest throughout the biochemical and biophysical communities. The affinity and specificity of macromolecular interactions are the result of both structural and electrostatic factors. Significant advances have been made in characterizing structural features of stable protein-protein interfaces through the techniques of modern structural biology, but much less is understood about how electrostatic factors promote and stabilize specific functional macromolecular interactions over all possible choices presented to a given molecule in a crowded environment. In this Feature Article, we describe how vibrational Stark effect (VSE) spectroscopy is being applied to measure electrostatic fields at protein-protein interfaces, focusing on measurements of guanosine triphosphate (GTP)-binding proteins of the Ras superfamily binding with structurally related but functionally distinct downstream effector proteins. In VSE spectroscopy, spectral shifts of a probe oscillator's energy are related directly to that probe's local electrostatic environment. By performing this experiment repeatedly throughout a protein-protein interface, an experimental map of measured electrostatic fields generated at that interface is determined. These data can be used to rationalize selective binding of similarly structured proteins in both in vitro and in vivo environments. Furthermore, these data can be used to compare to computational predictions of electrostatic fields to explore the level of simulation detail that is necessary to accurately predict our experimental findings.

Why a diaminopyrrolic tripodal receptor binds mannosides in acetonitrile but not in water?

Intermolecular interactions involving carbohydrates and their natural receptors play important roles in several biological processes. The development of synthetic receptors is very useful to study these recognition processes. Recently, it was synthetized a diaminopyrrolic tripodal receptor that is selective for mannosides, which are obtained from mannose, a sugar with significant relevance in living systems. However, this receptor is significantly more active in acetonitrile than in water. In this work, we performed several molecular dynamics and constant-pH molecular dynamics simulations in acetonitrile and water to evaluate the conformational space of the receptor and to understand the molecular detail of the receptor–mannoside interaction. The protonation states sampled by the receptor show that the positive charges are always as distant as possible in order to avoid large intramolecular repulsions. Moreover, the conformational space of the receptor is very similar in water above pH 4.0 and in acetonitrile. From the simulations with the mannoside, we observe that the interactions are more specific in acetonitrile (mainly hydrogen bonds) than in water (mainly hydrophobic). Our results suggest that the readiness of the receptor to bind mannoside is not significantly affected in water (above pH 4.0). Probably, the hydrogen bond network that is formed in acetonitrile (which is weaker in water) is the main reason for the higher activity in this solvent. This work also presents a new implementation of the stochastic titration constant-pH molecular dynamics method to a synthetic receptor of sugars and attests its ability to describe the protonation/conformation coupling in these molecules.

A simulated intermediate state for folding and aggregation provides insights into ΔN6 β2-microglobulin amyloidogenic behavior

A major component of ex vivo amyloid plaques of patients with dialysis-related amyloidosis (DRA) is a cleaved variant of β₂-microglobulin (ΔN6) lacking the first six N-terminal residues. Here we perform a computational study on ΔN6, which provides clues to understand the amyloidogenicity of the full-length β₂-microglobulin. Contrary to the wild-type form, ΔN6 is able to efficiently nucleate fibrillogenesis in vitro at physiological pH. This behavior is enhanced by a mild acidification of the medium such as that occurring in the synovial fluid of DRA patients. Results reported in this work, based on molecular simulations, indicate that deletion of the N-terminal hexapeptide triggers the formation of an intermediate state for folding and aggregation with an unstructured strand A and a native-like core. Strand A plays a pivotal role in aggregation by acting as a sticky hook in dimer assembly. This study further predicts that the detachment of strand A from the core is maximized at pH 6.2 resulting into higher aggregation efficiency. The structural mapping of the dimerization interface suggests that Tyr10, His13, Phe30 and His84 are hot-spot residues in ΔN6 amyloidogenesis.

Dialysis-related amyloidosis (DRA) is a conformational disease that affects individuals undergoing long-term haemodialysis. In DRA the progressive accumulation of protein human β₂-microglobulin (Hβ₂m) in the osteoarticular system, followed by its assembly into amyloid fibrils, eventually leads to tissue erosion and destruction. Disclosing the aggregation mechanism of Hβ₂m under physiologically relevant conditions represents a major challenge due to the inability of the protein to efficiently nucleate fibrillogenesis in vitro at physiological pH. On the other hand, ΔN6, a truncated variant of Hβ₂m, which is also a major component of ex vivo amyloid deposits extracted from DRA patients, is able to efficiently form amyloid fibrils de novo in physiological conditions. This amyloidogenic behavior is dramatically enhanced in a slightly more acidic pH (6.2) compatible with the mild acidification that occurs in the synovial fluid of DRA patients. In this work, an innovative three-stage methodological approach, relying on an array of molecular simulations, spanning different levels of resolution is used to investigate the initial stage of the de novo aggregation mechanism of ΔN6 in a physiologically relevant pH range. We identify an intermediate state for folding and aggregation, whose potential to dimerize is enhanced at pH 6.2. Our results provide rationalizations for previous experimental observations and new insights into the molecular basis of DRA.