Fast, approximate algorithm for detection of solvent-inaccessible atoms

Microsecond dynamics of H10N7 influenza neuraminidase reveals the plasticity of loop regions and drug resistance due to the R292K mutation

At the beginning of the last century, multiple pandemics caused by influenza (flu) viruses severely impacted public health. Despite the development of vaccinations and antiviral medications to prevent and control impending flu outbreaks, unforeseen novel strains and continuously evolving old strains continue to represent a serious threat to human life. Therefore, the recently identified H10N7, for which not much data is available for rational structure-based drug design, needs to be further explored. Here, we investigated the structural dynamics of neuraminidase N7 upon binding of inhibitors, and the drug resistance mechanisms against the oseltamivir (OTV) and laninamivir (LNV) antivirals due to the crucial R292K mutation on the N7 using the computational microscope, molecular dynamics (MD) simulations. In this study, each system underwent long 2 × 1 μs MD simulations to answer the conformational changes and drug resistance mechanisms. These long time-scale dynamics simulations and free energy landscapes demonstrated that the mutant systems showed a high degree of conformational variation compared to their wildtype (WT) counterparts, and the LNV-bound mutant exhibited an extended 150-loop conformation. Further, the molecular mechanics Poisson-Boltzmann surface area (MM/PBSA) calculation and MM/GBSA free energy decomposition were used to characterize the binding of OTV and LNV with WT, and R292K mutated N7, revealing the R292K mutation as drug-resistant, facilitated by a decline in binding interaction and a reduction in the dehydration penalty. Due to the broader binding pocket cavity of the smaller K292 mutant residue relative to the wildtype, the drug carboxylate to K292 hydrogen bonding was lost, and the area surrounding the K292 residue was more accessible to water molecules. This implies that drug resistance could be reduced by strengthening the hydrogen bond contacts between N7 inhibitors and altered N7, creating inhibitors that can form a hydrogen bond to the mutant K292, or preserving the closed cavity conformations.

Decoding the Host-Parasite Protein Interactions Involved in Cerebral Malaria Through Glares of Molecular Dynamics Simulations

Malaria causes millions of deaths every year. The malaria parasite spends a substantial part of its life cycle inside human erythrocytes. Inside erythrocytes, it synthesizes and displays various proteins onto the erythrocyte surface, such as Plasmodium falciparum erythrocytic membrane protein-1 (PfEMP1). This protein contains cysteine-rich interdomain region (CIDR) domains which have many subtypes based on sequence diversity and can cross-talk with host molecules. The CIDRα1.4 subtype can attach host endothelial protein C receptor (EPCR). This interaction facilitates infected erythrocyte adherence to brain endothelium and subsequent development of cerebral malaria. Through molecular dynamics simulations in conjunction with the molecular mechanics Poisson-Boltzmann surface area (MM/PBSA) method, we explored the mechanism of interaction in the CIDRα1-EPCR complex. We examined the structural behavior of two CIDRα1 molecules (encoded by HB3-isolate var03-gene and IT4-isolate var07-gene) with EPCR unbound and bound (complex) forms. ^HB3var03CIDRα1 in apo and complexed with EPCR was comparatively more stable than ^IT4var07CIDRα1. Both of the complexes adopted two distinct conformational energy states. The hydrophobic residues played a crucial role in the binding of both complexes. For ^HB3var03CIDRα1-EPCR, the dominant energetic components were total polar interactions, while in ^IT4var07CIDRα1-EPCR, the primary interaction was van der Waals and nonpolar solvation energy. The study also revealed details such as correlated conformational motions and secondary structure evolution. Further, it elucidated various hotspot residues involved in protein-protein recognition. Overall, our study provides additional information on the structural behavior of CIDR molecules in unbound and receptor-bound states, which will help to design potent inhibitors.

Molecular Mechanism of Inhibiting WNK Binding to OSR1 by Targeting the Allosteric Pocket of the OSR1-CCT Domain with Potential Antihypertensive Inhibitors: An In Silico Study

