Long-range electrostatic interactions in molecular dynamics: an endothelin-1 case study

Identification of potential inhibitors for drug resistance in acute lymphoblastic leukemia through differentially expressed gene analysis and in silico screening

Acute lymphoblastic leukemia (ALL) is a disease of lymphocyte origin predominantly diagnosed in children. While its 5-year survival rate is high, resistance to chemotherapy drugs is still an obstacle. Our aim is to determine differentially expressed genes (DEGs) related to Asparaginase, Daunorubicin, Prednisolone, and Vincristine resistance and identify potential inhibitors via docking. Three datasets were accessed from the Gene Expression Omnibus database; GSE635, GSE19143, and GSE22529. The microarray data was analyzed using R4.2.0 and Bioconductor packages, and pathway and protein-protein interaction analysis were performed. We identified 1294 upregulated DEGs, with 12 genes consistently upregulated in all four resistant groups. KEGG analysis revealed an association with the PI3K-Akt pathway. Among DEGs, 33 hub genes including MDM2 and USP7 were pinpointed. Within common genes, CLDN9 and HS3ST3A1 were subjected to molecular docking against 3556 molecules. Following ADMET analysis, three drugs emerged as potential inhibitors: Flunarizine, Talniflumate, and Eltrombopag. Molecular dynamics analysis for HS3ST3A1 indicated all candidates had the potential to overcome drug resistance, Eltrombopag displaying particularly promising results. This study promotes a further understanding of drug resistance in ALL, introducing novel genes for consideration in diagnostic screening. It also presents potential inhibitor candidates to tackle drug resistance through repurposing.

Analysis of Glycan Recognition by Concanavalin A Using Absolute Binding Free Energy Calculations

Carbohydrates are key biological mediators of molecular recognition and signaling processes. In this case study, we explore the ability of absolute binding free energy (ABFE) calculations to predict the affinities of a set of five related carbohydrate ligands for the lectin protein, concanavalin A, ranging from 27-atom monosaccharides to a 120-atom complex-type N-linked glycan core pentasaccharide. ABFE calculations quantitatively rank and estimate the affinity of the ligands in relation to microcalorimetry, with a mean signed error in the binding free energy of -0.63 ± 0.04 kcal/mol. Consequently, the diminished binding efficiencies of the larger carbohydrate ligands are closely reproduced: the ligand efficiency values from isothermal titration calorimetry for the glycan core pentasaccharide and its constituent trisaccharide and monosaccharide compounds are respectively -0.14, -0.22, and -0.41 kcal/mol per heavy atom. ABFE calculations predict these ligand efficiencies to be -0.14 ± 0.02, -0.24 ± 0.03, and -0.46 ± 0.06 kcal/mol per heavy atom, respectively. Consequently, the ABFE method correctly identifies the high affinity of the key anchoring mannose residue and the negligible contribution to binding of both β-GlcNAc arms of the pentasaccharide. While challenges remain in sampling the conformation and interactions of these polar, flexible, and weakly bound ligands, we nevertheless find that the ABFE method performs well for this lectin system. The approach shows promise as a quantitative tool for predicting and deconvoluting carbohydrate-protein interactions, with potential application to design of therapeutics, vaccines, and diagnostics.

DFT and molecular simulation validation of the binding activity of PDEδ inhibitors for repression of oncogenic k-Ras

The development of effective drugs targeting the K-Ras oncogene product is a significant focus in anticancer drug development. Despite the lack of successful Ras signaling inhibitors, recent research has identified PDEδ, a KRAS transporter, as a potential target for inhibiting the oncogenic KRAS signaling pathway. This study aims to investigate the interactions between eight K-Ras inhibitors (deltarazine, deltaflexin 1 and 2, and its analogues) and PDEδ to understand their binding modes. The research will utilize computational techniques such as density functional theory (DFT) and molecular electrostatic surface potential (MESP), molecular docking, binding site analyses, molecular dynamic (MD) simulations, electronic structure computations, and predictions of the binding free energy. Molecular dynamic simulations (MD) will be used to predict the binding conformations and pharmacophoric features in the active site of PDEδ for the examined structures. The binding free energies determined using the MMPB(GB)SA method will be compared with the observed potency values of the tested compounds. This computational approach aims to enhance understanding of the PDEδ selective mechanism, which could contribute to the development of novel selective inhibitors for K-Ras signaling.

