Structural and energetic effects of truncating long ranged interactions in ionic and polar fluids

Analytic Gradients for the Electrostatic Embedding QM/MM Model in Periodic Boundary Conditions Using Particle-Mesh Ewald Sums and Electrostatic Potential Fitted Charge Operators

Long-range electrostatic effects are fundamental for describing chemical reactivity in the condensed phase. Here, we present the methodology of an efficient quantum mechanical/molecular mechanical (QM/MM) model in periodic boundary conditions (PBC) compatible with QM/MM boundaries at chemical bonds. The method combines electrostatic potential fitted charge operators and electrostatic potentials derived from the smooth particle-mesh Ewald (PME) sum approach. The total energy and its analytic first derivatives with respect to QM, MM, and lattice vectors allow QM/MM molecular dynamics (MD) in the most common thermodynamic ensembles. We demonstrate the robustness of the method by performing a QM/MM MD equilibration of methanol in water. We simulate the cis/trans isomerization free-energy profiles in water of proline amino acid and a proline-containing oligopeptide, showing a correct description of the reaction barrier. Our PBC-compatible QM/MM model can efficiently be used to study the chemical reactivity in the condensed phase and enzymatic catalysis.

Absolute local conformation of poly(methyl methacrylate) chains adsorbed on a quartz surface

Polymer chains at a buried interface with an inorganic solid play a critical role in the performance of polymer nanocomposites and adhesives. Sum frequency generation (SFG) vibrational spectroscopy with a sub-nanometer depth resolution provides valuable information regarding the orientation angle of functional groups at interfaces. However, in the case of conventional SFG, since the signal intensity is proportional to the square of the second-order nonlinear optical susceptibility and thereby loses phase information, it cannot be unambiguously determined whether the functional groups face upward or downward. This problem can be solved by phase-sensitive SFG (ps-SFG). We here applied ps-SFG to poly(methyl methacrylate) (PMMA) chains in direct contact with a quartz surface, shedding light on the local conformation of chains adsorbed onto the solid surface. The measurements made it possible to determine the absolute orientation of the ester methyl groups of PMMA, which were oriented toward the quartz interface. Combining ps-SFG with all-atomistic molecular dynamics simulation, the distribution of the local conformation and the driving force are also discussed.

Mechanisms of dihydropyridine agonists and antagonists in view of cryo-EM structures of calcium and sodium channels

Opposite effects of 1,4-dihydropyridine (DHP) agonists and antagonists on the L-type calcium channels are a challenging problem. Cryo-EM structures visualized DHPs between the pore-lining helices S6III and S6IV in agreement with published mutational data. However, the channel conformations in the presence of DHP agonists and antagonists are virtually the same, and the mechanisms of the ligands' action remain unclear. We docked the DHP agonist S-Bay k 8644 and antagonist R-Bay k 8644 in Cav1.1 channel models with or without π-bulges in helices S6III and S6IV. Cryo-EM structures of the DHP-bound Cav1.1 channel show a π-bulge in helix S6III but not in S6IV. The antagonist's hydrophobic group fits into the hydrophobic pocket formed by residues in S6IV. The agonists' polar NO2 group is too small to fill up the pocket. A water molecule could sterically fit into the void space, but its contacts with isoleucine in helix S6IV (motif INLF) would be unfavorable. In a model with π-bulged S6IV, this isoleucine turns away from the DHP molecule and its position is occupied by the asparagine from the same motif INLF. The asparagine provides favorable contacts for the water molecule at the agonist's NO2 group but unfavorable contacts for the antagonist's methoxy group. In our models, the DHP antagonist stabilizes entirely α-helical S6IV. In contrast, the DHP agonist stabilizes π-bulged helix S6IV whose C-terminal part turned and rearranged the activation-gate region. This would stabilize the open channel. Thus, agonists, but not antagonists, would promote channel opening by stabilizing π-bulged helix S6IV.

Molecular Modeling of Cardiac Sodium Channel with Mexiletine

A sodium channel blocker mexiletine (MEX) is used to treat chronic pain, myotonia and some arrhythmias. Mutations in the pore domain (PD) of voltage-gated sodium channels differently affect tonic block (TB) and use-dependent block (UDB) by MEX. Previous studies identified several MEX-sensing residues in the hNav1.5 channel and demonstrated that the channel block by MEX increases with activation of the voltage-sensing domain III (VSDIII), whereas MEX stabilizes the activated state of VSDIII. Structural rationales for these observations are unclear. Here, Monte Carlo (MC) energy minimizations were used to dock MEX and its more potent analog, Thio-Me2, into the hNav1.5 cryo-EM structure with activated VSDs and presumably inactivated PD. Computations yielded two ensembles of ligand binding poses in close contacts with known MEX-sensing residues in helices S6III, S6IV and P1IV. In both ensembles, the ligand NH3 group approached the cation-attractive site between backbone carbonyls at the outer-pore bottom, while the aromatic ring protruded ether into the inner pore (putative UDB pose) or into the III/IV fenestration (putative TB pose). In silico deactivation of VSDIII shifted helices S4-S5III, S5III, S6III and S6IV and tightened the TB site. In a model with activated VSDIII and three resting VSDs, MC-minimized energy profile of MEX pulled from the TB site towards lipids shows a deep local minimum due to interactions with 11 residues in S5III, P1III, S6III and S6IV. The minimum may correspond to an interim binding site for MEX in the hydrophobic path to the TB site along the lipid-exposed sides of repeats III and IV where 15 polar and aromatic residues would attract cationic blockers. The study explains numerous experimental data and suggests the mechanism of allosteric modification of the MEX binding site by VSDIII.

Non-Ewald methods for evaluating the electrostatic interactions of charge systems: similarity and difference

In molecular simulations, it is essential to properly calculate the electrostatic interactions of particles in the physical system of interest. Here we consider a method called the non-Ewald method, which does not rely on the standard Ewald method with periodic boundary conditions, but instead relies on the cutoff-based techniques. We focus on the physicochemical and mathematical conceptual aspects of the method in order to gain a deeper understanding of the simulation methodology. In particular, we take into account the reaction field (RF) method, the isotropic periodic sum (IPS) method, and the zero-multipole summation method (ZMM). These cutoff-based methods are based on different physical ideas and are completely distinguishable in their underlying concepts. The RF and IPS methods are "additive" methods that incorporate information outside the cutoff region, via dielectric medium and isotropic boundary condition, respectively. In contrast, the ZMM is a "subtraction" method that tries to remove the artificial effects, generated near the boundary, from the cutoff sphere. Nonetheless, we find physical and/or mathematical similarities between these methods. In particular, the modified RF method can be derived by the principle of neutralization utilized in the ZMM, and we also found a direct relationship between IPS and ZMM.

