Journal of Computational Chemistry

Lee, Hwayoung; Im, Wonpil

doi: 10.1002/jcc.70108pmid: 40251888

Diacylglycerol O‐acyltransferase 1 (DGAT1) is an integral membrane protein that uses acyl‐coenzyme A (acyl‐CoA) and diacylglycerol (DAG) to catalyze the formation of triacylglycerides (TAGs). The acyl transfer reaction occurs between the activated carboxylate group of the fatty acid and the free hydroxyl group on the glycerol backbone of DAG. However, how the two substrates enter DGAT1's catalytic reaction chamber and interact with DGAT1 remains elusive. This study aims to explore the structural basis of DGAT1's substrate recognition by investigating each substrate's pathway to the reaction chamber. Using a human DGAT1 cryo‐EM structure in complex with an oleoyl‐CoA substrate, we designed two different all‐atom molecular dynamics (MD) simulation systems: DGAT1away (both acyl‐CoA and DAG away from the reaction chamber) and DGAT1bound (acyl‐CoA bound in and DAG away from the reaction chamber). Our DGAT1away simulations reveal that acyl‐CoA approaches the reaction chamber via interactions with positively charged residues in transmembrane helix 7. DGAT1bound simulations show DAGs entering into the reaction chamber from the cytosol leaflet. The bound acyl‐CoA's fatty acid lines up with the headgroup of DAG, which appears to be competent to TAG formation. We then converted them into TAG and coenzyme (CoA) and used adaptive biasing force (ABF) simulations to explore the egress pathways of the products. We identify their escape routes, which are aligned with their respective entry pathways. Visualization of the substrate and product pathways and their interactions with DGAT1 is expected to guide future experimental design to better understand DGAT1 structure and function.

Guvench, Olgun; Straffin, Andrew L.

doi: 10.1002/jcc.70110pmid: 40251812

The carbon and oxygen atoms of tetrahydropyran form the common substructure of pyranose monosaccharides in vertebrate glycans. This substructure can assume various ring puckering chair and skew‐boat conformations, and can thereby impact glycan conformations relevant for biomolecular structure and signaling. The Protein Data Bank (PDB) provides a wealth of experimental glycan structural biology data that can be useful in the development and validation of molecular mechanics force fields for these molecules. However, these experimental data are typically from solvent‐depleted crystalline environments at very low temperatures, in contrast to biological conditions that are aqueous and near ambient temperature, which is the regime targeted by biomolecular force fields. To determine if these PDB X‐ray crystal data can be of utility as references for carbohydrate force fields, we compared ring puckering conformations from these experimental data to both vacuum and explicit aqueous solvent puckering free energy data from extended‐system adaptive biasing force (eABF) molecular dynamics simulations using the previously validated CHARMM36 force field. We found that, for monosaccharides that are not charged (glucose, N‐acetylglucosamine, galactose, N‐acetylgalactosamine, mannose, xylose, and fucose), both the vacuum and aqueous simulation puckering preferences strongly correlate with PDB data, and therefore with each other. In contrast, all charged monosaccharides that were considered (the conjugate bases of N‐acetylneuraminic acid, glucuronic acid, and iduronic acid) had puckering preferences correlating with PDB data only in aqueous simulations and not in vacuum simulations. These results suggest that comparing puckering preferences from aqueous simulations to PDB X‐ray crystal puckering conformation data can be a valid and useful component of carbohydrate force field development and validation.

Vyboishchikov, Sergei F.

doi: 10.1002/jcc.70104pmid: 40251880

This paper introduces AtomicESE, an artificial neural network for calculating solvation‐free energies ΔG°solv of molecules in organic solvents. AtomicESE calculates ΔG°solv by summing atomic contributions, each evaluated by a dense neural network. This atomic network uses 13 physically relevant input features, comprising six local atomic features, two global charge‐related molecular properties, and five solvent‐specific properties. For neutral solutes, AtomicESE achieves an average RMSE below 0.6 kcal/mol, demonstrating strong performance across all solvent classes, with particularly high accuracy for aromatic, haloaromatic, alkane, and ketone solvents. AtomicESE also works reliably for ionic solutes.

Wittmann, Lukas; Garcia‐Ratés, Miquel; Riplinger, Christoph

doi: 10.1002/jcc.70099pmid: 40272182

