Molecular dynamics simulations of peptides and proteins with a continuum electrostatic model based on screened Coulomb potentials

An unusual dual sugar-binding lectin domain controls the substrate specificity of a mucin-type O-glycosyltransferase

N-acetylgalactosaminyl-transferases (GalNAc-Ts) initiate mucin-type O-glycosylation, an abundant and complex posttranslational modification that regulates host-microbe interactions, tissue development, and metabolism. GalNAc-Ts contain a lectin domain consisting of three homologous repeats (α, β, and γ), where α and β can potentially interact with O-GalNAc on substrates to enhance activity toward a nearby acceptor Thr/Ser. The ubiquitous isoenzyme GalNAc-T1 modulates heart development, immunity, and SARS-CoV-2 infectivity, but its substrates are largely unknown. Here, we show that both α and β in GalNAc-T1 uniquely orchestrate the O-glycosylation of various glycopeptide substrates. The α repeat directs O-glycosylation to acceptor sites carboxyl-terminal to an existing GalNAc, while the β repeat directs O-glycosylation to amino-terminal sites. In addition, GalNAc-T1 incorporates α and β into various substrate binding modes to cooperatively increase the specificity toward an acceptor site located between two existing O-glycans. Our studies highlight a unique mechanism by which dual lectin repeats expand substrate specificity and provide crucial information for identifying the biological substrates of GalNAc-T1.

Probing electrophysiological activity of amphiphilic Dynorphin A in planar neutral membranes reveals both ion channel-like activity and neuropeptide translocation

Dynorphin A (DynA) is an endogenous neuropeptide that besides acting as a ligand of the κ-opioid receptor, presents some non-opioid pathophysiological properties associated to its ability to induce cell permeability similarly to cell-penetrating peptides (CPPs). Here, we use electrophysiology experiments to show that amphiphilic DynA generates aqueous pores in neutral membranes similar to those reported previously in charged membranes, but we also find other events thermodynamically incompatible with voltage-driven ion channel activity (i.e. non-zero currents with no applied voltage in symmetric salt conditions, reversal potentials that exceed the theoretical limit for a given salt concentration gradient). By comparison with current traces generated by other amphiphilic molecule known to spontaneously cross membranes, we hypothesize that DynA could directly translocate across neutral bilayers, a feature never observed in charged membranes following the same electrophysiological protocol. Our findings suggest that DynA interaction with the cellular membrane is modulated by the lipid charge distribution, enabling either passive ionic transport via membrane remodeling and pore formation or by peptide direct internalization independent of cellular transduction pathways.

Computational Study of the Allosteric Effects of p5 on CDK5-p25 Hyperactivity as Alternative Inhibitory Mechanisms in Neurodegeneration

The cyclin-dependent kinase (CDK5) forms a stable complex with its activator p25, leading to the hyperphosphorylation of tau proteins and to the formation of plaques and tangles that are considered to be one of the typical causes of Alzheimer's disease (AD). Hence, the pathological CDK5-p25 complex is a promising therapeutic target for AD. Small peptides, obtained from the truncation of CDK5 physiological activator p35, have shown promise in inhibiting the pathological complex effectively while also crossing the blood-brain barrier. One such small 24-residue peptide, p5, has shown selective inhibition toward the pathological complex in vivo. Our previous research focused on the characterization of a computationally predicted CDK5-p5 binding mode and of its pharmacophore, which was consistent with competitive inhibition. In continuation of our previous work, herein, we investigate four additional binding modes to explore other possible mechanisms of interaction between CDK5 and p5. The quantitative description of the pharmacophore is consistent with both competitive and allosteric p5-induced inhibition mechanisms of CDK5-p25 pathology. The gained insights can direct further in vivo/in vitro tests and help design small peptides, linear or cyclic, or peptidomimetic compounds as adjuvants of orthosteric inhibitors or as part of a cocktail of drugs with enhanced effectiveness and lower side effects.

Modulation of free energy landscapes as a strategy for the design of antimicrobial peptides

Computational design of antimicrobial peptides (AMPs) is a promising area of research for developing novel agents against drug-resistant bacteria. AMPs are present naturally in many organisms, from bacteria to humans, a time-tested mechanism that makes them attractive as effective antibiotics. Depending on the environment, AMPs can exhibit α-helical or β-sheet conformations, a mix of both, or lack secondary structure; they can be linear or cyclic. Prediction of their structures is challenging but critical for rational design. Promising AMP leads can be developed using essentially two approaches: traditional modeling of the physicochemical mechanisms that determine peptide behavior in aqueous and membrane environments and knowledge-based, e.g., machine learning (ML) techniques, that exploit ever-growing AMP databases. Here, we explore the conformational landscapes of two recently ML-designed AMPs, characterize the dependence of these landscapes on the medium conditions, and identify features in peptide and membrane landscapes that mediate protein-membrane association. For both peptides, we observe greater conformational diversity in an aqueous solvent than in a less polar solvent, and one peptide is seen to alter its conformation more dramatically than the other upon the change of solvent. Our results support the view that structural rearrangement in response to environmental changes is central to the mechanism of membrane-structure disruption by linear peptides. We expect that the design of AMPs by ML will benefit from the incorporation of peptide conformational substates as quantified here with molecular simulations.

Blockade of GRTH/DDX25 Phosphorylation by Cyclic Peptides Provides an Avenue for Developing a Nonhormonal Male Contraceptive

Gonadotropin-regulated testicular RNA helicase (GRTH)/DDX25 is a DEAD-box RNA helicase essential for the completion of spermatogenesis. Our previous studies indicated that blocking the GRTH phospho-site or perturbing the GRTH/protein kinase A (PKA) interface could provide an avenue for developing a nonhormonal male contraceptive. In this study, cyclic peptides were rationally designed and synthesized as promising therapeutic agents. The peptides showed effective delivery into COS-1 and germ cells and a dose-dependent inhibitory effect on GRTH phosphorylation. The peptides inhibit GRTH phosphorylation in the presence of PKA, and binding to the helicase resulted in thermal stabilization of non-phospho GRTH. Increased efficiency in fluorescence resonance energy transfer (FRET) assay revealed their interaction with GRTH. Cyclic peptide exposure of cultures from mice seminiferous tubules resulted in significant inhibition of phospho GRTH. These peptides did not exhibit toxicity. Effective delivery and targeted decrease of in vitro expression of phospho GRTH by cyclic peptides provide a promising angle to develop effective compounds as a nonhormonal male contraceptive.

Adsorption of lysozyme on gold surfaces in the presence of an external electric potential