DFT-based QM/MM with Particle-Mesh Ewald for Direct, Long-Range Electrostatic Embedding

We present a DFT-based, QM/MM implementation with long-range electrostatic embedding achieved by direct real-space integration of the particle mesh Ewald (PME) computed electrostatic potential. The key transformation is the interpolation of the electrostatic potential from the PME grid to the DFT quadrature grid, from which integrals are easily evaluated utilizing standard DFT machinery. We provide benchmarks of the numerical accuracy with choice of grid size and real-space corrections, and demonstrate that good convergence is achieved while introducing nominal computational overhead. Furthermore, the approach requires only small modification to existing software packages, as is demonstrated with our implementation in the OpenMM and Psi4 software. After presenting convergence benchmarks, we evaluate the importance of long-range electrostatic embedding in three solute/solvent systems modeled with QM/MM. Water and BMIM/BF4 ionic liquid were considered as ``simple' and ``complex' solvents respectively, with water and p-phenylenediamine (PPD) solute molecules treated at QM level of theory. While electrostatic embedding with standard real-space truncation may introduce negligible error for simple systems such as water solute in water solvent, errors become more significant when QM/MM is applied to complex solvents such as ionic liquids. An extreme example is the electrostatic embedding energy for oxidized PPD in BMIM/BF4 for which real-space truncation produces severe error even at 2-3 nm cutoff distances. This latter example illustrates that utilization of QM/MM to compute redox potentials within concentrated electrolytes/ionic media requires carefully chosen long-range electrostatic embedding algorithms, with our presented algorithm providing a general and robust approach.

Intersegment Contacts of Potentially Damaging Variants of Cardiac Sodium Channel

Over 1,500 missense variants of sodium channel hNav1.5, which are reported in the ClinVar database, are associated with cardiac diseases. For most of the variants, the clinical significance is uncertain (VUS), not provided (NP), or has conflicting interpretations of pathogenicity (CIP). Reclassifying these variants as pathogenic/likely pathogenic (P/LP) variants is important for diagnosing genotyped patients. In our earlier work, several bioinformatics tools and paralogue annotation method consensually predicted that 74 VUS/NP/CIP variants of 54 wild type residues (set w54) are potentially damaging variants (PDVs). Atomic mechanisms underlying dysfunction of the PDVs are unknown. Here we employed a recent cryo-EM structure of the hNav1.5 channel with likely inactivated pore domain (PD) and activated voltage-sensing domains (VSDs), and ad hoc models of the closed and open PD and resting VSDs to explore intersegment contacts of w54 residues. We found that 44 residues from set w54 contact 84 residues with 118 disease missense variants. These include 104 VUS/NP/CIP variants, most of which are associated with the loss-of-function Brugada syndrome (BrS1) or gain-of-function long QT syndrome (LQT3). Matrix representation of the PDVs and their contact variants facilitated recognition of coupled mutations associated with the same disease. In particular, BrS1-associated coupled mutations, which disturb the P-loops region with the selectivity filter slow inactivation gate, would cause the channel dysfunction. Other likely causes of the channel dysfunction include coupled BrS1-associated variants within VSDs that would destabilize their activated states and coupled LQT3-associated variants, which would stabilize the open PD or activated VSDs. Our study proposes mechanisms of channel dysfunction for scores of BrS1- and LQT3-associated variants, confirms status for 82% of PDVs, and suggests damaging status for their contact variants, which are currently categorized as VUS/NP/CIP variants.

Possible Interactions of Extracellular Loop IVP2-S6 With Voltage-Sensing Domain III in Cardiac Sodium Channel

Motion transmission from voltage sensors to inactivation gates is an important problem in the general physiology of ion channels. In a cryo-EM structure of channel hNa_v1.5, residues N1736 and R1739 in the extracellular loop IVP2-S6 approach glutamates E1225 and E1295, respectively, in the voltage-sensing domain III (VSD-III). ClinVar-reported variants E1230K, E1295K, and R1739W/Q and other variants in loops IVP2-S6, IIIS1-S2, and IIIS3-S4 are associated with cardiac arrhythmias, highlighting the interface between IVP2-S6 and VSD-III as a hot spot of disease mutations. Atomic mechanisms of the channel dysfunction caused by these mutations are unknown. Here, we generated mutants E1295R, R1739E, E1295R/R1739E, and N1736R, expressed them in HEK-293T cells, and explored biophysical properties. Mutation E1295R reduced steady-state fast inactivation and enhanced steady-state slow inactivation. In contrast, mutation R1739E slightly enhanced fast inactivation and attenuated slow inactivation. Characteristics of the double mutant E1295R/R1739E were rather similar to those of the wild-type channel. Mutation N1736R attenuated slow inactivation. Molecular modeling predicted salt bridging of R1739E with the outermost lysine in the activated voltage-sensing helix IIIS4. In contrast, the loss-of-function substitution E1295R repelled R1739, thus destabilizing the activated VSD-III in agreement with our data that E1295R caused a depolarizing shift of the G-V curve. In silico deactivation of VSD-III with constraint-maintained salt bridge E1295-R1739 resulted in the following changes: 1) contacts between IIIS4 and IVS5 were switched; 2) contacts of the linker-helix IIIS4-S5 with IVS5, IVS6, and fast inactivation tripeptide IFM were modified; 3) contacts of the IFM tripeptide with helices IVS5 and IVS6 were altered; 4) mobile loop IVP2-S6 shifted helix IVP2 that contributes to the slow inactivation gate and helix IVS6 that contributes to the fast inactivation gate. The likelihood of salt bridge E1295-R1739 in deactivated VSD-III is supported by Poisson–Boltzmann calculations and state-dependent energetics of loop IVP2-S6. Taken together, our results suggest that loop IVP2-S6 is involved in motion transmission from VSD-III to the inactivation gates.

Improving the surface hydrophobicity by the solvent effect to reduce the water erosion of the CL-20/TNT cocrystal explosive