AMOEBA Force Field Predicts Accurate Hydrogen Bond Counts of Nitriles in SNase by Revealing Water-Protein Interaction in Vibrational Absorption Frequencies

Precisely quantifying the magnitude and direction of electric fields in proteins has long been an outstanding challenge in understanding biological functions. Nitrile vibrational Stark effect probes have been shown to be minimally disruptive to the protein structure and can be better direct reporters of local electrostatic field in the native state of a protein than other measures such as pKa shifts of titratable residues. However, interpretations of the connection between measured vibrational energy and electric field rely on the accurate molecular understanding of interactions of the nitrile group and its environment, particularly from hydrogen bonding. In this work, we compared the extent of hydrogen bonding calculated in two common force fields, the fixed charge force field Amber03 and polarizable force field AMOEBA, at 10 locations of cyanocysteine (CNC) in staphylococcal nuclease (SNase) against the experimental nitrile absorption frequency in terms of full width at half-maximum (FWHM) and frequency temperature line slope (FTLS). We observed that the number of hydrogen bonds correlated well in AMOEBA trajectories with respect to both the FWHM (r = 0.88) and the FTLS (r = -0.85), whereas the correlation of Amber03 trajectories was less reliable because the Amber03 force field predicted more hydrogen bonds in some mutants. Moreover, we demonstrated that contributions from the interactions between CNC and nearby water molecules were significant in AMOEBA trajectories but were not predicted by Amber03. We conclude that although the nitrile absorption peak shape could be qualitatively predicted by the fixed charge Amber03 force field, the detailed electrostatic environment measured by the nitrile probe in terms of the extent of hydrogen bonding could only be accurately observed in the AMOEBA trajectories, where the permanent dipole, quadrupole, and dipole-induced-dipole polarizable interactions were all taken into account. The significance of this finding to the goal of accurately predicting electric fields in complex biomolecular environments is discussed.

AMOEBA Force Field Trajectories Improve Predictions of Accurate pKa Values of the GFP Fluorophore: The Importance of Polarizability and Water Interactions

Precisely quantifying the magnitude, direction, and biological functions of electric fields in proteins has long been an outstanding challenge in the field. The most widely implemented experimental method to measure such electric fields at a particular residue in a protein has been through changes in pKa of titratable residues. While many computational strategies exist to predict these values, it has been difficult to do this accurately or connect predicted results to key structural or mechanistic features of the molecule. Here, we used experimentally determined pKa values of the fluorophore in superfolder green fluorescent protein (GFP) with amino acid mutations made at position Thr 203 to evaluate the pKa prediction ability of molecular dynamics (MD) simulations using a polarizable force field, AMOEBA. Structure ensembles from AMOEBA were used to calculate pKa values of the GFP fluorophore. The calculated pKa values were then compared to trajectories using a conventional fixed charge force field (Amber03 ff). We found that the position of water molecules included in the pKa calculation had opposite effects on the pKa values between the trajectories from AMOEBA and Amber03 force fields. In AMOEBA trajectories, the inclusion of water molecules within 35 Å of the fluorophore decreased the difference between the predicted and experimental values, resulting in calculated pKa values that were within an average of 0.8 pKa unit from the experimental results. On the other hand, in Amber03 trajectories, including water molecules that were more than 5 Å from the fluorophore increased the differences between the calculated and experimental pKa values. The inaccuracy of pKa predictions determined from Amber03 trajectories was caused by a significant stabilization of the deprotonated chromophore's free energy compared to the result in AMOEBA. We rationalize the cutoffs for explicit water molecules when calculating pKa to better predict the electrostatic environment surrounding the fluorophore buried in GFP. We discuss how the results from this work will assist the prospective prediction of pKa values or other electrostatic effects in a wide variety of folded proteins.

The impact of cycleanine in cancer research: a computational study