Adsorbed protein films consist of essential building blocks of many biotechnological and biomedical devices. The electrostatic potential may significantly modulate the protein behaviour on surfaces, affecting their structure and biological activity. In this study, lysozyme was used to investigate the effects of applied electric potentials on adsorption and the protein structure. The pH and the surface charge determine the amount and secondary structure of adsorbed lysozyme on a gold surface. In-situ measurements using polarization modulation infrared reflection absorption spectroscopy indicated that the concentration of both the adsorbed anions and the lysozyme led to conformational changes in the protein film, which was demonstrated by a greater amount of aggregated β-sheets in films fabricated at net positive charges of the Au electrode (E_ads > E_pzc). The changes in secondary structure involved two parallel processes. One comprised changes in the hydration/hydrogen-bond network at helices, leading to diverse helical structures: α-, 3₁₀- and/or π-helices. In the second process β-turns, β-sheets, and random coils displayed an ability to form aggregated β-sheet structures. The study illuminates the understanding of electrical potential-dependent changes involved in the protein misfolding process.

Reformulation of the self-guided molecular simulation method

Self-guided molecular/Langevin dynamics (SGMD/SGLD) simulation methods were developed to enhance conformational sampling through promoting low frequency motion of molecular systems and have been successfully applied in many simulation studies. Quantitative understanding of conformational distribution in SGLD has been achieved by separating microscopic properties according to frequency. However, a missing link between the guiding factors and conformational distributions makes it highly empirical and system dependent when choosing the values of the guiding parameters. Based on the understanding that molecular interactions are the source of energy barriers and diffusion friction, this work reformulates the equation of the low frequency motion to resemble Langevin dynamics. This reformulation leads to new forms of guiding forces and establishes a relation between the guiding factors and conformational distributions. We call simulations with these new guiding forces the generalized self-guided molecular/Langevin dynamics (SGMDg/SGLDg). In addition, we present a new way to calculate low frequency properties and an efficient algorithm to implement SGMDg/SGLDg that minimizes memory usage and inter-processor communication. Through example simulations with a skewed double well system, an argon fluid, and a cryo-EM map flexible fitting case, we demonstrate the guiding effects on conformational distributions and conformational searching.

The carbon monoxide dehydrogenase accessory protein CooJ is a histidine-rich multidomain dimer containing an unexpected Ni(II)-binding site

Activation of nickel enzymes requires specific accessory proteins organized in multiprotein complexes controlling metal transfer to the active site. Histidine-rich clusters are generally present in at least one of the metallochaperones involved in nickel delivery. The maturation of carbon monoxide dehydrogenase in the proteobacterium Rhodospirillum rubrum requires three accessory proteins, CooC, CooT, and CooJ, dedicated to nickel insertion into the active site, a distorted [NiFe₃S₄] cluster coordinated to an iron site. Previously, CooJ from R. rubrum (RrCooJ) has been described as a nickel chaperone with 16 histidines and 2 cysteines at its C terminus. Here, the X-ray structure of a truncated version of RrCooJ, combined with small-angle X-ray scattering data and a modeling study of the full-length protein, revealed a homodimer comprising a coiled coil with two independent and highly flexible His tails. Using isothermal calorimetry, we characterized several metal-binding sites (four per dimer) involving the His-rich motifs and having similar metal affinity (K_D = 1.6 μm). Remarkably, biophysical approaches, site-directed mutagenesis, and X-ray crystallography uncovered an additional nickel-binding site at the dimer interface, which binds Ni(II) with an affinity of 380 nm Although RrCooJ was initially thought to be a unique protein, a proteome database search identified at least 46 bacterial CooJ homologs. These homologs all possess two spatially separated nickel-binding motifs: a variable C-terminal histidine tail and a strictly conserved H(W/F)X ₂HX ₃H motif, identified in this study, suggesting a dual function for CooJ both as a nickel chaperone and as a nickel storage protein.

Peptide and Protein Structure Prediction with a Simplified Continuum Solvent Model

A continuum solvent model based on screened Coulomb potentials has been simplified and parametrized to sample native-like structures in replica-exchange simulations of each of six different peptides and miniproteins. Low-energy, native, and non-native structures were used to iteratively refine 11 parameter values. The centroid of the largest cluster of structures sampled in simulations initiated from an extended conformation represents the predicted structure. The main-chain rms deviation of this prediction from the experimental structure was 0.47 Å for the 12-residue Trp-zip2, 0.86 Å for the 14-residue MBH12, 2.53 Å for the 17-residue U(1-17)T9D, 2.03 Å for the 20-residue BS1, 1.08 Å for the 20-residue Trp-cage, and 3.64 Å for the 35-residue villin headpiece subdomain HP35. The centroid of the sixth largest cluster sampled for HP35 deviated by 0.91 Å. The CHARMM22/CMAP force field was used, with an additional ψ torsion term for residues other than glycine and proline. Six parameters govern the dielectric response of the continuum solvent, and four values of surface tension approximate nonpolar effects. An atom's self-energy and interaction energies are screened independently, each depending on whether the atom is part of a charged group, a neutral hydrogen-bonding main-chain group, or any other neutral group. The parameters inferred result in strong main-chain hydrogen bonds, consistent with the view that protein folding is dominated by the formation of these bonds. (1,2) Conformations of MBH12 and BS1 were excluded from the energy-function refinement, suggesting the parameters, referred to as SCP18, are transferable. An efficient estimate of solvent-accessible surface area is also described.

Incorporation of Solvent Effect into Multi-Objective Evolutionary Algorithm for Improved Protein Structure Prediction

The problem of predicting the three-dimensional (3-D) structure of a protein from its one-dimensional sequence has been called the "holy grail of molecular biology", and it has become an important part of structural genomics projects. Despite the rapid developments in computer technology and computational intelligence, it remains challenging and fascinating. In this paper, to solve it we propose a multi-objective evolutionary algorithm. We decompose the protein energy function Chemistry at HARvard Macromolecular Mechanics force fields into bond and non-bond energies as the first and second objectives. Considering the effect of solvent, we innovatively adopt a solvent-accessible surface area as the third objective. We use 66 benchmark proteins to verify the proposed method and obtain better or competitive results in comparison with the existing methods. The results suggest the necessity to incorporate the effect of solvent into a multi-objective evolutionary algorithm to improve protein structure prediction in terms of accuracy and efficiency.

Computational study of the inhibitory mechanism of the kinase CDK5 hyperactivity by peptide p5 and derivation of a pharmacophore