Explosives should be isolated from the air to reduce the surface erosion by water vapor in the storage and transport processes. The CL-20/TNT cocrystal was chosen as the research object on account of its structural arrangement and weak intermolecular interactions. A relatively extreme assumption that water molecules and the CL-20/TNT cocrystal existed as a solution system on the interface was proposed to investigate the influence mechanism. CL-20/TNT-water interfacial models were constructed, based on four important stable surfaces predicted using crystal morphology theory. The roughness and the electrostatic potential of each surface were analyzed to judge the strength of interactions between the water layer and the cocrystal surface. The effects of the water layer on the crystal surfaces were quantified in terms of binding energy and the radial distribution function using molecular dynamics simulations. Based on the analysis above, the (0 0 2) face was identified as the least affected by water erosion and its growth should be promoted. The crystal morphologies obtained by crystallization in different solvents are significantly different on account of the solvent effect, which can be used to search for a suitable solvent for crystallization to obtain a cocrystal with a larger hydrophobic surface.

Mapping the interaction surface of scorpion β-toxins with an insect sodium channel

The interaction of insect-selective scorpion depressant β-toxins (LqhIT2 and Lqh-dprIT3 from Leiurus quinquestriatus hebraeus) with the Blattella germanica sodium channel, BgNav1-1a, was investigated using site-directed mutagenesis, electrophysiological analyses, and structural modeling. Focusing on the pharmacologically defined binding site-4 of scorpion β-toxins at the voltage-sensing domain II (VSD-II), we found that charge neutralization of D802 in VSD-II greatly enhanced the channel sensitivity to Lqh-dprIT3. This was consistent with the high sensitivity of the splice variant BgNav2-1, bearing G802, to Lqh-dprIT3, and low sensitivity of BgNav2-1 mutant, G802D, to the toxin. Further mutational and electrophysiological analyses revealed that the sensitivity of the WT = D802E < D802G < D802A < D802K channel mutants to Lqh-dprIT3 correlated with the depolarizing shifts of activation in toxin-free channels. However, the sensitivity of single mutants involving IIS4 basic residues (K4E = WT << R1E < R2E < R3E) or double mutants (D802K = K4E/D802K = R3E/D802K > R2E/D802K > R1E/D802K > WT) did not correlate with the activation shifts. Using the cryo-EM structure of the Periplaneta americana channel, NavPaS, as a template and the crystal structure of LqhIT2, we constructed structural models of LqhIT2 and Lqh-dprIT3-c in complex with BgNav1-1a. These models along with the mutational analysis suggest that depressant toxins approach the salt-bridge between R1 and D802 at VSD-II to form contacts with linkers IIS1-S2, IIS3-S4, IIIP5-P1 and IIIP2-S6. Elimination of this salt-bridge enables deeper penetration of the toxin into a VSD-II gorge to form new contacts with the channel, leading to increased channel sensitivity to Lqh-dprIT3.

Theoretical study on the adsorption mechanism of PbCl2/CdCl2 by kaolinite during municipal solid waste pyrolysis

In the process of municipal solid waste (MSW) pyrolysis, kaolinite possesses an outstanding trapping effect on semi-volatile metal vapors (Pb, Cd) through physical and chemical adsorption. In this paper, the microscopic mechanism of PbCl₂ and CdCl₂ adsorption on the surface of Al rings and Si rings of kaolinite was investigated by combining Monte Carlo method with density functional theory (DFT). The calculations indicate that the continuously enriched pore structure in the process of dehydroxylation indirectly influences the adsorption of PbCl₂/CdCl₂ by kaolinite. Under the non-bond interaction and electron transfer induction, PbCl₂ molecules are more conveniently adsorbed on the Al-(001) surface than CdCl₂, while the adsorption sites of CdCl₂ molecules are more widely distributed on the Si-(001) surface. Moreover, the transform in the Al-coordination and the exposed active oxygen atoms significantly affect the adsorption activity of kaolinite (the capability to gain and lose electrons). Considering the energy barrier and electrophilic nucleophilic sensitivity, it is more feasible for PbCl₂/CdCl₂ to be adsorbed near IV/V-coordinated Al and active O under Van der Waals action. Subsequently, IV/V-coordinated Al will act as an electron acceptor, and the active oxygen atoms after dehydrogenation will serve as an electron donor. Under the induction of the energy difference of frontier orbitals, the electrons transfer will encourage the formation of more stable adsorption states. The results shed new light on strengthening the adsorption activity of kaolinite to PbCl₂/CdCl₂ in the process of MSW pyrolysis.

Solvation Free Energies and Adsorption Energies at the Metal/Water Interface from Hybrid Quantum-Mechanical/Molecular Mechanics Simulations

Modeling adsorption at metal/water interfaces is a cornerstone toward an improved understanding in a variety of fields from heterogeneous catalysis to corrosion. We propose and validate a hybrid scheme that combines the adsorption free energies obtained in the gas phase at the density functional theory level with the variation in solvation from the bulk phase to the interface evaluated using a MM-based alchemical transformation, denoted MMsolv. Using the GAL17 force field for the platinum/water interaction, we retrieve a qualitatively correct interaction energy of the water solvent at the interface. This interaction is of near chemisorption character and thus challenging, both for the alchemical transformation and also for the fixed point-charge electrostatics. Our scheme passes through a state characterized by a well-behaved physisorption potential for the Pt(111)/H₂O interaction to converge the free energy difference. The workflow is implemented in the freely available SolvHybrid package. We first assess the adsorption of a water molecule at the Pt/water interface, which turns out to be a stringent test. The intrinsic error of our quantum-mechanical/molecular mechanics (QM/MM) hybrid scheme is limited to 6 kcal mol^-1 through the introduction of a correction term to attenuate the electrostatic interaction between near-chemisorbed water molecules and the underlying Pt atoms. Next, we show that the MMsolv solvation free energy of Pt (-0.46 J m^-2) is in good agreement with the experimental estimate (-0.32 J m^-2). Furthermore, we show that the entropy contribution at room temperature is roughly of equal magnitude as the free energy but with an opposite sign. Finally, we compute the adsorption energy of benzene and phenol at the Pt(111)/water interface, one of the rare systems for which experimental data are available. In qualitative agreement with the experiment, but in stark contrast with a standard implicit solvent model, the adsorption of these aromatic molecules is strongly reduced (i.e., less exothermic by ∼30 and 40 kcal mol^-1 for our QM/MM hybrid scheme and experiment, respectively, but ∼0 with the implicit solvent) at the solid/liquid interface compared to the solid/gas interface. This reduction occurs mainly because of the competition between the organic adsorbate and the solvent for adsorption on the metallic surface. The semiquantitative agreement with experimental estimates for the adsorption energy of aromatic molecules thus validates the soundness of our hybrid QM/MM scheme.