The oxidative-stress-responsive kinase 1 (OSR1) and the STE20/SPS1-related proline-alanine-rich kinase (SPAK) are physiological substrates of the with-no-lysine (WNK) kinase. They are the master regulators of cation Cl^- cotransporters that could be targeted for discovering novel antihypertensive agents. Both kinases have a conserved carboxy-terminal (CCT) domain that recognizes a unique peptide motif (Arg-Phe-Xaa-Val) present in their upstream kinases and downstream substrates. Here, we have combined molecular docking with molecular dynamics simulations and free-energy calculations to identify potential inhibitors that can bind to the allosteric pocket of the OSR1-CCT domain and impede its interaction with the WNK peptide. Our study revealed that STOCK1S-14279 and Closantel bound strongly to the allosteric pocket of OSR1 and displaced the WNK peptide from the primary pocket of OSR1. We showed that primarily Arg1004 and Gln1006 of the WNK4-peptide motif were involved in strong H-bond interactions with Glu453 and Arg451 of OSR1. Besides, our study revealed that atoms of Arg1004 were solvent-exposed in cases of STOCK1S-14279 and Closantel, implying that the WNK4 peptide was moved out of the pocket. Overall, the predicted potential inhibitors altogether abolish the OSR1-WNK4-peptide interaction, suggesting their potency as a prospective allosteric inhibitor against OSR1.

Design of Thioether Cyclic Peptide Scaffolds with Passive Permeability and Oral Exposure

Advances in the design of permeable peptides and in the synthesis of large arrays of macrocyclic peptides with diverse amino acids have evolved on parallel but independent tracks. Less precedent combines their respective attributes, thereby limiting the potential to identify permeable peptide ligands for key targets. Herein, we present novel 6-, 7-, and 8-mer cyclic peptides (MW 774-1076 g·mol-1) with passive permeability and oral exposure that feature the amino acids and thioether ring-closing common to large array formats, including DNA- and RNA-templated synthesis. Each oral peptide herein, selected from virtual libraries of partially N-methylated peptides using in silico methods, reflects the subset consistent with low energy conformations, low desolvation penalties, and passive permeability. We envision that, by retaining the backbone N-methylation pattern and consequent bias toward permeability, one can generate large peptide arrays with sufficient side chain diversity to identify permeability-biased ligands to a variety of protein targets.

Computational Investigation of Structural Dynamics of SARS-CoV-2 Methyltransferase-Stimulatory Factor Heterodimer nsp16/nsp10 Bound to the Cofactor SAM

Recently, a highly contagious novel coronavirus disease 2019 (COVID-19), caused by SARS-CoV-2, has emerged, posing a global threat to public health. Identifying a potential target and developing vaccines or antiviral drugs is an urgent demand in the absence of approved therapeutic agents. The 5'-capping mechanism of eukaryotic mRNA and some viruses such as coronaviruses (CoVs) are essential for maintaining the RNA stability and protein translation in the virus. SARS-CoV-2 encodes S-adenosyl-L-methionine (SAM) dependent methyltransferase (MTase) enzyme characterized by nsp16 (2'-O-MTase) for generating the capped structure. The present study highlights the binding mechanism of nsp16 and nsp10 to identify the role of nsp10 in MTase activity. Furthermore, we investigated the conformational dynamics and energetics behind the binding of SAM to nsp16 and nsp16/nsp10 heterodimer by employing molecular dynamics simulations in conjunction with the Molecular Mechanics Poisson-Boltzmann Surface Area (MM/PBSA) method. We observed from our simulations that the presence of nsp10 increases the favorable van der Waals and electrostatic interactions between SAM and nsp16. Thus, nsp10 acts as a stimulator for the strong binding of SAM to nsp16. The hydrophobic interactions were predominately identified for the nsp16-nsp10 interactions. Also, the stable hydrogen bonds between Ala83 (nsp16) and Tyr96 (nsp10), and between Gln87 (nsp16) and Leu45 (nsp10) play a vital role in the dimerization of nsp16 and nsp10. Besides, Computational Alanine Scanning (CAS) mutagenesis was performed, which revealed hotspot mutants, namely I40A, V104A, and R86A for the dimer association. Hence, the dimer interface of nsp16/nsp10 could also be a potential target in retarding the 2'-O-MTase activity in SARS-CoV-2. Overall, our study provides a comprehensive understanding of the dynamic and thermodynamic process of binding nsp16 and nsp10 that will contribute to the novel design of peptide inhibitors based on nsp16.