The hyperactivity of the cyclic dependent kinase 5 (CDK5) induced by the activator protein p25 has been linked to a number of pathologies of the brain. The CDK5-p25 complex has thus emerged as a major therapeutic target for Alzheimer's disease (AD) and other neurodegenerative conditions. Experiments have shown that the peptide p5 reduces the CDK5-p25 activity without affecting the endogenous CDK5-p35 activity, whereas the peptide TFP5, obtained from p5, elicits similar inhibition, crosses the blood-brain barrier, and exhibits behavioral rescue of AD mice models with no toxic side effects. The molecular basis of the kinase inhibition is not currently known, and is here investigated by computer simulations. It is shown that p5 binds the kinase at the same CDK5/p25 and CDK5/p35 interfaces, and is thus a non-selective competitor of both activators, in agreement with available experimental data in vitro. Binding of p5 is enthalpically driven with an affinity estimated in the low µM range. A quantitative description of the binding site and pharmacophore is presented, and options are discussed to increase the binding affinity and selectivity in the design of drug-like compounds against AD.

Altered secondary structure of Dynorphin A associates with loss of opioid signalling and NMDA-mediated excitotoxicity in SCA23

Spinocerebellar ataxia type 23 (SCA23) is caused by missense mutations in prodynorphin, encoding the precursor protein for the opioid neuropeptides α-neoendorphin, Dynorphin (Dyn) A and Dyn B, leading to neurotoxic elevated mutant Dyn A levels. Dyn A acts on opioid receptors to reduce pain in the spinal cord, but its cerebellar function remains largely unknown. Increased concentration of or prolonged exposure to Dyn A is neurotoxic and these deleterious effects are very likely caused by an N-methyl-d-aspartate-mediated non-opioid mechanism as Dyn A peptides were shown to bind NMDA receptors and potentiate their glutamate-evoked currents. In the present study, we investigated the cellular mechanisms underlying SCA23-mutant Dyn A neurotoxicity. We show that SCA23 mutations in the Dyn A-coding region disrupted peptide secondary structure leading to a loss of the N-terminal α-helix associated with decreased κ-opioid receptor affinity. Additionally, the altered secondary structure led to increased peptide stability of R6W and R9C Dyn A, as these peptides showed marked degradation resistance, which coincided with decreased peptide solubility. Notably, L5S Dyn A displayed increased degradation and no aggregation. R6W and wt Dyn A peptides were most toxic to primary cerebellar neurons. For R6W Dyn A, this is likely because of a switch from opioid to NMDA- receptor signalling, while for wt Dyn A, this switch was not observed. We propose that the pathology of SCA23 results from converging mechanisms of loss of opioid-mediated neuroprotection and NMDA-mediated excitotoxicity.

Evaluating the dynamics and electrostatic interactions of folded proteins in implicit solvents

Three implicit solvent models, namely GBMVII, FACTS, and SCPISM, were evaluated for their abilities to emulate an explicit solvent environment by comparing the simulated conformational ensembles, dynamics, and electrostatic interactions of the Src SH2 domain and the Lyn kinase domain. This assessment in terms of structural features in folded proteins expands upon the use of hydration energy as a metric for comparison. All-against-all rms coordinate deviation, average positional fluctuations, and ion-pair distance distribution were used to compare the implicit solvent models with the TIP3P explicit solvent model. Our study shows that the Src SH2 domains solvated with TIP3P, GBMVII, and FACTS sample similar global conformations. Additionally, the Src SH2 ion-pair distance distributions of solvent-exposed side chains corresponding to TIP3P, GBMVII, and FACTS do not differ substantially, indicating that GBMVII and FACTS are capable of modeling these electrostatic interactions. The ion-pair distance distributions of SCPISM are distinct from others, demonstrating that these electrostatic interactions are not adequately reproduced with the SCPISM model. On the other hand, for the Lyn kinase domain, a non-globular protein with bilobal structure and a large concavity on the surface, implicit solvent does not accurately model solvation to faithfully reproduce partially buried electrostatic interactions and lobe-lobe conformations. Our work reveals that local structure and dynamics of small, globular proteins are modeled well using FACTS and GBMVII. Nonetheless, global conformations and electrostatic interactions in concavities of multi-lobal proteins resulting from simulations with implicit solvent models do not match those obtained from explicit water simulations.

Detection and characterization of nonspecific, sparsely populated binding modes in the early stages of complexation

A method is proposed to study protein-ligand binding in a system governed by specific and nonspecific interactions. Strong associations lead to narrow distributions in the proteins configuration space; weak and ultraweak associations lead instead to broader distributions, a manifestation of nonspecific, sparsely populated binding modes with multiple interfaces. The method is based on the notion that a discrete set of preferential first-encounter modes are metastable states from which stable (prerelaxation) complexes at equilibrium evolve. The method can be used to explore alternative pathways of complexation with statistical significance and can be integrated into a general algorithm to study protein interaction networks. The method is applied to a peptide-protein complex. The peptide adopts several low-population conformers and binds in a variety of modes with a broad range of affinities. The system is thus well suited to analyze general features of binding, including conformational selection, multiplicity of binding modes, and nonspecific interactions, and to illustrate how the method can be applied to study these problems systematically. The equilibrium distributions can be used to generate biasing functions for simulations of multiprotein systems from which bulk thermodynamic quantities can be calculated.

Implicit treatment of solvent dispersion forces in protein simulations

A model is proposed for the evaluation of dispersive forces in a continuum solvent representation for use in large-scale computer simulations. The model captures the short- and long-range effects of water-exclusion in conditions of partial and anisotropic hydration. The model introduces three parameters, one of which represents the degree of hydration (water occupancy) at any point in the system, which depends on the solute conformation, and two that represent the strength of water-water and water-solute dispersive interactions. The model is optimized for proteins, using hydration data of a suboptimally hydrated binding site and results from dynamics simulations in explicit water. The model is applied to a series of aliphatic-alcohol/protein complexes and a set of binary and ternary complexes of various sizes. Implications for weak and ultra-weak protein-protein association and for simulation in crowded media are discussed.

Specific and non-specific protein association in solution: computation of solvent effects and prediction of first-encounter modes for efficient configurational bias Monte Carlo simulations

Weak and ultraweak protein-protein association play a role in molecular recognition and can drive spontaneous self-assembly and aggregation. Such interactions are difficult to detect experimentally, and are a challenge to the force field and sampling technique. A method is proposed to identify low-population protein-protein binding modes in aqueous solution. The method is designed to identify preferential first-encounter complexes from which the final complex(es) at equilibrium evolve. A continuum model is used to represent the effects of the solvent, which accounts for short- and long-range effects of water exclusion and for liquid-structure forces at protein/liquid interfaces. These effects control the behavior of proteins in close proximity and are optimized on the basis of binding enthalpy data and simulations. An algorithm is described to construct a biasing function for self-adaptive configurational-bias Monte Carlo of a set of interacting proteins. The function allows mixing large and local changes in the spatial distribution of proteins, thereby enhancing sampling of relevant microstates. The method is applied to three binary systems. Generalization to multiprotein complexes is discussed.