Cancer is imposing a global health burden because of the steady increase in new cases. Moreover, current anticancer therapeutics are associated with many drawbacks, mainly the emergence of resistance and the severe adverse effects. Therefore, there is a continuous need for developing new anticancer agents with novel mechanisms of action and lower side effects. Natural products have been a rich source of anticancer medication. Cycleanine, a natural product, was reported to exert an antiproliferative effect on ovarian cancer cells by causing apoptosis through activation of caspases 3/7 and cleavage of poly (ADP-ribose) polymerase to form poly (ADP-ribose) polymerase-1 (PARP1). It is well-established that PARP1 is associated with carcinogenesis, and different PARP1 inhibitors are approved as anticancer drugs. In this study, the cytotoxic activity of cycleanine was computationally investigated to determine whether it is a PARP1 inhibitor or a caspase activator. Molecular docking and molecular dynamics (MD) simulations were utilized for this purpose. The results showed that cycleanine has a good binding affinity to PARP1; moreover, MD simulation showed that it forms a stable complex with the enzyme. Consequently, the results showed that cycleanine is a potential inhibitor of the PARP1 enzyme.

1,2,3-Triazolyl-tetrahydropyrimidine Conjugates as Potential Sterol Carrier Protein-2 Inhibitors: Larvicidal Activity against the Malaria Vector Anopheles arabiensis and In Silico Molecular Docking Study

Alteration of insect growth regulators by the action of inhibitors is becoming an attractive strategy to combat disease-transmitting insects. In the present study, we investigated the larvicidal effect of 1,2,3-triazolyl-pyrimidinone derivatives against the larvae of the mosquito Anopheles arabiensis, a vector of malaria. All compounds demonstrated insecticidal activity against mosquito larvae in a dose-dependent fashion. A preliminary study of the structure-activity relationship indicated that the electron-withdrawing substituent in the para position of the 4-phenyl-pyrimidinone moiety enhanced the molecules' potency. A docking study of these derivatives revealed favorable binding affinity for the sterol carrier protein-2 receptor, a protein present in the intestine of the mosquito larvae. Being effective insecticides against the malaria-transmitting Anopheles arabiensis, 1,2,3-triazole-based pyrimidinones represent a starting point to develop novel inhibitors of insect growth regulators.

The design of TOPK inhibitors using structure-based pharmacophore modeling and molecular docking based on an MD-refined homology model

The TOPK enzyme (also known as PBK) is a serine-threonine protein kinase that is rarely detected in normal tissues yet is found to be overexpressed and activated in a variety of cancers such as lung, colorectal, breast, and esophageal cancer. Its prevalence in cancerous cells is associated with their poor prognosis and responsiveness to treatment. This enzyme plays a vital role in cell division, specifically in regulating cytokinesis. Unlike drugs targeting early phases in mitosis, inhibition of cytokinesis by targeting biomolecules that are unique to multiplying cells poses no threat to the normal function of non-multiplying cells. Studies have shown that inhibition of cytokinesis is promising in suppressing the growth of proliferating cancerous cells as exemplified by the complete tumor regression seen with the suppression of TOPK. Herein, we report the identification of potent TOPK inhibitors with anticancer potential using a structure-based drug design approach. The only available crystal structure of TOPK corresponds to a double mutant (T9E and T198E) dimer with a distorted N-lobe conformation, thus 3D homology modeling was implemented to rebuild the enzyme's native conformation. The resulting refined model was used to generate 3D pharmacophore models for the virtual screening of small molecules databases. Retrieved hits were filtered, docked into the ATP binding site of the enzyme, rescored, and the binding free energies for the top consensually scoring hits were calculated. Consequently, 45 compounds were selected and their in vitro inhibitory activity against TOPK was tested revealing four potential hits with the most active compound having an IC₅₀ of 3.85 µM. This compound will be chosen as a lead compound to synthesize analogs aiming to enhance potency and drug-like properties and to enrich the SAR data.

In vitro anti-TB properties, in silico target validation, molecular docking and dynamics studies of substituted 1,2,4-oxadiazole analogues against Mycobacterium tuberculosis

