Applications of MMPBSA to Membrane Proteins I: Efficient Numerical Solutions of Periodic Poisson-Boltzmann Equation

Inhibitory Potential of the Truncated Isoforms on Glutamate Transporter Oligomerization Identified by Computational Analysis of Gene-Centric Isoform Maps

OBJECTIVE

Glutamate transporters play a key role in central nervous system physiology by maintaining excitatory neurotransmitter homeostasis. Biological assemblies of the transporters, consisting of cyclic homotrimers, emerge as a crucial aspect of glutamate transporter modulation. Hence targeting heteromerization promises an effective approach for modulator design. On the other hand, the dynamic nature of transcription allows for the generation of transporter isoforms in structurally distinct manners.

METHODS

The potential isoforms were identified through the analysis of computationally generated gene-centric isoform maps. The conserved features of isoform sequences were revealed by computational chemistry methods and subsequent structural analysis of AlphaFold2 predictions. Truncated isoforms were further subjected to a wide range of docking analyses, 50ns molecular dynamics simulations, and evolutionary coupling analyses.

RESULTS

Energetic landscapes of isoform-canonical transporter complexes suggested an inhibitory potential of truncated isoforms on glutamate transporter bio-assembly. Moreover, isoforms that mimic the trimerization domain (in particular, TM2 helices) exhibited stronger interactions with canonical transporters, underscoring the role of transmembrane helices in isoform interactions. Additionally, self-assembly dynamics observed in truncated isoforms mimicking canonical TM5 helices indicate a potential protective role against unwanted interactions with canonical transporters.

CONCLUSION

Our computational studies on glutamate transporters offer insights into the roles of alternative splicing on protein interactions and identifies potential drug targets for physiological or pathological processes.

Identification of Gα12-vs-Gα13-coupling determinants and development of a Gα12/13-coupled designer GPCR

G-protein-coupled receptors (GPCRs) transduce diverse signals into the cell by coupling to one or several Gα subtypes. Of the 16 Gα subtypes in human cells, Gα12 and Gα13 belong to the G12 subfamily and are reported to be functionally different. Notably, certain GPCRs display selective coupling to either Gα12 or Gα13, highlighting their significance in various cellular contexts. However, the structural basis underlying this selectivity remains unclear. Here, using a Gα12-coupled designer receptor exclusively activated by designer drugs (DREADD; G12D) as a model system, we identified residues in the α5 helix and the receptor that collaboratively determine Gα12-vs-Gα13 selectivity. Residue-swapping experiments showed that G12D distinguishes differences between Gα12 and Gα13 in the positions G.H5.09 and G.H5.23 in the α5 helix. Molecular dynamics simulations observed that I378G.H5.23 in Gα12 interacts with N1032.39, S1693.53 and Y17634.53 in G12D, while H364G.H5.09 in Gα12 interact with Q2645.71 in G12D. Screening of mutations at these positions in G12D identified G12D mutants that enhanced coupling with Gα12 and to an even greater extent with Gα13. Combined mutations, most notably the dual Y17634.53H and Q2645.71R mutant, further enhanced Gα12/13 coupling, thereby serving as a potential Gα12/13-DREADD. Such novel Gα12/13-DREADD may be useful in future efforts to develop drugs that target Gα12/13 signaling as well as to identify their therapeutic indications.

Inhibition of polymorphic MexXY-OprM efflux system in Pseudomonas aeruginosa clinical isolates by Berberine derivatives

The MexXY-OprM multidrug efflux pump (EP) in aminoglycosides resistant Pseudomonas aeruginosa is one of the major resistance mechanisms, which is often overexpressed in strains isolated from pulmonary chronic disease such as cystic fibrosis.[1-3] In this research, we focused on the design of potential efflux pumps inhibitors, targeting MexY, the inner membrane component, in an allosteric site. Berberine[4] has been considered as lead molecule since we previously demonstrated its effectiveness in targeting MexY in laboratory reference strains.[5,6] Since this protein is often present in polymorphic variants in clinical strains, we sequenced and modeled all the mutated forms and we synthesized and evaluated by computational techniques, some berberine derivatives carrying an aromatic functionalization in its 13-C ring position. These compounds were tested in vitro against clinical P. aeruginosa strains for antimicrobial and antibiofilm activity. In conclusion, the results demonstrated the importance of the aromatic moiety functionalization in exerting the EP inhibitory activity in synergy with the aminoglycoside tobramycin. More, we found that aminoacidic composition of MexY in different strains must be considered for predicting potential binding site and affects the different activity of berberine derivatives. Finally, the antibiofilm effect of these new EPIs is promising, particularly for o-CH3-berberine derivative.