Weak intra-ring allosteric communications of the archaeal chaperonin thermosome revealed by normal mode analysis

Chaperonins are molecular machines that use ATP-driven cycles to assist misfolded substrate proteins to reach the native state. During the functional cycle, these machines adopt distinct nucleotide-dependent conformational states, which reflect large-scale allosteric changes in individual subunits. Distinct allosteric kinetics has been described for the two chaperonin classes. Bacterial (group I) chaperonins, such as GroEL, undergo concerted subunit motions within each ring, whereas archaeal and eukaryotic chaperonins (group II) undergo sequential subunit motions. We study these distinct mechanisms through a comparative normal mode analysis of monomer and double-ring structures of the archaeal chaperonin thermosome and GroEL. We find that thermosome monomers of each type exhibit common low-frequency behavior of normal modes. The observed distinct higher-frequency modes are attributed to functional specialization of these subunit types. The thermosome double-ring structure has larger contribution from higher-frequency modes, as it is found in the GroEL case. We find that long-range intersubunit correlation of amino-acid pairs is weaker in the thermosome ring than in GroEL. Overall, our results indicate that distinct allosteric behavior of the two chaperonin classes originates from different wiring of individual subunits as well as of the intersubunit communications.

Generating reservoir conformations for replica exchange through the use of the conformational space annealing method

Temperature replica exchange molecular dynamics (T-REM) has been successfully used to improve the conformational search for model peptides and small proteins. However, for larger and more complicated systems, the use of T-REM is computationally intensive since the complexity of the free energy landscape and number of replicas required increase with system size. Achieving convergence of systems with slow transition kinetics is often difficult. Several methods have been proposed to overcome the size and convergence speed issues of standard T-REM. One of these is the Reservoir Replica Exchange Method (R-REM), in which the conformational search and temperature equilibration are separated by exchanging with a pre-existing reservoir of structures. This approach allows the integration of computationally efficient search algorithms with replica exchange. The Conformational Space Annealing (CSA) method has been shown to be able to determine the global energy minimum of proteins efficiently and has been used in structure prediction successfully. CSA uses a genetic algorithm to generate a diverse set of conformations to determine the minimum energy structure. We combine these methods by using conformations generated by the CSA method to build a reservoir. R-REM is then used to seed the top replica with the structures from the reservoir; fast convergence at every temperature is observed. The efficiency of this method is then demonstrated with model peptides and small proteins, and significant improvement of efficiency is observed while maintaining the overall shape of the free energy landscape.

Replica exchanging self-guided Langevin dynamics for efficient and accurate conformational sampling

This work presents a replica exchanging self-guided Langevin dynamics (RXSGLD) simulation method for efficient conformational searching and sampling. Unlike temperature-based replica exchanging simulations, which use high temperatures to accelerate conformational motion, this method uses self-guided Langevin dynamics (SGLD) to enhance conformational searching without the need to elevate temperatures. A RXSGLD simulation includes a series of SGLD simulations, with simulation conditions differing in the guiding effect and/or temperature. These simulation conditions are called stages and the base stage is one with no guiding effect. Replicas of a simulation system are simulated at the stages and are exchanged according to the replica exchanging probability derived from the SGLD partition function. Because SGLD causes less perturbation on conformational distribution than high temperatures, exchanges between SGLD stages have much higher probabilities than those between different temperatures. Therefore, RXSGLD simulations have higher conformational searching ability than temperature based replica exchange simulations. Through three example systems, we demonstrate that RXSGLD can generate target canonical ensemble distribution at the base stage and achieve accelerated conformational searching. Especially for large systems, RXSGLD has remarkable advantages in terms of replica exchange efficiency, conformational searching ability, and system size extensiveness.

GPU/CPU Algorithm for Generalized Born/Solvent-Accessible Surface Area Implicit Solvent Calculations

Molecular dynamics methodologies comprise a vital research tool for structural biology. Molecular dynamics has benefited from technological advances in computing, such as multi-core CPUs and graphics processing units (GPUs), but harnessing the full power of hybrid GPU/CPU computers remains difficult. The generalized Born/solvent-accessible surface area implicit solvent model (GB/SA) stands to benefit from hybrid GPU/CPU computers, employing the GPU for the GB calculation and the CPU for the SA calculation. Here, we explore the computational challenges facing GB/SA calculations on hybrid GPU/CPU computers and demonstrate how NAMD, a parallel molecular dynamics program, is able to efficiently utilize GPUs and CPUs simultaneously for fast GB/SA simulations. The hybrid computation principles demonstrated here are generally applicable to parallel applications employing hybrid GPU/CPU calculations.

Water-exclusion and liquid-structure forces in implicit solvation

A continuum model of solvation is proposed to describe (i) long-range electrostatic effects of water exclusion resulting from incomplete and anisotropic hydration in crowded environments and (ii) short-range effects of liquid-structure forces on the hydrogen-bond interactions at solute/water interfaces. The model is an extension of the phenomenological screened coulomb potential-based implicit model of solvation. The developments reported here allow a more realistic representation of highly crowded and spatially heterogeneous environments, such as those in the interior of a living cell. Only the solvent is treated as a continuum medium. It is shown that the electrostatic effects of long-range water-exclusion can strongly affect protein-protein binding energies and are then related to the thermodynamics of complex formation. Hydrogen-bond interactions modulated by the liquid structure at interfaces are calibrated based on systematic calculations of potentials of mean force in explicit water. The electrostatic component of the model is parametrized for monovalent, divalent and trivalent ions. The conceptual and practical aspects of the model are discussed based on simulations of protein complexation and peptide folding. The current implementation is ~1.5 times slower than the gas-phase force field and exhibits good parallel performance.

Parallel Generalized Born Implicit Solvent Calculations with NAMD

Accurate electrostatic descriptions of aqueous solvent are critical for simulation studies of bio-molecules, but the computational cost of explicit treatment of solvent is very high. A computationally more feasible alternative is a generalized Born implicit solvent description which models polar solvent as a dielectric continuum. Unfortunately, the attainable simulation speedup does not transfer to the massive parallel computers often employed for simulation of large structures. Longer cutoff distances, spatially heterogenous distribution of atoms and the necessary three-fold iteration over atom-pairs in each timestep combine to challenge efficient parallel performance of generalized Born implicit solvent algorithms. Here we report how NAMD, a parallel molecular dynamics program, meets the challenge through a unique parallelization strategy. NAMD now permits efficient simulation of large systems whose slow conformational motions benefit most from implicit solvent descriptions due to the inherent low viscosity. NAMD's implicit solvent performance is benchmarked and then illustrated in simulating the ratcheting Escherichia coli ribosome involving ~250,000 atoms.

