Deep Learning Total Energies and Orbital Energies of Large Organic Molecules Using Hybridization of Molecular Fingerprints.

The ability to predict material properties without the need for resource-consuming experimental efforts can immensely accelerate material and drug discovery. Although ab initio methods can be reliable and accurate in making such predictions, they are computationally too expensive on a large scale. The recent advancements in artificial intelligence and machine learning as well as the availability of large quantum mechanics derived datasets enable us to train models on these datasets as a benchmark and to make fast predictions on much larger datasets. The success of these machine learning models highly depends on the machine-readable fingerprints of the molecules that capture their chemical properties as well as topological information. In this work, we propose a common deep learning-based framework to combine different types of molecular fingerprints to enhance prediction accuracy. A graph neural network (GNN), many-body tensor representation (MBTR), and a set of simple molecular descriptors (MD) were used to predict the total energies, highest occupied molecular orbital (HOMO) energies, and lowest unoccupied molecular orbital (LUMO) energies of a dataset containing ∼62k large organic molecules with complex aromatic rings and remarkably diverse functional groups. The results demonstrate that a combination of best performing molecular fingerprints can produce better results than the individual ones. The simple and flexible deep learning framework developed in this work can be easily adapted to incorporate other types of molecular fingerprints.

Configurational Disorder of Water Hydrogen-Bond Network at the Protein Dynamical Transition.

We introduce a novel strategy to quantify the disorder of extended water-water hydrogen-bond (HB) networks sampled in particle-based computer simulations. The method relies on the conformational clustering of the HB connectivity states. We successfully applied it to unveil the fine relationship among the protein dynamical transition in hydrated powder, which marks the activation of protein flexibility at Td ≈ 240 K, and the sudden increase in the configurational disorder of the water HB network enveloping the proteins. Our finding links, in the spirit of the Adam-Gibbs relationship, the diffusivity of protein atoms, as quantified by the hydrogen mean-square displacements, and the thermodynamic solvent configurational entropy.

A structural insight into mechanical strength of graphene-like carbon and carbon nitride networks.

Graphene, one of the strongest materials ever discovered, triggered the exploration of many 2D materials in the last decade. However, the successful synthesis of a stable nanomaterial requires a rudimentary understanding of the relationship between its structure and strength. In the present study, we investigate the mechanical properties of eight different carbon-based 2D nanomaterials by performing extensive density functional theory calculations. The considered structures were just recently either experimentally synthesized or theoretically predicted. The corresponding stress-strain curves and elastic moduli are reported. They can be useful in training force field parameters for large scale simulations. A comparative analysis of these results revealed a direct relationship between atomic density per area and elastic modulus. Furthermore, for the networks that have an armchair and a zigzag orientation, we observed that they were more stretchable in the zigzag direction than the armchair direction. A critical analysis of the angular distributions and radial distribution functions suggested that it could be due to the higher ability of the networks to suppress the elongations of the bonds in the zigzag direction by deforming the bond angles. The structural interpretations provided in this work not only improve the general understanding of a 2D material's strength but also enables us to rationally design them for higher qualities.

Mechanical responses of borophene sheets: a first-principles study.

Recent experimental advances for the fabrication of various borophene sheets introduced new structures with a wide range of applications. Borophene is the boron atom analogue of graphene. Borophene exhibits various structural polymorphs all of which are metallic. In this work, we employed first-principles density functional theory calculations to investigate the mechanical properties of five different single-layer borophene sheets. In particular, we analyzed the effect of the loading direction and point vacancy on the mechanical response of borophene. Moreover, we compared the thermal stabilities of the considered borophene systems. Based on the results of our modelling, borophene films depending on the atomic configurations and the loading direction can yield a remarkable elastic modulus in the range of 163-382 GPa nm and a high ultimate tensile strength from 13.5 GPa nm to around 22.8 GPa nm at the corresponding strain from 0.1 to 0.21. Our study reveals the remarkable mechanical characteristics of borophene films.

Role of Internal Water on Protein Thermal Stability: The Case of Homologous G Domains.