Exploring the potency of currently used drugs against HIV-1 protease of subtype D variant by using multiscale simulations

Acquired immune deficiency syndrome (AIDS) is caused by the human immunodeficiency virus (HIV), type 1 and 2. Further, the diversity in HIV-1 has given rise to many serotypes and recombinant strains. The currently used protease inhibitors have been developed for subtype B, although non-B subtype strains account for ∼ 90% of the global HIV infections. Subtype D is spreading rapidly and infecting a large population in North Africa and the Middle East. In the current study, molecular dynamics simulations in conjunction with the molecular mechanics/Poisson-Boltzmann surface area (MM-PBSA) scheme was used to investigate the potency of four drugs, namely atazanavir (ATV), darunavir (DRV), lopinavir (LPV) and tipranavir (TPV) against the subtype D variant. Our calculations predicted that the potency of the inhibitors decreased in the order TPV > ATV > DRV > LPV. TPV was found to be the most potent against subtype D due to an increase in van der Waals and electrostatic interactions and reduction in the desolvation energy compared to other inhibitors. This result is further supported by the hydrogen bond interactions between inhibitors and protease. Furthermore, our analyses suggested that the binding of TPV induced a more closed conformation of the flap compared to apo or other complexes. It was observed that TPV/PR^D has a lower cavity volume relative to the other three complexes leading to a tighter binding. The open conformation of the flap was observed for LPV/PR^D. We expect that this study might be useful for designing more potent inhibitors against HIV-1 subtype D. Communicated by Ramaswamy H. Sarma.

Investigating Phosphorylation-Induced Conformational Changes in WNK1 Kinase by Molecular Dynamics Simulations

The With-No-Lysine (WNK) kinase is considered to be a master regulator for various cation-chloride cotransporters involved in maintaining cell-volume and ion homeostasis. Here, we have investigated the phosphorylation-induced structural dynamics of the WNK1 kinase bound to an inhibitor via atomistic molecular dynamics simulations. Results from our simulations show that the phosphorylation at Ser³⁸² could stabilize the otherwise flexible activation loop (A-loop). The intrahelix salt-bridge formed between Arg²⁶⁴ and Glu²⁶⁸ in the unphosphorylated system is disengaged after the phosphorylation, and Glu²⁶⁸ reorients itself and forms a stable salt-bridge with Arg³⁴⁸. The dynamic cross-correlation analysis shows that phosphorylation diminishes anticorrelated motions and increases correlated motions between different domains. Structural network analysis reveals that the phosphorylation causes structural rearrangements and shortens the communication path between the αC-helix and catalytic loop, making the binding pocket more suitable for accommodating the ligand. Overall, we have characterized the structural changes in the WNK kinase because of phosphorylation in the A-loop, which might help in designing rational drugs.

Recognition and interaction of CphA from Aeromonas hydrophila with imipenem and biapenem by spectroscopic analysis in combination with molecular docking

The molecular recognition and interaction of CphA from Aeromonas hydrophila with imipenem (Imip) and biapenem (Biap) were studied by means of the combined use of fluorescence spectra and molecular docking. The results showed that both the fluorescence quenching of CphA by Imip and Biap were caused through the combined dynamic and static quenching, and the latter was dominating in the process; the microenvironment and conformational of CphA were altered upon the addition of Imip and Biap from synchronous and three-dimensional fluorescence. The binding of CphA with Imip or Biap caused a conformational change in the loop of CphA, and through the conformational change, the loop opened the binding pocket of CphA to allow for an induced fit of the newly introduced ligand. In the binding of CphA with Imip, the whole molecule entered into the active pocket of CphA. The binding was driven by enthalpy change, and the binding force between them was mainly hydrogen bonding and Van der Waals force; whereas in the binding of CphA with Biap, only the beta-lactam ring of Biap entered into the binding pocket of CphA while the side chain was located outside the active pocket. The binding was driven by the enthalpy change and entropy change together, and the binding force between them was mainly electrostatic interaction. This study provided an insight into the recognition and binding of CphA with antibiotics, which may be helpful for designing new substrate for beta-lactamase and developing new antibiotics resistant to superbugs.