The alarming increase in multi- and extensively drug-resistant (MDR and XDR) strains of Mycobacterium tuberculosis (MTB) has triggered the scientific community to search for novel, effective, and safer therapeutics. To this end, a series of 3,5-disubstituted-1,2,4-oxadiazole derivatives (3a-3i) were tested against H37Rv, MDR and XDR strains of MTB. Of which, compound 3a with para-trifluorophenyl substituted oxadiazole showed excellent activity against the susceptible H37Rv and MDR-MTB strain with a MIC values of 8 and 16 µg/ml, respectively.To understand the mechanism of action of these compounds (3a-3i) and identify their putative drug target, molecular docking and dynamics studies were employed against a panel of 20 mycobacterial enzymes reported to be essential for mycobacterial growth and survival. These computational studies revealed polyketide synthase (Pks13) enzyme as the putative target. Moreover, in silico ADMET predictions showed satisfactory properties for these compounds, collectively, making them, particularly compound 3a, promising leads worthy of further optimisation.

Structural characterization of cassava linamarase-linamarin enzyme complex: an integrated computational approach

Cassava linamarase is a hydrolyzing enzyme that belongs to a glycoside hydrolase family 1 (GH1). It is responsible for breaking down linamarin to toxic cyanide. The enzyme provides a defensive mechanism for plants against herbivores and has various applications in many fields. Understanding the structure of linamarase at the molecular level is a key to avail its reaction mechanism. In this study, the three-dimensional (3D) structure of linamarase was built for the first time using homology modelling and used to study its interaction with linamarin. Molecular docking calculations established the binding and orientation nature of linamarin, while molecular dynamics (MD) simulation established protein-ligand complexes' stability. Binding-free energy based on MM/PBSA was further used to rescore the docking results. An ensemble structure was found to be relatively stable compared to the modelled structure. This study sheds light on the exploration of linamarase towards understanding its reaction mechanisms.Communicated by Ramaswamy H. Sarma.

CHK1 kinase inhibition: identification of allosteric hits using MD simulations, pharmacophore modeling, docking and MM-PBSA calculations

The CHK1 kinase plays a pivotal role in the DNA damage response pathway. Hence, inhibition of CHK1 appeared as a promising strategy to overcome the resistance problem of chemotherapeutic agents resulting from the overexpression of CHK1 that enables cancerous cells to repair their chemotherapy-induced DNA damage. In this study, different computational drug design techniques were employed to identify new CHK1 inhibitors targeting its allosteric pocket. A 1 μs MD simulation of the apo form of the enzyme was run to study its native dynamics. The resulting trajectory was analyzed to select a frame where the ATP binding pocket is most occluded while its allosteric counterpart is most exposed to be used in the design of potential allosteric inhibitors that could trap the enzyme in such nearly inactive state. Besides the selected frame, another three crystal structures of CHK1 complexed with allosteric inhibitors were utilized to generate structure-based pharmacophore models. Seven pharmacophores were generated and utilized in virtual screening of different databases. The retrieved hits were filtered and then docked into the allosteric pocket. Finally, the binding energies of the top-ranked docked hits were calculated. Twenty compounds were selected as candidates for biological evaluation against CHK1 enzyme. The biological screening results showed moderate activities where the percentage of CHK1 inhibition ranged from zero to 28.26%. Four of the tested compounds showed percentage of CHK1 inhibition greater than 20%, of which, two compounds were identified as allosteric hits that upon further optimization could be converted into lead-like compounds.

Tackling COVID-19: identification of potential main protease inhibitors via structural analysis, virtual screening, molecular docking and MM-PBSA calculations