Study of SQ109 analogs binding to mycobacterium MmpL3 transporter using MD simulations and alchemical relative binding free energy calculations

N-geranyl-N΄-(2-adamantyl)ethane-1,2-diamine (SQ109) is a tuberculosis drug that has high potency against Mycobacterium tuberculosis (Mtb) and may function by blocking cell wall biosynthesis. After the crystal structure of MmpL3 from Mycobacterium smegmatis in complex with SQ109 became available, it was suggested that SQ109 inhibits Mmpl3 mycolic acid transporter. Here, we showed using molecular dynamics (MD) simulations that the binding profile of nine SQ109 analogs with inhibitory potency against Mtb and alkyl or aryl adducts at C-2 or C-1 adamantyl carbon to MmpL3 was consistent with the X-ray structure of MmpL3 - SQ109 complex. We showed that rotation of SQ109 around carbon-carbon bond in the monoprotonated ethylenediamine unit favors two gauche conformations as minima in water and lipophilic solvent using DFT calculations as well as inside the transporter's binding area using MD simulations. The binding assays in micelles suggested that the binding affinity of the SQ109 analogs was increased for the larger, more hydrophobic adducts, which was consistent with our results from MD simulations of the SQ109 analogues suggesting that sizeable C-2 adamantyl adducts of SQ109 can fill a lipophilic region between Y257, Y646, F260 and F649 in MmpL3. This was confirmed quantitatively by our calculations of the relative binding free energies using the thermodynamic integration coupled with MD simulations method with a mean assigned error of 0.74 kcal mol^-1 compared to the experimental values.

SARS-CoV-2 spike evolutionary behaviors; simulation of N501Y mutation outcomes in terms of immunogenicity and structural characteristic

Since the emergence of severe acute respiratory syndrome coronavirus 2 (SARS‐CoV‐2), a large number of mutations in its genome have been reported. Some of the mutations occur in noncoding regions without affecting the pathobiology of the virus, while mutations in coding regions are significant. One of the regions where a mutation can occur, affecting the function of the virus is at the receptor‐binding domain (RBD) of the spike protein. RBD interacts with angiotensin‐converting enzyme 2 (ACE2) and facilitates the entry of the virus into the host cells. There is a lot of focus on RBD mutations, especially the displacement of N501Y which is observed in the UK/Kent, South Africa, and Brazilian lineages of SARS‐CoV‐2. Our group utilizes computational biology approaches such as immunoinformatics, protein–protein interaction analysis, molecular dynamics, free energy computation, and tertiary structure analysis to disclose the consequences of N501Y mutation at the molecular level. Surprisingly, we discovered that this mutation reduces the immunogenicity of the spike protein; also, displacement of Asn with Tyr reduces protein compactness and significantly increases the stability of the spike protein and its affinity to ACE2. Moreover, following the N501Y mutation secondary structure and folding of the spike protein changed dramatically.

Estimating the Roles of Protonation and Electronic Polarization in Absolute Binding Affinity Simulations