In this work, we address the question of whether the enhanced stability of thermophilic proteins has a direct connection with internal hydration. Our model systems are two homologous G domains of different stability: the mesophilic G domain of the elongation factor thermal unstable protein from E. coli and the hyperthermophilic G domain of the EF-1α protein from S. solfataricus. Using molecular dynamics simulation at the microsecond time scale, we show that both proteins host water molecules in internal cavities and that these molecules exchange with the external solution in the nanosecond time scale. The hydration free energy of these sites evaluated via extensive calculations is found to be favorable for both systems, with the hyperthermophilic protein offering a slightly more favorable environment to host water molecules. We estimate that, under ambient conditions, the free energy gain due to internal hydration is about 1.3 kcal/mol in favor of the hyperthermophilic variant. However, we also find that, at the high working temperature of the hyperthermophile, the cavities are rather dehydrated, meaning that under extreme conditions other molecular factors secure the stability of the protein. Interestingly, we detect a clear correlation between the hydration of internal cavities and the protein conformational landscape. The emerging picture is that internal hydration is an effective observable to probe the conformational landscape of proteins. In the specific context of our investigation, the analysis confirms that the hyperthermophilic G domain is characterized by multiple states and it has a more flexible structure than its mesophilic homologue.

How conformational flexibility stabilizes the hyperthermophilic elongation factor G-domain.

Proteins from thermophilic organisms are stable and functional well above ambient temperature. Understanding the molecular mechanism underlying such a resistance is of crucial interest for many technological applications. For some time, thermal stability has been assumed to correlate with high mechanical rigidity of the protein matrix. In this work we address this common belief by carefully studying a pair of homologous G-domain proteins, with their melting temperatures differing by 40 K. To probe the thermal-stability content of the two proteins we use extensive simulations covering the microsecond time range and employ several different indicators to assess the salient features of the conformational landscape and the role of internal fluctuations at ambient condition. At the atomistic level, while the magnitude of fluctuations is comparable, the distribution of flexible and rigid stretches of amino-acids is more regular in the thermophilic protein causing a cage-like correlation of amplitudes along the sequence. This caging effect is suggested to favor stability at high T by confining the mechanical excitations. Moreover, it is found that the thermophilic protein, when folded, visits a higher number of conformational substates than the mesophilic homologue. The entropy associated with the occupation of the different substates and the thermal resilience of the protein intrinsic compressibility provide a qualitative insight on the thermal stability of the thermophilic protein as compared to its mesophilic homologue. Our findings potentially open the route to new strategies in the design of thermostable proteins.

Evaluation of several two-step scoring functions based on linear interaction energy, effective ligand size, and empirical pair potentials for prediction of protein-ligand binding geometry and free energy.

The performances of several two-step scoring approaches for molecular docking were assessed for their ability to predict binding geometries and free energies. Two new scoring functions designed for "step 2 discrimination" were proposed and compared to our CHARMM implementation of the linear interaction energy (LIE) approach using the Generalized-Born with Molecular Volume (GBMV) implicit solvation model. A scoring function S1 was proposed by considering only "interacting" ligand atoms as the "effective size" of the ligand and extended to an empirical regression-based pair potential S2. The S1 and S2 scoring schemes were trained and 5-fold cross-validated on a diverse set of 259 protein-ligand complexes from the Ligand Protein Database (LPDB). The regression-based parameters for S1 and S2 also demonstrated reasonable transferability in the CSARdock 2010 benchmark using a new data set (NRC HiQ) of diverse protein-ligand complexes. The ability of the scoring functions to accurately predict ligand geometry was evaluated by calculating the discriminative power (DP) of the scoring functions to identify native poses. The parameters for the LIE scoring function with the optimal discriminative power (DP) for geometry (step 1 discrimination) were found to be very similar to the best-fit parameters for binding free energy over a large number of protein-ligand complexes (step 2 discrimination). Reasonable performance of the scoring functions in enrichment of active compounds in four different protein target classes established that the parameters for S1 and S2 provided reasonable accuracy and transferability. Additional analysis was performed to definitively separate scoring function performance from molecular weight effects. This analysis included the prediction of ligand binding efficiencies for a subset of the CSARdock NRC HiQ data set where the number of ligand heavy atoms ranged from 17 to 35. This range of ligand heavy atoms is where improved accuracy of predicted ligand efficiencies is most relevant to real-world drug design efforts.

Development of a ReaxFF reactive force field for glycine and application to solvent effect and tautomerization.