Structural and dynamic basis of acid amido synthetase GH3.1: an investigation of substrate selectivity and major active site access channels

Auxin/IAA (Indole-3-acetic acid) plays critical roles in many aspects of plant growth and development. Gretchen Hagen 3.1 (GH3.1) enzyme from grapevine (Vitis vinifera) catalyzes the ATP-dependent conjugation of aspartate to IAA suggesting a significant way to modulate levels of cellular auxins/IAA. It is reported that VvGH3.1 prefers IAA as the substrate than 1-naphthaleneacetic acid (NAA) and benzothiazole-2-oxyacetic acid (BTOA), whereas the detailed interaction mechanism of these substrates remains unclear. In this study, based on the recently reported crystal structure of VvGH3.1 and AtGH3.12, the open form of VvGH3.1 was built. Then combined computational techniques including molecular docking, molecular dynamics (MD) simulations, Molecular Mechanics Generalized Born Surface Area and Poisson-Boltzmann Surface Area (MM-GB/PBSA) methods and active site access channel analysis were utilized to investigate the binding mechanism of three substrates, IAA, NAA and BTOA, and dynamic behaviors of GH3.1 induced by substrate binding. The predicted binding free energy is in agreement with the experimental work with an order of IAA < NAA < BTOA. Key residues interacting with three substrates and residues specifically binding to one of them were identified. The strong hydrogen bond interaction formed between IAA and Ser108, Ser339, Gln560 makes it the best substrate over NAA and BTOA. The dynamic behaviors, especially in the dominant access channels of IAA and NAA bound GH3.1, were found to be different from those in the BTOA-GH3.1 system with larger bottleneck radius and shorter length. This study provided novel insight to understand the substrate selectivity mechanism of GH3.1 and also for the rational design of novel auxin-based herbicides and growth regulators.

CUSA and CUDE: GPU-accelerated methods for estimating solvent accessible surface area and desolvation

It is well-established that a linear correlation exists between accessible surface areas and experimentally measured solvation energies. Combining this knowledge with an analytic formula for calculation of solvent accessible surfaces, we derive a simple model of desolvation energy as a differentiable function of atomic positions. Additionally, we find that this algorithm is particularly well suited for hardware acceleration on graphics processing units (GPUs), outperforming the CPU by up to two orders of magnitude. We explore the scaling of this desolvation algorithm and provide implementation details applicable to general pairwise algorithms.

Molecular simulations enlighten the binding mode of quercetin to lipoxygenase-3

Inhibition of lipoxygenases (LOXs) by flavonoid compounds is now well documented, but the description of the associated mechanism remains controversial due to a lack of information at the molecular level. For instance, X-ray determination of quercetin/LOX-3 system has led to a structure where the enzyme was cocrystallized with a degradation product of the substrate, which rendered the interpretation of the reported interactions between this flavonoid compound and the enzyme difficult. Molecular modeling simulations can in principle allow obtaining precious insights that could fill this lack of structural information. Thus, in this study, we have investigated various binding modes of quercetin to LOX-3 enzyme in order to understand the first step of the inhibition process, that is the association of the two entities. Molecular dynamics simulations and free energy calculations suggest that quercetin binds the metal center via its 3-hydroxychromone function. Moreover, enzyme/substrate interactions within the cavity impose steric hindrances to quercetin that may activate a direct dioxygen addition on the substrate.

Conformational analysis of trimeric maleimide substituted 1,5,9-triazacyclododecane HIV fusion scaffolds

An analysis of the conformational preferences of three trimeric maleimide substituted 1,5,9-triazacyclododecane derivatives, proposed as cross linking reagents for HIV-1 fusion inhibitors, is presented. Exhaustive sampling was performed using the mixed Low Mode Monte Carlo conformational searching technique on the corresponding OPLS2005/GBSA(water) potential energy surface. Geometric structure, molecular length, and hydrogen bonding patterns of the compounds are analyzed. Global minimum energy structures were verified as minima using B3LYP/6-31G * geometry optimization. All structures adopt a crown-like 12-membered ring conformation; however, the system with the shortest maleimide arms (1a) can also adopt alternative ring orientations. Overall, derivatives with longer maleimide arms were more flexible and resulted in ensembles with a larger number of low energy structures. Comparison with biological inhibition data indicates that there is very little relationship between molecular size and the ability of the scaffold to orient CD4M9 miniproteins for optimal inhibition; however hydrophobicity may play a role.