Accurate prediction of binding free energies is critical to streamlining the drug development and protein design process. With the advent of GPU acceleration, absolute alchemical methods, which simulate the removal of ligand electrostatics and van der Waals interactions with the protein, have become routinely accessible and provide a physically rigorous approach that enables full consideration of flexibility and solvent interaction. However, standard explicit solvent simulations are unable to model protonation or electronic polarization changes upon ligand transfer from water to the protein interior, leading to inaccurate prediction of binding affinities for charged molecules. Here, we perform extensive simulation totaling ∼540 μs to benchmark the impact of modeling conditions on predictive accuracy for absolute alchemical simulations. Binding to urokinase plasminogen activator (UPA), a protein frequently overexpressed in metastatic tumors, is evaluated for a set of 10 inhibitors with extended flexibility, highly charged character, and titratable properties. We demonstrate that the alchemical simulations can be adapted to utilize the MBAR/PBSA method to improve the accuracy upon incorporating electronic polarization, highlighting the importance of polarization in alchemical simulations of binding affinities. Comparison of binding energy prediction at various protonation states indicates that proper electrostatic setup is also crucial in binding affinity prediction of charged systems, prompting us to propose an alternative binding mode with protonated ligand phenol and Hid-46 at the binding site, a testable hypothesis for future experimental validation.

Elucidation of WW domain ligand binding specificities in the Hippo pathway reveals STXBP4 as YAP inhibitor

The Hippo pathway, which plays a critical role in organ size control and cancer, features numerous WW domain-based protein-protein interactions. However, ~100 WW domains and 2,000 PY motif-containing peptide ligands are found in the human proteome, raising a "WW-PY" binding specificity issue in the Hippo pathway. In this study, we have established the WW domain binding specificity for Hippo pathway components and uncovered a unique amino acid sequence required for it. By using this criterion, we have identified a WW domain-containing protein, STXBP4, as a negative regulator of YAP. Mechanistically, STXBP4 assembles a protein complex comprising α-catenin and a group of Hippo PY motif-containing components/regulators to inhibit YAP, a process that is regulated by actin cytoskeleton tension. Interestingly, STXBP4 is a potential tumor suppressor for human kidney cancer, whose downregulation is correlated with YAP activation in clear cell renal cell carcinoma. Taken together, our study not only elucidates the WW domain binding specificity for the Hippo pathway, but also reveals STXBP4 as a player in actin cytoskeleton tension-mediated Hippo pathway regulation.

Improved Poisson-Boltzmann Methods for High-Performance Computing

Implicit solvent models based on the Poisson-Boltzmann equation (PBE) have been widely used to study electrostatic interactions in biophysical processes. These models often treat the solvent and solute regions as high and low dielectric continua, leading to a large jump in dielectrics across the molecular surface which is difficult to handle. Higher order interface schemes are often needed to seek higher accuracy for PBE applications. However, these methods are usually very liberal in the use of grid points nearby the molecular surface, making them difficult to use on high-performance computing platforms. Alternatively, the harmonic average (HA) method has been used to approximate dielectric interface conditions near the molecular surface with surprisingly good convergence and is well suited for high-performance computing. By adopting a 7-point stencil, the HA method is advantageous in generating simple 7-banded coefficient matrices, which greatly facilitate linear system solution with dense data parallelism, on high-performance computing platforms such as a graphics processing unit (GPU). However, the HA method is limited due to its lower accuracy. Therefore, it would be of great interest for high-performance applications to develop more accurate methods while retaining the simplicity and effectiveness of the 7-point stencil discretization scheme. In this study, we have developed two new algorithms based on the spirit of the HA method by introducing more physical interface relations and imposing the discretized Poisson's equation to the second order, respectively. Our testing shows that, for typical biomolecules, the new methods significantly improve the numerical accuracy to that comparable to the second-order solvers and with ∼65% overall efficiency gain on widely available high-performance GPU platforms.

Understanding G Protein Selectivity of Muscarinic Acetylcholine Receptors Using Computational Methods