Probing the structural determinants for the function of intracellular loop 2 in structurally cognate G-protein-coupled receptors

Intracellular loop 2 (IL2) in G-protein-coupled receptors (GPCRs) is functionally important, e.g., in binding to G-protein and β-arrestin. Differences in secondary structure of IL2 in the crystal structures of the very similar β(1)- and β(2)-adrenergic receptors (β(1)AR and β(2)AR, respectively), i.e., an α-helix and an L-shaped strand, respectively, emphasize the need to understand the structural basis for IL2 functionality. We studied the properties of IL2 in the context of experimental data using a Monte Carlo-based ab initio method. The procedure was validated first by verifying that the IL2 structures in β(1)AR and β(2)AR crystals were correctly reproduced, even after conformational ensemble searches at >1200 K where most secondary structure had been lost. We found that IL2 in β(1)AR and β(2)AR sampled each other's conformation but adopted different energetically preferred conformations, consistent with the crystal structures. The results indicate a persistent contextual preference for the structure of IL2, which was conserved when the IL2 sequences were interchanged between the receptors. We conclude that the protein environment, more than the IL2 sequence, regulates the IL2 structures. We extended the approach to the molecular model of 5-HT(2A)R for which no crystal structure is available and found that IL2 is predominantly helical, similar to IL2 in β(1)AR. Because the P3.57A mutation in IL2 had been shown to decrease β-arrestin binding and internalization, we predicted the effects of the mutation and found that it decreased the propensity of IL2 to form helix, identifying the helical IL2 as a component of the GPCR active form.

CHARMM: the biomolecular simulation program

CHARMM (Chemistry at HARvard Molecular Mechanics) is a highly versatile and widely used molecular simulation program. It has been developed over the last three decades with a primary focus on molecules of biological interest, including proteins, peptides, lipids, nucleic acids, carbohydrates, and small molecule ligands, as they occur in solution, crystals, and membrane environments. For the study of such systems, the program provides a large suite of computational tools that include numerous conformational and path sampling methods, free energy estimators, molecular minimization, dynamics, and analysis techniques, and model-building capabilities. The CHARMM program is applicable to problems involving a much broader class of many-particle systems. Calculations with CHARMM can be performed using a number of different energy functions and models, from mixed quantum mechanical-molecular mechanical force fields, to all-atom classical potential energy functions with explicit solvent and various boundary conditions, to implicit solvent and membrane models. The program has been ported to numerous platforms in both serial and parallel architectures. This article provides an overview of the program as it exists today with an emphasis on developments since the publication of the original CHARMM article in 1983.

Intramolecular disulfide bonds of the prolactin receptor short form are required for its inhibitory action on the function of the long form of the receptor

The short form (S1b) of the prolactin receptor (PRLR) silences prolactin-induced activation of gene transcription by the PRLR long form (LF). The functional and structural contributions of two intramolecular disulfide (S-S) bonds within the extracellular subdomain 1 (D1) of S1b to its inhibitory function on the LF were investigated. Mutagenesis of the paired cysteines eliminated the inhibitory action of S1b. The expression of the mutated S1b (S1bx) on the cell surface was not affected, indicating native-like folding of the receptor. The constitutive JAK2 phosphorylation observed in S1b was not present in cells expressing S1bx, and JAK2 association was disrupted. BRET(50) (BRET(50) represents the relative affinity as acceptor/donor ratio required to reach half-maximal BRET [bioluminescence resonance energy transfer] values) showed decreased LF/S1bx heterodimeric-association and increased affinity in S1bx homodimerization, thus favoring LF homodimerization and prolactin-induced signaling. Computer modeling based on the PRLR crystal structure showed that minor changes in the tertiary structure of D1 upon S-S bond disruption propagated to the quaternary structure of the homodimer, affecting the dimerization interface. These changes explain the higher homodimerization affinity of S1bx and provide a structural basis for its lack of inhibitory function. The PRLR conformation as stabilized by S-S bonds is required for the inhibitory action of S1b on prolactin-induced LF-mediated function and JAK2 association.

Normal mode partitioning of Langevin dynamics for biomolecules

We propose a novel normal mode multiple time stepping Langevin dynamics integrator called NML. The aim is to approximate the kinetics or thermodynamics of a biomolecule by a reduced model based on a normal mode decomposition of the dynamical space. Our basis set uses the eigenvectors of a mass reweighted Hessian matrix calculated with a biomolecular force field. This particular choice has the advantage of an ordering according to the eigenvalues, which have a physical meaning of being the square of the mode frequency. Low frequency eigenvalues correspond to more collective motions, whereas the highest frequency eigenvalues are the limiting factor for the stability of the integrator. In NML, the higher frequency modes are overdamped and relaxed near their energy minimum while respecting the subspace of low frequency dynamical modes. Our numerical results confirm that both sampling and rates are conserved for an implicitly solvated alanine dipeptide model, with only 30% of the modes propagated, when compared to the full model. For implicitly solvated systems, NML gives a twofold improvement in efficiency over plain Langevin dynamics for sampling a small 22 atom (alanine dipeptide) model and in excess of an order of magnitude for sampling an 882 atom (bovine pancreatic trypsin inhibitor) model, with good scaling with system size subject to the number of modes propagated. NML has been implemented in the open source software PROTOMOL.

FACTS: Fast analytical continuum treatment of solvation

An efficient method for calculating the free energy of solvation of a (macro)molecule embedded in a continuum solvent is presented. It is based on the fully analytical evaluation of the volume and spatial symmetry of the solvent that is displaced from around a solute atom by its neighboring atoms. The two measures of solvent displacement are combined in empirical equations to approximate the atomic (or self) electrostatic solvation energy and the solvent accessible surface area. The former directly yields the effective Born radius, which is used in the generalized Born (GB) formula to calculate the solvent-screened electrostatic interaction energy. A comparison with finite-difference Poisson data shows that atomic solvation energies, pair interaction energies, and their sums are evaluated with a precision comparable to the most accurate GB implementations. Furthermore, solvation energies of a large set of protein conformations have an error of only 1.5%. The solvent accessible surface area is used to approximate the nonpolar contribution to solvation. The empirical approach, called FACTS (Fast Analytical Continuum Treatment of Solvation), is only four times slower than using the vacuum energy in molecular dynamics simulations of proteins. Notably, the folded state of structured peptides and proteins is stable at room temperature in 100-ns molecular dynamics simulations using FACTS and the CHARMM force field.

Evidence that proline focuses movement of the floppy loop of arylalkylamine N-acetyltransferase (EC 2.3.1.87)