The widespread of the COVID-19 disease, caused by the novel severe acute respiratory syndrome coronavirus (SARS-CoV-2), had severely affected the entire world. Unfortunately, no successful vaccines or antiviral drugs are currently available which leaves the scientific community under huge pressure to tackle this pandemic. Among the identified promising druggable targets, specific to this virus, is the main protease (M^pro) enzyme, which is vital for viral replication, transcription and packaging within the host cells. In this study, selective inhibition of the M^pro was sought via thorough analysis of its available structural data in the Protein Data Bank. To this end, COVID-19 M^pro crystal complexes were explored and the key interacting residues (KIRs) within its active site, that are expected to be vital for effective ligand binding, were identified. Based on these KIRs, 3D pharmacophore models were generated and used in virtual screening of different databases. Retrieved hits were docked into the active site of the enzyme and their MM-PBSA based free binding energies were calculated. Finally, ADMET descriptors were calculated to aid the selection of top scoring hits with best ADMET properties. Nine compounds with different chemotypes were identified as potential M^pro inhibitors. Further, MD simulations of a virtual complex of M^pro with one of the promising hits revealed stable binding which is indicative of good inhibitory potential. The identified compounds in this study are expected to support the global drug discovery efforts in fighting against this highly contagious virus by narrowing the searchable chemical space for potential effective therapeutics.Communicated by Ramaswamy H. Sarma.

Scale-Space Splatting: Reforming Spacetime for Cross-Scale Exploration of Integral Measures in Molecular Dynamics

Understanding large amounts of spatiotemporal data from particle-based simulations, such as molecular dynamics, often relies on the computation and analysis of aggregate measures. These, however, by virtue of aggregation, hide structural information about the space/time localization of the studied phenomena. This leads to degenerate cases where the measures fail to capture distinct behaviour. In order to drill into these aggregate values, we propose a multi-scale visual exploration technique. Our novel representation, based on partial domain aggregation, enables the construction of a continuous scale-space for discrete datasets and the simultaneous exploration of scales in both space and time. We link these two scale-spaces in a scale-space space-time cube and model linked views as orthogonal slices through this cube, thus enabling the rapid identification of spatio-temporal patterns at multiple scales. To demonstrate the effectiveness of our approach, we showcase an advanced exploration of a protein-ligand simulation.

Sparsification of long range force networks for molecular dynamics simulations

Atomic interactions in solid materials are described using network theory. The tools of network theory focus on understanding the properties of a system based upon the underlying interactions which govern their dynamics. While the full atomistic network is dense, we apply a spectral sparsification technique to construct a sparse interaction network model that reduces the computational complexity while preserving macroscopic conservation properties. This sparse network is compared to a reduced network created using a cut-off radius (threshold method) that is commonly used to speed-up computations while approximating interatomic forces. The approximations used to estimate the total forces on each atom are quantified to assess how local interatomic force errors propagate errors at the global or continuum scale by comparing spectral sparsification to thresholding. In particular, we quantify the performance of the spectral sparsification algorithm for the short-range Lennard-Jones potential and the long-range Coulomb potential. Spectral sparsification of the Lennard–Jones potential yields comparable results to thresholding while spectral sparsification yields improvements when considering a long-range Coulomb potential. The present network-theoretic formulation is implemented on two sample problems: relaxation of atoms near a surface and a tensile test of a solid with a circular hole.

Enhancement of levodopa stability when complexed with β-cyclodextrin in transdermal patches

Levodopa is a promising candidate for administration via the transdermal route because it exhibits a short plasma half-life and has a small window of absorption in the upper section of the small intestine. The aim of this study was to prepare stable levodopa transdermal patches. Both xanthan gum and Carbopol 971 polymers were selected with ethylcellulose constituting the backing layer of the prepared patches. The effect of adding β-cyclodextrin on the prepared patches was investigated. The uniformity in thickness, weight and content of the studied patches was acceptable. Physicochemical characterization revealed that there was no interaction between levodopa and the applied polymer. The results proved that levodopa precipitated as an amorphous form in carbopol patches. Controlled drug release was achieved for all the tested patches over a 6 h period. However, increased permeation was achieved for the carbopol patches. Although cyclodextrin did not enhance levodopa permeation, the stability study confirmed that levodopa stability was enhanced when complexed with β-cyclodextrin. The cumulative amount of drug released from carbopol patches is slightly higher than that of xanthan patches. The optimal stability was achieved in the carbopol/levodopa:β-cyclodextrin patch. The levodopa-β-cyclodextrin complex was successfully characterized using X-ray diffraction, NMR analysis and molecular dynamics simulations. In conclusion, carbopol/levodopa:β-cyclodextrin patches can be considered as a promising stable and effective transdermal drug-delivery system.

Explaining the autoinhibition of the SMYD enzyme family: A theoretical study