The neurotransmitter molecule acetylcholine is capable of activating five muscarinic acetylcholine receptors, M1 through M5, which belong to the superfamily of G-protein-coupled receptors (GPCRs). These five receptors share high sequence and structure homology; however, the M1, M3, and M5 receptor subtypes signal preferentially through the Gαq/11 subset of G proteins, whereas the M2 and M4 receptor subtypes signal through the Gαi/o subset of G proteins, resulting in very different intracellular signaling cascades and physiological effects. The structural basis for this innate ability of the M1/M3/M5 set of receptors and the highly homologous M2/M4 set of receptors to couple to different G proteins is poorly understood. In this study, we used molecular dynamics (MD) simulations coupled with thermodynamic analyses of M1 and M2 receptors coupled to both Gαi and Gαq to understand the structural basis of the M1 receptor’s preference for the Gαq protein and the M2 receptor’s preference for the Gαi protein. The MD studies showed that the M1 and M2 receptors can couple to both Gα proteins such that the M1 receptor engages with the two Gα proteins in slightly different orientations and the M2 receptor engages with the two Gα proteins in the same orientation. Thermodynamic studies of the free energy of binding of the receptors to the Gα proteins showed that the M1 and M2 receptors bind more strongly to their cognate Gα proteins compared to their non-cognate ones, which is in line with previous experimental studies on the M3 receptor. A detailed analysis of receptor–G protein interactions showed some cognate-complex-specific interactions for the M2:Gαi complex; however, G protein selectivity determinants are spread over a large overlapping subset of residues. Conserved interaction between transmembrane helices 5 and 6 far away from the G-protein-binding receptor interface was found only in the two cognate complexes and not in the non-cognate complexes. An analysis of residues implicated previously in G protein selectivity, in light of the cognate and non-cognate structures, shaded a more nuanced role of those residues in affecting G protein selectivity. The simulation of both cognate and non-cognate receptor–G protein complexes fills a structural gap due to difficulties in determining non-cognate complex structures and provides an enhanced framework to probe the mechanisms of G protein selectivity exhibited by most GPCRs.

NFBTA: A Potent Cytotoxic Agent against Glioblastoma

Piplartine (PPL), also known as piperlongumine, is a biologically active alkaloid extracted from the Piper genus which has been found to have highly effective anticancer activity against several tumor cell lines. This study investigates in detail the antitumoral potential of a PPL analogue; (E)-N-(4-fluorobenzyl)-3-(3,4,5-trimethoxyphenyl) acrylamide (NFBTA). The anticancer potential of NFBTA on the glioblastoma multiforme (GBM) cell line (U87MG) was determined by 3-(4,5-dimethyl-2-thia-zolyl)-2, 5-diphenyl-2H-tetrazolium bromide (MTT), and lactate dehydrogenase (LDH) release analysis, and the selectivity index (SI) was calculated. To detect cell apoptosis, fluorescent staining via flow cytometry and Hoechst 33258 staining were performed. Oxidative alterations were assessed via colorimetric measurement methods. Alterations in expressions of key genes related to carcinogenesis were determined. Additionally, in terms of NFBTA cytotoxic, oxidative, and genotoxic damage potential, the biosafety of this novel agent was evaluated in cultured human whole blood cells. Cell viability analyses revealed that NFBTA exhibited strong cytotoxic activity in cultured U87MG cells, with high selectivity and inhibitory activity in apoptotic processes, as well as potential for altering the principal molecular genetic responses in U87MG cell growth. Molecular docking studies strongly suggested a plausible anti-proliferative mechanism for NBFTA. The results of the experimental in vitro human glioblastoma model and computational approach revealed promising cytotoxic activity for NFBTA, helping to orient further studies evaluating its antitumor profile for safe and effective therapeutic applications.

Heterogeneous Dielectric Implicit Membrane Model for the Calculation of MMPBSA Binding Free Energies

Membrane-bound protein receptors are a primary biological drug target, but the computational analysis of membrane proteins has been limited. In order to improve molecular mechanics Poisson-Boltzmann surface area (MMPBSA) binding free energy calculations for membrane protein-ligand systems, we have optimized a new heterogeneous dielectric implicit membrane model, with respect to free energy simulations in explicit membrane and explicit water, and implemented it into the Amber software suite. This new model supersedes our previous uniform, single dielectric implicit membrane model by allowing the dielectric constant to vary with depth within the membrane. We calculated MMPBSA binding free energies for the human purinergic platelet receptor (P2Y₁₂R) and two of the muscarinic acetylcholine receptors (M2R and M3R) bound to various antagonist ligands using both membrane models, and we found that the heterogeneous dielectric membrane model has a stronger correlation with experimental binding affinities compared to the older model under otherwise identical conditions. This improved membrane model increases the utility of MMPBSA calculations for the rational design and improvement of future drug candidates.