Within this work, we present the derivation and implementation of analytical gradients for the Gaussian‐switching (SwiG) Conductor‐like Polarizable Continuum Model (CPCM) with general nuclear coordinate‐dependent non‐static radii used for the creation of van der Waals‐type cavities. This is done using the recently presented dynamic radii adjustment for continuum solvation (Draco) scheme. This allows for efficient geometry optimization and reasonable numerical Hessian calculations. The derived gradient is implemented in ORCA, and therefore is easily applicable. The derivation and implementation is validated by comparing analytical and numerical gradients and testing geometry optimizations on a diverse test set, including small organic compounds, metal‐organic complexes, and highly charged species. We additionally test the continuity of the potential energy surface using an example where very strong changes in the radii occur. The computational efficiency of the derived gradient is investigated.

Chahal, Pooja; Dhilip Kumar, T. J.

doi: 10.1002/jcc.70103pmid: 40251868

Cyanopolyyne and protonated‐dicyanopolyyne molecules always get special attention for their detection in the interstellar medium. The rotational quantum dynamics for the collision of recently detected protonated dicyanoacetylene (NC4NH+$$ {\mathrm{NC}}_4{\mathrm{NH}}^{+} $$) with He is studied to get the inelastic rate coefficients till temperature range of 100 K. An accurate potential energy surface (PES), computed using ab initio methods, has been developed for the NC4NH+$$ {\mathrm{NC}}_4{\mathrm{NH}}^{+} $$–He collision system. The PES was developed with the coupled cluster, that is, the CCSD(T)‐F12b method in combination with the aug‐cc‐pVTZ basis set. The 2D PES has a global minimum with a value of −239.19 cm−1$$ {\mathrm{cm}}^{-1} $$. The analytical fitting of this 2D PES is done to obtain the radial coefficients, that give cross‐sections for NC4NH+$$ {\mathrm{NC}}_4{\mathrm{NH}}^{+} $$ molecule till collisional energy range of 300 cm−1$$ {\mathrm{cm}}^{-1} $$. The rate coefficients are achieved for the first 20 rotational transitions. An important trend is observed when comparing the de‐excitation rate coefficients at different temperatures. For transitions below Δj=10$$ \Delta j=10 $$, a preference for odd Δj$$ \Delta j $$ values is evident, which can be attributed to the anisotropy in the PES of the NC4NH+$$ {\mathrm{NC}}_4{\mathrm{NH}}^{+} $$–He collision. This similar behavior is observed for HC3NH+$$ {\mathrm{HC}}_3{\mathrm{NH}}^{+} $$–He collision. However, for higher transitions, a strong propensity for even Δj$$ \Delta j $$ transitions emerges. The results obtained in the present work will enable us to estimate the abundance of NC4NH+$$ {\mathrm{NC}}_4{\mathrm{NH}}^{+} $$ in the ISM under non‐local thermal equilibrium conditions.

Paz, Holden; Beck, Silvan; Lee, Richmond; Ho, Junming; Yu, Haibo

doi: 10.1002/jcc.70109pmid: 40251879

Molecular properties in combined quantum mechanics and molecular mechanics (QM/MM) simulations have been shown to be dependent on the size of the quantum mechanical (QM) region and the amount of conformational sampling. Previous studies have largely focused on enzymatic systems, which have made it difficult to distinguish the effects of QM region size and conformational sampling from other factors including QM‐MM boundary artifacts and the boundary effects. This study uses the difference‐based adaptive solvation QM/MM method to investigate the tautomerization reactions of alanine and aspartate in explicit solvent. The choice of computationally tractable systems enables the decoupling of QM region size effects from other factors and a direct comparison of free energy surfaces with potential energy surfaces (PESs). The results show that (1) it is crucial to properly account for thermal fluctuations along the reaction pathways, and (2) free energy surfaces converge rapidly with increasing QM region size, whereas charge transfer requires a slightly larger QM region to achieve convergence. These findings are expected to guide future studies of enzymatic systems and other complex systems where QM/MM methods are applied.

Zhang, Runhua; Zhang, Xin; Zhao, Shulin; Zou, Quan; Ding, Yijie; Guo, Xiaoyi; Wu, Hongjie

doi: 10.1002/jcc.70111pmid: 40271912

In pursuit of unraveling novel structural inhibitors for treating monkeypox virus, targeting the VP37 protein, which is bioactive in response to ST‐246, to discern pharmaceutical molecules specifically tailored to combat monkeypox virus. We employed a semi‐flexible molecular docking, molecular dynamic simulation, and ADME screening methodology, which are based on structure, to screen compounds from CMNPD and TCM in silico. These methodologies allowed us to find potential candidates depending on their binding values and interactions with the binding site of main protease. To further evaluate the stability of these interactions, we conducted molecular dynamics simulations and calculated binding energies. Herein, employing methods such as binding energy calculations, comparative analyses, and molecular dynamics simulations for activity computations, the six top hits of the compounds were validated as five kinds of good inhibitors, surpassing its reference compound ST‐246, for better in vitro drug candidates against MPXV.

