Molecular scale theories and simulations of fluid systems: from microscopic to macroscopic effects

Amphetamine-type stimulants (ATS) are a group of incitation and psychedelic drugs affecting the central nervous system. Physicochemical data for these compounds are essential for understanding the stimulating mechanism, for assessing their environmental impacts, and for developing new drug detection methods. However, experimental data are scarce due to tight regulation of such illicit drugs, yet conventional methods to estimate their properties are often unreliable. Here we introduce a tailor-made multiscale procedure for predicting the hydration free energies and the solvation structures of ATS molecules by a combination of first principles calculations and the classical density functional theory. We demonstrate that the multiscale procedure performs well for a training set with similar molecular characteristics and yields good agreement with a testing set not used in the training. The theoretical predictions serve as a benchmark for the missing experimental data and, importantly, provide microscopic insights into manipulating the hydrophobicity of ATS compounds by chemical modifications.

Implicit solvent models offer an attractive way to estimate the effects of a solvent environment on the properties of small or large solutes without the complications of explicit simulations. One common test of accuracy is to compute the free energy of transfer from gas to liquid for a variety of small molecules, since many of these values have been measured. Studies of the temperature dependence of these values (i.e. solvation enthalpies and entropies) can provide additional insights into the performance of implicit solvent models. Here, we show how to compute temperature derivatives of hydration free energies for the 3D-RISM integral equation approach. We have computed hydration free energies of 1123 small drug-like molecules (both neutral and charged). Temperature derivatives were also used to calculate hydration energies and entropies of 74 of these molecules (both neutral and charged) for which experimental data is available. While direct results have rather poor agreement with experiment, we have found that several previously proposed linear hydration free energy correction schemes give good agreement with experiment. These corrections also provide good agreement for hydration energies and entropies though simple extensions are required in some cases.

For neutral hard-sphere solutes, we compare the reduced density profile of water around a solute g(r), solvation free energy μ, energy U, and entropy S under the isochoric condition predicted by the two theories: dielectrically consistent reference interaction site model (DRISM) and angle-dependent integral equation (ADIE) theories. A molecular model for water pertinent to each theory is adopted. The hypernetted-chain (HNC) closure is employed in the ADIE theory, and the HNC and Kovalenko–Hirata (K–H) closures are tested in the DRISM theory. We also calculate g(r), U, S, and μ of the same solute in a hard-sphere solvent whose molecular diameter and number density are set at those of water, in which case the radial-symmetric integral equation (RSIE) theory is employed. The dependences of μ, U, and S on the excluded volume and solvent-accessible surface area are analyzed using the morphometric approach (MA). The results from the ADIE theory are in by far better agreement with those from computer simulations available for g(r), U, and μ. For the DRISM theory, g(r) in the vicinity of the solute is quite high and becomes progressively higher as the solute diameter d_U increases. By contrast, for the ADIE theory, it is much lower and becomes further lower as d_U increases. Due to unphysically positive U and significantly larger |S|, μ from the DRISM theory becomes too high. It is interesting that μ, U, and S from the K–H closure are worse than those from the HNC closure. Overall, the results from the DRISM theory with a molecular model for water are quite similar to those from the RSIE theory with the hard-sphere solvent. Based on the results of the MA analysis, we comparatively discuss the different theoretical methods for cases where they are applied to studies on the solvation of a protein.

Complex formation between molecules in solution is the key process by which molecular interactions are translated into functional systems. These processes are governed by the binding or free energy of association which depends on both direct molecular interactions and the solvation contribution. A design goal frequently addressed in pharmaceutical sciences is the optimization of chemical properties of the complex partners in the sense of minimizing their binding free energy with respect to a change in chemical structure. Here, we demonstrate that liquid-state theory in the form of the solute–solute equation of the reference interaction site model provides all necessary information for such a task with high efficiency. In particular, computing derivatives of the potential of mean force (PMF), which defines the free-energy surface of complex formation, with respect to potential parameters can be viewed as a means to define a direction in chemical space toward better binders. We illustrate the methodology in the benchmark case of alkali ion binding to the crown ether 18-crown-6 in aqueous solution. In order to examine the validity of the underlying solute–solute theory, we first compare PMFs computed by different approaches, including explicit free-energy molecular dynamics simulations as a reference. Predictions of an optimally binding ion radius based on free-energy derivatives are then shown to yield consistent results for different ion parameter sets and to compare well with earlier, orders-of-magnitude more costly explicit simulation results. This proof-of-principle study, therefore, demonstrates the potential of liquid-state theory for molecular design problems.

The molecular recognition process of the carbohydrate-binding module family 36 (CBM36) was examined theoretically. The mechanism of xylan binding by CBM36 and the role of Ca²⁺ were investigated by the combined use of molecular dynamics simulations and the 3D reference interaction site model method. The CBM36 showed affinity for xylan after Ca²⁺ binding, but not after Mg²⁺ binding. Free-energy component analysis of the xylan-binding process revealed that the major factor for xylan-binding affinity is the electrostatic interaction between the Ca²⁺ and the hydroxyl oxygens of xylan. The van der Waals interaction between the hydrophobic side chain of CBM36 and the glucopyranose ring of xylan also contributes to the stabilization of the xylan-binding state. Dehydration on the formation of the complex has the opposite effect on these interactions. The affinity of CBM36 for xylan results from a balance of the interactions between the binding ion and solvents, hydrophilic residues around xylan, and the hydroxyl oxygens of xylan. When CBM binds Ca²⁺, these interactions are well balanced; in contrast, when CBM binds Mg²⁺, the dehydration penalty is excessively large.