Robustness and Efficiency of Poisson-Boltzmann Modeling on Graphics Processing Units

Poisson-Boltzmann equation (PBE) based continuum electrostatics models have been widely used in modeling electrostatic interactions in biochemical processes, particularly in estimating protein-ligand binding affinities. Fast convergence of PBE solvers is crucial in binding affinity computations as numerous snapshots need to be processed. Efforts have been reported to develop PBE solvers on graphics processing units (GPUs) for efficient modeling of biomolecules, though only relatively simple successive over-relaxation and conjugate gradient methods were implemented. However, neither convergence nor scaling properties of the two methods are optimal for large biomolecules. On the other hand, geometric multigrid (MG) has been shown to be an optimal solver on CPUs, though no MG have been reported for biomolecular applications on GPUs. This is not a surprise as it is a more complex method and depends on simpler but limited iterative methods such as Gauss-Seidel in its core relaxation procedure. The robustness and efficiency of MG on GPUs are also unclear. Here we present an implementation and a thorough analysis of MG on GPUs. Our analysis shows that robustness is a more pronounced issue than efficiency for both MG and other tested solvers when the single precision is used for complex biomolecules. We further show how to balance robustness and efficiency utilizing MG's overall efficiency and conjugate gradient's robustness, pointing to a hybrid GPU solver with a good balance of efficiency and accuracy. The new PBE solver will significantly improve the computational throughput for a range of biomolecular applications on the GPU platforms.

Computational Analysis for the Rational Design of Anti-Amyloid Beta (Aβ) Antibodies

Alzheimer's disease (AD) is a neurodegenerative disorder that lacks effective treatment options. Anti-amyloid beta (Aβ) antibodies are the leading drug candidates to treat AD, but the results of clinical trials have been disappointing. Introducing rational mutations into anti-Aβ antibodies to increase their effectiveness is a way forward, but the path to take is unclear. In this study, we demonstrate the use of computational fragment-based docking and MMPBSA binding free energy calculations in the analysis of anti-Aβ antibodies for rational drug design efforts. Our fragment-based docking method successfully predicts the emergence of the common EFRH epitope. MD simulations coupled with MMPBSA binding free energy calculations are used to analyze scenarios described in prior studies, and we computationally introduce rational mutations into PFA1 to predict mutations that can improve its binding affinity toward the pE3-Aβ_3-8 form of Aβ. Two out of our four proposed mutations are predicted to stabilize binding. Our study demonstrates that a computational approach may lead to an improved drug candidate for AD in the future.

Molecular insights into substrate binding mechanism of undecaprenyl pyrophosphate with membrane integrated phosphatidyl glycerophosphate phosphatase B (PgpB) using molecular dynamics simulation approach

Undecaprenyl phosphate (C55-P) acts as carrier lipid in the synthesis of peptidoglycan, which is de novo synthesized from dephosphorylation of undecaprenyl pyrophosphate (C55-PP). The phosphatidylglycerol phosphate phosphatase B (PgpB) catalyzes the dephosphorylation of C55-PP and forms C55-P. As no structural study has been made regarding the binding of C55-PP to PgpB, in the current study, in silico molecular docking, followed by 150 ns molecular dynamics simulation of the putative binding complex in membrane/solvent environment has been performed to understand conformational dynamics. Results are compared with simulated apo form and PE inhibitor-bound form. Analysis of correlated residual fluctuation network in apo form, C55-PP bound and PE inhibitor-bound form suggests that difference in dynamic coupling between TM domain and α2 and α3 helix of periplasmic domain provides ligand binding to facilitate catalysis or to show inhibitory activity. Distance distribution in catalytic residual pair, H207-R104; H207-R201 and H207-D211 which stabilizes phosphate-enzyme intermediate shows a narrow peak in 2.4-3.6 Å in substrate-bound compared to apo form. Binding interactions and binding free energy analyses complement the partial inhibition of PE where PE has less binding free energy compared to the C55-PP substrate as well as the difference in binding interaction with catalytic pocket. Thus, the present study provides how substrate binding couples the movement in TM domain and periplasmic domain which might help in the understanding of active site communication in PgpB. C55-PP phosphatase interactions with a catalytic pocket of PgpB provide new insight for designing drugs against bacterial infection.