Arylalkylamine N-acetyltransferase (AANAT) catalyzes the N-acetylation of serotonin, the penultimate step in the synthesis of melatonin. Pineal AANAT activity increases at night in all vertebrates, resulting in increased melatonin production. This increases circulating levels of melatonin, thereby providing a hormonal signal of darkness. Kinetic and structural analysis of AANAT has determined that one element is floppy. This element, termed Loop 1, is one of three loops that comprise the arylalkylamine binding pocket. During the course of chordate evolution, Loop 1 acquired the tripeptide CPL, and the enzyme became highly active. Here we focused on the functional importance of the CPL tripeptide and found that activity was markedly reduced when it was absent. Moreover, increasing the local flexibility of this tripeptide region by P64G and P64A mutations had the counterintuitive effect of reducing activity and reducing the overall movement of Loop 1, as estimated from Langevin dynamics simulations. Binding studies indicate that these mutations increased the off-rate constant of a model substrate without altering the dissociation constant. The structural kink and local rigidity imposed by Pro-64 may enhance activity by favoring configurations of Loop 1 that facilitate catalysis and do not become immobilized by intramolecular interactions.

Balancing an accurate representation of the molecular surface in generalized born formalisms with integrator stability in molecular dynamics simulations

Different integrator time steps in NVT and NVE simulations of protein and nucleic acid systems are tested with the GBMV (Generalized Born using Molecular Volume) and GBSW (Generalized Born with simple SWitching) methods. The simulation stability and energy conservation is investigated in relation to the agreement with the Poisson theory. It is found that very close agreement between generalized Born methods and the Poisson theory based on the commonly used sharp molecular surface definition results in energy drift and simulation artifacts in molecular dynamics simulation protocols with standard 2-fs time steps. New parameters are proposed for the GBMV method, which maintains very good agreement with the Poisson theory while providing energy conservation and stable simulations at time steps of 1 to 1.5 fs.

Liquid-structure forces and electrostatic modulation of biomolecular interactions in solution

Molecular interactions in solution are controlled by the bulk medium and by the forces originating in the structured region of the solvent close to the solutes. In this paper, a model of electrostatic and liquid-structure forces for dynamics simulations of biomolecules is presented. The model introduces information on the microscopic nature of the liquid in the vicinity of polar and charged groups and the associated non-pairwise character of the forces, thus improving upon conventional continuum representations. The solvent is treated as a polar and polarizable medium, with dielectric properties described by an inhomogeneous version of the Onsager theory. This treatment leads to an effective position-dependent dielectric permittivity that incorporates saturation effects of the electric field and the spatial variation of the liquid density. The non-pairwise additivity of the liquid-structure forces is represented by centers of force located at specific points in the liquid phase. These out-of-the-solute centers are positioned at the peaks of liquid density and exert local, external forces on the atoms of the solute. The density is calculated from a barometric law, using a Lennard-Jones-type solute-liquid effective interaction potential. The conceptual aspects of the model and its exact numerical solutions are discussed for single alkali and halide ions and for ion-pair interactions. The practical aspects of the model and the simplifications introduced for efficient computation of forces in molecular solutes are discussed in the context of polar and charged amino acid dimers. The model reproduces the contact and solvent-separated minima and the desolvation barriers of intermolecular potentials of mean force of amino acid dimers, as observed in atomistic dynamics simulations. Possible refinements based on an improved treatment of molecular correlations are discussed.

Development and application of hybrid structure based method for efficient screening of ligands binding to G-protein coupled receptors

G-protein coupled receptors (GPCRs) comprise a large superfamily of proteins that are targets for nearly 50% of drugs in clinical use today. In the past, the use of structure-based drug design strategies to develop better drug candidates has been severely hampered due to the absence of the receptor's three-dimensional structure. However, with recent advances in molecular modeling techniques and better computing power, atomic level details of these receptors can be derived from computationally derived molecular models. Using information from these models coupled with experimental evidence, it has become feasible to build receptor pharmacophores. In this study, we demonstrate the use of the Hybrid Structure Based (HSB) method that can be used effectively to screen and identify prospective ligands that bind to GPCRs. Essentially; this multi-step method combines ligand-based methods for building enriched libraries of small molecules and structure-based methods for screening molecules against the GPCR target. The HSB method was validated to identify retinal and its analogues from a random dataset of approximately 300,000 molecules. The results from this study showed that the 9 top-ranking molecules are indeed analogues of retinal. The method was also tested to identify analogues of dopamine binding to the dopamine D2 receptor. Six of the ten top-ranking molecules are known analogues of dopamine including a prodrug, while the other thirty-four molecules are currently being tested for their activity against all dopamine receptors. The results from both these test cases have proved that the HSB method provides a realistic solution to bridge the gap between the ever-increasing demand for new drugs to treat psychiatric disorders and the lack of efficient screening methods for GPCRs.

Ab initio computational modeling of loops in G-protein-coupled receptors: lessons from the crystal structure of rhodopsin

With the help of the crystal structure of rhodopsin an ab initio method has been developed to calculate the three-dimensional structure of the loops that connect the transmembrane helices (TMHs). The goal of this procedure is to calculate the loop structures in other G-protein coupled receptors (GPCRs) for which only model coordinates of the TMHs are available. To mimic this situation a construct of rhodopsin was used that only includes the experimental coordinates of the TMHs while the rest of the structure, including the terminal domains, has been removed. To calculate the structure of the loops a method was designed based on Monte Carlo (MC) simulations which use a temperature annealing protocol, and a scaled collective variables (SCV) technique with proper structural constraints. Because only part of the protein is used in the calculations the usual approach of modeling loops, which consists of finding a single, lowest energy conformation of the system, is abandoned because such a single structure may not be a representative member of the native ensemble. Instead, the method was designed to generate structural ensembles from which the single lowest free energy ensemble is identified as representative of the native folding of the loop. To find the native ensemble a successive series of SCV-MC simulations are carried out to allow the loops to undergo structural changes in a controlled manner. To increase the chances of finding the native funnel for the loop, some of the SCV-MC simulations are carried out at elevated temperatures. The native ensemble can be identified by an MC search starting from any conformation already in the native funnel. The hypothesis is that native structures are trapped in the conformational space because of the high-energy barriers that surround the native funnel. The existence of such ensembles is demonstrated by generating multiple copies of the loops from their crystal structures in rhodopsin and carrying out an extended SCV-MC search. For the extracellular loops e1 and e3, and the intracellular loop i1 that were used in this work, the procedure resulted in dense clusters of structures with Calpha-RMSD approximately 0.5 angstroms. To test the predictive power of the method the crystal structure of each loop was replaced by its extended conformations. For e1 and i1 the procedure identifies native clusters with Calpha-RMSD approximately 0.5 angstroms and good structural overlap of the side chains; for e3, two clusters were found with Calpha-RMSD approximately 1.1 angstroms each, but with poor overlap of the side chains. Further searching led to a single cluster with lower Calpha-RMSD but higher energy than the two previous clusters. This discrepancy was found to be due to the missing elements in the constructs available from experiment for use in the calculations. Because this problem will likely appear whenever parts of the structural information are missing, possible solutions are discussed.