In silico profiling of tyrosine kinases binding specificity and drug resistance using Monte Carlo simulations with the ensembles of protein kinase crystal structures

The molecular basis of the tyrosine kinases binding specificity and drug resistance against cancer drugs Imatinib and Dasatinib is elucidated using Monte Carlo simulations of the inhibitor-receptor binding with the ensembles of protein kinase crystal structures. In silico proteomics analysis unravels mechanisms by which structural plasticity of the tyrosine kinases is linked with the conformational preferences of Imatinib and Dasatinib in achieving effective drug binding with a distinct spectrum of the tyrosine kinome. The differences in the inhibitor sensitivities to the ABL kinase mutants are rationalized based on variations in the binding free energy profiles with the conformational states of the ABL kinase. While Imatinib binding is highly sensitive to the activation state of the enzyme, the computed binding profile of Dasatinib is remarkably tolerant to the conformational state of ABL. A comparative analysis of the inhibitor binding profiles with the clinically important ABL kinase mutants has revealed an excellent agreement with the biochemical and proteomics data. We have found that conformational adaptability of the kinase inhibitors to structurally different conformational states of the tyrosine kinases may have pharmacological relevance in acquiring a specific array of potent activities and regulating a scope of the inhibitor resistance mutations. This study outlines a useful approach for understanding and predicting the molecular basis of the inhibitor sensitivity against potential kinase targets and drug resistance.

Computational proteomics of biomolecular interactions in the sequence and structure space of the tyrosine kinome: Deciphering the molecular basis of the kinase inhibitors selectivity

Understanding and predicting the molecular basis of protein kinases specificity against existing therapeutic agents remains highly challenging and deciphering this complexity presents an important problem in discovery and development of effective cancer drugs. We explore a recently introduced computational approach for in silico profiling of the tyrosine kinases binding specificity with a class of the pyrido-[2,3-d]pyrimidine kinase inhibitors. Computational proteomics analysis of the ligand-protein interactions using parallel simulated tempering with an ensemble of the tyrosine kinases crystal structures reveals an important molecular determinant of the kinase specificity. The pyrido-[2,3-d]pyrimidine inhibitors are capable of dynamically interacting with both active and inactive forms of the tyrosine kinases, accommodating structurally different kinase conformations with a similar binding affinity. Conformational tolerance of the protein tyrosine kinases binding with the pyrido[2,3-d]pyrimidine inhibitors provides the molecular basis for the broad spectrum of potent activities and agrees with the experimental inhibition profiles. The analysis of the pyrido[2,3-d]pyrimidine sensitivities against a number of clinically relevant ABL kinase mutants suggests an important role of conformational adaptability of multitargeted kinase inhibitors in developing drug resistance mechanisms. The presented computational approach may be useful in complementing proteomics technologies to characterize activity signatures of small molecules against a large number of potential kinase targets.

Cages, Baskets, Ladders, and Tubes: Conformational Studies of Polyhedral Oligomeric Silsesquioxanes

The conformational flexibility of a series of cage, basket, ladder, and tube polyhedral oligomeric silsesquioxanes (POSS) has been examined using the Low Mode:Monte Carlo conformational search method in conjunction with the MM3/GBSA(CHCl3) surface. An ensemble of low energy structures was generated and used to explore the molecular shape and flexibility of each system. The results indicate that, except for the ladder molecule, the incompletely condensed systems that are studied are relatively rigid. Even in cases where the molecule is able to adopt numerous low energy conformations, the overall shape remains cage-like and the conformations differ only by small angles or substituent orientations. The ladder molecule is the most flexible and this ensemble clusters into two families: one that is cage-like and the other that is more open and ladder-like. The conformational flexibilities in the gas and solvent phases, as approximated using the GBSA continuum solvent model, are very similar.