In this work we have assessed the ability of a recently proposed three-dimensional integral equation approach to describe the explicit spatial distribution of molecular hydrogen confined in a crystal formed by short-capped nanotubes of C₅₀ H₁₀. To that aim we have resorted to extensive molecular simulation calculations whose results have been compared with our three-dimensional integral equation approximation. We have first tested the ability of a single C₅₀ H₁₀ nanocage for the encapsulation of H₂ by means of molecular dynamics simulations, in particular using targeted molecular dynamics to estimate the binding Gibbs energy of a host hydrogen molecule inside the nanocage. Then, we have investigated the adsorption isotherm of the nanocage crystal using grand canonical Monte Carlo simulations in order to evaluate the maximum load of molecular hydrogen. For a packing close to the maximum load explicit hydrogen density maps and density profiles have been determined using molecular dynamics simulations and the three-dimensional Ornstein–Zernike equation with a hypernetted chain closure. In these conditions of extremely tight confinement the theoretical approach has shown to be able to reproduce the three-dimensional structure of the adsorbed fluid with accuracy down to the finest details.

The interaction between any two biological molecules must compete with their interaction with water molecules. This makes water the most important molecule in medicine, as it controls the interactions of every therapeutic with its target. A small molecule binding to a protein is able to recognize a unique binding site on a protein by displacing bound water molecules from specific hydration sites. Quantifying the interactions of these water molecules allows us to estimate the potential of the protein to bind a small molecule. This is referred to as ligandability. In the study, we describe a method to predict ligandability by performing a search of all possible combinations of hydration sites on protein surfaces. We predict ligandability as the summed binding free energy for each of the constituent hydration sites, computed using inhomogeneous fluid solvation theory. We compared the predicted ligandability with the maximum observed binding affinity for 20 proteins in the human bromodomain family. Based on this comparison, it was determined that effective inhibitors have been developed for the majority of bromodomains, in the range from 10 to 100 nM. However, we predict that more potent inhibitors can be developed for the bromodomains BPTF and BRD7 with relative ease, but that further efforts to develop inhibitors for ATAD2 will be extremely challenging. We have also made predictions for the 14 bromodomains with no reported small molecule K_d values by isothermal titration calorimetry. The calculations predict that PBRM1(1) will be a challenging target, while others such as TAF1L(2), PBRM1(4) and TAF1(2), should be highly ligandable. As an outcome of this work, we assembled a database of experimental maximal K_d that can serve as a community resource assisting medicinal chemistry efforts focused on BRDs. Effective prediction of ligandability would be a very useful tool in the drug discovery process.

As an example of charged, dipolar soft matter, the ionic liquid 1-ethyl-3-methyl-imidazolium dicyanamide is studied by coarse-grained molecular dynamics simulations. We focus on the link between microscopic and mesoscopic properties for both structure and dynamics. Thereby, the generalized Kirkwood g_K-factor plays a central role in establishing this link which is not possible on the basis of molecular hydrodynamics. The decoupling between translational and rotational motion is indicative of the dynamical heterogeneity in ionic liquids.

A model is proposed for liquid water that is a continuum between the ordered state with predominantly tetrahedral coordination, linear hydrogen bonds and activated dynamics and a disordered state with a continuous distribution of multiple coordinations, multiple types of hydrogen bond, and diffusive dynamics, similar to that of normal liquids. Central to water's heterogeneous structure is the ability of hydrogen to donate to either one acceptor in a conventional linear hydrogen bond or to multiple acceptors as a furcated hydrogen. Linear hydrogen bonds are marked by slow, activated kinetics for hydrogen-bond switching to more crowded acceptors and sharp first peaks in the hydrogen-oxygen radial distribution function. Furcated hydrogens, equivalent to free, broken, dangling or distorted hydrogens, have barrierless, rapid kinetics and poorly defined first peaks in their hydrogen-oxygen radial distribution function. They involve the weakest donor in a local excess of donors, such that barrierless whole-molecule vibration rapidly swaps them between the linear and furcated forms. Despite the low number of furcated hydrogens and their transient existence, they are readily created in a single hydrogen-bond switch and free up the dynamics of numerous surrounding molecules, bringing about the disordered state. Hydrogens in the ordered state switch with activated dynamics to make the non-tetrahedral coordinations of the disordered state, which can also combine to make the ordered state. Consequently, the ordered and disordered states are both connected by diffusive dynamics and differentiated by activated dynamics, bringing about water's continuous heterogeneity.