Krishna, ; Saini, Lalit Kumar; Pandey, Mukesh

doi: 10.1002/jcc.70114pmid: 40251885

Ab initio and density functional theory (DFT) calculations are employed to investigate LiH:nNH3 (n = 1–4) cluster complexes. The nature of the interactions is analyzed using molecular electrostatic potential maps, quantum theory of atoms in molecules, delocalization indices, and electron density difference maps. In the presence of LiH, NH3 molecules engage in several types of noncovalent interactions, namely, hydrogen bonding (HB), lithium bonding (LB), and dihydrogen bonding (DHB). The LiH:NH3 dimer is stabilized primarily through Li···N interactions. The role of these noncovalent interactions in complexes having more than one NH3 molecule, for example, hetero‐trimer, tetramer, and pentamer structures, is also examined. Increasing the number of NH3 molecules enhances the number of HB sites. Additionally, the strengths of LB and DHB interactions associated with HB‐bonded NH3 molecules increase. Interaction energy estimates and many‐body energy decomposition analysis suggest that increasing NH3 molecules increases cooperativity, approaching ~10% of the total interaction's energy in the case of tetramers and pentamers.

Liao, Jun; Wu, Mincong; Meng, Fanjun; Chen, Changjun

doi: 10.1002/jcc.70113pmid: 40249089

Residue interaction networks determine various characteristics of proteins, such as the folding rate, thermostability, and allosteric process. The interactions between residues can be described by distances or energies. The former is simple but less rigorous. The latter is complicated but more precise, especially when considering the solvent effect. In this work, we apply an existing energy decomposition method based on the Poisson–Boltzmann equation solver. The calculation is especially accelerated on GPU for higher performance. In four formal applications, the constructed interaction energy (IE) network shows good results. First, it is found that the protein folding rate has a stronger correlation with the energy‐based contact order than the distance‐based contact order. The Pearson correlation coefficient (PCC) is 0.839 versus 0.784 on a dataset of non‐two‐state proteins. Second, we find that most thermophilic proteins have lower IEs than mesophilic proteins. The IE in solvent acts as an indicator to evaluate the thermostabilities of proteins. Third, we use the IE network to predict the key residues in the formation of the insulin dimer. Most key residues are in agreement with the findings in previous alanine‐scanning experiments. Lastly, we propose a novel method (called APFN) to predict the allosteric pathway based on the IE network. The method gives the same allosteric pathway for CheY protein as in previous nuclear magnetic resonance spectroscopy experiments. On the whole, the IE network in the solvent has been demonstrated to be reliable in describing the characteristics embedded in protein structures.

Stanley, Katherine; Givhan, R. Houston; Turney, Justin M.; Schaefer, Henry F.

doi: 10.1002/jcc.70105pmid: 40251860

The molecule 5,6‐diaminobenzene‐1,2,3,4‐tetracarbonnitrile (MOI) was first synthesized by Müllen and coworkers in 2016 and boasts an ultrastrong dipole moment of 14.1±0.7$$ 14.1\pm 0.7 $$ Debye in THF. Gas phase DFT computations do not fully reflect this ultrastrong dipole moment, demonstrating the role of solvent in increasing this dipole moment. Here, we investigate the effect of solvent molecule position on the dipole moment of this species, computationally examining systems with giant dipole moments. These systems are optimized in the gas phase with the B3LYP functional, employing the aug‐cc‐pVTZ and def2‐TZVP basis sets, as well as the B3LYP‐D3BJ/aug‐cc‐pVTZ functional in Orca. Single point DLPNO‐CCSD/aug‐cc‐pVDZ results were obtained from Orca and Psi4, as well as DLPNO‐CCSD(T)/CBS information from Psi4. Additionally, these are compared to the dipole moments of di‐ and tri‐hydrated systems, and the SMD models for THF and water at the B3LYP/aug‐cc‐pVTZ level of theory. The dissociation energies, HOMO‐LUMO energy gaps, and dipole moments are presented. These metrics show the nh1nh1′ THF system boasts the largest dissociation energy and dipole moment of the singly solvated systems, due to its strong hydrogen bonding. The importance of solvent placement is highlighted and may guide the synthesis of macromolecules or organic frameworks incorporating the MOI or MOI‐like subunits. Remarkably, a single solvent molecule provides a good model for the difference between the gas phase and solvated species. The predicted gas phase dipole moments computed with B3LYP/aug‐cc‐pVTZ for the MOI, its monohydrated complex, dihydrated complex, and its trihydrated complex are 9.6, 14.2, 16.0, and 16.8 Debye, respectively.