Tautomerization of amino acids between the neutral form (NF) and the zwitterionic form (ZW) in water has been extensively studied, often using glycine as a model to understand this fundamental process. In spite of many advanced studies, the tautomerization reaction remains poorly understood because of the intrinsic complexities of the system, including multiple accessible reaction pathways, charge transfer, and variations of solvation structure. To establish an accurate model that can be used for molecular dynamics simulations, a ReaxFF reactive force field has been developed for glycine. A training set for the ReaxFF hydrocarbon potential was augmented with several glycine conformers and glycine-water complexes. The force field parameters were optimized to reproduce the quantum mechanically derived energies of the species in the training set. The optimized potential could accurately describe the properties of gas-phase glycine. It was applied to investigate the effect of solvation on the conformational distribution of glycine. Molecular dynamics simulations indicated significant differences in the dominant conformers in the gas phase and in water. This suggests that the tautomerization of glycine occurs through a conformational isomerization followed by the proton transfer event. The direct reaction mechanism of the NF → ZW proton transfer reaction in water, as well as mechanisms mediated by one or two water molecules, were investigated using molecular dynamics simulations. The results suggest that the proton transfer reaction is most likely mediated by a single water molecule. The ReaxFF potential developed in this work provides an accurate description of proton transfer in glycine and thus provides a useful methodology for simulating proton transfer reactions in organic molecules in the aqueous environment.

Development and validation of a ReaxFF reactive force field for Cu cation/water interactions and copper metal/metal oxide/metal hydroxide condensed phases.

To enable large-scale reactive dynamic simulations of copper oxide/water and copper ion/water interactions we have extended the ReaxFF reactive force field framework to Cu/O/H interactions. To this end, we employed a multistage force field development strategy, where the initial training set (containing metal/metal oxide/metal hydroxide condensed phase data and [Cu(H(2)O)(n)](2+) cluster structures and energies) is augmented by single-point quantum mechanices (QM) energies from [Cu(H(2)O)(n)](2+) clusters abstracted from a ReaxFF molecular dynamics simulation. This provides a convenient strategy to both enrich the training set and to validate the final force field. To further validate the force field description we performed molecular dynamics simulations on Cu(2+)/water systems. We found good agreement between our results and earlier experimental and QM-based molecular dynamics work for the average Cu/water coordination, Jahn-Teller distortion, and inversion in [Cu(H(2)O)(6)](2+) clusters and first- and second-shell O-Cu-O angular distributions, indicating that this force field gives a satisfactory description of the Cu-cation/water interactions. We believe that this force field provides a computationally convenient method for studying the solution and surface chemistry of metal cations and metal oxides and, as such, has applications for studying protein/metal cation complexes, pH-dependent crystal growth/dissolution, and surface catalysis.

Development of a ReaxFF reactive force field for aqueous chloride and copper chloride.

Copper ions play crucial roles in many enzymatic and aqueous processes. A critical analysis of the fundamental properties of copper complexes is essential to understand their impact on a wide range of chemical interactions. However the study of copper complexes is complicated by the presence of strong polarization and charge transfer effects, multiple oxidation states, and quantum effects like Jahn-Teller distortions. These complications make the experimental observations difficult to interpret. In order to provide a computationally inexpensive yet reliable method for simulation of aqueous-phase copper chemistry, ReaxFF reactive force field parameters have been developed. The force field parameters have been trained against a large set of DFT-derived energies for condensed-phase copper-chloride clusters as well as chloride/water and copper-chloride/water clusters sampled from molecular dynamics (MD) simulations. The parameters were optimized by iteratively training them against configurations generated from ReaxFF MD simulations that are performed multiple times with improved sets of parameters. This cycle was repeated until the ReaxFF results were in accordance with the DFT-derived values. We have performed MD simulations on chloride/water and copper-chloride/water systems to validate the optimized force field. The structural properties of the chloride/water system are in accord with previous experimental and computational studies. The properties of copper-chloride/water agreed with the experimental observations including evidence of the Jahn-Teller distortion. The results of this study demonstrate the applicability of ReaxFF for the precise characterization of aqueous copper chloride. This force field provides a base for the design of a computationally inexpensive tool for the investigation of various properties and functions of metal ions in industrial, environmental, and biological environments.

Molecular dynamics simulations of a DMPC bilayer using nonadditive interaction models.