The general anesthetic etomidate and fenamate mefenamic acid oppositely affect GABAAR and GlyR: a structural explanation

GABA and glycine act as inhibitory neurotransmitters in the CNS. Inhibitory neurotransmission is mediated via activation of ionotropic GABA_A and glycine receptors. We used a modeling approach to explain the opposite effects of the general anesthetic etomidate (ETM) and fenamate mefenamic acid (MFA) on GABA- and glycine-activated currents recorded in isolated cerebellar Purkinje cells and hippocampal pyramidal neurons, respectively. These drugs potentiated GABA_ARs but blocked GlyRs. We built a homology model of α1β GlyR based on the cryo-EM structure of open α1 GlyR, used the α1β3γ2 GABA_AR structure from the PDB, and applied Monte-Carlo energy minimization to optimize models of receptors and ligand-receptor complexes. In silico docking suggests that ETM/MFA bind at the transmembrane β( +)/α( -) intersubunit interface in GABA_AR. Our models predict that the bulky side chain of the highly conserved Arg^19' residue at the plus interface side wedges the interface and maintains the conducting receptor state. We hypothesized that MFA/ETM binding at the β( +)/α( -) interface leads to prolongation of receptor life-time in the open state. Having analyzed different GABA_AR and GlyR structures available in the PDB, we found that mutual arrangement of the Arg^19' and Gln^-26' side chains at the plus and minus interface sides, respectively, plays an important role when the receptor switches from the open to closed state. We show that this process is accompanied by narrowing of the intersubunit interfaces, leading to extrusion of the Arg^19' side chain from the interface. Our models allow us to explain the lack of GlyR potentiation in our electrophysiological experiments.

Phenylalkylamines in calcium channels: computational analysis of experimental structures

Experimental 3D structures of calcium channels with phenylalkylamines (PAAs) provide basis for further analysis of atomic mechanisms of these important cardiovascular drugs. In the crystal structure of the engineered calcium channel CavAb with Br-verapamil and in the cryo-EM structure of the Cav1.1 channel with verapamil, the ligands bind in the inner pore. However, there are significant differences between these structures. In the crystal structure the ligand ammonium group is much closer to the ion in the selectivity-filter region Site 3, which is most proximal to the inner pore, than in the cryo-EM structure. Here we used Monte Carlo energy minimizations to dock PAAs in calcium channels. Our computations suggest that in the crystal structure Site 3 is occupied by a water molecule rather than by a calcium ion. Analysis of the published electron density map does not rule out this possibility. In the cryo-EM structures the ammonium group of verapamil is shifted from the calcium ion in Site 3 either along the pore axis, towards the cytoplasm or away from the axis. Our unbiased docking reproduced these binding modes. However, in the cryo-EM structures detergent and lipid molecules interact with verapamil. When we removed these molecules, the nitrile group of verapamil bound to the calcium ion in Site 3. Models of Cav1.2 with different PAAs suggest similar binding modes and direct contacts of the ligands electronegative atoms with the calcium ion in Site 3. Such interactions explain paradoxes in structure-activity relationships of PAAs.

Comparative investigation on the thermostability, sensitivity, and mechanical performance of RDX/HMX energetic cocrystal and its mixture

Molecular mechanics (MM) and molecular dynamics (MD) simulation method were applied to explore the impact of temperature (220-380 K) on the thermostability, sensitivity, and mechanical performance of RDX (1,3,5-trinitro-1,3,5-triazacyco-hexane)/HMX (1,3,5,7-tetranitro-1,3,5,7-tetrazocane) energetic cocrystal and mixture models. The mechanical property, the maximum trigger bond length ([Formula: see text]), binding energy, and cohesive energy density (CED) of the pure RDX, β-HMX crystal, the cocrystal, and mixture models were acquired and compared. The results manifest that temperature has an important impact on the binding capacity between the components of the cocrystal and mixture. The binding energies decrease as the temperature rises, and the cocrystal has larger values than those of mixture. For all the models, the [Formula: see text] increases and the CEDs decrease with the rising temperature, implying that the sensitivity of the explosives increases, while the [Formula: see text] values of the cocrystal are smaller than those of HMX and the CED values are between those of RDX and β-HMX, indicating that the sensitivity has been enhanced through co-crystallization. As the temperature increases, the shear modulus (G), bulk modulus (K), and tensile modulus (E) values of all models have an evident downtrend. Simultaneously, G, K, and E values of the cocrystal model are less than those of RDX and β-HMX, while the K/G ratio and Cauchy pressure (C₁₂-C₄₄) are larger, signifying that co-crystallization can weaken the brittleness and enhance the ductility of the pure crystals. Compared with the mixture, the cocrystal has better ductility and stability.

Distributed and multicore QuBiLS-MIDAS software v2.0: Computing chiral, fuzzy, weighted and truncated geometrical molecular descriptors based on tensor algebra