Computational analysis of ligand binding dynamics at the intermolecular hot spots with the aid of simulated tempering and binding free energy calculations

Equilibrium binding dynamics is studied for a panel of benzimidazole-containing compounds at the remodeled interface between human growth hormone (hGH) and the extracellular domain of its receptor (hGHbp), engineered by mutating to glycine hot spot residues T175 from the hormone and W104 from the receptor. Binding energetics is predicted in a good agreement with the experimental data for a panel of these small molecules that complement the engineered defect and restore the binding affinity of the wild-type hGH-hGHbp complex. The results of simulated tempering ligand dynamics at the protein-protein interface reveals a diversity of ligand binding modes that is consistent with the structural orientation of the benzimidazole ring which closely mimics the position of the mutated W104 hot spot residue in the wild-type hGH-hGHbp complex. This structural positioning of the benzimidazole core motif is shown to be a critical feature of the low-energy ligand conformations binding in the engineered cavity. The binding free energy analysis provides a plausible rationale behind the experimental dissociation constants and suggests a key role of ligand-protein van der Waals interactions in restoring binding affinity.

Computational detection of the binding-site hot spot at the remodeled human growth hormone-receptor interface

A hierarchical computational approach is used to identify the engineered binding-site cavity at the remodeled intermolecular interface between the mutants of human growth hormone (hGH) and the extracellular domain of its receptor (hGHbp). Multiple docking simulations are conducted with the remodeled hGH-hGHbp complex for a panel of potent benzimidazole-containing inhibitors that can restore the binding affinity of the wild-type complex, and for a set of known nonactive small molecules that contain different heterocyclic motifs. Structural clustering of ligand-bound conformations and binding free-energy calculations, using the AMBER force field and a continuum solvation model, can rapidly locate and screen numerous ligand-binding modes on the protein surface and detect the binding-site hot spot at the intermolecular interface. Structural orientation of the benzimidazole motif in the binding-site cavity closely mimics the position of the hot spot residue W104 in the crystal structure of the wild-type complex, which is recognized as an important structural requirement for restoring binding affinity. Despite numerous pockets on the protein surface of the mutant hGH-hGHbp complex, the binding-site cavity presents the energetically favorable hot spot for the benzimidazole-containing inhibitors, whereas for a set of nonactive molecules, the lowest energy ligand conformations do not necessarily bind in the engineered cavity. The results reveal a dominant role of the intermolecular van der Waals interactions in providing favorable ligand-protein energetics in the redesigned interface, in agreement with the experimental and computational alanine scanning of the hGH-hGHbp complex.

A comparison of the Low Mode and Monte Carlo conformational search methods

The Low Mode (LM) and Monte Carlo (MC) conformational search methods were compared on three diverse molecular systems; (4R, 5S, 6S, 7R)-hexahydro-5,6-dihydroxy-1,3,4,7-tetrakis(phenylmethyl)-2H-1,3-diazapin-2-one (1), 2-methoxy-2-phenyl-2-triflouromethyl-N-alpha-methyl benzyl propanamide (2) and a trimeric 39-membered polyazamacrolide (3). We find that either method, or a combination of the methods, is equally efficient at searching the conformational space of the smaller molecular systems while a 50:50 hybrid of Low Mode and Monte Carlo is most efficient at searching the space of the larger molecular system.

Monte Carlo simulations of the peptide recognition at the consensus binding site of the constant fragment of human immunoglobulin G: the energy landscape analysis of a hot spot at the intermolecular interface