We present a polarizable force field based on the charge-equilibration formalism for molecular dynamics simulations of phospholipid bilayers. We discuss refinement of headgroup dihedral potential parameters to reproduce ab initio conformational energies of dimethylphosphate calculated at the MP2/cc-pVTZ level of theory. We also address the refinement of electrostatic and Lennard-Jones (van der Waals) parameters to reproduce ab initio polarizabilities and water interaction energies of dimethylphosphate and tetramethylammonium. We present results of molecular dynamics simulations of a solvated dimyristoylphosphatidylcholine bilayer using this polarizable force field as well as the nonpolarizable, fixed-charge CHARMM27 and CHARMM27r force fields for comparison. Calculated atomic and electron-density profiles, deuterium order parameters, and headgroup orientations are found to be consistent with previous simulations and with experiment. Polarizable interaction models for solvent and lipid exhibit greater water penetration into the lipid interior; this is due to the variation of water molecular dipole moment from a bulk value of 2.6 Debye to a value of 1.9 Debye in the membrane interior. The reduction in the electrostatic component of the desolvation free-energy penalty allows for greater water density. The surface dipole potential predicted by the polarizable model is 0.95 V compared to the value of 0.8 V based on nonpolarizable force-field calculations. Effects of inclusion of explicit polarization are discussed in relation to water dipole moment and varying charge distributions. Dielectric permittivity profiles for polarizable and nonpolarizable interactions exhibit subtle differences arising from the nature of the individual component parameterizations; for the polarizable force field, we obtain a bulk dielectric permittivity of 79, whereas the nonpolarizable force field plateaus at 97 (the value for pure TIP3P water). In the membrane interior, both models predict unit permittivities, with the polarizable models contributing from one to two more units due to the optical dielectric (high-frequency dipole fluctuations). This contribution is a step toward the continuing development of a CHARMM (Chemistry at Harvard Molecular Mechanics) polarizable force field for simulations of biomacromolecular systems.

Unfolding pathways of human serum albumin: evidence for sequential unfolding and folding of its three domains.

Human serum albumin (HSA) contains three alpha-helical domains (I-III). The unfolding process of these domains was monitored using covalently bound fluorescence probes; domain I was monitored by N-(1-pyrene)maleimide (PM) conjugated with cys-34, domain II was monitored by the lone tryptophan residue and domain III was followed by p-nitrophenyl anthranilate (NPA) conjugated with Tyrosine-411 (Tyr-411). Using domain-specific probes, we found that guanidium hydrochloride-induced unfolding of HSA occurred sequentially. The unfolding of domain II preceded that of domain I and the unfolding of domain III followed that of domain I. In addition, the domains I and III refolded within the dead time of the fluorescence recovery experiment while the refolding of domain II occurred slowly. The results suggest that individual domain of a multi-domain protein can fold and unfold sequentially.

Multiple-probe analysis of folding and unfolding pathways of human serum albumin. Evidence for a framework mechanism of folding.

The changes in the far-UV CD signal, intrinsic tryptophan fluorescence and bilirubin absorbance showed that the guanidine hydrochloride (GdnHCl)-induced unfolding of a multidomain protein, human serum albumin (HSA), followed a two-state process. However, using environment sensitive Nile red fluorescence, the unfolding and folding pathways of HSA were found to follow a three-state process and an intermediate was detected in the range 0.25-1.5 m GdnHCl. The intermediate state displayed 45% higher fluorescence intensity than that of the native state. The increase in the Nile red fluorescence was found to be due to an increase in the quantum yield of the HSA-bound Nile red. Low concentrations of GdnHCl neither altered the binding affinity of Nile red to HSA nor induced the aggregation of HSA. In addition, the secondary structure of HSA was not perturbed during the first unfolding transition (<1.5 m GdnHCl); however, the secondary structure was completely lost during the second transition. The data together showed that the half maximal loss of the tertiary structure occurred at a lower GdnHCl concentration than the loss of the secondary structure. Further kinetic studies of the refolding process of HSA using multiple spectroscopic techniques showed that the folding occurred in two phases, a burst phase followed by a slow phase. An intermediate with native-like secondary structure but only a partial tertiary structure was found to form in the burst phase of refolding. Then, the intermediate slowly folded into the native state. An analysis of the refolding data suggested that the folding of HSA could be best explained by the framework model.