The SMYD enzymes (SMYD1-5) are lysine methyltransferases that have diverse biological functions including gene expression and regulation of skeletal and cardiac muscle development and function. Recently, they have gained more attention as potential drug targets because of their involvement in cardiovascular diseases and in the progression of different cancer types. Their activity has been suggested to be regulated by a posttranslational mechanism and by autoinhibition. The later relies on a hinge-like movement of the N- and C-lobes to adopt an open or closed conformation, consequently, determining the accessibility of the active site and substrate specificity. In this study we aim to investigate and explain the possibility of the regulatory autoinhibition process of the SMYD enzymes by a thorough computational exploration of their dynamic, energetic, and structural changes by using extended molecular dynamics simulations; normal mode analysis (NMA); and energy correlations. Three SMYD models (SMYD1-3) were used in this study. Our results showed an obvious hinge-like motion between the N- and C-lobes. Also, we identified interaction energy pathways within the 3D structures of the proteins, and hot spots on their surfaces that could be of particular importance for the regulation of their activities via allosteric means. These results can help in a better understanding of the nature of these promising drug targets; and in designing selective drugs that can interfere with (inhibit) the function of a specific SMYD member by disrupting its dynamical and conformational behaviour without disrupting the function of the entire SMYD proteins.

Nonbond interactions between graphene nanosheets and polymers: a computational study

Based on the geometries from molecular dynamics simulations and a package compiled by us, the interactions between graphene nanosheet (GNS) and nine types of flexible polymers have been investigated with force field. Both the van der Waals (vdW) interaction and the electrostatic interaction (EI) for two same polymer chains and between a polymer and a GNS were calculated and compared. The effect of cut-off distance was explored. It was found that the cut-off distance plays a significant role in EI energy, but a less important role in vdW energy when the cut-off distance is over 9.5 Å. The reasonable cut-off distances for EI and vdW interactions for simulation are suggested.

The role of electrostatic energy in prediction of obligate protein-protein interactions

BACKGROUND

Prediction and analysis of protein-protein interactions (PPI) and specifically types of PPIs is an important problem in life science research because of the fundamental roles of PPIs in many biological processes in living cells. In addition, electrostatic interactions are important in understanding inter-molecular interactions, since they are long-range, and because of their influence in charged molecules. This is the main motivation for using electrostatic energy for prediction of PPI types.

RESULTS

We propose a prediction model to analyze protein interaction types, namely obligate and non-obligate, using electrostatic energy values as properties. The prediction approach uses electrostatic energy values for pairs of atoms and amino acids present in interfaces where the interaction occurs. The main features of the complexes are found and then the prediction is performed via several state-of-the-art classification techniques, including linear dimensionality reduction (LDR), support vector machine (SVM), naive Bayes (NB) and k-nearest neighbor (k-NN). For an in-depth analysis of classification results, some other experiments were performed by varying the distance cutoffs between atom pairs of interacting chains, ranging from 5Å to 13Å. Moreover, several feature selection algorithms including gain ratio (GR), information gain (IG), chi-square (Chi2) and minimum redundancy maximum relevance (mRMR) are applied on the available datasets to obtain more discriminative pairs of atom types and amino acid types as features for prediction.

CONCLUSIONS

Our results on two well-known datasets of obligate and non-obligate complexes confirm that electrostatic energy is an important property to predict obligate and non-obligate protein interaction types on the basis of all the experimental results, achieving accuracies of over 98%. Furthermore, a comparison performed by changing the distance cutoff demonstrates that the best values for prediction of PPI types using electrostatic energy range from 9Å to 12Å, which show that electrostatic interactions are long-range and cover a broader area in the interface. In addition, the results on using feature selection before prediction confirm that (a) a few pairs of atoms and amino acids are appropriate for prediction, and (b) prediction performance can be improved by eliminating irrelevant and noisy features and selecting the most discriminative ones.

Potential human cholesterol esterase inhibitor design: benefits from the molecular dynamics simulations and pharmacophore modeling studies