Advances to the distributed, multi-core and fully cross-platform QuBiLS-MIDAS software v2.0 (http://tomocomd.com/qubils-midas) are reported in this article since the v1.0 release. The QuBiLS-MIDAS software is the only one that computes atom-pair and alignment-free geometrical MDs (3D-MDs) from several distance metrics other than the Euclidean distance, as well as alignment-free 3D-MDs that codify structural information regarding the relations among three and four atoms of a molecule. The most recent features added to the QuBiLS-MIDAS software v2.0 are related (a) to the calculation of atomic weightings from indices based on the vertex-degree invariant (e.g., Alikhanidi index); (b) to consider central chirality during the molecular encoding; (c) to use measures based on clustering methods and statistical functions to codify structural information among more than two atoms; (d) to the use of a novel method based on fuzzy membership functions to spherically truncate inter-atomic relations; and (e) to the use of weighted and fuzzy aggregation operators to compute global 3D-MDs according to the importance and/or interrelation of the atoms of a molecule during the molecular encoding. Moreover, a novel module to compute QuBiLS-MIDAS 3D-MDs from their headings was also developed. This module can be used either by the graphical user interface or by means of the software library. By using the library, both the predictive models built with the QuBiLS-MIDAS 3D-MDs and the QuBiLS-MIDAS 3D-MDs calculation can be embedded in other tools. A set of predefined QuBiLS-MIDAS 3D-MDs with high information content and low redundancy on a set comprised of 20,469 compounds is also provided to be employed in further cheminformatics tasks. This set of predefined 3D-MDs evidenced better performance than all the universe of Dragon (v5.5) and PaDEL 0D-to-3D MDs in variability studies, whereas a linear independence study proved that these QuBiLS-MIDAS 3D-MDs codify chemical information orthogonal to the Dragon 0D-to-3D MDs. This set of predefined 3D-MDs would be periodically updated as long as new results be achieved. In general, this report highlights our continued efforts to provide a better tool for a most suitable characterization of compounds, and in this way, to contribute to obtaining better outcomes in future applications.

The mechanisms of potentiation and inhibition of GABAA receptors by non-steroidal anti-inflammatory drugs, mefenamic and niflumic acids

Fenamates mefanamic and niflumic acids (MFA and NFA) induced dual potentiating and inhibitory effects on GABA currents recorded in isolated cerebellar Purkinje cells using the whole-cell patch-clamp and fast-application techniques. Regardless of the concentration, both drugs induced a pronounced prolongation of the current response. We demonstrated that the same concentration of drugs can produce both potentiating and inhibitory effects, depending on the GABA concentration, which indicates that both processes take place simultaneously and the net effect depends on the concentrations of both the agonist and fenamate. We found that the NFA-induced block is strongly voltage-dependent. The Woodhull analysis of the block suggests that NFA has two binding sites in the pore - shallow and deep. We built a homology model of the open GABA_AR based on the cryo-EM structure of the open α1 GlyR and applied Monte-Carlo energy minimization to optimize the ligand-receptor complexes. A systematic search for MFA/NFA binding sites in the GABA_AR pore revealed the existence of two sites, the location of which coincides well with predictions of the Woodhull model. In silico docking suggests that two fenamate molecules are necessary to occlude the pore. We showed that MFA, acting as a PAM, competes with an intravenous anesthetic etomidate for a common binding site. We built structural models of MFA and NFA binding at the transmembrane β(+)/α(-) intersubunit interface. We suggested a hypothesis on the molecular mechanism underlying the prolongation of the receptor lifetime in open state after MFA/NFA binding and β subunit specificity of the fenamate potentiation.

Smoothed Spherical Truncation based on Fuzzy Membership Functions: Application to the Molecular Encoding

A novel spherical truncation method, based on fuzzy membership functions, is introduced to truncate interatomic (or interaminoacid) relations according to smoothing values computed from fuzzy membership degrees. In this method, the molecules are circumscribed into a sphere, so that the geometric centers of the molecules are the centers of the spheres. The fuzzy membership degree of each atom (or aminoacid) is computed from its distance with respect to the geometric center of the molecule, by using a fuzzy membership function. So, the smoothing value to be applied in the truncation of a relation (or interaction) is computed by averaging the fuzzy membership degrees of the atoms (or aminoacids) involved in the relation. This truncation method is rather different from the existing ones, at considering the geometric center for the whole molecule and not only for atom-groups, as well as for using fuzzy membership functions to compute the smoothing values. A variability study on a set comprised of 20,469 compounds (15,050 drug-like compounds, 2994 drugs approved, 880 natural products from African sources, and 1545 plant-derived natural compounds exhibiting anti-cancerous activity) demonstrated that the truncation method proposed allows to determine molecular encodings with better ability for discriminating among structurally different molecules than the encodings obtained without applying truncation or applying non-fuzzy truncation functions. Moreover, a principal component analysis revealed that orthogonal chemical information of the molecules is encoded by using the method proposed. Lastly, a modeling study proved that the truncation method improves the modeling ability of existing geometric molecular descriptors, at allowing to develop more robust models than the ones built only using non-truncated descriptors. In this sense, a comparison and statistical assessment were performed on eight chemical datasets. As a result, the models based on the truncated molecular encodings yielded statistically better results than 12 procedures considered from the literature. It can thus be stated that the proposed truncation method is a relevant strategy for obtaining better molecular encodings, which will be ultimately useful in enhancing the modeling ability of existing encodings both on small-to-medium size molecules and biomacromolecules. © 2019 Wiley Periodicals, Inc.

Development of 1,3-thiazole analogues of imidazopyridines as potent positive allosteric modulators of GABAA receptors

Structure-activity relationship studies were conducted in the search for 1,3-thiazole isosteric analogs of imidazopyridine drugs (Zolpidem, Alpidem). Three series of novel γ-aminobutyric acid receptor (GABA_AR) ligands belonging to imidazo[2,1-b]thiazoles, imidazo[2,1-b][1,3,4]thiadiazoles, and benzo[d]imidazo[2,1-b]thiazoles were synthesized and characterized as active agents against GABA_AR benzodiazepine-binding site. In each of these series, potent compounds were discovered using a radioligand competition binding assay. The functional properties of highest-affinity compounds 28 and 37 as GABA_AR positive allosteric modulators (PAMs) were determined by electrophysiological measurements. In vivo studies on zebrafish demonstrated their potential for the further development of anxiolytics. Using the OECD "Fish, Acute Toxicity Test" active compounds were found safe and non-toxic. Structural bases for activity of benzo[d]imidazo[2,1-b]thiazoles were proposed using molecular docking studies. The isosteric replacement of the pyridine nuclei by 1,3-thiazole, 1,3,4-thiadiazole, or 1,3-benzothiazole in the ring-fused imidazole class of GABA_AR PAMs was shown to be promising for the development of novel hypnotics, anxiolytics, anticonvulsants, and sedatives drug-candidates.

Atomic Mechanisms of Timothy Syndrome-Associated Mutations in Calcium Channel Cav1.2

Timothy syndrome (TS) is a very rare multisystem disorder almost exclusively associated with mutations G402S and G406R in helix IS6 of Cav1.2. Recently, mutations R518C/H in helix IIS0 of the voltage sensing domain II (VSD-II) were described as a cause of cardiac-only TS. The three mutations are known to decelerate voltage-dependent inactivation (VDI). Here, we report a case of cardiac-only TS caused by mutation R518C. To explore possible impact of the three mutations on interdomain contacts, we modeled channel Cav1.2 using as templates Class Ia and Class II cryo-EM structures of presumably inactivated channel Cav1.1. In both models, R518 and several other residues in VSD-II donated H-bonds to the IS6-linked α1-interaction domain (AID). We further employed steered Monte Carlo energy minimizations to move helices S4–S5, S5, and S6 from the inactivated-state positions to those seen in the X-ray structures of the open and closed NavAb channel. In the open-state models, positions of AID and VSD-II were similar to those in Cav1.1. In the closed-state models, AID moved along the β subunit (Cavβ) toward the pore axis and shifted AID-bound VSD-II. In all the models R518 retained strong contacts with AID. Our calculations suggest that conformational changes in VSD-II upon its deactivation would shift AID along Cavβ toward the pore axis. The AID-linked IS6 would bend at flexible G402 and G406, facilitating the activation gate closure. Mutations R518C/H weakened the IIS0-AID contacts and would retard the AID shift. Mutations G406R and G402S stabilized the open state and would resist the pore closure. Several Cav1.2 mutations associated with long QT syndromes are consistent with this proposition. Our results provide a mechanistic rationale for the VDI deceleration caused by TS-associated mutations and suggest targets for further studies of calcium channelopathies.

A molecular modeling study for miscibility of polyimide/polythene mixing systems with/without compatibilizer

Molecular models were established to predict the miscibility of polyimide/polythene mixing systems and the enhancing effects of compatibilizer addition of maleic anhydride grafted polythene (MAH-g-PE). Molecular dynamics simulations were applied to investigate radial distribution functions and Flory-Huggins parameters of the mixing systems. Results show that polyimide/polythene is miscible to a certain degree, and the miscibility gets better after adding MAH-g-PE. Dissipative particle dynamics (DPD) simulations display that micro-phase separation occurs in the polyimide/polythene mixing systems, however, effective interfaces appear between polyimide and polythene phases after adding MAH-g-PE. The results of molecular mechanics simulations indicate that the ability of mixing systems to resist stretch, compression and shear deformation increases after adding MAH-g-PE. This work offers a promising technique to predict miscibility properties for polyimide/polythene system prior to actual production and attempt to find a suitable compatibilizer for that system.

Sensitivity of electrospray molecular dynamics simulations to long-range Coulomb interaction models

Molecular dynamics (MD) electrospray simulations of 1-ethyl-3-methylimidazolium tetrafluoroborate (EMIM-BF_{4}) ion liquid were performed with the goal of evaluating the influence of long-range Coulomb models on ion emission characteristics. The direct Coulomb (DC), shifted force Coulomb sum (SFCS), and particle-particle particle-mesh (PPPM) long-range Coulomb models were considered in this work. The DC method with a sufficiently large cutoff radius was found to be the most accurate approach for modeling electrosprays, but, it is computationally expensive. The Coulomb potential energy modeled by the DC method in combination with the radial electric fields were found to be necessary to generate the Taylor cone. The differences observed between the SFCS and the DC in terms of predicting the total ion emission suggest that the former should not be used in MD electrospray simulations. Furthermore, the common assumption of domain periodicity was observed to be detrimental to the accuracy of the capillary-based electrospray simulations.

Mechanism of sodium channel block by local anesthetics, antiarrhythmics, and anticonvulsants

Local anesthetics, antiarrhythmics, and anticonvulsants include both charged and electroneutral compounds that block voltage-gated sodium channels. Prior studies have revealed a common drug-binding region within the pore, but details about the binding sites and mechanism of block remain unclear. Here, we use the x-ray structure of a prokaryotic sodium channel, NavMs, to model a eukaryotic channel and dock representative ligands. These include lidocaine, QX-314, cocaine, quinidine, lamotrigine, carbamazepine (CMZ), phenytoin, lacosamide, sipatrigine, and bisphenol A. Preliminary calculations demonstrated that a sodium ion near the selectivity filter attracts electroneutral CMZ but repels cationic lidocaine. Therefore, we further docked electroneutral and cationic drugs with and without a sodium ion, respectively. In our models, all the drugs interact with a phenylalanine in helix IVS6. Electroneutral drugs trap a sodium ion in the proximity of the selectivity filter, and this same site attracts the charged group of cationic ligands. At this position, even small drugs can block the permeation pathway by an electrostatic or steric mechanism. Our study proposes a common pharmacophore for these diverse drugs. It includes a cationic moiety and an aromatic moiety, which are usually linked by four bonds.

Graphene Nanoribbons on Highly Porous 3D Graphene for High-Capacity and Ultrastable Al-Ion Batteries

Graphene nanoribbons on highly porous 3D graphene foam as the binder-free cathode for flexible Al-ion batteries exhibit low charge voltage, high capacity, excellent cycling ability (even after 10 000 cycles there is no capacity decay), and fast charging and slow discharging performance (the battery can be fully charged in 80 s and discharged in more than 3100 s).

Side chain flexibility and the pore dimensions in the GABAA receptor

Permeation of ions through open channels and their accessibility to pore-targeting drugs depend on the pore cross-sectional dimensions, which are known only for static X-ray and cryo-EM structures. Here, we have built homology models of the closed, open and desensitized α1β2γ2 GABAA receptor (GABAAR). The models are based, respectively, on the X-ray structure of α3 glycine receptor (α3 GlyR), cryo-EM structure of α1 GlyR and X-ray structure of β3 GABAAR. We employed Monte Carlo energy minimizations to explore how the pore lumen may increase due to repulsions of flexible side chains from a variable-diameter electroneutral atom (an expanding sphere) pulled through the pore. The expanding sphere computations predicted that the pore diameter averaged along the permeation pathway is larger by approximately 3 Å than that computed for the models with fixed sidechains. Our models predict three major pore constrictions located at the levels of -2', 9' and 20' residues. Residues around the -2' and 9' rings are known to form the desensitization and activation gates of GABAAR. Our computations predict that the 20' ring may also serve as GABAAR gate whose physiological role is unclear. The side chain flexibility of residues -2', 9' and 20' and hence the dimensions of the constrictions depend on the GABAAR functional state.

Long-range electrostatic corrections in multipolar/polarizable QM/MM simulations

Taking long-range electrostatic effects into account in classical and hybrid quantum mechanics-molecular mechanics (QM/MM) simulations is necessary for an accurate description of the system under study. We have recently developed a method, termed long-range electrostatic corrections (LREC), for monopolar QM/MM calculations. Here, we present an extension of LREC for multipolar/polarizable QM/MM simulations within the LICHEM software package. Reaction barriers and QM-MM interaction energies converge with a LREC cutoff between 20 and 25 Å, in agreement with our previous results. Additionally, the LREC approach for the QM-MM interactions can be smoothly combined with standard shifting or Ewald summation methods in the MM calculations. We recommend the use of QM(LREC)/MM(PME), where the QM region is treated with LREC and the MM region is treated with particle mesh Ewald (PME) summation. This combination is an excellent compromise between simplicity, speed, and accuracy for large QM/MM simulations.

Direct simulation of electron transfer in the cobalt hexammine(ii/iii) self-exchange reaction

We present an atomistic simulation of electron transfer in a transition metal complex system using path integral methods.

Efficient Multistate Reactive Molecular Dynamics Approach Based on Short-Range Effective Potentials

Nonbonded interactions between molecules usually include the van der Waals force and computationally expensive long-range electrostatic interactions. This article develops a more efficient approach: the effective-interaction multistate empirical-valence-bond (EI-MS-EVB) model. The EI-MS-EVB method relies on a mapping of all interactions onto a short-range and thus, computationally efficient effective potential. The effective potential is tabulated by matching its force to known trajectories obtained from the full-potential empirical multistate empirical-valence-bond (MS-EVB) model. The effective pairwise interaction depends on and is uniquely determined by the atomic configuration of the system, varying only with respect to the hydrogen-bonding topology. By comparing the EI-MS-EVB and full MS-EVB calculations of several equilibrium and dynamic properties important to hydrated excess proton solvation and transport, we show that the EI-MS-EVB model produces very accurate results for the specific system in which the tabulated potentials were generated. The EI-MS-EVB potential also transfers reasonably well to similar systems with different temperatures and box sizes. The EI-MS-EVB method also reduces the computational cost of the nonbonded interactions by about 1 order of magnitude in comparison with the full algorithm.

A new smoothing function to introduce long-range electrostatic effects in QM/MM calculations

A new method to account for long range electrostatic contributions is proposed and implemented for quantum mechanics/molecular mechanics long range electrostatic correction (QM/MM-LREC) calculations. This method involves the use of the minimum image convention under periodic boundary conditions and a new smoothing function for energies and forces at the cutoff boundary for the Coulomb interactions. Compared to conventional QM/MM calculations without long-range electrostatic corrections, the new method effectively includes effects on the MM environment in the primary image from its replicas in the neighborhood. QM/MM-LREC offers three useful features including the avoidance of calculations in reciprocal space (k-space), with the concomitant avoidance of having to reproduce (analytically or approximately) the QM charge density in k-space, and the straightforward availability of analytical Hessians. The new method is tested and compared with results from smooth particle mesh Ewald (PME) for three systems including a box of neat water, a double proton transfer reaction, and the geometry optimization of the critical point structures for the rate limiting step of the DNA dealkylase AlkB. As with other smoothing or shifting functions, relatively large cutoffs are necessary to achieve comparable accuracy with PME. For the double-proton transfer reaction, the use of a 22 Å cutoff shows a close reaction energy profile and geometries of stationary structures with QM/MM-LREC compared to conventional QM/MM with no truncation. Geometry optimization of stationary structures for the hydrogen abstraction step by AlkB shows some differences between QM/MM-LREC and the conventional QM/MM. These differences underscore the necessity of the inclusion of the long-range electrostatic contribution.

Assessing the accuracy of the isotropic periodic sum method through Madelung energy computation

We tested the isotropic periodic sum (IPS) method for computing Madelung energies of ionic crystals. The performance of the method, both in its nonpolar (IPSn) and polar (IPSp) forms, was compared with that of the zero-charge and Wolf potentials [D. Wolf, P. Keblinski, S. R. Phillpot, and J. Eggebrecht, J. Chem. Phys. 110, 8254 (1999)]. The results show that the IPSn and IPSp methods converge the Madelung energy to its reference value with an average deviation of ∼10(-4) and ∼10(-7) energy units, respectively, for a cutoff range of 18-24a (a/2 being the nearest-neighbor ion separation). However, minor oscillations were detected for the IPS methods when deviations of the computed Madelung energies were plotted on a logarithmic scale as a function of the cutoff distance. To remove such oscillations, we introduced a modified IPSn potential in which both the local-region and long-range electrostatic terms are damped, in analogy to the Wolf potential. With the damped-IPSn potential, a smoother convergence was achieved. In addition, we observed a better agreement between the damped-IPSn and IPSp methods, which suggests that damping the IPSn potential is in effect similar to adding a screening potential in IPSp.

Fluid infiltration pressure for hydrophobic nanochannels

We investigate water infiltration pressure for hydrophobic nanochannels through molecular dynamics simulations. It is found that the entrance energy barrier significantly raises the infiltration pressure, which makes the classic Young-Laplace equation invalid for nanochannels. As the channel surface is tuned from superhydrophobic to hydrophobic, the infiltration pressure is greatly reduced mainly due to the decrease of the capillary pressure (Young-Laplace equation) caused by the contact angle change, while the contribution of the entrance energy barrier to the infiltration pressure, which is termed entrance barrier pressure, increases from 25% to 60%.

Toward the correction of effective electrostatic forces in explicit-solvent molecular dynamics simulations: restraints on solvent-generated electrostatic potential and solvent polarization

Despite considerable advances in computing power, atomistic simulations under nonperiodic boundary conditions, with Coulombic electrostatic interactions and in systems large enough to reduce finite-size associated errors in thermodynamic quantities to within the thermal energy, are still not affordable. As a result, periodic boundary conditions, systems of microscopic size and effective electrostatic interaction functions are frequently resorted to. Ensuing artifacts in thermodynamic quantities are nowadays routinely corrected a posteriori, but the underlying configurational sampling still descends from spurious forces. The present study addresses this problem through the introduction of on-the-fly corrections to the physical forces during an atomistic molecular dynamics simulation. Two different approaches are suggested, where the force corrections are derived from special potential energy terms. In the first approach, the solvent-generated electrostatic potential sampled at a given atom site is restrained to a target value involving corrections for electrostatic artifacts. In the second approach, the long-range regime of the solvent polarization around a given atom site is restrained to the Born polarization, i.e., the solvent polarization corresponding to the ideal situation of a macroscopic system under nonperiodic boundary conditions and governed by Coulombic electrostatic interactions. The restraints are applied to the explicit-water simulation of a hydrated sodium ion, and the effect of the restraints on the structural and energetic properties of the solvent is illustrated. Furthermore, by means of the calculation of the charging free energy of a hydrated sodium ion, it is shown how the electrostatic potential restraint translates into the on-the-fly consideration of the corresponding free-energy correction terms. It is discussed how the restraints can be generalized to situations involving several solute particles. Although the present study considers a very simple system only, it is an important step toward the on-the-fly elimination of finite-size and approximate-electrostatic artifacts during atomistic molecular dynamics simulations.

Simulating chemical reactions in ionic liquids using QM/MM methodology

The use of ionic liquids as a reaction medium for chemical reactions has dramatically increased in recent years due in large part to the numerous reported advances in catalysis and organic synthesis. In some extreme cases, ionic liquids have been shown to induce mechanistic changes relative to conventional solvents. Despite the large interest in the solvents, a clear understanding of the molecular factors behind their chemical impact is largely unknown. This feature article reviews our efforts developing and applying mixed quantum and molecular mechanical (QM/MM) methodology to elucidate the microscopic details of how these solvents operate to enhance rates and alter mechanisms for industrially and academically important reactions, e.g., Diels-Alder, Kemp eliminations, nucleophilic aromatic substitutions, and β-eliminations. Explicit solvent representation provided the medium dependence of the activation barriers and atomic-level characterization of the solute-solvent interactions responsible for the experimentally observed "ionic liquid effects". Technical advances are also discussed, including a linear-scaling pairwise electrostatic interaction alternative to Ewald sums, an efficient polynomial fitting method for modeling proton transfers, and the development of a custom ionic liquid OPLS-AA force field.

Block of GABA(A) receptor ion channel by penicillin: electrophysiological and modeling insights toward the mechanism

GABA(A) receptors (GABA(A)R) mainly mediate fast inhibitory neurotransmission in the central nervous system. Different classes of modulators target GABA(A)R properties. Penicillin G (PNG) belongs to the class of noncompetitive antagonists blocking the open GABA(A)R and is a prototype of β-lactam antibiotics. In this study, we combined electrophysiological and modeling approaches to investigate the peculiarities of PNG blockade of GABA-activated currents recorded from isolated rat Purkinje cells and to predict the PNG binding site. Whole-cell patch-сlamp recording and fast application system was used in the electrophysiological experiments. PNG block developed after channel activation and increased with membrane depolarization suggesting that the ligand binds within the open channel pore. PNG blocked stationary component of GABA-activated currents in a concentration-dependent manner with IC50 value of 1.12mM at -70mV. The termination of GABA and PNG co-application was followed by a transient tail current. Protection of the tail current from bicuculline block and dependence of its kinetic parameters on agonist affinity suggest that PNG acts as a sequential open channel blocker that prevents agonist dissociation while the channel remains blocked. We built the GABA(A)R models based on nAChR and GLIC structures and performed an unbiased systematic search of the PNG binding site. Monte-Carlo energy minimization was used to find the lowest energy binding modes. We have shown that PNG binds close to the intracellular vestibule. In both models the maximum contribution to the energy of ligand-receptor interactions revealed residues located on the level of 2', 6' and 9' rings formed by a bundle of M2 transmembrane segments, indicating that these residues most likely participate in PNG binding. The predicted structural models support the described mechanism of PNG block.

Quantum and Molecular Mechanical (QM/MM) Monte Carlo Techniques for Modeling Condensed-Phase Reactions

A recent review (Acc. Chem. Res. 2010, 43:142-151) examined our use and development of a combined quantum and molecular mechanical (QM/MM) technique for modelling organic and enzymatic reactions. Advances included the PDDG/PM3 semiempirical QM (SQM) method, computation of multi-dimensional potentials of mean force (PMF), incorporation of on-the-fly QM in Monte Carlo simulations, and a polynomial quadrature method for rapidly treating proton-transfer reactions. The current article serves as a follow up on our progress. Highlights include new reactions, alternative SQM methods, a polarizable OPLS force field, and novel solvent environments, e.g., "on water" and room temperature ionic liquids. The methodology is strikingly accurate across a wide range of condensed-phase and antibody-catalyzed reactions including substitution, decarboxylation, elimination, isomerization, and pericyclic classes. Comparisons are made to systems treated with continuum-based solvents and ab initio or density functional theory (DFT) methods. Overall, the QM/MM methodology provides detailed characterization of reaction paths, proper configurational sampling, several advantages over implicit solvent models, and a reasonable computational cost.

Isotropic Periodic Sum Treatment of Long-Range Electrostatic Interactions in Combined Quantum Mechanical and Molecular Mechanical Calculations

The isotropic periodic sum (IPS) method was extended to describe long-range electrostatic interactions in combined quantum mechanical and molecular mechanical (QM/MM) calculations. The resulting method, designated QM/MM-IPS, was tested for two ion association processes and a model SN2 reaction in aqueous solution. Potential of mean force (PMF) profiles and radial distribution functions computed from the QM/MM-IPS simulations were compared with those obtained by using the existing QM/MM-Ewald sum and cutoff (QM/MM-Cutoff) methods. In contrast to the QM/MM-Cutoff method, with which PMFs of ion separation tend to display a spurious linear drift, the QM/MM-IPS method successfully eliminates such artifacts, in excellent agreement with the QM/MM-Ewald results. The PMF obtained with the QM/MM-IPS method for the SN2 reaction that transfers an NH3 group between two chloride anions closely resembles that from the QM/MM-Ewald simulations. Compared with QM/MM-Ewald, the QM/MM-IPS method reduces the computational cost by 60-70% when a local region of 12 to 14 Å is used. These results suggest that the QM/MM-IPS method can be used as a reliable and efficient alternative to the QM/MM-Ewald method to incorporate long-range electrostatic effects in simulating solution-phase chemical reactions.

Molecular dynamics simulations of double-stranded DNA in an explicit solvent model with the zero-dipole summation method

Molecular dynamics (MD) simulations of a double-stranded DNA with explicit water and small ions were performed with the zero-dipole summation (ZD) method, which was recently developed as one of the non-Ewald methods. Double-stranded DNA is highly charged and polar, with phosphate groups in its backbone and their counterions, and thus precise treatment for the long-range electrostatic interactions is always required to maintain the stable and native double-stranded form. A simple truncation method deforms it profoundly. On the contrary, the ZD method, which considers the neutralities of charges and dipoles in a truncated subset, well reproduced the electrostatic energies of the DNA system calculated by the Ewald method. The MD simulations using the ZD method provided a stable DNA system, with similar structures and dynamic properties to those produced by the conventional Particle mesh Ewald method.