Ab initio computational modeling of long loops in G-protein coupled receptors

A newly developed approach for predicting the structure of segments that connect known elements of secondary structure in proteins has been applied to some of the longer loops in the G-protein coupled receptors (GPCRs) rhodopsin and the dopamine receptor D2R. The algorithm uses Monte Carlo (MC) simulation in a temperature annealing protocol combined with a scaled collective variables (SCV) technique to search conformation space for loop structures that could belong to the native ensemble. Except for rhodopsin, structural information is only available for the transmembrane helices (TMHs), and therefore the usual approach of finding a single conformation of lowest energy has to be abandoned. Instead the MC search aims to find the ensemble located at the absolute minimum free energy, i.e., the native ensemble. It is assumed that structures in the native ensemble can be found by an MC search starting from any conformation in the native funnel. The hypothesis is that native structures are trapped in this part of conformational space because of the high-energy barriers that surround the native funnel. In this work it is shown that the crystal structure of the second extracellular loop (e2) of rhodopsin is a member of this loop's native ensemble. In contrast, the crystal structure of the third intracellular loop is quite different in the different crystal structures that have been reported. Our calculations indicate, that of three crystal structures examined, two show features characteristic of native ensembles while the other one does not. Finally the protocol is used to calculate the structure of the e2 loop in D2R. Here, the crystal structure is not known, but it is shown that several side chains that are involved in interaction with a class of substituted benzamides assume conformations that point into the active site. Thus, they are poised to interact with the incoming ligand.

Dynamic models of G-protein coupled receptor dimers: indications of asymmetry in the rhodopsin dimer from molecular dynamics simulations in a POPC bilayer

Based on the growing evidence that G-protein coupled receptors (GPCRs) form homo- and hetero-oligomers, models of GPCR signaling are now considering macromolecular assemblies rather than monomers, with the homo-dimer regarded as the minimal oligomeric arrangement required for functional coupling to the G-protein. The dynamic mechanisms of such signaling assemblies are unknown. To gain some insight into properties of GPCR dimers that may be relevant to functional mechanisms, we study their current structural prototype, rhodopsin. We have carried out nanosecond time-scale molecular dynamics (MD) simulations of a rhodopsin dimer and compared the results to the monomer simulated in the same type of bilayer membrane model composed of an equilibrated unit cell of hydrated palmitoyl-oleoyl-phosphatidyl choline (POPC). The dynamic representation of the homo-dimer reveals the location of structural changes in several regions of the monomeric subunits. These changes appear to be more pronounced at the dimerization interface that had been shown to be involved in the activation process [Proc Natl Acad Sci USA 102:17495, 2005]. The results are consistent with a model of GPCR activation that involves allosteric modulation through a single GPCR subunit per dimer.

Factors influencing the energetics of electron and proton transfers in proteins. What can be learned from calculations

A protein structure should provide the information needed to understand its observed properties. Significant progress has been made in developing accurate calculations of acid/base and oxidation/reduction reactions in proteins. Current methods and their strengths and weaknesses are discussed. The distribution and calculated ionization states in a survey of proteins is described, showing that a significant minority of acidic and basic residues are buried in the protein and that most of these remain ionized. The electrochemistry of heme and quinones are considered. Proton transfers in bacteriorhodopsin and coupled electron and proton transfers in photosynthetic reaction centers, 5-coordinate heme binding proteins and cytochrome c oxidase are highlighted as systems where calculations have provided insight into the reaction mechanism.

Exploring peptide energy landscapes: a test of force fields and implicit solvent models

A biased Monte Carlo-minimization/annealing conformational search was used to characterize five descriptions of the energy landscape for each of three model systems: the 20-residue "trp-cage" miniprotein, the 20-residue "BS1" peptide, and the 17-residue "U(1-17)T9D" peptide. The EEF1 and SASA energy landscapes were studied as well as those defined by using the GB/ACE implicit water model with one of three protein force fields: CHARMM19, CHARMM22, and CHARMM22/CMAP. The lowest-energy structures of the trp-cage and BS1 peptides found for the EEF1 landscape have main-chain root-mean-square deviations (rmsds) from the respective NMR structures of less than 2 A; for U(1-17)T9D, the deviation is less than 3 A using EEF1. The main-chain rmsd of the minimum-energy trp-cage conformation obtained for the GB/ACE/CHARMM22/CMAP landscape is less than 1 A. However, this energy function strongly favored helical structures for the two peptides shown by NMR to form beta-sheet structures. Brief annealing of the system following main-chain conformational changes was found to enhance the exploration of low-energy states. The thousands of simulations reported here suggest that the prediction of protein structure might be improved by the simultaneous use of a CMAP-like description of the main chain and an EEF1-like description of the solvent.

Hallucinogen actions on 5-HT receptors reveal distinct mechanisms of activation and signaling by G protein-coupled receptors

We review the effect of some key advances in the characterization of molecular mechanisms of signaling by G protein-coupled receptors (GPCRs) on our current understanding of mechanisms of drugs of abuse. These advances are illustrated by results from our ongoing work on the actions of hallucinogens on serotonin (5-HT) receptors. We show how a combined computational and experimental approach can reveal specific modes of receptor activation underlying the difference in properties of hallucinogens compared with nonhallucinogenic congeners. These modes of activation-that can produce distinct ligand-dependent receptor states-are identified in terms of structural motifs (SM) in molecular models of the receptors, which were shown to constitute conserved functional microdomains (FM). The role of several SM/FMs in the activation mechanism of the GPCRs is presented in detail to illustrate how this mechanism can lead to ligand-dependent modes of signaling by the receptors. Novel bioinformatics tools are described that were designed to support the quantitative mathematical modeling of ligand-specific signaling pathways activated by the 5-HT receptors targeted by hallucinogens. The approaches for mathematical modeling of signaling pathways activated by 5-HT receptors are described briefly in the context of ongoing work on detailed biochemical models of 5-HT2A, and combined 5-HT2A/5-HT1A, receptor-mediated activation of the MAPK 1,2 pathway. The continuing need for increasingly more realistic representation of signaling in dynamic compartments within the cell, endowed with spatio-temporal characteristics obtained from experiment, is emphasized. Such developments are essential for attaining a quantitative understanding of how the multiple functions of a cell are coordinated and regulated, and to evaluate the specifics of the perturbations caused by the drugs of abuse that target GPCRs.

Long dynamics simulations of proteins using atomistic force fields and a continuum representation of solvent effects: calculation of structural and dynamic properties