Monte Carlo simulations of molecular recognition at the consensus binding site of the constant fragment (Fc) of human immunoglobulin G (Ig) protein have been performed to analyze structural and thermodynamic aspects of binding for the 13-residue cyclic peptide DCAWHLGELVWCT. The energy landscape analysis of a hot spot at the intermolecular interface using alanine scanning and equilibrium-simulated tempering dynamics with the simplified, knowledge-based energy function has enabled the role of the protein hot spot residues in providing the thermodynamic stability of the native structure to be determined. We have found that hydrophobic interactions between the peptide and the Met-252, Ile-253, His-433, and His-435 protein residues are critical to guarantee the thermodynamic stability of the crystallographic binding mode of the complex. Binding free energy calculations, using a molecular mechanics force field and a solvation energy model, combined with alanine scanning have been conducted to determine the energetic contribution of the protein hot spot residues in binding affinity. The conserved Asn-434, Ser-254, and Tyr-436 protein residues contribute significantly to the binding affinity of the peptide-protein complex, serving as an energetic hot spot at the intermolecular interface. The results suggest that evolutionary conserved hot spot protein residues at the intermolecular interface may be partitioned in fulfilling thermodynamic stability of the native binding mode and contributing to the binding affinity of the complex.

Application of the frozen atom approximation to the GB/SA continuum model for solvation free energy

The generalized Born/surface area (GB/SA) continuum model for solvation free energy is a fast and accurate alternative to using discrete water molecules in molecular simulations of solvated systems. However, computational studies of large solvated molecular systems such as enzyme-ligand complexes can still be computationally expensive even with continuum solvation methods simply because of the large number of atoms in the solute molecules. Because in such systems often only a relatively small portion of the system such as the ligand binding site is under study, it becomes less attractive to calculate energies and derivatives for all atoms in the system. To curtail computation while still maintaining high energetic accuracy, atoms distant from the site of interest are often frozen; that is, their coordinates are made invariant. Such frozen atoms do not require energetic and derivative updates during the course of a simulation. Herein we describe methodology and results for applying the frozen atom approach to both the generalized Born (GB) and the solvent accessible surface area (SASA) parts of the GB/SA continuum model for solvation free energy. For strictly pairwise energetic terms, such as the Coulombic and van-der-Waals energies, contributions from pairs of frozen atoms can be ignored. This leaves energetic differences unaffected for conformations that vary only in the positions of nonfrozen atoms. Due to the nonlocal nature of the GB analytical form, however, excluding such pairs from a GB calculation leads to unacceptable inaccuracies. To apply a frozen-atom scheme to GB calculations, a buffer region within the frozen-atom zone is generated based on a user-definable cutoff distance from the nonfrozen atoms. Certain pairwise interactions between frozen atoms in the buffer region are retained in the GB computation. This allows high accuracy in conformational GB comparisons to be maintained while achieving significant savings in computational time compared to the full (nonfrozen) calculation. A similar approach for using a buffer region of frozen atoms is taken for the SASA calculation. The SASA calculation is local in nature, and thus exact SASA energies are maintained. With a buffer region of 8 A for the frozen-atom cases, excellent agreement in differences in energies for three different conformations of cytochrome P450 with a bound camphor ligand are obtained with respect to the nonfrozen cases. For various minimization protocols, simulations run 2 to 10.5 times faster and memory usage is reduced by a factor of 1.5 to 5. Application of the frozen atom method for GB/SA calculations thus can render computationally tractable biologically and medically important simulations such as those used to study ligand-receptor binding conformations and energies in a solvated environment.

Approximate solvent-accessible surface areas from tetrahedrally directed neighbor densities

A fast analytical formula (TDND) has been derived for the calculation of approximate atomic and molecular solvent-accessible surface areas (SASA), as well as the first and second derivatives of these quantities with respect to atomic coordinates. Extending the work of Stouten et al. (Molecular Simulation, 1993, Vol. 10, pp. 97-120), as well as our own (Journal of Computational Chemistry, 1999, Vol. 20, pp. 586-596), the method makes use of a Gaussian function to calculate the neighbor density in four tetrahedral directions in three-dimensional space, sometimes twice with different orientations. SASA and first derivatives of the 2366 heavy atoms of penicillopepsin are computed in 0.13 s on an SGI R10000/194 MHz processor. When second derivatives are computed as well, the total time is 0.23 s. This is considerably faster than timings reported previously for other algorithms. Based on a parameterization set of nineteen compounds of different size (11-4346 atoms) and class (organics, proteins, DNA, and various complexes) consisting of a total 23,197 atoms, the method exhibits relative errors in the range 0.2-12.6% for total molecular surface areas and average absolute atomic surface area deviations in the range 1.7 to 3.6 A(2).