A Finite Element Solution of Lateral Periodic Poisson-Boltzmann Model for Membrane Channel Proteins

Membrane channel proteins control the diffusion of ions across biological membranes. They are closely related to the processes of various organizational mechanisms, such as: cardiac impulse, muscle contraction and hormone secretion. Introducing a membrane region into implicit solvation models extends the ability of the Poisson–Boltzmann (PB) equation to handle membrane proteins. The use of lateral periodic boundary conditions can properly simulate the discrete distribution of membrane proteins on the membrane plane and avoid boundary effects, which are caused by the finite box size in the traditional PB calculations. In this work, we: (1) develop a first finite element solver (FEPB) to solve the PB equation with a two-dimensional periodicity for membrane channel proteins, with different numerical treatments of the singular charges distributions in the channel protein; (2) add the membrane as a dielectric slab in the PB model, and use an improved mesh construction method to automatically identify the membrane channel/pore region even with a tilt angle relative to the z-axis; and (3) add a non-polar solvation energy term to complete the estimation of the total solvation energy of a membrane protein. A mesh resolution of about 0.25 Å (cubic grid space)/0.36 Å (tetrahedron edge length) is found to be most accurate in linear finite element calculation of the PB solvation energy. Computational studies are performed on a few exemplary molecules. The results indicate that all factors, the membrane thickness, the length of periodic box, membrane dielectric constant, pore region dielectric constant, and ionic strength, have individually considerable influence on the solvation energy of a channel protein. This demonstrates the necessity to treat all of those effects in the PB model for membrane protein simulations.

The Model Structures of the Complement Component 5a Receptor (C5aR) Bound to the Native and Engineered hC5a

The interaction of ^hC5a with C5aR, previously hypothesized to involve a “two-site” binding, (i) recognition of the bulk of ^hC5a by the N-terminus (NT) of C5aR (“site1”), and (ii) recognition of C-terminus (CT) of ^hC5a by the extra cellular surface (ECS) of the C5aR (“site2”). However, the pharmacological landscapes of such recognition sites are yet to be illuminated at atomistic resolution. In the context, unique model complexes of C5aR, harboring pharmacophores of diverse functionality at the “site2” has recently been described. The current study provides a rational illustration of the “two-site” binding paradigm in C5aR, by recruiting the native agonist ^hC5a and engineered antagonist ^hC5a(A8). The ^hC5a-C5aR and ^hC5a(A8)-C5aR complexes studied over 250 ns of molecular dynamics (MD) each in POPC bilayer illuminate the hallmark of activation mechanism in C5aR. The intermolecular interactions in the model complexes are well supported by the molecular mechanics Poisson–Boltzmann surface area (MM-PBSA) based binding free energy calculation, strongly correlating with the reported mutational studies. Exemplified in two unique and contrasting molecular complexes, the study provides an exceptional understanding of the pharmacological divergence observed in C5aR, which will certainly be useful for search and optimization of new generation “neutraligands” targeting the ^hC5a-C5aR interaction.

Recent Developments and Applications of the MMPBSA Method

The Molecular Mechanics Poisson-Boltzmann Surface Area (MMPBSA) approach has been widely applied as an efficient and reliable free energy simulation method to model molecular recognition, such as for protein-ligand binding interactions. In this review, we focus on recent developments and applications of the MMPBSA method. The methodology review covers solvation terms, the entropy term, extensions to membrane proteins and high-speed screening, and new automation toolkits. Recent applications in various important biomedical and chemical fields are also reviewed. We conclude with a few future directions aimed at making MMPBSA a more robust and efficient method.

Ionic Solution: What Goes Right and Wrong with Continuum Solvation Modeling