Human pancreatic cholesterol esterase (hCEase) is one of the lipases found to involve in the digestion of large and broad spectrum of substrates including triglycerides, phospholipids, cholesteryl esters, etc. The presence of bile salts is found to be very important for the activation of hCEase. Molecular dynamic simulations were performed for the apoform and bile salt complexed form of hCEase using the co-ordinates of two bile salts from bovine CEase. The stability of the systems throughout the simulation time was checked and two representative structures from the highly populated regions were selected using cluster analysis. These two representative structures were used in pharmacophore model generation. The generated pharmacophore models were validated and used in database screening. The screened hits were refined for their drug-like properties based on Lipinski's rule of five and ADMET properties. The drug-like compounds were further refined by molecular docking simulation using GOLD program based on the GOLD fitness score, mode of binding, and molecular interactions with the active site amino acids. Finally, three hits of novel scaffolds were selected as potential leads to be used in novel and potent hCEase inhibitor design. The stability of binding modes and molecular interactions of these final hits were re-assured by molecular dynamics simulations.

Molecular dynamics simulations of a reversibly folding beta-heptapeptide in methanol: influence of the treatment of long-range electrostatic interactions

Eight 100-ns molecular dynamics simulations of a beta-heptapeptide in methanol at 340 K (within cubic periodic computational boxes of about 6-nm edge) are reported and compared. These simulations were performed with three different charge-state combinations at the peptide termini, one of them with or without a neutralizing chloride counterion, and using either the lattice-sum (LS) or reaction-field (RF) scheme to handle electrostatic interactions. The choice of the electrostatic scheme has essentially no influence on the folding-unfolding equilibrium when the peptide termini are uncharged and only a small influence when the peptide is positively charged at its N-terminus (with or without inclusion of a neutralizing chloride counterion). However, when the peptide is zwitterionic, the LS scheme leads to preferential sampling of the high-dipole folded helical state, whereas the RF scheme leads to preferential sampling of a low-dipole unfolded salt-bridged state. A continuum electrostatics analysis based on the sampled configurations (zwitterionic case) suggests that the LS scheme stabilizes the helical state through artificial periodicity, but that the magnitude of this perturbation is essentially negligible (compared to the thermal energy) for the large box size and relatively polar solvent considered. The results thus provide clear evidence (continuum electrostatics analysis) for the absence of LS artifacts and some indications (still not definitive because of the limited sampling of the folding-unfolding transition) for the presence of RF artifacts in this specific system.

Liquid water simulation: A critical examination of cutoff length

Cutoff treatment is the simplest approach for evaluating intermolecular interactions in molecular dynamics simulations. It has been believed that increasing cutoff length makes simulation results better. On the contrary, our results of the bulk water simulations studied within the range of cutoff lengths, 9-18 A, showed an opposite tendency: the artifact was enhanced by increasing the cutoff length. Especially, in terms of the distance dependent Kirkwood factor GK(r), it was clearly shown that the orientational behavior of water molecules becomes gradually worse as the cutoff length becomes longer. The artifact enhanced by the increased cutoff length led to a reported spurious artifact, i.e., phase transition [Y. Yonetani, Chem. Phys. Lett. 406, 49 (2005)]. Though the cutoff artifact was largely reduced by adopting a force switching technique, it did not completely remove the anomalous cutoff length dependence of the artifact. These results suggest that increasing the cutoff should not be attempted regardless of whether the switching force is adopted or not.

Development of a lattice-sum method emulating nonperiodic boundary conditions for the treatment of electrostatic interactions in molecular simulations: A continuum-electrostatics study

Artifacts induced by the application of periodic boundary conditions and lattice-sum methods in explicit-solvent simulations of (bio-)molecular systems are nowadays a major concern in the computer-simulation community. The present article reports a first step toward the design of a modified lattice-sum algorithm emulating nonperiodic boundary conditions, and therefore exempt of such periodicity-induced artifacts. This result is achieved here in the (more simple) context of continuum electrostatics. It is shown that an appropriate modification of the periodic Poisson equation and of its boundary conditions leads to a continuum-electrostatics scheme, which, although applied under periodic boundary conditions, exactly mimics the nonperiodic situation. The possible extension of this scheme to explicit-solvent simulations is outlined and its practical implementation will be described in more details in a forthcoming article.