Long dynamics simulations were carried out on the B1 immunoglobulin-binding domain of streptococcal protein G (ProtG) and bovine pancreatic trypsin inhibitor (BPTI) using atomistic descriptions of the proteins and a continuum representation of solvent effects. To mimic frictional and random collision effects, Langevin dynamics (LD) were used. The main goal of the calculations was to explore the stability of tens-of-nanosecond trajectories as generated by this molecular mechanics approximation and to analyze in detail structural and dynamical properties. Conformational fluctuations, order parameters, cross correlation matrices, residue solvent accessibilities, pKa values of titratable groups, and hydrogen-bonding (HB) patterns were calculated from all of the trajectories and compared with available experimental data. The simulations comprised over 40 ns per trajectory for ProtG and over 30 ns per trajectory for BPTI. For comparison, explicit water molecular dynamics simulations (EW/MD) of 3 ns and 4 ns, respectively, were also carried out. Two continuum simulations were performed on each protein using the CHARMM program, one with the all-atom PAR22 representation of the protein force field (here referred to as PAR22/LD simulations) and the other with the modifications introduced by the recently developed CMAP potential (CMAP/LD simulations). The explicit solvent simulations were performed with PAR22 only. Solvent effects are described by a continuum model based on screened Coulomb potentials (SCP) reported earlier, i.e., the SCP-based implicit solvent model (SCP-ISM). For ProtG, both the PAR22/LD and the CMAP/LD 40-ns trajectories were stable, yielding C(alpha) root mean square deviations (RMSD) of about 1.0 and 0.8 A respectively along the entire simulation time, compared to 0.8 A for the EW/MD simulation. For BPTI, only the CMAP/LD trajectory was stable for the entire 30-ns simulation, with a C(alpha) RMSD of approximately 1.4 A, while the PAR22/LD trajectory became unstable early in the simulation, reaching a C(alpha) RMSD of about 2.7 A and remaining at this value until the end of the simulation; the C(alpha) RMSD of the EW/MD simulation was about 1.5 A. The source of the instabilities of the BPTI trajectories in the PAR22/LD simulations was explored by an analysis of the backbone torsion angles. To further validate the findings from this analysis of BPTI, a 35-ns SCP-ISM simulation of Ubiquitin (Ubq) was carried out. For this protein, the CMAP/LD simulation was stable for the entire simulation time (C(alpha) RMSD of approximately 1.0 A), while the PAR22/LD trajectory showed a trend similar to that in BPTI, reaching a C(alpha) RMSD of approximately 1.5 A at 7 ns. All the calculated properties were found to be in agreement with the corresponding experimental values, although local deviations were also observed. HB patterns were also well reproduced by all the continuum solvent simulations with the exception of solvent-exposed side chain-side chain (sc-sc) HB in ProtG, where several of the HB interactions observed in the crystal structure and in the EW/MD simulation were lost. The overall analysis reported in this work suggests that the combination of an atomistic representation of a protein with a CMAP/CHARMM force field and a continuum representation of solvent effects such as the SCP-ISM provides a good description of structural and dynamic properties obtained from long computer simulations. Although the SCP-ISM simulations (CMAP/LD) reported here were shown to be stable and the properties well reproduced, further refinement is needed to attain a level of accuracy suitable for more challenging biological applications, particularly the study of protein-protein interactions.

Recent advances in the development and application of implicit solvent models in biomolecule simulations

Advances have recently been made in the development of implicit solvent methodologies and their application to the modeling of biomolecules, particularly with regard to generalized Born approaches, dielectric screening function formulations and models based on solvent-accessible surface areas. Interesting new developments include more refined non-polar solvation energy estimators, and implicit methods for modeling low-dielectric and heterogeneous environments such as membrane systems. These have been successfully applied to molecular dynamics simulations, the scoring of protein conformations, and the calculation of binding affinities and folding free energy landscapes.

Molecular determinants for the complex binding specificity of the PDZ domain in PICK1

PICK1 (protein interacting with C kinase 1) contains a single PDZ domain known to mediate interaction with the C termini of several receptors, transporters, ion channels, and kinases. In contrast to most PDZ domains, the PICK1 PDZ domain interacts with binding sequences classifiable as type I (terminating in (S/T)XPhi; X, any residue) as well as type II (PhiXPhi; Phi, any hydrophobic residue). To enable direct assessment of the affinity of the PICK1 PDZ domain for its binding partners we developed a purification scheme for PICK1 and a novel quantitative binding assay based on fluorescence polarization. Our results showed that the PICK1 PDZ domain binds the type II sequence presented by the human dopamine transporter (-WLKV) with an almost 15-fold and >100-fold higher affinity than the type I sequences presented by protein kinase Calpha (-QSAV) and the beta(2)-adrenergic receptor (-DSLL), respectively. Mutational analysis of Lys(83) in the alphaB1 position of the PDZ domain suggested that this residue mimics the function of hydrophobic residues present in this position in regular type II PDZ domains. The PICK1 PDZ domain was moreover found to prefer small hydrophobic residues in the C-terminal P(0) position of the ligand. Molecular modeling predicted a rank order of (Val > Ile > Leu) that was verified experimentally with up to a approximately 16-fold difference in binding affinity between a valine and a leucine in P(0). The results define the structural basis for the unusual binding pattern of the PICK1 PDZ domain by substantiating the critical role of the alphaB1 position (Lys(83)) and of discrete side chain differences in position P(0) of the ligands.

Nonuniform charge scaling (NUCS): A practical approximation of solvent electrostatic screening in proteins

In molecular mechanics calculations, electrostatic interactions between chemical groups are usually represented by a Coulomb potential between the partial atomic charges of the groups. In aqueous solution these interactions are modified by the polarizable solvent. Although the electrostatic effects of the polarized solvent on the protein are well described by the Poisson--Boltzmann equation, its numerical solution is computationally expensive for large molecules such as proteins. The procedure of nonuniform charge scaling (NUCS) is a pragmatic approach to implicit solvation that approximates the solvent screening effect by individually scaling the partial charges on the explicit atoms of the macromolecule so as to reproduce electrostatic interaction energies obtained from an initial Poisson--Boltzmann analysis. Once the screening factors have been determined for a protein the scaled charges can be easily used in any molecular mechanics program that implements a Coulomb term. The approach is particularly suitable for minimization-based simulations, such as normal mode analysis, certain conformational reaction path or ligand binding techniques for which bulk solvent cannot be included explicitly, and for combined quantum mechanical/molecular mechanical calculations when the interface to more elaborate continuum solvent models is lacking. The method is illustrated using reaction path calculations of the Tyr 35 ring flip in the bovine pancreatic trypsin inhibitor.