Solvent-mediated electrostatic interactions were well recognized to be important in the structure and function of molecular systems. Ionic interaction is an important component in electrostatic interactions, especially in highly charged molecules, such as nucleic acids. Here, we focus on the quality of the widely used Poisson-Boltzmann surface area (PBSA) continuum models in modeling ionic interactions by comparing with both explicit solvent simulations and the experiment. In this work, the molality-dependent chemical potentials for sodium chloride (NaCl) electrolyte were first simulated in the SPC/E explicit solvent. Our high-quality simulation agrees well with both the previous study and the experiment. Given the free-energy simulations in SPC/E as the benchmark, we used the same sets of snapshots collected in the SPC/E solvent model for PBSA free-energy calculations in the hope to achieve the maximum consistency between the two solvent models. Our comparative analysis shows that the molality-dependent chemical potentials of NaCl were reproduced well with both linear PB and nonlinear PB methods, although nonlinear PB agrees better with SPC/E and the experiment. Our free-energy simulations also show that the presence of salt increases the hydrophobic effect in a nonlinear fashion, in qualitative agreement with previous theoretical studies of Onsager and Samaras. However, the lack of molality-dependency in the nonelectrostatics continuum models dramatically reduces the overall quality of PBSA methods in modeling salt-dependent energetics. These analyses point to further improvements needed for more robust modeling of solvent-mediated interactions by the continuum solvation frameworks.

A Continuum Poisson-Boltzmann Model for Membrane Channel Proteins

Membrane proteins constitute a large portion of the human proteome and perform a variety of important functions as membrane receptors, transport proteins, enzymes, signaling proteins, and more. Computational studies of membrane proteins are usually much more complicated than those of globular proteins. Here, we propose a new continuum model for Poisson-Boltzmann calculations of membrane channel proteins. Major improvements over the existing continuum slab model are as follows: (1) The location and thickness of the slab model are fine-tuned based on explicit-solvent MD simulations. (2) The highly different accessibilities in the membrane and water regions are addressed with a two-step, two-probe grid-labeling procedure. (3) The water pores/channels are automatically identified. The new continuum membrane model is optimized (by adjusting the membrane probe, as well as the slab thickness and center) to best reproduce the distributions of buried water molecules in the membrane region as sampled in explicit water simulations. Our optimization also shows that the widely adopted water probe of 1.4 Å for globular proteins is a very reasonable default value for membrane protein simulations. It gives the best compromise in reproducing the explicit water distributions in membrane channel proteins, at least in the water accessible pore/channel regions. Finally, we validate the new membrane model by carrying out binding affinity calculations for a potassium channel, and we observe good agreement with the experimental results.

Acceleration of Linear Finite-Difference Poisson-Boltzmann Methods on Graphics Processing Units

Electrostatic interactions play crucial roles in biophysical processes such as protein folding and molecular recognition. Poisson-Boltzmann equation (PBE)-based models have emerged as widely used in modeling these important processes. Though great efforts have been put into developing efficient PBE numerical models, challenges still remain due to the high dimensionality of typical biomolecular systems. In this study, we implemented and analyzed commonly used linear PBE solvers for the ever-improving graphics processing units (GPU) for biomolecular simulations, including both standard and preconditioned conjugate gradient (CG) solvers with several alternative preconditioners. Our implementation utilizes the standard Nvidia CUDA libraries cuSPARSE, cuBLAS, and CUSP. Extensive tests show that good numerical accuracy can be achieved given that the single precision is often used for numerical applications on GPU platforms. The optimal GPU performance was observed with the Jacobi-preconditioned CG solver, with a significant speedup over standard CG solver on CPU in our diversified test cases. Our analysis further shows that different matrix storage formats also considerably affect the efficiency of different linear PBE solvers on GPU, with the diagonal format best suited for our standard finite-difference linear systems. Further efficiency may be possible with matrix-free operations and integrated grid stencil setup specifically tailored for the banded matrices in PBE-specific linear systems.

Numerical interpretation of molecular surface field in dielectric modeling of solvation