A protocol to calculate the binding free energy of a host–guest system is proposed based on the MM/3D-RISM method, taking cyclodextrin derivatives and their ligands as model systems. The protocol involves the procedure to identify the most probable binding mode (MPBM) of receptors and ligands by means of the umbrella sampling method. The binding free energies calculated by the MM/3D-RISM method for the complexes of the seven ligands with the MPBM of the cyclodextrin, and with the fluctuated structures around it, are in agreement with the corresponding experimental data in a semi-quantitative manner. It suggests that the protocol proposed here is promising for predicting the binding affinity of a small ligand to a relatively rigid receptor such as cyclodextrin.

The atom-scale characterisation of interfaces between transition metal oxides and liquid water is fundamental to our mechanistic understanding of diverse phenomena ranging from crystal growth to biogeochemical transformations to solar fuel production. Here we report on the results of large-scale hybrid density functional theory-based molecular dynamics simulations for the hematite(001)-liquid water interface. A specific focus is placed on understanding how different terminations of the same surface influence surface solvation. We find that the two dominant terminations for the hematite(001) surface exhibit strong differences both in terms of the active species formed on the surface and the strength of surface solvation. According to present simulations, we find that charged oxyanions (-O⁻) and doubly protonated oxygens (-OH$_{2}^{+}$) can be formed on the iron terminated layer via autoionization of neutral -OH groups. No such charged species are found for the oxygen terminated surface. In addition, the missing iron sublayer in the iron terminated surface strongly influences the solvation structure, which becomes less well ordered in the vicinity of the interface. These pronounced differences are likely to affect the reactivity of the two surface terminations, and in particular the energetics of excess charge carriers at the surface.

Natural asphaltenes are defined as polyaromatic compounds whose chemical composition and structure are dependent on their geological origin and production history, hence are regarded as complex molecules with aromatic cores and aliphatic tails that occur in the heaviest fraction of crude oil. The aggregation of asphaltenes presents a range of technical challenges to the production and processing of oil. In this work we study the behaviour of the model asphaltene-like molecule hexa-tert-butylhexa-peri-hexabenzocoronene (HTBHBC) using molecular dynamics simulation. It was found that the regular arrangement of the tert-butyl side chains prevents the formation of strongly-bound dimers by severely restricting the configurational space of the aggregation pathway. In contrast, a modified molecule with only 3 side chains is readily able to form dimers. This work therefore confirms the influence of the molecular structure of polyaromatic compounds on their aggregation mechanism, and reveals the unexpected design rules required for model systems that can mimic the behavior of asphaltenes.

Molecular dynamics simulations of the liquid–vapour interface of dimethyl sulphoxide (DMSO)-water mixtures of 11 different compositions, including two neat systems are performed on the canonical (N, V, T) ensemble at 298 K. The molecules constituting the surface layer of these systems are selected by means of the identification of the truly interfacial molecules (ITIM) method, and their local lateral environment at the liquid surface is investigated by performing Voronoi analysis. The obtained results reveal that both molecules prefer to be in a mixed local environment, consisting of both kinds of molecules, at the liquid surface, and this preference is even stronger here than in the bulk liquid phase. Neat-like patches, in which a molecule is surrounded by like neighbours, are not found. However, vacancies that are surrounded solely by water molecules are observed at the liquid surface. Our results show that strongly hydrogen bonded DMSO·H₂O complexes, known to exist in the bulk phase of these mixtures, are absent from the liquid surface.

We propose a closure to the integral equations of molecular theory of solvation, particularly suitable for polar and charged macromolecules in electrolyte solution. This includes such systems as oligomeric polyelectrolytes at a finite concentration in aqueous and various non-aqueous solutions, as well as drug-like compounds in solution. The new closure by Kobryn, Gusarov, and Kovalenko (KGK closure) imposes the mean spherical approximation (MSA) almost everywhere in the solvation shell but levels out the density distribution function to zero (with the continuity at joint boundaries) inside the repulsive core and in the spatial regions of strong density depletion emerging due to molecular associative interactions. Similarly to MSA, the KGK closure reduces the problem to a linear equation for the direct correlation function which is predefined analytically on most of the solvation shells and has to be determined numerically on a relatively small (three-dimensional) domain of strong depletion, typically within the repulsive core. The KGK closure leads to the solvation free energy in the form of the Gaussian fluctuation (GF) functional. We first test the performance of the KGK closure coupled to the reference interaction site model (RISM) integral equations on the examples of Lennard-Jones liquids, polar and nonpolar molecular solvents, including water, and aqueous solutions of simple ions. The solvation structure, solvation chemical potential, and compressibility obtained from RISM with the KGK closure favorably compare to the results of the hypernetted chain (HNC) and Kovalenko–Hirata (KH) closures, including their combination with the GF solvation free energy. We then use the KGK closure coupled to RISM to obtain the solvation structure and thermodynamics of oligomeric polyelectrolytes and drug-like compounds at a finite concentration in electrolyte solution, for which no convergence is obtained with other closures. For comparison, we calculate their solvation structure from molecular dynamics (MD) simulations. We further couple the 3D-RISM integral equation with the 3D-version of the KGK closure, and solve it for molecular mixtures as well as oligomeric polyelectrolytes and drug-like molecules in electrolyte solutions.