Continuum solvent models, particularly those based on the Poisson-Boltzmann equation (PBE), are widely used in the studies of biomolecular structures and functions. Existing PBE developments have been mainly focused on how to obtain more accurate and/or more efficient numerical potentials and energies. However to adopt the PBE models for molecular dynamics simulations, a difficulty is how to interpret dielectric boundary forces accurately and efficiently for robust dynamics simulations. This study documents the implementation and analysis of a range of standard fitting schemes, including both one-sided and two-sided methods with both first-order and second-order Taylor expansions, to calculate molecular surface electric fields to facilitate the numerical calculation of dielectric boundary forces. These efforts prompted us to develop an efficient approximated one-dimensional method, which is to fit the surface field one dimension at a time, for biomolecular applications without much compromise in accuracy. We also developed a surface-to-atom force partition scheme given a level set representation of analytical molecular surfaces to facilitate their applications to molecular simulations. Testing of these fitting methods in the dielectric boundary force calculations shows that the second-order methods, including the one-dimensional method, consistently perform among the best in the molecular test cases. Finally, the timing analysis shows the approximated one-dimensional method is far more efficient than standard second-order methods in the PBE force calculations. © 2017 Wiley Periodicals, Inc.

Crystal Structure of StnA for the Biosynthesis of Antitumor Drug Streptonigrin Reveals a Unique Substrate Binding Mode

Streptonigrin methylesterase A (StnA) is one of the tailoring enzymes that modify the aminoquinone skeleton in the biosynthesis pathway of Streptomyces species. Although StnA has no significant sequence homology with the reported α/β-fold hydrolases, it shows typical hydrolytic activity in vivo and in vitro. In order to reveal its functional characteristics, the crystal structures of the selenomethionine substituted StnA (SeMet-StnA) and the complex (S185A mutant) with its substrate were resolved to the resolution of 2.71 Å and 2.90 Å, respectively. The overall structure of StnA can be described as an α-helix cap domain on top of a common α/β hydrolase domain. The substrate methyl ester of 10'-demethoxystreptonigrin binds in a hydrophobic pocket that mainly consists of cap domain residues and is close to the catalytic triad Ser185-His349-Asp308. The transition state is stabilized by an oxyanion hole formed by the backbone amides of Ala102 and Leu186. The substrate binding appears to be dominated by interactions with several specific hydrophobic contacts and hydrogen bonds in the cap domain. The molecular dynamics simulation and site-directed mutagenesis confirmed the important roles of the key interacting residues in the cap domain. Structural alignment and phylogenetic tree analysis indicate that StnA represents a new subfamily of lipolytic enzymes with the specific binding pocket located at the cap domain instead of the interface between the two domains.

Modeling Membrane Protein-Ligand Binding Interactions: The Human Purinergic Platelet Receptor

Membrane proteins, due to their roles as cell receptors and signaling mediators, make prime candidates for drug targets. The computational analysis of protein-ligand binding affinities has been widely employed as a tool in rational drug design efforts. Although efficient implicit solvent-based methods for modeling globular protein-ligand binding have been around for many years, the extension of such methods to membrane protein-ligand binding is still in its infancy. In this study, we extended the widely used Amber/MMPBSA method to model membrane protein-ligand systems, and we used it to analyze protein-ligand binding for the human purinergic platelet receptor (P2Y12R), a prominent drug target in the inhibition of platelet aggregation for the prevention of myocardial infarction and stroke. The binding affinities, computed by the Amber/MMPBSA method using standard parameters, correlate well with experiment. A detailed investigation of these parameters was conducted to assess their impact on the accuracy of the method. These analyses show the importance of properly treating the nonpolar solvation interactions and the electrostatic polarization in the binding of nucleotide agonists and non-nucleotide antagonists to P2Y12R. On the basis of the crystal structures and the experimental conditions in the binding assay, we further hypothesized that the nucleotide agonists lose their bound magnesium ion upon binding to P2Y12R, and our computational study supports this hypothesis. Ultimately, this work illustrates the value of computational analysis in the interpretation of experimental binding reactions.